[R] Cholesky Decomposition in R

[Manli Yan]

Hi everyone:
I try to use r to do the Cholesky Decomposition,which is A=LDL',so far I
only found how to decomposite A in to LL' by using chol(A),the function
Cholesky(A) doesnt work,any one know other command to decomposte A in to
LDL'

My r code is:
library(Matrix)
A=matrix(c(1,1,1,1,5,5,1,5,14),nrow=3)

> chol(A)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 0 2 2
[3,] 0 0 3

> Cholesky(A)
Error in function (classes, fdef, mtable) :
unable to find an inherited method for function "Cholesky", for signature
"matrix"

whatz wrong???
thanks~

[[alternative HTML version deleted]]

______________________________________________
R-h...@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.

[Ravi Varadhan]

You can use singular value decomposition:

?svd

> svd(A)
$d
[1] 16.3405917 2.8996176 0.7597907 # This is your diagonal matrix "D"

$u
[,1] [,2] [,3]
[1,] 0.08585595 -0.2420411 0.96645997
[2,] 0.40826313 -0.8763116 -0.25573252
[3,] 0.90881790 0.4165261 0.02357989

$v
[,1] [,2] [,3]
[1,] 0.08585595 -0.2420411 0.96645997
[2,] 0.40826313 -0.8763116 -0.25573252
[3,] 0.90881790 0.4165261 0.02357989

>

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvar...@jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

----------------------------------------------------------------------------
--------

[Albyn Jones]

try Cholesky() in package Matrix.

albyn

[Douglas Bates]

On Tue, Mar 10, 2009 at 4:33 PM, Manli Yan <manliy...@gmail.com> wrote:
>  Hi everyone:
>  I try to use r to do the Cholesky Decomposition,which is A=LDL',so far I
> only found how to decomposite A in to  LL' by using chol(A),the function
> Cholesky(A) doesnt work,any one know other command to decomposte A in to
> LDL'
>
>  My r code is:
> library(Matrix)
> A=matrix(c(1,1,1,1,5,5,1,5,14),nrow=3)
>
>> chol(A)
>     [,1] [,2] [,3]
> [1,]    1    1    1
> [2,]    0    2    2
> [3,]    0    0    3
>
>> Cholesky(A)
> Error in function (classes, fdef, mtable)  :
>  unable to find an inherited method for function "Cholesky", for signature
> "matrix"
>
> whatz wrong???

The answer, surprisingly, is in the documentation, accessible as

?Cholesky

which says that the first argument has to be a sparse, symmetric
matrix. Because the Cholesky function is intended for sparse matrices
it is not the best approach. The object returned is rather obscure

> Cholesky(as(A, "dsCMatrix"), LDL = TRUE)
'MatrixFactorization' of Formal class 'dCHMsimpl' [package "Matrix"]
with 10 slots
..@ x : num [1:6] 1 1 1 4 1 9
..@ p : int [1:4] 0 3 5 6
..@ i : int [1:6] 0 1 2 1 2 2
..@ nz : int [1:3] 3 2 1
..@ nxt : int [1:5] 1 2 3 -1 0
..@ prv : int [1:5] 4 0 1 2 -1
..@ colcount: int [1:3] 3 2 1
..@ perm : int [1:3] 0 1 2
..@ type : int [1:4] 2 0 0 1
..@ Dim : int [1:2] 3 3

It turns out that the factorization you want is encoded in the 'x'
slot but not in an obvious way. Even if you ask for an expansion

> expand(Cholesky(as(A, "dsCMatrix"), LDL = TRUE))
$P

[1,] | . .
[2,] . | .
[3,] . . |

$L

[1,] 1 . .
[2,] 1 2 .
[3,] 1 2 3

the result is converted from the LDL' factor to the LL' factor.

A better approach is to consider how the LDL' factorization is related
to the R'R form of the factorization returned by chol()

> ch <- chol(A)
> dd <- diag(ch)
> L <- t(ch/dd)
> DD <- dd^2
> L
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 1 1 0
[3,] 1 1 1
> DD
[1] 1 4 9

This is all rather backwards in that the whole purpose of the LDL'
form of the factorization is to avoid taking square roots to get the
diagonal elements and to contend with positive semidefinite matrices.
In other words, the LDL' form avoids some of the possible problems of
the LL' form but not if you go through the LL' form to get to it.

I think the underlying reason that an LDL' form is not directly
available in R is because there is no Lapack subroutine for it.

Let me know what our grade on the homework is.

[megh]

A=matrix(c(1,1,1,1,5,5,1,5,14),nrow=3)
t(eigen(A)$vector) %*% A %*% eigen(A)$vector
mat1 = t(eigen(A)$vector)
mat2 = diag(eigen(A)$values) # this is your diagonal matrix

--
View this message in context: http://www.nabble.com/Cholesky-Decomposition-in-R-tp22444083p22448498.html
Sent from the R help mailing list archive at Nabble.com.

[Ravi Varadhan]

You got an "A+" on the homework, Doug!

I got a "C-" for suggesting svd(), which clearly doesn't yield a lower (or
upper) triangular factorization.

Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvar...@jhmi.edu

Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

----------------------------------------------------------------------------
--------

-----Original Message-----
From: r-help-...@r-project.org [mailto:r-help-...@r-project.org] On

[Simon Wood]

A <- matrix(runif(25),5,5); A <- A%*%t(A) ## Example +ve def matrix

U <- chol(A) ## start from factor given by chol
D <- diag(U) ## extract sqrt(D)
L <- t(U/D) ## get unit lower triangular factor
D <- diag(D^2) ## and diagonal
## So now A = LDL'
range(A-L%*%D%*%t(L))

#best,
#Simon

--
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603 www.maths.bath.ac.uk/~sw283

[Manli Yan]

Thank you guys so much~~~~
Regards~

2009/3/11 Simon Wood <s.w...@bath.ac.uk>

> > http://www.R-project.org/posting-guide.html<http://www.r-project.org/posting-guide.html>and provide commented, minimal,