PJSR: 1185 Win7 Math.svd bug?

[mschuster]

Hi Juan,

Here is a 6x6 matrix that causes 1185 Win7 Math.svd to generate NaN values in its result.

The matrix is high dynamic range, but FYI both Matlab and Mathematica give reasonable results with no NaNs. It would be helpful if PI did likewise.

The priority of this is low, I am guessing that I will be able to work around this by zeroing the tiny values.

Thanks,
Mike

Code: [Select]

var matrix = new Matrix([
   4844580310818187,
   4844570391466517,
   46134247463213.92,
   8.207893678085278e-146,
   8.207893678041447e-146,
   -1.8615598742947518e-147,
   4844570391466517,
   4844564048327026,
   45898970675508.98,
   8.207893678041447e-146,
   8.207893677997616e-146,
   -1.861559874284811e-147,
   46134247463213.92,
   45898970675508.98,
   25823630441271668,
   -1.8615598742947518e-147,
   -1.8615598742848108e-147,
   4.222039540834665e-149,
   8.207893678085278e-146,
   8.207893678041447e-146,
   -1.8615598742947518e-147,
   1.7014839556157924e-300,
   1.7014839556067063e-300,
   -3.8589854873334084e-302,
   8.207893678041447e-146,
   8.207893677997616e-146,
   -1.8615598742848108e-147,
   1.7014839556067063e-300,
   1.7014839555976202e-300,
   -3.858985487312801e-302,
   -1.8615598742947518e-147,
   -1.861559874284811e-147,
   4.222039540834665e-149,
   -3.8589854873334084e-302,
   -3.858985487312801e-302,
   8.752224164265075e-304
], 6, 6);

var svd = Math.svd(matrix);
var u = svd[0];
var w = svd[1];
var v = svd[2];

console.writeln("u: ", u.toArray());
console.writeln("w: ", w.toArray());
console.writeln("v: ", v.toArray());


[mschuster]

Hi Juan,

FYI: update on my Math.svd(matrix) problem.

Maybe you are using the Numerical Recipes SVD implementation here, or something similar.

If so, please read on.

This implementation in several places assumes that for double precision v if v != 0 then isFinite(1 / v) holds. For small values this is not true.

Three places in the code do something like this: if (z != 0) {... 1 / z ...}

I change this to: if (z != 0 && isFinite(1 / z)) { ... 1 / z ...}

This solves the problem.

So my workaround now is to use my own PJSR SVD implementation.

Thanks,
Mike

Edit: Here are some small floating point values that print as non-zero, yet their inverses print as infinite. Output and code below.

Code: [Select]

g: 1.29892e-318, 1 / g: Infinity
g: -2.06070878658906e-310, 1 / g: -Infinity
g: 3.422704e-318, 1 / g: Infinity
g: 3.69996e-319, 1 / g: Infinity


Code: [Select]

g = 1.29892e-318;
console.writeln("g: ", g, ", 1 / g: ", 1 / g);

g = -2.06070878658906e-310;
console.writeln("g: ", g, ", 1 / g: ", 1 / g);

g = 3.422704e-318;
console.writeln("g: ", g, ", 1 / g: ", 1 / g);

g = 3.69996e-319;
console.writeln("g: ", g, ", 1 / g: ", 1 / g);

Edit: These numbers are "denormalized", i.e., smaller than the smallest IEEE "normalized" double precision number. This appears to be the reason why their inverses are not representable.

[Juan Conejero]

Hi Mike,

Yes, the current SVD routine in the core application is an adaptation of the NR3 implementation. I have implemented your solution and it seems to work well in all of my tests. For your example 6x6 matrix, the output of your test script is now:

Code: [Select]

u: 0.7071062183462742,0.7071015740954208,-0.0028565517339196995,-3.3015825279789923e-161,0,0,-0.7071073439970429,0.7071004851896847,-0.0028474391933923254,5.698060336481597e-161,0,0,-0.000006440381235082195,-0.004033330669572706,-0.9999918660680353,-1.442544086577126e-163,0,0,-4.499243768466547e-161,1.198112770725562e-161,5.395651329327801e-164,-0.7070158669675028,0,0,-4.499243768476779e-161,1.198112770719164e-161,5.395651329299175e-164,-0.7070158669637272,0,0,1.0204337425167995e-162,-2.7173337598387453e-163,-1.223740024457242e-165,0.016035202453346702,0,0
w: 1787028157.7260127,9688880090589428,25823892922799290,3.403840351129283e-300,0,0
v: 0.7071062183462746,0.7071015740954211,-0.0028565517339195212,0,0,0,-0.7071073439970428,0.7071004851896849,-0.0028474391933924607,2.3964725526739963e-161,-3.064165658794164e-171,8.985720755944261e-173,-0.00000644038123510659,-0.004033330669572709,-0.9999918660680357,-1.4455511948398978e-163,1.848292220693647e-173,-5.420170265107162e-175,-1.1980831043289655e-161,1.198112770725431e-161,5.395651328686843e-164,-0.7070158670579629,-0.7070158668733328,0.016035202450446904,-1.1980831043225685e-161,1.1981127707190339e-161,5.395651328658029e-164,-0.7070158668733328,0.7071664970096749,0.0066414980550424775,2.717266476104717e-163,-2.7173337598384507e-163,-1.2237400243117817e-165,0.016035202450446907,0.0066414980550424775,0.9998493700482879
where [25823892922799290, 9688880090589428, 1787028157.7260127, 3.403840351129283e-300, 0, 0] look like good singular values. As you see, the new routine replaces NaNs with zeros in the output matrices. Let me know if the above results are valid.

