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 );
} |