[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.
--------------06648C7828389FFEA9FD7F0C
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
this.
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
http://www.maths.dur.ac.uk/stats/people/jcr/jcr.html
--------------06648C7828389FFEA9FD7F0C
Content-Type: text/plain; charset=us-ascii;
name="chol.R"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="chol.R"
# 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,
as.integer(p),
as.integer(p),
double(p),
piv = as.integer(rep(0, p)),
as.integer(pivot),
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
}
robj
}
chol2inv <- function(x, size=ncol(x))
{
if(!is.numeric(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",
x=x,
nr,
size,
v=matrix(0, nr=size, nc=size),
info=integer(1),
DUP=FALSE, PACKAGE="base")
if(z$info)
stop("singular matrix in chol2inv")
z$v
}
--------------06648C7828389FFEA9FD7F0C
Content-Type: text/plain; charset=us-ascii;
name="chol.Rd"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="chol.Rd"
\name{chol}
\alias{chol}
\title{The Choleski Decomposition}
\description{
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}. }
\usage{
chol(x, pivot = FALSE)
}
\arguments{
\item{x}{a symmetric, positive definite matrix}
\item{pivot}{Should pivoting be used?}
}
\value{
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
information).
}
\details{
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
symmetric.
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.
}
\section{Warning}{
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.
}
\references{
Dongarra, J. J., Bunch, J. R., Moler, C. B. and Stewart, G. W. (1978)
\emph{LINPACK Users Guide.} Philadelphia: SIAM Publications.
}
\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{
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!
}
\keyword{algebra}
\keyword{array}
--------------06648C7828389FFEA9FD7F0C--
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._