[R] about the Choleski factorization
Ravi Varadhan
rvaradhan at jhmi.edu
Fri Mar 27 18:34:46 CET 2009
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.
More information about the R-help
mailing list