OptiVec: MF_SVdecompose

MF_SVdecompose

MD_SVdecompose

ME_SVdecompose

Function

Singular Value Decomposition

Syntax C/C++	#include <MFstd.h> int MF_SVdecompose( fMatrix U, fMatrix V, fVector W, fMatrix MA, ui htA, ui lenA );
C++ MatObj	#include <OptiVec.h> int matrix<T>::SVdecompose( matrix<T> MV, vector<T> W, const matrix<T>& MA ); int matrix<T>::SVdecompose( matrix<T>* MV, vector<T>* W, const matrix<T>& MA );
Pascal/Delphi	uses MFstd; function MF_SVdecompose( MU, MV:fMatrix; W:fVector; MA:fMatrix; htA, lenA:UIntSize ): IntBool;

Description

The matrix MA of dimensions [ht, len] is decomposed into a product
MA = MU * MW * MV^T,
where MU has the dimensions [max(ht,len), len]. MW is a diagonal matrix with all elements positive or zero. Actually, only the diagonal of this matrix is stored in the vector W of size [len], MV, finally, is a square matrix [len, len]. Both MU and MV are orthogonal:
MU^T * MU = MV^T * MV = (1). Due to these orthogonality relations, the solution of a linear system MA * X = B is straightforward, once MA has been singular-value decomposed:
X = MV * W^-1 * MU^T
(This equation must be evaluated from right to left.) As MU and MV are orthogonal, only W must explicitly be inverted. This, however, is easy, as W is diagonal, and the inverse of a diagonal matrix consists of the inverse of the diagonal elements. This is the important point which makes SVD so useful, because it provides both a diagnose and a cure for singularities. If any element of W is very small compared to the largest element, this corresponds to a singularity (at least in the numerical sense, where dividing by extremely small numbers is almost as bad as dividing by zero). A singularity means that the underlying linear system is under-determined. This, in turn, means that one can arbitrarily choose one particular solution. Specifically, one may choose W^-1_ii = 0 for very small W_ii, which is one of the rare cases where it makes sense to set the result of a division by (near-)zero to zero instead of infinity. This reduction of the dimension by "Singular Value editing" is the crucial point in all functions relying on SVD, like MF_SVsolve, MF_safeSolve), VF_linfit and others.

If you wish to change the default editing threshold, you can do so by calling MF_SVDsetEdit. As MF_SVDsetEdit is not thread-safe, this function should not be used to set different editing thresholds for different calls to MF_SVsolve etc. Rather than repeatedly changing the default value then, use the "wEdit" variant of MF_SVsolve etc., namely MF_SVsolvewEdit, MF_solveBySVDwEdit and MF_safeSolvewEdit.

The singular values are stored in W in arbitrary order. For some applications, like visual inspection of the singular values, dimension reduction, calculation of the rank of a matrix, etc., it is desirable to bring the singular values into descending order, re-ordering also the right and left singular vectors (stored in the matrices MU and MV) accordingly. This is done by calling MF_SVsort.

If you compare results of MF_SVdecompose with SVD routines by other vendors, please note the following points:

MF_SVdecompose itself does not perform a dimension reduction (replacement of small singular values by 0), but leaves this step to following functions or to manual inspection by the user.
Singular values can be permutated, if the same permutation is also applied to the corresponding columns of both MU and MV. Often, this property in of SVD is exploited to bring the singular values into a certain (for example, into descending) order. As the price for this increased beauty is an increased work-load, MF_SVdecompose itself keeps the singular values in arbitrary order. Only if really needed, the task of sorting may be performed by calling the separate function MF_SVsort after MF_SVdecompose, as described above.
The sign of a singular value may be exchanged if the corresponding column of MV is also multiplied by -1. MF_SVdecompose uses this property to make all non-vanishing singular values positive.

Verifying SVD results

The following routine demonstrates SV decomposition and the verification of its results.

#include <VDstd.h>
#include <MDstd.h>
#include <conio.h>
void SVDtest( void )
{
    unsigned ht=386, len=52,            // choose any values
             maxhtlen, sizeA, sizeV;
    dMatrix MA, MU, MV, MOne, MA2, MUV2, MDiffA, MDiffOne;
    dVector W;
    double err;

    maxhtlen = (ht >= len ? ht : len);
    sizeA = ht*len;
    sizeV = len*len;

    MA   = MD_matrix( ht, len );   // input matrix
    MU   = MD_matrix( maxhtlen, len );
    MV   = MD_matrix( len, len );
    W    = VD_vector( len );
    MOne = MD_matrix( len, len );   // unity matrix for comparison
    MD_equ1( MOne, len );

    MA2  = MD_matrix( maxhtlen, len );
      // will be filled with MU*W*MVT to compare with input matrix
      // leading dimension maxhtlen is needed for under-determined matrices
    MUV2 = MD_matrix( len, len );
      // will be filled with MUT*MU, MVT*MV, MV*MVT for comparison with MOne
    MDiffA = MD_matrix( ht, len );  // will hold temporary result
    MDiffOne = MD_matrix( len, len );  // will hold temporary result

    MD_random( MA, ht, len, 1, -1., +1. );
      // fill MA with random numbers between -1 and +1
      // alternatively generate or read in any other test matrix
    MD_SVdecompose( MU, MV, W, MA, ht, len );

      /* check if MA = MU*W*MVT, evaluating from right to left : */
    MDdia_mulMT( MUV2, W, MV, len, len ); // use MUV2 to store W * MVT
    MD_mulM( MA2, MU, MUV2, maxhtlen, len, len );
    MD_subM( MDiffA, MA, MA2, ht, len );
      // for under-determined matrices, just ignore the rows from ht to maxhtlen-1
    err = VD_absmax( MDiffA[0], sizeA );
    printf( "SVDtest: check the SVD routine of OptiVec\n\n"
            "In the following tests in double precision,\n"
            "errors of up to a few times 1.e-14 are okay\n"
            "for matrices on the order of 100x100 to 1000x1000 elements:\n\n"
            "max. error of matrix reconstruction MA = MU*W*MVT: %lg", err );

    /* check if MU is orthogonal, i.e. MUT*MU = (1): */
    MD_TmulM( MUV2, MU, MU, maxhtlen, len, len );
    MD_subM( MDiffOne, MUV2, MOne, len, len );
    err = VD_absmax( MDiffOne[0], sizeV );
    printf( "\nmax. error of MUT*MU = (1): %lg", err );

    /* check if MV is orthogonal, i.e. MVT*MV = (1): */
    MD_TmulM( MUV2, MV, MV, len, len, len );
    MD_subM( MDiffOne, MUV2, MOne, len, len );
    err = VD_absmax( MDiffOne[0], sizeV );
    printf( "\nmax. error of MVT*MV = (1): %lg", err );

    /* check if MV is also orthonormal, i.e. MV*MVT = (1): */
    MD_mulMT( MUV2, MV, MV, len, len, len );
    MD_subM( MDiffOne, MUV2, MOne, len, len );
    err = VD_absmax( MDiffOne[0], sizeV );
    printf( "\nmax. error of MV*MVT = (1): %lg", err );
    printf( "\n\nHit any key to end SVDtest" ); _getch();

    M_nfree( 8, MA, MU, MV, MOne, MA2, MUV2, MDiffA, MDiffOne );
    V_free( W );
}