ericjang / svd3

Fast singular value decomposition, diagonalization, QR decomposition of 3x3 matrices.
MIT License
149 stars 19 forks source link

Does not pass even simple test cases. #7

Open changyoonpark opened 7 years ago

changyoonpark commented 7 years ago

For anybody who is looking for a 3x3 SVD implementation in c++, be aware when using this snippet. This implementation even fails simple test cases, which may be the reason why the "snow simulation" video linked on the documentation does not look anything like the original simulation by Stomakhin.

ericjang commented 7 years ago

That's harsh. I can take another look though, there was a recent bug found in the last few months. Can you please provide an example test case?

changyoonpark commented 7 years ago

Sorry, didn't mean to be a jerk about it.

I wasn't trying be harsh on you, I actually think your implementation is compact and much more readable compared to the madly-optimized original implementation by Sifakis (2011), but I wasn't able to produce desirable results (three identity matricies) using a identity matrix as the input. I'm writing a visco-elasto-plasticity solver.

Here is what I'm doing to test your code.

    float u11, u12, u13;
    float u21, u22, u23;
    float u31, u32, u33;

    float s11, s12, s13;
    float s21, s22, s23;
    float s31, s32, s33;

    float v11, v12, v13;
    float v21, v22, v23;
    float v31, v32, v33;

  svd(1.0,0.0,0.0,
        0.0,1.0,0.0,
        0.0,0.0,1.0,

        u11,u12,u13,
        u21,u22,u23,
        u31,u32,u33,

        s11,s12,s13,
        s21,s22,s23,
        s31,s32,s33,

        v11,v12,v13,
        v21,v22,v23,
        v31,v32,v33        

      );

  printf("%f %f %f \n",u11,u12,u13);
  printf("%f %f %f \n",u21,u22,u23);
  printf("%f %f %f \n\n",u31,u32,u33);

  printf("%f %f %f \n",s11,s12,s13);
  printf("%f %f %f \n",s21,s22,s23);
  printf("%f %f %f \n\n",s31,s32,s33);

  printf("%f %f %f \n",v11,v12,v13);
  printf("%f %f %f \n",v21,v22,v23);
  printf("%f %f %f \n\n",v31,v32,v33);

The output U,S,V I get are

U : 
0.000497 -0.000000 0.996568 
0.000000 -0.997067 -0.000000 
0.999499 0.000000 -0.000495 

S : 
0.999499 -0.000000 0.000472 
0.000000 0.948047 -0.000000 
-0.000495 0.000000 0.947571 

V:
0.000000 -0.000000 0.950835 
-0.000000 -0.950835 -0.000000 
1.000000 0.000000 -0.000000 

U * S * V.T() : 
8.97891719e-01   0.00000000e+00   3.44984300e-06
0.00000000e+00   8.98792357e-01   0.00000000e+00
2.58235661e-06   0.00000000e+00   9.98998496e-01

For continuum mechanical purposes, the singular values are not close enough to 1. This means that the continuum model will exhibit stress even in zero-deformation states. If you were to use these singular values as a measure of strain in continuum mechanical simulations with materials with huge work function gradients, it will deteriorate the simulation very quickly.

I am currently using the code written by Sifakis (http://pages.cs.wisc.edu/~sifakis/project_pages/svd.html), which gives me more desirable results (almost exact identity matricies for U,S,V).

If you think I'm somehow not using your SVD function properly please let me know, I would be actually very happy to replace my SVD solver to what you have here.

Thanks for your reply.

xchern commented 6 years ago

I experienced exactly the same problem while implementing the paper of snow. The actual problem of this code is that it uses the famous fast inverse square algorithm (in the rsqrt, rsqrt1 functions, https://github.com/ericjang/svd3/blob/master/svd3.h#L36).

Adding more iterations to the newton iterations solves the problem.

l0calh05t commented 6 years ago

We replaced those by actual sqrt calls (and divisions). I have to agree that the original is unusable as-is.

xchern commented 6 years ago

@l0calh05t But I'm fine with my snow simulator now... Perhaps you could add more iterations steps (from 4 to 18 or more) additionally at https://github.com/ericjang/svd3/blob/master/svd3.h#L231, according to the other issue of this repo.

evs625 commented 6 years ago

If you add more iterations (from 4 to 8 on line https://github.com/ericjang/svd3/blob/master/svd3.h#L231) it makes things worse...