svd.cpp

Go to the documentation of this file.

//===========================================================================
//===========================================================================

#include <cmath>
#include <cstdlib>
#include <SharkDefs.h>
#include <LinAlg/LinAlg.h>

static double SIGN(double a, double b)
{
    return b > 0 ? fabs(a) : -fabs(a);
}

//===========================================================================
void svd
(
    const Array2D< double >& amatA,
    Array2D< double >& umatA,
    Array2D< double >& vmatA,
    Array  < double >& w, // singular values
    unsigned maxIterations,
    bool ignoreThreshold
)
{
    unsigned m = amatA.rows(); /* rows */
    unsigned n = amatA.cols(); /* cols */
//  double** const a = amatA.ptrArr(); /* input matrix */
//  double **u = umatA.ptrArr(); /* left vectors */
//  double **v = vmatA.ptrArr(); /* right vectors */

    int flag;
    unsigned i, its, j, jj, k, l, nm(0);
    double anorm, c, f, g, h, s, scale, x, y, z;

    Array< double > rv1(n);

    /* copy A to U */
    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++) {
            umatA(i, j) = amatA(i, j);
        }
    }

    /* householder reduction to bidiagonal form */
    g = scale = anorm = 0.0;

    for (i = 0; i < n; i++) {
        l = i + 1;
        rv1(i) = scale * g;
        g = s = scale = 0.0;

        if (i < m) {
            for (k = i; k < m; k++) {
                scale += fabs(umatA(k, i));
            }

            if (scale != 0.0) {
                for (k = i; k < m; k++) {
                    umatA(k, i) /= scale;
                    s += umatA(k, i) * umatA(k, i);
                }

                f = umatA(i, i);
                g = -SIGN(sqrt(s), f);
                h = f * g - s;
                umatA(i, i) = f - g;

                for (j = l; j < n; j++) {
                    s = 0.0;
                    for (k = i; k < m; k++) {
                        s += umatA(k, i) * umatA(k, j);
                    }

                    f = s / h;
                    for (k = i; k < m; k++) {
                        umatA(k, j) += f * umatA(k, i);
                    }
                }

                for (k = i; k < m; k++) {
                    umatA(k, i) *= scale;
                }
            }
        }

        w(i) = scale * g;
        g = s = scale = 0.0;

        if (i < m && i != n - 1) {
            for (k = l; k < n; k++) {
                scale += fabs(umatA(i, k));
            }

            if (scale != 0.0) {
                for (k = l; k < n; k++) {
                    umatA(i, k) /= scale;
                    s += umatA(i, k) * umatA(i, k);
                }

                f = umatA(i, l);
                g = -SIGN(sqrt(s), f);
                h = f * g - s;
                umatA(i, l) = f - g;

                for (k = l; k < n; k++) {
                    rv1(k) = umatA(i, k) / h;
                }

                for (j = l; j < m; j++) {
                    s = 0.0;
                    for (k = l; k < n; k++) {
                        s += umatA(j, k) * umatA(i, k);
                    }

                    for (k = l; k < n; k++) {
                        umatA(j, k) += s * rv1(k);
                    }
                }

                for (k = l; k < n; k++) {
                    umatA(i, k) *= scale;
                }
            }
        }

        anorm = Shark::max(anorm, fabs(w(i)) + fabs(rv1(i)));
    }

    /* accumulation of right-hand transformations */
    for (l = i = n; i--; l--) {
        if (l < n) {
            if (g != 0.0) {
                for (j = l; j < n; j++) {
                    /* double division avoids possible underflow */
                    vmatA(j, i) = (umatA(i, j) / umatA(i, l)) / g;
                }

                for (j = l; j < n; j++) {
                    s = 0.0;
                    for (k = l; k < n; k++) {
                        s += umatA(i, k) * vmatA(k, j);
                    }

                    for (k = l; k < n; k++) {
                        vmatA(k, j) += s * vmatA(k, i);
                    }
                }
            }

            for (j = l; j < n; j++) {
                vmatA(i, j) = vmatA(j, i) = 0.0;
            }
        }

        vmatA(i, i) = 1.0;
        g = rv1(i);
    }

    /* accumulation of left-hand transformations */
    for (l = i = Shark::min(m, n); i--; l--) {
        g = w(i);

        for (j = l; j < n; j++) {
            umatA(i, j) = 0.0;
        }

        if (g != 0.0) {
            g = 1.0 / g;

            for (j = l; j < n; j++) {
                s = 0.0;
                for (k = l; k < m; k++) {
                    s += umatA(k, i) * umatA(k, j);
                }

                /* double division avoids possible underflow */
                f = (s / umatA(i, i)) * g;

                for (k = i; k < m; k++) {
                    umatA(k, j) += f * umatA(k, i);
                }
            }

            for (j = i; j < m; j++) {
                umatA(j, i) *= g;
            }
        }
        else {
            for (j = i; j < m; j++) {
                umatA(j, i) = 0.0;
            }
        }

        umatA(i, i)++;
    }

    /* diagonalization of the bidiagonal form */
    for (k = n; k--;) {
        for (its = 1; its <= maxIterations; its++) {
            flag = 1;

            /* test for splitting */
//          for (l = k + 1; l--; )
            /* Thomas Buecher:
                    nderung, die fehlerhaft sein kann, aber Absturz verhindert:
                l kann 0 werden --> nm (zwei Zeilen tiefer) wird riesig, da 0-1 auf
                unsigned typ berechnet wird
            */
            for (l = k; l > 0; l--) {
                /* rv1(0) is always zero, so there is no exit */
                nm = l - 1;

                if (fabs(rv1(l)) + anorm == anorm) {
                    flag = 0;
                    break;
                }

                if (fabs(w(nm)) + anorm == anorm) {
                    break;
                }
            }

            if (flag) {
                /* cancellation of rv1(l) if l greater than 0 */
                c = 0.0;
                s = 1.0;

                for (i = l; i <= k; i++) {
                    f = s * rv1(i);
                    rv1(i) *= c;

                    if (fabs(f) + anorm == anorm) {
                        break;
                    }

                    g = w(i);
                    h = hypot(f, g);
                    w(i) = h;
                    h = 1.0 / h;
                    c = g * h;
                    s = -f * h;

                    for (j = 0; j < m; j++) {
                        y = umatA(j, nm);
                        z = umatA(j, i);
                        umatA(j, nm) = y * c + z * s;
                        umatA(j, i ) = z * c - y * s;
                    }
                }
            }

            /* test for convergence */
            z = w(k);

            if (l == k) {
                if (z < 0.0) {
                    w(k) = -z;
                    for (j = 0; j < n; j++) {
                        vmatA(j, k) = -vmatA(j, k);
                    }
                }
                break;
            }

            if( its == maxIterations ) { 
                if(ignoreThreshold)
                    break ;
                else
                    throw SHARKEXCEPTION("too many iterations for diagonalization of the bidiagonal form during SVD");
            }


            /* shift from bottom 2 by 2 minor */
            x = w(l);
            nm = k - 1;
            y = w(nm);
            g = rv1(nm);
            h = rv1(k);
            f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
            g = hypot(f, 1.0);
            f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x;

            /* next qr transformation */
            c = s = 1.0;

            for (j = l; j < k; j++) {
                i = j + 1;
                g = rv1(i);
                y = w(i);
                h = s * g;
                g *= c;
                z = hypot(f, h);
                rv1(j) = z;
                c = f / z;
                s = h / z;
                f = x * c + g * s;
                g = g * c - x * s;
                h = y * s;
                y *= c;

                for (jj = 0; jj < n; jj++) {
                    x = vmatA(jj, j);
                    z = vmatA(jj, i);
                    vmatA(jj, j) = x * c + z * s;
                    vmatA(jj, i) = z * c - x * s;
                }

                z = hypot(f, h);
                w(j) = z;

                /* rotation can be arbitrary if z is zero */
                if (z != 0.0) {
                    z = 1.0 / z;
                    c = f * z;
                    s = h * z;
                }

                f = c * g + s * y;
                x = c * y - s * g;

                for (jj = 0; jj < m; jj++) {
                    y = umatA(jj, j);
                    z = umatA(jj, i);
                    umatA(jj, j) = y * c + z * s;
                    umatA(jj, i) = z * c - y * s;
                }
            }

            rv1(l) = 0.0;
            rv1(k) = f;
            w(k) = x;
        }
    }
}


//===========================================================================
// void svd
// (
//  const Array< double >& A,
//  Array< double >& U,
//  Array< double >& V,
//  Array< double >& W
// )
// {
//  SIZE_CHECK(A.ndim() == 2 && A.dim(0) == A.dim(1))
// 
//  Array2DReference< double > amatL(const_cast< Array< double >& >(A));
//  Array2DReference< double > umatL(U);
//  Array2DReference< double > vmatL(V);
// 
//  U.resize(A.dim(0), A.dim(1), false);
//  V.resize(A.dim(1), A.dim(1), false);
//  W.resize(A.dim(1), false);
// 
//  svd(amatL, umatL, vmatL, W);
//  svdsort(umatL, vmatL, W);
// }