Inverted Matrix Values are NaN and Infinity


(Ayhan Tekgul) #1

Hi,
I am solving a least squares adjustment problem, and trying to invert a 3x3 matrix, but the result matrix values are getting NaN or Infinity. I have tried to invert this 3x3 matrix in Matlab14 with inv() function, there is no problem. I have tried to use another matrix library to get inversion of same 3x3 matrix, there is no problem.

This is the situation I have:
N is a square matrix with 39610x39610 size. It is symmetric and sparse.
N has a special design which is
| N11 N12 |
| N21 N22 |
N11 and N22 are hyperdiagonal matrixes like
| xxxx… |
| xxxx… |
| xxxx… |
| xxxx… |
| …xxxx… |
| …xxxx… |
| …xxxx… |
| …xxxx… |
| …xxxx… |
| …xxxx… |
| …xxxx… |
| …xxxx… |

and N21 = N12.transpose().

In somewhere in calculation I need to invert N22. Since it is hyperdiagonal, I can invert sub square diagonal matrixes individually and fill N22Inv matrix. For this purpose, I am getting submatrixes by tmp = N.SubMatrix(x, 3, x, 3), inverting them tmpi = tmp.Inversion() and writing to a new matrix by N22Inv.SetSubMatrix(x, y, tmpi). But tmpi matrixes are coming filled with NaN or Infinity values.

One of the matrixes is;
| 0.888889 0.0 2.29907 |
| 0.0 0.888889 -0.0609424 |
| 2.29907 -0.0609424 5.9506 |

Is there a bug or am I doing a wrong process?

Best regards…


(Christoph Rüegg) #2

You seem to have managed to pin down the issue to a small 3x3 matrix which should simplify things.

I tried to repro with the specific 3x3 matrix you posted. While it is very ill-conditioned (3e+6) I don’t seem to have problems with NaN or infinity. In LINQPad:

var m = Matrix<double>.Build.DenseOfArray(
new[,] { { 0.888889, 0.0, 2.29907 },
         { 0.0, 0.888889, -0.0609424 },
         { 2.29907, -0.0609424, 5.9506} });
m.ConditionNumber().Dump();
m.Inverse().Dump();

returns:

3351477.20148778
DenseMatrix 3x3-Double
-426034   11293.1    164718
11293.1  -298.226  -4366.25
 164718  -4366.25  -63684.9

Do you get a different result here?


(Ayhan Tekgul) #3

Hi Christoph,

The inverse matrix you have written is same as the Matlab result.

The condition number I got from double ic = tmp.ConditionNumber(); is 1.8146075638064666E+17.

And inversion matrix is still
tmpi.ToMatrixString();
" NaN NaN NaN\r\nInfinity -Infinity -Infinity\r\nInfinity -Infinity -Infinity\r\n"

I changed the matrix creation as
tmp = Matrix.Build.DenseOfArray(NMat.SubMatrix(0, 3, 0, 3).ToArray());

But results are still same as expected.

I wonder how are the results different between us?

best regards…


(Christoph Rüegg) #4

What version are you using? My test was with v3.8.0.


(Ayhan Tekgul) #5

I am using 3.7.0.

Now I will install the 3.8.0 and try again. Also I will inform you about the results.

Best regards…


(Ayhan Tekgul) #6

Hi Christoph,
I have installed 3.8.0, but the results are same.

tmp.ToMatrixString()
“0.888889 0 2.29907\r\n 0 0.888889 -0.0609424\r\n 2.29907 -0.0609424 5.9506\r\n"
tmp.ConditionNumber();
1.8146075638064666E+17
tmpi.ToMatrixString()
” NaN NaN NaN\r\nSonsuz -Sonsuz -Sonsuz\r\nSonsuz -Sonsuz -Sonsuz\r\n"

“Sonsuz” is “Infinity” in Turkish.

here is my code;

     private void Solvex()
     {
         int imgunk = Images.Count * 6;
         Matrix<double> N11 = SparseMatrix.Create(imgunk, imgunk, 0.0d);

         Matrix<double> tmp = null;
         for (int i = 0; i < Images.Count; i++)
         {
             int offset = i * 6;
             tmp = NMat.SubMatrix(offset, 6, offset, 6);
             N11.SetSubMatrix(offset, offset, tmp);
         }

         int tieunk = TIEs.Count * 3;

         Matrix<double> N22I = SparseMatrix.Create(tieunk, tieunk, 0.0d);
         Matrix<double> tmpi = null;
         for (int i = 0; i < TIEs.Count; i++)
         {
             int offset = i * 3;
             int offset2 = offset + imgunk;

             //tmp = NMat.SubMatrix(offset2, 3, offset2, 3);
             tmp = Matrix<double>.Build.DenseOfArray(NMat.SubMatrix(offset2, 3, offset2, 3).ToArray());

             double ic = tmp.ConditionNumber();
             tmpi = tmp.Inverse();
             N22I.SetSubMatrix(offset, offset, tmpi);
         }

         Matrix<double> N12 = NMat.SubMatrix(0, imgunk, imgunk, tieunk);
         Matrix<double> N21 = N12.Transpose();

         Matrix<double> n1 = nMat.SubMatrix(0, imgunk, 0, 1);
         Matrix<double> n2 = nMat.SubMatrix(imgunk, tieunk, 0, 1);

         Matrix<double> A = N12 * N22I;
         Matrix<double> B = (N11 - A * N21);
         Matrix<double> C = B.Inverse();
         Matrix<double> D = n1 - A * n2;

         Matrix<double> x1 = C * D;
         Matrix<double> x2 = N22I * (n2 - N21 * x1);
         xMat.SetSubMatrix(0, 0, x1);
         xMat.SetSubMatrix(imgunk, 0, x2);
     }

(Ayhan Tekgul) #7

Hi Christoph,
I have tried the code you have sent, and took the same result with yours.

mi.ToMatrixString()
"-426034 11293.1 164718\r\n11293.1 -298.226 -4366.25\r\n 164718 -4366.25 -63684.9\r\n"
m.ConditionNumber()
3351477.2004373926