\HeaderA{chol}{The Choleski Decomposition}{chol}
\keyword{algebra}{chol}
\keyword{array}{chol}
\begin{Description}\relax
Compute the Choleski factorization of a real symmetric
positive-definite square matrix.
\end{Description}
\begin{Usage}
\begin{verbatim}
chol(x, pivot = FALSE,  LINPACK = pivot)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{x}] a real symmetric, positive-definite matrix
\item[\code{pivot}] Should pivoting be used?
\item[\code{LINPACK}] logical.  Should LINPACK be used in the non-pivoting
case (for compatibility with \R{} < 1.7.0)?
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
This is an interface to the LAPACK routine DPOTRF and the
LINPACK routines DPOFA and DCHDC.

Note that only the upper triangular part of \code{x} is used, so
that \eqn{R'R = x}{} when \code{x} is symmetric.

If \code{pivot = FALSE} and \code{x} is not non-negative definite an
error occurs.  If \code{x} is positive semi-definite (i.e., some zero
eigenvalues) an error will also occur, as a numerical tolerance is used.

If \code{pivot = TRUE}, then the Choleski decomposition of a positive
semi-definite \code{x} can be computed.  The rank of \code{x} is
returned as \code{attr(Q, "rank")}, subject to numerical errors.
The pivot is returned as \code{attr(Q, "pivot")}.  It is no longer
the case that \code{t(Q) \%*\% Q} equals \code{x}.  However, setting
\code{pivot <- attr(Q, "pivot")} and \code{oo <- order(pivot)}, it
is true that \code{t(Q[, oo]) \%*\% Q[, oo]} equals \code{x},
or, alternatively, \code{t(Q) \%*\% Q} equals \code{x[pivot,
  pivot]}.  See the examples.
\end{Details}
\begin{Value}
The upper triangular factor of the Choleski decomposition, i.e., the
matrix \eqn{R}{} such that \eqn{R'R = x}{} (see example).

If pivoting is used, then two additional attributes
\code{"pivot"} and \code{"rank"} are also returned.
\end{Value}
\begin{Section}{Warning}
The code does not check for symmetry.

If \code{pivot = TRUE} and \code{x} is not non-negative
definite then there will be a warning message but a meaningless
result will occur.  So only use \code{pivot = TRUE} when \code{x} is
non-negative definite by construction.
\end{Section}
\begin{References}\relax
Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988)
\emph{The New S Language}.
Wadsworth \& Brooks/Cole.

Dongarra, J. J., Bunch, J. R., Moler, C. B. and Stewart, G. W. (1978)
\emph{LINPACK Users Guide.}  Philadelphia: SIAM Publications.

Anderson. E. and ten others (1999)
\emph{LAPACK Users' Guide}. Third Edition. SIAM.\\
Available on-line at
\url{http://www.netlib.org/lapack/lug/lapack_lug.html}.
\end{References}
\begin{SeeAlso}\relax
\code{\LinkA{chol2inv}{chol2inv}} for its \emph{inverse} (without pivoting),
\code{\LinkA{backsolve}{backsolve}} for solving linear systems with upper
triangular left sides.

\code{\LinkA{qr}{qr}}, \code{\LinkA{svd}{svd}} for related matrix factorizations.
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
( m <- matrix(c(5,1,1,3),2,2) )
( cm <- chol(m) )
t(cm) %*% cm  #-- = 'm'
crossprod(cm)  #-- = 'm'

# now for something positive semi-definite
x <- matrix(c(1:5, (1:5)^2), 5, 2)
x <- cbind(x, x[, 1] + 3*x[, 2])
m <- crossprod(x)
qr(m)$rank # is 2, as it should be

# chol() may fail, depending on numerical rounding:
# chol() unlike qr() does not use a tolerance.
try(chol(m))

(Q <- chol(m, pivot = TRUE)) # NB wrong rank here ... see Warning section.
## we can use this by
pivot <- attr(Q, "pivot")
oo <- order(pivot)
t(Q[, oo]) %*% Q[, oo] # recover m

## now for a non-positive-definite matrix
( m <- matrix(c(5,-5,-5,3),2,2) )
try(chol(m))  # fails
try(chol(m, LINPACK=TRUE))   # fails
(Q <- chol(m, pivot = TRUE)) # warning
crossprod(Q)  # not equal to m
\end{ExampleCode}
\end{Examples}

