Very nice, Duncan.
Here is a little function called loch() that implements your idea for the
Lochesky factorization:
loch <- function(mat) {
n <- ncol(mat)
rev <- diag(1, n)[, n: 1]
rev %*% chol(rev %*% mat %*% rev) %*% rev
}
x=matrix(c(5,1,2,1,3,1,2,1,4),3,3)
L <- loch(x)
all.equal(x, t(L) %*% L)
A <- matrix(rnorm(36), 6, 6)
A <- A %*% t(A)
L <- loch(x)
all.equal(x, t(L) %*% L)
Ravi.
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
----- Original Message -----
From: 93354504 <93354...@nccu.edu.tw>
Date: Friday, March 27, 2009 11:58 am
Subject: [R] about the Choleski factorization
To: r-help <r-help@r-project.org>
> Hi there,
>
> Given a positive definite symmetric matrix, I can use chol(x) to
> obtain U where U is upper triangular
> and x=U'U. For example,
>
> x=matrix(c(5,1,2,1,3,1,2,1,4),3,3)
> U=chol(x)
> U
> # [,1] [,2] [,3]
> #[1,] 2.236068 0.4472136 0.8944272
> #[2,] 0.000000 1.6733201 0.3585686
> #[3,] 0.000000 0.0000000 1.7525492
> t(U)%*%U # this is exactly x
>
> Does anyone know how to obtain L such that L is lower triangular and
> x=L'L? Thank you.
>
> Alex
>
> ______________________________________________
> R-help@r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
