Matrix Inverse NaN/Infinity/-Infinity

[koryakinp]

I am trying to Inverse 16x16 matrix, but the resulting values are NaN/Infinity/-Infinity.

var json = File.ReadAllText("data.json"); var arr = JsonConvert.DeserializeObject<double[]>(json);

var matrix = Matrix.Build.Dense(16,16, arr); var inverted = matrix.Inverse();

Here is data.json

I am using MathNet.Numerics 4.5.1

[gaston-e-nusimovich]

@koryakinp Pavel, Hi: did you calculate the determinant of said matrix?

If the determinant of any given matrix is either zero or too close to zero (from either the negative side or the positive side of the real axis), when you attempt to calculate the inverse matrix of said matrix, you may get -Infinity or +Infinity as result.

If you get NaN as a result, it could mean that at least one of the items in the array given as input is not a number (beware: JSON data is basically text data, so, it could easily "hide" a value that "almost" looks like a number, but actually, it is not!)

There is an old principle in Information Systems, a principle with a rather funny name: GIGO.

It is an acronym that says: "Garbage In, Garbage Out".

If your input contains Garbage data, the result will certainly be Garbage as well!

The only practical way to mitigate GIGO is with proper validation measures with ALL input.

Kind regards, GEN

[koryakinp]

@gaston-e-nusimovich, Python calculates inverted matrix with no problem, although the determinant is indeed close to 0.

import numpy as np
import pandas as pd
arr = np.array(pd.read_json("data.json")[0])
arr.shape = (16,16)
inversed = np.linalg.inv(arr)

[gaston-e-nusimovich]

Pavel, the point is that Python is designed to work on just about any computer, which means that it has to work with the least amount of resources.

Based on similar principles, it is designed to be easy to learn, so, maybe it is "too tolerant" to many things that are not consistent with how mathematics is consistent within itself.

Please, do not take me wrong, I am a huge fan of Python, but, at the very same time, I am fully aware all of Python's shortcomings.

If the problem we are facing is math based in such a way that it requires a numerical solution (that is a mathematical solution that requires one or more computers to be solved), that means that we have to know how robust and stable is the set of algorithms that we are using, and the implementation of said set of algorithms in whatever language/OS/Computer we are using, so, we have to be incredibly responsible with our choices of tools.

If I were in your shoes, I would choose C# and a set of libraries like MathNet to solve ALL THINGS MATHEMATICS, but that choice may require from you to be more EXIGENT WITH YOUR MATH: C# and MathNet are as exigent as the math behind them, which happens to be HOW IT SHOULD BE!

The fact that Python is SO tolerant is not necessarily A DREAM COME TRUE.

Kind regards, GEN

[koryakinp]

@gaston-e-nusimovich, when using a library, specifically designed for numerical computation, I would expect it to be able to solve such fundamental operations as matrix inversion. Provided a valid input I would expect to receive back a sensible answer. NaN, Infinity or -Infinity is not a sensible answer. It is confusing and not intuitive. I would much prefer an exception to be thrown.

Stepping all that aside, I found a workaround simply to limit the precision of the matrix to 3 decimals:

var json = File.ReadAllText("data.json");
var arr = JsonConvert.DeserializeObject<double[]>(json);
var matrix = Matrix.Build.Dense(16,16, arr);
matrix.MapInplace(q => Math.Round(q,3))
var inverted = matrix.Inverse();

[gaston-e-nusimovich]

It is interesting how the rounding to 3 positions works the way it does, in this case in particular.

I wonder if the condition number of the matrix without rounding to 3 positions is much higher than the value of the condition number for the matrix after rounding to 3 positions.

If that is so, then, that system of equations represented by the matrix (that is, the original matrix) unluckily happens to be ill conditioned.

[koryakinp]

@gaston-e-nusimovich here you are:

var json = File.ReadAllText("data.json");
var arr = JsonConvert.DeserializeObject<double[]>(json);

