Index: trunk/doc/pslib/psLibADD.tex
===================================================================
--- trunk/doc/pslib/psLibADD.tex	(revision 301)
+++ trunk/doc/pslib/psLibADD.tex	(revision 301)
@@ -0,0 +1,254 @@
+\documentclass[panstarrs]{panstarrs}
+%\documentclass[panstarrs]{panstarrs}
+
+% basic document variables
+\title{Pan-STARRS IPP Algorithm Design Description}
+\author{Eugene Magnier, Paul Price, Robert Lupton}
+\shorttitle{PSLib ADD}
+\group{Pan-STARRS Algorithm Group}
+\project{Pan-STARRS Image Processing Pipeline}
+\organization{Institute for Astronomy}
+\version{DR}
+\docnumber{PSDC-430-006}
+
+\begin{document}
+\maketitle
+
+% -- Revision History --
+% provide explicit values for the old versions
+% use '\theversion' for the current version (set above)
+% use \hline between each table row
+\RevisionsStart
+% version     Date         Description
+\theversion & 2003 Mar 11 & Hacking \\
+\RevisionsEnd
+
+\pagebreak
+\tableofcontents
+
+\pagebreak 
+\pagenumbering{arabic}
+
+\section{PanSTARRS Library PSLib}
+
+\subsection{Data Manipulation}
+
+\subsubsection{Matrix Operations}
+
+In this section, we define the linear algebra operations performed on
+matrices.  We have defined APIs to implement the following matrix
+functions:
+
+\begin{itemize}
+\item Invert a matrix;
+\item Calculate a matrix determinant;
+\item Perform matrix addition, subtraction and multiplication;
+\item Transpose a matrix; and
+\item Convert a matrix to a vector.
+\end{itemize}
+
+Many of these operations are implemented using LU decomposition.  We
+define LU decomposition in the following paragraph.  Implementation of
+LU decomposition shall make use of the GSL function
+\code{gsl_linalg_LU_decomp}.
+
+\paragraph{LU Decomposition}
+\label{LUdecomp}
+
+We wish to decompose the matrix $A$ with elements $a_{ij}$ into
+diagonal matrices that satisfy the relationship $A = L U$ where $L$ is
+a lower-diagonal matrix of the form:
+\[ L = \left(
+\begin{matrix}
+\alpha_{11} & 0           & 0           & 0 \\
+\alpha_{21} & \alpha_{22} & 0           & 0 \\
+\alpha_{31} & \alpha_{32} & \alpha_{33} & 0 \\
+\alpha_{41} & \alpha_{42} & \alpha_{43} & \alpha_{44} \\
+\end{matrix} \right)
+\]
+%
+(where all diagonal values $\alpha){ii} = 1$) and $U$ is an upper-diagonal matrix of the form:
+\[ U = \left(
+\begin{matrix}
+\beta_{11} & \beta_{12} & \beta_{13} & \beta_{14} \\
+0          & \beta_{22} & \beta_{23} & \beta_{24} \\
+0          & 0          & \beta_{33} & \beta_{34} \\
+0          & 0          & 0          & \beta_{44} \\
+\end{matrix} \right)
+\]
+%
+We can find the values of $\alpha_{ij}$ and $\beta_{ij}$ by following
+this procedure.  First, $\alpha_{ii} = 1$.  For all values of $j = 1,
+2, \dots, N$, follow the next steps: For each value of $i = 1, 2,
+\dots, j$, solve for $\beta_{ij}$ using the relationship:
+\[ \beta_{ij} = a_{ij} - \sum_{k=1}^{i-1} \alpha_{ik}\beta_{kj} \]
+For the values of $i = j+1, j+2, \dots, N$, solve for the values of
+$\alpha_{ij}$ with the following:
+%
+\[ \alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) \]
+
+\paragraph{Calculate a matrix determinant}
+
+The determinant $D$ of a matrix $a_{ij}$ is calculated from the
+product of the diagonal elements of the LU decomposition: 
+\[ D = P_{i=1}{N} U_{ii} \]
+
+This shall be calculated using the GSL function
+\code{gsl_linalg_LU_det}.  If the matrix is large or the magnitude of
+the elements is large, the determinant may overflow the data type.  In
+these cases, the (natural) logarithm of the determinant may be
+calculated instead (as the sum of the logarithms of the diagonal
+elements).  In this case, the GSL function \code{gsl_linalg_LU_lndet}
+shall be used.
+
+\paragraph{Solving a Linear Equation}
+
+The LU decomposition of a matrix may be used to solve the
+matrix equation $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $
+
+Given the LU decomposition of a matrix into $\alpha_{ij}$ and
+$\beta_{ij}$, the technique of back-substitution is used to solve the
+equation above. First solve the equation $L y = B$:
+\[ y_1 = \frac{b_1}{\alpha_{11}}\]
+\[ y_i = \frac{1}{\alpha_{ii}}\left(b_i - \sum_{j=1}^{i-1} \alpha_{ij} y_i\right) \]
+% 
+The values of $y$ may be used to solve for $x_i$ using the
+relationship $U x = b$:
+\[ x_N = \frac{y_N}{\beta_{NN}}\]
+\[ x_i = \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right) \]
+
+\paragraph{Invert a matrix}
+
+Inversion of a matrix using the LU decomposition is performed by
+performing back-subsitution to solve a series of linear equations of
+the form $\sum_{j=1}^{N} A_{i,j} \times x_j = B_j $ where the values
+of $B_j$ represent the columns of the identity matrix.  The solution
+vectors $x_j$ represent the columns of the inverse matrix.  This
+operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}.
+
+\paragraph{Perform matrix addition, subtraction and multiplication}
+
+Matrix binary arithmetic operations differ from image binary
+arithmetic operations in a fundamental way:  For image operations, the
+resulting pixel value is determined by performing the operation on the
+two corresponding input operand pixels.  For matrix operations, the
+resulting element value results from the operation on the
+corresponding pair of row and colum from the input matrix.  So, for
+example, given the input matrices $\alpha_{ij}$ and $\beta_{ij}$, the
+image and matrix multiplication correspond to:
+%
+\[ \mbox{image}_{ij} = \alpha_{ij} \times \beta_{ij} \]
+\[ \mbox{matrix}_{jk} = \sum_{i=1}^N \alpha_{ij} \times \beta_{ki} \]
+%
+The matrix math function \code{psMatrixOp} shall implement the matrix
+version of the binary arithmetical operations for the following
+operators: $+, -, \times, /$
+
+\paragraph{Transpose a matrix}
+
+The transpose of a matrix is simply the reorganization of the matrix
+elements by 'flipping' the matrix along a diagonal.  For non-square
+matrices with dimensions $N \times M$, the resulting matrix has
+dimensions $M \times N$.  The element of the output matrix $T_{ij}$ is
+given by
+%
+\[ T_{ij} = M_{ji} \]
+where $M_{ji}$ is the matrix to be transposed.
+
+\paragraph{Convert a matrix to a vector}
+\TBD{write me} 
+
+\subsubsection{General Polynomial Fitting}
+
+\TBD{define for use with psGetArrayPolynomial}
+
+\TBD{change psGetArrayPolynomial to return a psPolynomial}
+
+Given a set of data values $y_i$ with errors $\sigma_i$, related to
+independent data values $x_i$, we would like to fit a polynomial model
+of $N_{par}$ free parameters of the form $y = \sum_{j=0}^{N_{par}}
+\alpha_j x^j$.  The model is determined by minimizing the value of
+$\chi^2$ as defined above (\ref{chisq}), and the minimization is
+determined by solving for the set of parameters $\alpha_j$ for which
+the partial derivatives of the $\chi^2$ with respect to those
+parameters are zero.  The partial derivatives are
+
+\[ \frac{\partial \chi^2}{\partial \alpha_k} = -2 \sum_i (y_i - \sum_j x_i^j \alpha_j) \frac{x_i^k}{\sigma_i^2} \]
+
+Setting the derivatives to zero, this may be reduced to the following
+matrix equation:
+
+\[ \sum_j \alpha_j \sum_i \frac{x_i^k x_i^j}{\sigma_i^2} = \sum_i \frac{x_i^k y_i}{\sigma_i^2} \]
+
+This matrix equation may be solved with LU Decomposition
+(section~\ref{LUdecomp}).
+
+\subsubsection{Non-linear Fitting: Levenberg-Marquardt Method}
+
+\TBD{describe LMM for psMinimize and psMinimizeChi2}
+
+For models in which the system of equations defined by the partial
+derivatives cannot be solved with the linear technique, other options
+are necessary.  The Levenberg-Marquardt Method (LMM; see Numerical
+Recipes for a discusion) may be used for these situations.  In LMM, we
+make a guess at the input parameters, measure the $\chi^2$, vary the
+parameters by a particular choice based on the gradient, measure the
+$\chi^2$ again, and adjust the parameters and the parameter varient
+based on the results.
+
+Given a set of $N$ data values $y_i$ with errors $\sigma_i$, dependent
+on values $x_i$, we would like to find the parameters $a_k$ of the
+function $f(x_i; a_k)$ which minimize the $\chi^2$, defined in the
+usual manner (\ref{chisq}).  We start with a set of parameter guesses,
+$a_k$.  We calculate the gradient $\beta_k$ and the Hessian matrix
+$\alpha_{j,k}$ at this parameter selection as follows:
+
+\[ \beta_k = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \] 
+\[ \alpha_{j,k} = \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2} \] 
+%
+We now define the new parameter guess for $a_k$ based on the gradient
+and Hessian by defining $A_{j,k}$ as a variant on $\alpha_{j,k}$ as
+follows:
+
+\[ A_{j,k} = \alpha_{j,k} ~ (j \ne k) \]
+\[ A_{j,k} = (1 + \lambda) \alpha_{j,k} ~ (j = k) \]
+%
+and solve the system of equations represented by:
+\[ A_{j,k} \alpha^\prime_k = \beta_j \]
+%
+where $\alpha^\prime_k$ represents our new attempt at a parameter
+guess. We use this parameter set to calculate $\chi^2$.  If the new
+value of $\chi^2$ is lower than the previous guess, we accept this new
+set of parameters and decrease $\lambda$ by a factor of 10, otherwise
+we keep the old set, and increase the value of $\lambda$ by a factor
+of 10.  We repeat this process until the value of the reduced $\chi^2$
+changes by much less than 1.0.  The resulting values of $a_k$ are the
+best-fit parameters for the system.  If the errors are normally
+distributed, the formal errors on the parameters are then calculated
+by setting $\lambda = 0$ and calculating the covarience matrix
+$C_{i,j}$, the inverse of the matrix $\alpha_{j,k}$.  The independent
+68.3\% confidence limit on parameter $a_k$ is then $\sqrt{C_{k,k}}$.
+Confidence contours for sets of parameters may be defined as well by
+the function $\Delta = \delta\bar{a} P_{j,k}^{-1} \delta\bar{a}$ where
+$P_{j,k}$ is the projected matrix of $C_{j,k}$, ie those rows and
+columns of $C_{j,k}$ associated with the parameters of interest, the
+vector $\delta\bar{a}$.  The value of $\Delta$ is given by the table
+below for specific confidence limits and numbers of parameters.  
+Note that it is necessary to be able to calculate both the function as
+well as its derivative for any combination of parameters and dependent
+variables.
+
+\begin{center}
+\begin{tabular}{|l|r|r|r|}
+\hline
+{\bf P} & \multicolumn{3}{c|}{\bf $N_{par}$} \\
+        & 1    & 2    & 3    \\
+\hline
+68.3\%  & 1.00 & 2.30 & 3.53 \\
+95.4\%  & 4.00 & 6.17 & 8.02 \\
+99.73\% & 9.00 & 11.8 & 14.2 \\
+\hline
+\end{tabular}
+\end{center}
+
+\end{document}