about the Choleski factorization
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: rvaradhan at jhmi.edu
----- Original Message -----
From: Duncan Murdoch <murdoch at stats.uwo.ca>
Date: Friday, March 27, 2009 1:29 pm
Subject: Re: [R] about the Choleski factorization
To: 93354504 at nccu.edu.tw
Cc: r-help <r-help at r-project.org>
On 3/27/2009 11:46 AM, 93354504 wrote:
> 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 at r-project.org mailing list > > PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
> rev <- matrix(c(0,0,1,0,1,0,1,0,0),3,3) > rev
[,1] [,2] [,3] [1,] 0 0 1 [2,] 0 1 0 [3,] 1 0 0 (the matrix that reverses the row and column order when you pre and post multiply it). Then L <- rev %*% chol(rev %*% x %*% rev) %*% rev is what you want, i.e. you reverse the row and column order of the Choleski square root of the reversed x:
> x
[,1] [,2] [,3] [1,] 5 1 2 [2,] 1 3 1 [3,] 2 1 4
> L <- rev %*% chol(rev %*% x %*% rev) %*% rev > L
[,1] [,2] [,3] [1,] 1.9771421 0.000000 0 [2,] 0.3015113 1.658312 0 [3,] 1.0000000 0.500000 2
> t(L) %*% L
[,1] [,2] [,3] [1,] 5 1 2 [2,] 1 3 1 [3,] 2 1 4 Duncan Murdoch
______________________________________________ R-help at r-project.org mailing list PLEASE do read the posting guide and provide commented, minimal, self-contained, reproducible code.