Matrix decomposition



Donna and Dan Stanger  writes:

> Does anyone know a decomposition which will decompose a positive
> definite matrix Q into factors transpose(H)*H?

(c1) load("eigen");
Warning - you are redefining the Macsyma function eigenvalues
Warning - you are redefining the Macsyma function eigenvectors
(d1)            /usr/src/maxima-newlines/share/matrix/eigen.mac
(c2) innerproduct(x,y):=ratsimp(matrix(x) . Q . transpose(matrix(y)));
(d2) innerproduct(x, y) := ratsimp(matrix(x) . (Q . transpose(matrix(y))))
(c3) Q:matrix([1,1/2],[1/2,1]);
                                   [    1 ]
                                   [ 1  - ]
                                   [    2 ]
(d3)                               [      ]
                                   [ 1    ]
                                   [ -  1 ]
                                   [ 2    ]
(c4) gramschmidt([[1,0],[0,1]]);
                                          1
(d4)                          [[1, 0], [- -, 1]]
                                          2
(c5) A:apply('matrix,%);
                                  [  1   0 ]
                                  [        ]
(d5)                              [   1    ]
                                  [ - -  1 ]
                                  [   2    ]
(c6) ratsimp(A . Q . transpose(A));
                                   [ 1  0 ]
                                   [      ]
(d6)                               [    3 ]
                                   [ 0  - ]
                                   [    4 ]
(c7) D:%^(1/2);
                                [ 1     0    ]
                                [            ]
(d7)                            [    sqrt(3) ]
                                [ 0  ------- ]
                                [       2    ]
(c8) H:ratsimp(D . transpose(A)^^(-1));
                                [       1    ]
                                [ 1     -    ]
                                [       2    ]
(d8)                            [            ]
                                [    sqrt(3) ]
                                [ 0  ------- ]
                                [       2    ]
(c9) ratsimp(transpose(H) . H);
                                   [    1 ]
                                   [ 1  - ]
                                   [    2 ]
(d9)                               [      ]
                                   [ 1    ]
                                   [ -  1 ]
                                   [ 2    ]


Wolfgang