\\% Singular value decomposition (SVD) and pseudo inverse of a matrix.
\\ Author: Joerg Arndt
\\ online at http://www.jjj.de/pari/
\\ version: 2005-October-30 (16:57)

matSVDcore(A)=
\\ Core routine for the singular value decomposition:
\\ Return [U, d, V] so that U*d*V~==A
\\ d is a diagonal matrix
\\ U, V are orthogonal
{
    local(U, d, V);  \\ returned quantities
    local(t, R, d1);

    R = conj(A~)*A;  \\ R==V*d^2*V~

    \\  qfjacobi(x): eigenvalues and orthogonal matrix of eigenvectors
    \\  of the real symmetric matrix x.  Vectors are the columns.
    t = qfjacobi( R );  \\ fails with eigenvalues==zero

    V = t[2];
    d = real(sqrt(t[1]));
    d1 = d;
    for (k=1, length(d1), t=d1[k]; if (abs(t)>1e-16, t=1/t, t=0); d1[k]=t );
    d1 = matdiagonal(d1);
    d = matdiagonal(d);
    U = (A*V*d1);
    return( [U, d, V] );
} /* ------- */


matSVD(A)=
\\ Singular value decomposition:
\\ Return [U, d, V] so that U*d*V~==A
\\ d is a diagonal matrix
\\ U, V are orthogonal
{
    local(tq, t, U, d, V);
    t = matsize(A);
    tq=0;  if ( t[1]<t[2],  tq=1; A=A~; );
    t = matSVDcore(A);
    d = t[2];
    if ( tq,
        U = t[3];  V=t[1];
    , /* else */
        U = t[1];  V=t[3];
    );
    return( [U, d, V] );
} /* ------- */


matpseudoinv(A)=
\\ Return pseudo-inverse of A
{
    local(t, x, U, d, V);
    t = matSVD(A);
    U = t[1];  d = t[2];  V = t[3];
    for (k=1, matsize(d)[1],
        x=d[k,k];  if (x>1e-15, d[k,k]=1/x, d[k,k]=0);
    );
    return( V*d*U~ );
} /* ------- */


\\ ==== end of file ====