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
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 > 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.> 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
You want a factorizzation of the form: A = L' L. Am I right (we may name this a "Lochesky factorization)? By convention, Cholesky factorization is of the form A = L L', where L is a lower triangular matrix, and L', its transpose, is upper traingular. So, all numerical routines compute L according to this definition. R gives you U = L', which is obviously upper triangular. If you want to use a different definition: A = L' L, that is fine mathematically. Although there is no easy way to transform the result of existing routines to get what you want, you can actually derive an algorithm to convert the standard factorization to the form you want. Rather than go to this trouble, you might as well just code it up from scratch. The big question of course is why do you want the Lochesky factorization? It doesn't do anything special that the traditional Cholesky factorization can do for a symmetric, positive-definite matrix (mainly, solve a system of equations). 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: 93354504 <93354504 at nccu.edu.tw> Date: Friday, March 27, 2009 11:58 am Subject: [R] about the Choleski factorization To: r-help <r-help at 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 at r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.
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: 93354504 <93354504 at nccu.edu.tw> Date: Friday, March 27, 2009 11:58 am Subject: [R] about the Choleski factorization To: r-help <r-help at 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 at r-project.org mailing list > > PLEASE do read the posting guide > and provide commented, minimal, self-contained, reproducible code.