\HeaderA{chol}{The Choleski Decomposition}{chol}
\methaliasA{chol.default}{chol}{chol.default}
\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, ...)

## Default S3 method:
chol(x, pivot = FALSE,  LINPACK = pivot, ...)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{x}] an object for which a method exists.  The default method
applies to real symmetric, positive-definite matrices.
\item[\code{...}] arguments to be based to or from methods.
\item[\code{pivot}] Should pivoting be used?
\item[\code{LINPACK}] logical.  Should LINPACK be used in the non-pivoting
case (for compatibility with \R{} \textless{} 1.7.0)?
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
\code{chol} is generic: the description here applies to the default
method.

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}

