The R Project SVN R

% File src/library/base/man/chol.Rd
% Part of the R package, https://www.R-project.org
% Copyright 1995-2025 R Core Team
% Distributed under GPL 2 or later

\name{chol}
\alias{chol}
\alias{chol.default}
\title{The Cholesky Decomposition}
\description{
  Compute the Cholesky factorization of a real symmetric
  positive-definite square matrix.
}
\usage{
chol(x, \dots)

\method{chol}{default}(x, pivot = FALSE,  LINPACK = FALSE, tol = -1, \dots)
}
\arguments{
  \item{x}{an object for which a method exists.  The default method
    applies to numeric (or logical) symmetric, positive-definite matrices.}
  \item{\dots}{arguments to be passed to or from methods.}
  \item{pivot}{logical: should pivoting be used?}
  \item{LINPACK}{logical.  Defunct and gives an error.}
  \item{tol}{a numeric tolerance for use with \code{pivot = TRUE}.}
}
\value{
  The upper triangular factor of the Cholesky 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.
}
\details{
  \code{chol} is generic: the description here applies to the default
  method.

  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 Cholesky decomposition of a (rearranged)
  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]}. Notice also, that \code{Q[,oo]} is typically not triangular.
   See the examples.

  NOTE: For versions of R up to 4.5.x, the above wasn't quite true, 
  since LAPACK does not zero out the rows of \code{Q} beyond the rank.
  This is now done in R, but if you want code to be compatible with 
  earlier versions of R, you may want to ensure that you only use 
  the first \code{rank} rows.
  
  The value of \code{tol} is passed to LAPACK, with negative values
  selecting the default tolerance of (usually) \code{nrow(x) *
  .Machine$double.neg.eps * max(diag(x))}.  The algorithm terminates once
  the pivot is less than \code{tol}.

  Unsuccessful results from the underlying LAPACK code will result in an
  error giving a positive error code: these can only be interpreted by
  detailed study of the FORTRAN code.
}

\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.
}

\source{
  This is an interface to the LAPACK routines \code{DPOTRF} and
  \code{DPSTRF},

  LAPACK is from \url{https://netlib.org/lapack/} and its guide is listed
  in the references.
}
\references{
  \bibshow{R:Anderson+Bai+Bischof:1999,
    R:Becker+Chambers+Wilks:1988}
}

\seealso{
  \code{\link{chol2inv}} for its \emph{inverse} (without pivoting),
  \code{\link{backsolve}} for solving linear systems with upper
  triangular left sides.

  \code{\link{qr}}, \code{\link{svd}} for related matrix factorizations.
}

\examples{
( 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], 2*x[, 1] + x[, 2])
colnames(x) <- letters[20:23]
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))
## we can use this by
pivot <- attr(Q, "pivot")
crossprod(Q[, order(pivot)]) # recover m

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

## another example, this time with positive, negative, and zero eigenvalues
## notice that this (and the previous example) gets the rank wrong - it
## returns the rank of an approximating positive semidefinite matrix,
## effectively discounting negative eigenvalues.
(m <- matrix(c(1,1,2, 1,1,2, 2,2,1), 3, 3))
(Q <- chol(m, pivot=TRUE))
eigen(m)$values

## conjecture: indefinite cases should be detectable by comparing the 
## diagonal of the reconstructed matrix to the original
piv <- attr(Q, "pivot")
(D <- diag(crossprod(Q))) ## faster: colSums(Q^2)
diag(m)[piv] - D

}

\keyword{algebra}
\keyword{array}