Comparing SVD algorithm in C (svdlibc) with MATLAB?

Question

I wanted to compare the SVD (singular value decomposition) algorithm C code from https://github.com/lucasmaystre/svdlibc to MATLAB SVD implementation.

C code is in fact SVDLIBC.

I built the C project without problems, and replaced the main() with my simple test:

#define M (4)
#define N (4)

int main(int argc, char *argv[])
{

    // Test array
    double test[M][N] = {
            { 1.0e-07* 0.609979111792474,  1.0e-07* 0.442691149480391,  1.0e-07* 0.063440080518362, 1.0e-07*-0.346712389824400,},
            { 1.0e-07* 0.442691149480391,  1.0e-07* 0.484281334855962,  1.0e-07* 0.055069632301392, 1.0e-07*-0.300966260855239,},
            { 1.0e-07* 0.063440080518362,  1.0e-07* 0.055069632301392,  1.0e-07* 0.107891781688921, 1.0e-07*-0.043130123212035,},
            { 1.0e-07*-0.346712389824400,  1.0e-07*-0.300966260855239,  1.0e-07*-0.043130123212035, 1.0e-07* 0.335714519434403 }
    };

    // Create the input matrix
    DMat D = svdNewDMat(M, N);
    // Initialize with data
    for(int i=0;i<M;i++)
        D->value[i]=test[i];

    // We need sparse S matrix as input to main algo
    SMat S = svdConvertDtoS(D);

    // Other inputs (values taken from original main() )
    int dimensions = M;
    int iterations = 0;
    double las2end[2] = {-1.0e-30, 1.0e-30};
    double kappa = 1e-6;

    // Calculate SVD algo
    SVDRec R = svdLAS2(S, dimensions, iterations, las2end, kappa);

    // Print results from R
    for(int i=0;i<M;i++){
        for(int j=0;j<N;j++){
            printf("U[%d][%d] = %+.15f  ", i,j, R->Ut->value[i][j]);
        }
        printf("\n");
    }
    printf("\n");

    for(int i=0;i<M;i++){
        printf("S[%d] = %+.15f  ", i, R->S[i]);
    }
    printf("\n\n");

    for(int i=0;i<M;i++){
        for(int j=0;j<N;j++){
            printf("V[%d][%d] = %+.15f  ", i,j, R->Vt->value[i][j]);
        }
        printf("\n");
    }

    return 0;
}

However, the result I get is:

U[0][0] = +0.669469029923469  U[0][1] = +0.581137555529811  U[0][2] = +0.083280213210335  U[0][3] = -0.455142577236855  
U[1][0] = +0.000000000000000  U[1][1] = +0.000000000000000  U[1][2] = +0.000000000000000  U[1][3] = +0.000000000000000  
U[2][0] = +0.000000000000000  U[2][1] = +0.000000000000000  U[2][2] = +0.000000000000000  U[2][3] = +0.000000000000000  
U[3][0] = +0.000000000000000  U[3][1] = +0.000000000000000  U[3][2] = +0.000000000000000  U[3][3] = +0.000000000000000  

Expected values calculated using MATLAB are (transposed):

  -0.669469029923469  -0.581137555529811  -0.083280213210335   0.455142577236855
   0.000000000000000   0.112110571326639  -0.992948073424445  -0.038540151524073
  -0.000000000000000  -0.612706096707861  -0.038540151524073  -0.789370569363666
   0.742839967942847  -0.523738102878401  -0.075054555430211   0.410187756191207

Seems like I get only the first row, and that one is OK (only the sign is inverted).
But why are the other values all zero?

Am I handling the input array properly? The input DMat D is declared here as:

struct dmat {
  long rows;
  long cols;
  double **value; /* Accessed by [row][col]. Free value[0] and value to free.*/
};

Answer 1

One problem, based on seeing how value is defined:

struct dmat {
  long rows;
  long cols;
  double **value; /* Accessed by [row][col]. Free value[0] and value to free.*/
};

In this call DMat D = svdNewDMat(M, N);, member value is created as a dynamically allocated 2D array.

This then is not properly initializing value with test

 // Create the input matrix
    DMat D = svdNewDMat(M, N);
    // Initialize with data
    for(int i=0;i<M;i++)
        D->value[i]=test[i];

Change to the initialization to assign all values in testtest[i][j] to corresponding locations within D->value[i][j]:

 // Create the input matrix
    DMat D = svdNewDMat(M, N);
    // Initialize with data
    for(int i=0;i<M;i++)
       for(int j=0; j<N;  j++)
           D->value[i][j]=test[i][j];