[R] about the Choleski factorization

[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>


> 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.