var matrix = Matrix<double>.Build.Dense(16,16, arr);
var cn1 = matrix.ConditionNumber(); //Infinity

matrix.MapInplace(q => Math.Round(q, 3));
var cn2 = matrix.ConditionNumber(); //87130

[gaston-e-nusimovich]

Pavel, now that you have all this quantitative evidence, you can validate what it really means from a mathematical standpoint, considering the fact that the concept of ill conditioned and Condition Number are mathematical tools that are strictly important ONLY within the realm of numerical solutions.

Kind regards, GEN

[wo80]

@koryakinp You haven't told us what problem you are actually trying to solve. Since the linear system is ill-conditioned, computing the inverse doesn't make sense. If you are trying to solve a linear system, your rounding approach will completely mess up the solutions.

My general advice: never use matrix.Inverse()

In Math.NET, inverting the matrix is actually done using the LU factorization. So, instead of inverting the matrix, use the LU or - even better - the QR factorization

var lu = matrix.LU();

Console.WriteLine("LU determinant: {0}", lu.Determinant);

var qr = matrix.QR();

Console.WriteLine("QR determinant: {0}", qr.Determinant);
Console.WriteLine("QR full rank: {0}", qr.IsFullRank);

// Using Linq to determine the numerical rank with threshold 1e-8.
Console.WriteLine("QR numerical rank: {0}", qr.R.Diagonal().Count(a => Math.Abs(a) > 1e-8));

Now, lu.Solve(vector) won't work, since the matrix is not invertible, but using qr.Solve(vector) you will get a least-squares solution for the system.

[hhaldn]

I agree with @koryakinp, that it is not a good idea for Matrix.Inverse to return a matrix containing entries that are NaN or Infinity. I void much prefer Inverse to throw an exception, if the matrix is not invertible and the inverse therefore cannot be computed.

This is in fact the behavior that you get, if you use the MKL-provider, In this case Inverse will throw a SingularUMatrixException, if the matrix is not invertible. In my opinion, the Inverse method should try to handle a non-invertible matrix consistently, regardless of the provider used.

I therefore think the behavior of Matrix.Inverse should be changed to throw an exception instead of returning a matrix with entries that may be NaN or Infinity.

I am using version 4.8.1 of MathNet.Numerics and version 2.3.0 of MathNet.Numerics.MKL.Win.

[GeorchW]

I do think, too, that inverting a matrix should either return a useful result or throw an exception - I've just run into the same issue. Otherwise the error just propagates through the program, getting harder and harder to debug the larger the program is.

But if there are any reasons against this, I would at least leave a note in the member documentation for Matrix.Invert(). The same goes for @wo80's advice, although I think that there are valid reasons to require an inverted matrix.

This is slightly different to the original issue, should I open a new one?

[RejectKid]

@wo80 hitting a similar thing to OP's comment and i did your implementation of least squares. I get results!! but oddly not exactly the same as numpy.lstsq. 12/20 of the numbers are the same tho interestingly. is the qr.Solve(vector) the same as the numpy.lstsq implementation? (sorry if im a little math dumb in this area)

[wo80]

@RejectKid numpy is probably using some native LAPACK routine under the hood, which there are several of for least squares - so you'd have to find out which routine is called and what parameters are used.

But even if exactly the same QR algorithm is used, there might be deviations in the solution just because of different compiler flags used to compile LAPACK and the core CLR, and thus different CPU features used.

Would you consider one of the solutions better?

[RejectKid]

@wo80 thanks so much for the reply, that was so useful!

Been reading on this problem space more, so it would seem that since the A's determinant is 0 0 or infinite solutions exist. thus your answer in mathnet with the QR.Solve(vector) works as well as python's lstsq's works, and what we do with the values after is what matters. In our case the interpolation of the coefficients received. And as long as thats consistent then we're okay regardless if the coefficients created here are different with python, mathnet, octave, or matlabb