Dan, here is a function to compute the Cholesky decomposition.
Maybe this is what you want. The result is a lower-triangular
matrix L s.t. L.transpose(L) equals the original matrix.
I wrote this function some time ago so it's possible that
I didn't know what I was doing. However, in some simple
tests it did what it's supposed to.
Here's a test case:
X: matrix ([a, b, c], [d, e, f], [g, h, i])$
XX: transpose (X) . X$ /* ensure XX is p.d. */
L: cholesky (XX)$
L . transpose (L) - XX; /* should yield all 0 */
hth,
Robert Dodier
cholesky(A):=
block(
n:length(A),
L:copymatrix(A),
array(p,n,1),
for i thru n do
block(print("i is",i),
for j:i thru n do
block(print("j is",j),
x:L[i,j],
x: x - sum( L[j,k]*L[i,k], k, 1, i-1 ),
if i=j then
block(print("i=j and x is",x),
p[i,1]:1/sqrt(x))
else
block(print("i!=j and x is",x,"and p[",i,",1] is",p[i,1]),
L[j,i]:x*p[i,1]))),
for i thru n do print("p[",i,",1] is",p[i,1]),
for i thru n do
L[i,i]:1/p[i,1],
for i thru n do
for j:i+1 thru n do
L[i,j]:0,
return(L))$
__________________________________
Do you Yahoo!?
Yahoo! Mail - Helps protect you from nasty viruses.
http://promotions.yahoo.com/new_mail