Thank you for detecting this problem, and for your good analysis. We haven't seen it before because the existing tools don't pose so ill-conditioned problems, but it's nice to have a more capable implementation. The improved SVD routine will be available in the incoming 1.8.4.1187 version.

[Juan Conejero]

Hi Mike,

With this modified version of your script and the new version 1.8.4.1187 of PixInsight:

Code: [Select]

var matrix = new Matrix([
   4844580310818187,
   4844570391466517,
   46134247463213.92,
   8.207893678085278e-146,
   8.207893678041447e-146,
   -1.8615598742947518e-147,
   4844570391466517,
   4844564048327026,
   45898970675508.98,
   8.207893678041447e-146,
   8.207893677997616e-146,
   -1.861559874284811e-147,
   46134247463213.92,
   45898970675508.98,
   25823630441271668,
   -1.8615598742947518e-147,
   -1.8615598742848108e-147,
   4.222039540834665e-149,
   8.207893678085278e-146,
   8.207893678041447e-146,
   -1.8615598742947518e-147,
   1.7014839556157924e-300,
   1.7014839556067063e-300,
   -3.8589854873334084e-302,
   8.207893678041447e-146,
   8.207893677997616e-146,
   -1.8615598742848108e-147,
   1.7014839556067063e-300,
   1.7014839555976202e-300,
   -3.858985487312801e-302,
   -1.8615598742947518e-147,
   -1.861559874284811e-147,
   4.222039540834665e-149,
   -3.8589854873334084e-302,
   -3.858985487312801e-302,
   8.752224164265075e-304
], 6, 6);

/*
var matrix = new Matrix( 6, 6 );
matrix.setRandom();
*/

function matrixToString( A )
{
   let s = '';
   for ( let i = 0; i < A.rows; ++i )
   {
      for ( let j = 0; j < A.cols; ++j )
         s += format( "%+20.10g", A.at( i, j ) );
      s += '\n';
   }
   return s;
}

// Perform SVD decomposition
var svd = Math.svd(matrix);
var u = svd[0];
var w = svd[1];
var v = svd[2];

// Diagonal matrix W
var W = new Matrix( w.length, w.length )
W.assign( 0 );
for ( let i = 0; i < w.length; ++i )
   W.at( i, i%w.length, w.at( i ) );

// Reconstructed matrix A = u*W*vt
var A = u.mul( W ).mul( v.transpose() );

// Show the singular values sorted in descending order.
w.sort();
w.reverse()

console.writeln( "M:\n", matrixToString( matrix ) );
console.writeln( "u:\n", matrixToString( u ) );
console.writeln( "w:\n", w.toArray(), '\n' );
console.writeln( "v:\n", matrixToString( v ) );
console.writeln();
console.writeln( "W:\n", matrixToString( W ) );
console.writeln( "v<sup>T</sup>:\n", matrixToString( v.transpose() ) );
console.writeln( "A:\n", matrixToString( A ) );

I get the following output:

Code: [Select]

M:
    +4.844580311e+15    +4.844570391e+15    +4.613424746e+13   +8.207893678e-146   +8.207893678e-146   -1.861559874e-147
    +4.844570391e+15    +4.844564048e+15    +4.589897068e+13   +8.207893678e-146   +8.207893678e-146   -1.861559874e-147
    +4.613424746e+13    +4.589897068e+13    +2.582363044e+16   -1.861559874e-147   -1.861559874e-147   +4.222039541e-149
   +8.207893678e-146   +8.207893678e-146   -1.861559874e-147   +1.701483956e-300   +1.701483956e-300   -3.858985487e-302
   +8.207893678e-146   +8.207893678e-146   -1.861559874e-147   +1.701483956e-300   +1.701483956e-300   -3.858985487e-302
   -1.861559874e-147   -1.861559874e-147   +4.222039541e-149   -3.858985487e-302   -3.858985487e-302   +8.752224164e-304

u:
       +0.7071062183       +0.7071015741     -0.002856551734   -3.301582528e-161                  +0                  +0
        -0.707107344       +0.7071004852     -0.002847439193   +5.698060336e-161                  +0                  +0
    -6.440381235e-06      -0.00403333067       -0.9999918661   -1.442544087e-163                  +0                  +0
   -4.499243768e-161   +1.198112771e-161   +5.395651329e-164        -0.707015867                  +0                  +0
   -4.499243768e-161   +1.198112771e-161   +5.395651329e-164        -0.707015867                  +0                  +0
   +1.020433743e-162    -2.71733376e-163   -1.223740024e-165      +0.01603520245                  +0                  +0

