[Rd] Pivoting in chol

Jonathan Rougier J.C.Rougier@durham.ac.uk
Wed, 20 Feb 2002 10:34:27 +0000

This is a multi-part message in MIME format.
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Hi Everyone,

I have modified my version of R-1.4.1 to include choleski with pivoting
(like in Splus).  I thought R-core might consider including this in the
next version of R, so I give below the steps required to facilitate

1. Copied Linpack routine "dchdc.f" into src/appl

2. Inserted line F77_SUBROUTINE(dchdc) in src/appl/ROUTINES

3. Inserted "dchdc.f" into SOURCES_F in appl/Makefile.in

4. Modifed src/library/base/R/chol.R (attached) and 
src/library/base/man/chol.Rd (also attached)

5. Tidied up.  Removed "chol.f", "dpoco.f" and "dposl.f" from
src/appl/Makefile.in, src/appl/ROUTINES, src/include/R_ext/Applic.h,
src/include/R_ext/Linpack.h, src/unix/FFDecl.h and src/unix/FFTab.h.

6. Note: Cannot remove "dpofa.f" because it is used in "lbfgsb.c".

I was following my nose a bit during this process.  Everything seems to
be OK but if I have done something drastic I would appreciate it if
someone could tell me!

Cheers, Jonathan.
Jonathan Rougier                       Science Laboratories
Department of Mathematical Sciences    South Road
University of Durham                   Durham DH1 3LE
tel: +44 (0)191 374 2361, fax: +44 (0)191 374 7388
Content-Type: text/plain; charset=us-ascii;
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;

# modifed J.C.Rougier <J.C.Rougier@durham.ac.uk>, 19.02.02
# to allow cholesky decomposition with pivoting.
#chol <- function(x)
#    if(!is.numeric(x))
#	stop("non-numeric argument to chol")
#    if(is.matrix(x)) {
#	if(nrow(x) != ncol(x))
#	    stop("non-square matrix in chol")
#	n <- nrow(x)
#    }
#    else {
#	if(length(x) != 1)
#	    stop("non-matrix argument to chol")
#	n <- as.integer(1)
#    }
#    if(!is.double(x)) storage.mode(x) <- "double"
#    z <- .Fortran("chol",
#		  x=x,
#		  n,
#		  n,
#		  v=matrix(0, nr=n, nc=n),
#		  info=integer(1),
#		  DUP=FALSE, PACKAGE="base")
#    if(z$info)
#	stop("non-positive definite matrix in chol")
#    z$v

"chol" <- function(x, pivot = FALSE)
  if (is.complex(x))
    stop("complex matrices not permitted at present")
  else if (!is.numeric(x))
    stop("non-numeric argument")
  x <- as.matrix(x)
  if (length(p <- unique(dim(x))) > 1)
    stop("non-square matrix")

  x[lower.tri(x)] <- 0
  if (!is.double(x))
    storage.mode(x) <- "double"

  z <- .Fortran("dchdc",
    x = x,
    piv = as.integer(rep(0, p)),
    rank = integer(1),
    package = "BASE"
  )[c(1, 5, 7)]

  if (!pivot && z$rank<p)
    stop("matrix not positive definite")

  robj <- z$x
  if (pivot) {
    attr(robj, "pivot") <- z$piv
    attr(robj, "rank") <- z$rank

chol2inv <- function(x, size=ncol(x))
	stop("non-numeric argument to chol2inv")
    if (!is.null(attr(x, "pivot")))
        stop("not implemented for pivoted Choleski decomposition") # JCR added
    if(is.matrix(x)) {
	nr <- nrow(x)
	nc <- ncol(x)
    else {
	nr <- length(x)
	nc <- as.integer(1)
    size <- as.integer(size)
    if(size <= 0 || size > nr || size > nc)
	stop("invalid size argument in chol2inv")
    if(!is.double(x)) storage.mode(x) <- "double"
    z <- .Fortran("ch2inv",
		  v=matrix(0, nr=size, nc=size),
		  DUP=FALSE, PACKAGE="base")
	stop("singular matrix in chol2inv")

Content-Type: text/plain; charset=us-ascii;
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;



\title{The Choleski Decomposition}

  Compute the Choleski factorization of a symmetric real non-negative
  definite square matrix, i.e. the matrix \code{Q} for which
  \code{t(Q) \%*\% Q} equals \code{x}.  }

chol(x, pivot = FALSE)

  \item{x}{a symmetric, positive definite matrix}
  \item{pivot}{Should pivoting be used?}
  The Choleski decomposition of \code{x}, as an upper-triangular
  matrix.  If pivoting is used, then two additional attributes
  \code{pivot} and \code{rank} are also returned (see Details for more

  Note that only the upper triangular part of \code{x} is used, so
  that \code{t(Q) \%*\% Q} only equals \code{x} when \code{x} is

  If \code{pivot=FALSE} and \code{x} is not non-negative definite an
  error occurs.  If \code{x} is positive semi-definite (ie some zero
  eigenvalues) an error may or may not occur, depending on
  finite precision numerical errors.

  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 now 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 for more information.

  The code does not check for symmetry or for non-negative
  definiteness.  If \code{pivot=TRUE} and \code{x} is not non-negative
  definite then there will be no error message but a meaningless
  result will occur.  So only use \code{pivot=TRUE} when \code{x} is
  non-negative definite by construction.

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

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

x <- matrix(c(1:5, (1:5)^2), 5, 2)
m <- crossprod(x)
Q <- chol(m)
all.equal.numeric(t(Q) \%*\% Q, m) # TRUE

Q <- chol(m, pivot=TRUE)
pivot <- attr(Q, "pivot")
oo <- order(pivot)
all.equal.numeric(t(Q[, oo]) \%*\% Q[, oo], m) # TRUE
all.equal.numeric(t(Q) \%*\% Q, m[pivot, pivot]) # TRUE

# now for something postitive semi definite

x <- cbind(x, x[, 1]+3*x[, 2])
m <- crossprod(x)
qr(m)$rank # is 2, as it should be

test <- try(Q <- chol(m))
inherits(test, "try-error") # TRUE

Q <- chol(m, pivot=TRUE) # nb wrong rank here ... you have been warned!
pivot <- attr(Q, "pivot")
oo <- order(pivot)
all.equal.numeric(t(Q[, oo]) \%*\% Q[, oo], m) # TRUE
all.equal.numeric(t(Q) \%*\% Q, m[pivot, pivot]) # TRUE

# and something meaningless

m <- matrix(1:16, 4, 4)
m <- m + t(m) # not nnd
Q <- chol(m, pivot=TRUE)
pivot <- attr(Q, "pivot")
all.equal(t(Q) \%*\% Q, m[pivot, pivot]) # FALSE!



r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch