On 2012-09-13, Barton Willis <willisb at unk.edu> wrote:
> Th Cholesky code doesn't work for symbolic cases, I think.
Hmm ... there used to be a cholesky.mac in share which worked for
symbolic cases, but it got nuked in favor of the one in the
linearalgebra package. For the record here it is.
load ("cholesky.mac");
M : matrix ([1, p12, p13], [p12, 1, p23], [p13, p23, 1]);
M1 : cholesky (M1);
=> matrix([1,0,0],[p12,sqrt(1-p12^2),0],
[p13,(p23-p12*p13)/sqrt(1-p12^2),
sqrt(-(p23-p12*p13)^2/(1-p12^2)-p13^2+1)])
M1 . transpose (M1) - M;
=> matrix([0,0,0],[0,0,0],[0,0,0])
Probably there is a way to get the one version to punt to the other
(from symbolic to numerical, presumably -- that's from the general case
to the special).
All the best,
Robert Dodier
PS. Here's cholesky.mac.
/* Compute Cholesky decomposition of A,
* a lower-triangular matrix L such that L . transpose(L) = A
*
* Examples:
A : matrix ([a, b, c], [d, e, f], [g, h, i]);
A2 : transpose (A) . A;
B : cholesky (A2);
B . transpose (B) - A2;
A : matrix ([2, 3, 4], [-2, 2,- 3], [11, -2, 3]);
A2 : transpose (A) . A;
B : cholesky (A2);
B . transpose (B) - A2;
*
* copyright Robert Dodier, 2005/11/01
* Released under the terms of the GNU Public License
*
*/
cholesky (A):= block
([n : length (A),
L : copymatrix (A),
p : makelist (0, i, 1, length (A))],
for i thru n do
for j : i thru n do
(x : L[i, j],
x : x - sum (L[j, k] * L[i, k], k, 1, i - 1),
if i = j then
p[i] : 1 / sqrt(x)
else
L[j, i] : x * p[i]),
for i thru n do
L[i, i] : 1 / p[i],
for i thru n do
for j : i + 1 thru n do
L[i, j] : 0,
L);