Skip to content

about the Choleski factorization

6 messages · 93354504, Duncan Murdoch, Ravi Varadhan +1 more

93354504
# 
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
Duncan Murdoch
#
On 3/27/2009 11:46 AM, 93354504 wrote:
> 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
Ravi Varadhan
# 
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>
Ravi Varadhan
# 
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>
Ravi Varadhan
# 
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>
Peter Dalgaard
#
Duncan Murdoch wrote:
Or just

 > r<-3:1
 > chol(x[r,r])[r,r]
           [,1]     [,2] [,3]
[1,] 1.9771421 0.000000    0
[2,] 0.3015113 1.658312    0
[3,] 1.0000000 0.500000    2

(It is after all, just a matter of starting from the other end).