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.*/
};

What do you observe when removing all the 1.0 e-7 in the matrix ? The eigenvectors should be unchanged, a difference could indicate a particular rounding issue (or a test issue) in the algorithm. — Damien
– Damien, Commented Feb 3, 2022 at 13:19
What eigenvalues did you get in both cases, matlab and C? In the documentation, it is mentioned that This algorithm has the drawback that the low order singular values may be relatively imprecise. — Damien
– Damien, Commented Feb 3, 2022 at 13:32
What is the shape of D? i.e. 1D, 2D? your initialization with test (2 2D) suggests you are only populating part of D. — ryyker
– ryyker, Commented Feb 3, 2022 at 13:39
I am not familiar with these algorithms, but int iterations = 0; looks suspicious. — ryyker
– ryyker, Commented Feb 3, 2022 at 14:03

ryyker · Accepted Answer · 2022-02-03 15:06:57Z

0

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];

edited Feb 3, 2022 at 15:06

answered Feb 3, 2022 at 13:56

ryyker

23.2k3 gold badges46 silver badges96 bronze badges

Sign up to request clarification or add additional context in comments.

6 Comments

Danijel Over a year ago

Thanks. Already tried that, it gives the same result.

ryyker Over a year ago

@Danijel - I am not suggesting that you edit that section in your post,(that would be a bad idea.) but why did you choose to show a partial initialization if you knew it was incorrect?

Danijel Over a year ago

I didn't think it's incorrect. After all, it works the same.

ryyker Over a year ago

@Danijel - There is something else, in addition to the incomplete population of value problem. Have you seen my comment on iterations being set to zero? Is that a valid value to pass as an argument into svdLAS2()?

Danijel Over a year ago

Yes, iterations ends up set to dimensions if set to 0, so setting it to 0 seems OK.

|

Collectives™ on Stack Overflow

Comparing SVD algorithm in C (svdlibc) with MATLAB?

1 Answer 1

6 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

6 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related