w:
25823892922799290,9688880090589428,1787028157.7260127,3.403840351129283e-300,0,0

v:
       +0.7071062183       +0.7071015741     -0.002856551734                  +0                  +0                  +0
        -0.707107344       +0.7071004852     -0.002847439193   +2.396472553e-161   -3.064165659e-171   +8.985720756e-173
    -6.440381235e-06      -0.00403333067       -0.9999918661   -1.445551195e-163   +1.848292221e-173   -5.420170265e-175
   -1.198083104e-161   +1.198112771e-161   +5.395651329e-164       -0.7070158671       -0.7070158669      +0.01603520245
   -1.198083104e-161   +1.198112771e-161   +5.395651329e-164       -0.7070158669        +0.707166497     +0.006641498055
   +2.717266476e-163    -2.71733376e-163   -1.223740024e-165      +0.01603520245     +0.006641498055         +0.99984937


W:
         +1787028158                  +0                  +0                  +0                  +0                  +0
                  +0    +9.688880091e+15                  +0                  +0                  +0                  +0
                  +0                  +0    +2.582389292e+16                  +0                  +0                  +0
                  +0                  +0                  +0   +3.403840351e-300                  +0                  +0
                  +0                  +0                  +0                  +0                  +0                  +0
                  +0                  +0                  +0                  +0                  +0                  +0

vT:
       +0.7071062183        -0.707107344    -6.440381235e-06   -1.198083104e-161   -1.198083104e-161   +2.717266476e-163
       +0.7071015741       +0.7071004852      -0.00403333067   +1.198112771e-161   +1.198112771e-161    -2.71733376e-163
     -0.002856551734     -0.002847439193       -0.9999918661   +5.395651329e-164   +5.395651329e-164   -1.223740024e-165
                  +0   +2.396472553e-161   -1.445551195e-163       -0.7070158671       -0.7070158669      +0.01603520245
                  +0   -3.064165659e-171   +1.848292221e-173       -0.7070158669        +0.707166497     +0.006641498055
                  +0   +8.985720756e-173   -5.420170265e-175      +0.01603520245     +0.006641498055         +0.99984937

A:
    +4.844580311e+15    +4.844570391e+15    +4.613424746e+13   +8.207897849e-146   +8.207897849e-146    -1.86156082e-147
    +4.844570391e+15    +4.844564048e+15    +4.589897068e+13   +8.207889507e-146   +8.207889507e-146   -1.861558928e-147
    +4.613424746e+13    +4.589897068e+13    +2.582363044e+16   -1.861559875e-147   -1.861559874e-147   +4.222039541e-149
   +8.207893678e-146   +8.207893678e-146   -1.861559874e-147   +1.701483956e-300   +1.701483955e-300   -3.858985487e-302
   +8.207893678e-146   +8.207893678e-146   -1.861559874e-147   +1.701483956e-300   +1.701483955e-300   -3.858985487e-302
   -1.861559874e-147   -1.861559874e-147   +4.222039541e-149   -3.858985488e-302   -3.858985487e-302   +8.752224163e-304

As you see, the original matrix is now very well reconstructed, considering the difficulty of the problem. We have a much better SVD implementation now. Thank you for your help on these critical topics.

[mschuster]

Hi Juan,

Edited reply:

Thanks, this looks good. Comments:

1) My NR3SVD includes a final reorder(), that sort w and interchanges rows/cols of u and v to match. I was initially confused by your output above where w was sorted but u and v did not match.

2) NR3SVD throws an exception on non-convergence. How do you handle this? Is an exception thrown in PJSR and if so in what form? If not, something like a failure flag would be helpful: [u, w, v, failure] = Math.svd(matrix).

Thanks,
Mike

[Juan Conejero]

1) My NR3SVD includes a final reorder(), that sort w and interchanges rows/cols of u and v to match. I was initially confused by your output above where w was sorted but u and v did not match.

No, the PCL/PJSR SVD implementation does not reorder matrix and vector components. Writing a simple routine to detect the index of the largest or smallest singular value is very easy. Once you have the index, you can extract the corresponding eigenvector as a column vector from V. For example, this is what I do in StarAlignment to implement the eight-point algorithm. The SVD classes in PCL have IndexOfSmallestSingularValue() and IndexOfLargestSingularValue() member functions to facilitate these tasks.

2) NR3SVD throws an exception on non-convergence. How do you handle this? Is an exception thrown in PJSR and if so in what form? If not, something like a failure flag would be helpful: [u, w, v, failure] = Math.svd(matrix).

When there is no convergence after 30 iterations (as NR does) the Math.svd() routine throws an exception that you can (should) handle with a try{}catch{} construct. The exception is a normal Error object with the appropriate error message "No convergence".