NaN and -infinity when solving system of linear equations


#1

Hi there,

I’m trying to do something really simple but am not having any luck. I’m trying to solve a system of linear equations Ax=B

where A, x, and B are 4x4 matrices (x is a rotational and translational transform). I split up B into 4 column vectors and am trying to compute the columns of x separately (I couldn’t figure out how to do all at the same time from looking at the functions available in Math.Net). I’ve tried the following but get NaN and seemingly wrong answers when I compare it to the output from python

        double[][] ddp = new double[][] {
                    new double[] { 38, -144.5, 0, 1 },
                    new double[] { 52, -169.5, 0, 1 },
                    new double[] { 430, -144.5, 0, 1 },
                    new double[] { 66, -144.5, 0, 1 }
        };
        var cp1 = new double[] { 38, 52, 430, 66 };
        double[] q = MathNet.Numerics.LinearRegression.MultipleRegression.QR(ddp, cp1);

I should get an answer like [1, 0, 0, 0] but get [NaN, NaN, -infinity, 0] instead. I’ve also tried using Fit.LinearMultiDimFunc but get an error that the “number of matrix columns must be positive”.

Can someone give me tips on how to proceed?

Thanks in advance!
Robin


(Peter Vanderwaart) #2

The reason you are getting NaN values is that the regression arithmetic will only work if there is a unique solution. In your example, because of the column of zeros in ddp, the solution vector q could have any value in the third place. Multiple solutions leads to division by zero in the algorithm which leads to NaN value. The algebra class language is that the rows of ddp have to be linearly independent.

For your original problem, solving Ax=B where A and B are both 4x4, The easiest solution is to multiply both sides by the inverse of A giving x = A’B. In order for A’ to exist, the rows of A must be linearly independent for the same reasons as described above. The way to test for this is calculate the Rank which should be equal to the #rows = 4. B does not have to have rank = 4.

To program this in MathNet, the easy way is to use Matrices rather than arrays and use the built-in functions. Note: the MathNet matrix inverse does not give an error if the rows are not linearly independent and returns an incorrect answer. Check the rank explicitly. There is some helpful info on using MathNet matrices and vectors here: https://numerics.mathdotnet.com/Matrix.html.

Example below:

    // build matrices
    Matrix<double> A = Matrix<double>.Build.Dense(4, 4);
    Matrix<double> B = Matrix<double>.Build.Dense(4, 4);
    
    for (int r = 0; r < 4; r++)
    {
        for (int c = 0; c < 4; c++)
        {
            A[r, c] = Math.Tan(r + c);
            B[r, c] = Math.Cos(r - c);
        }
    }

    // calculate solution matrix x
    var x = A.Inverse() * B;
    // check that multiplication gives desired answer
    var z = A * x;

    Console.WriteLine("A Determinant = " + A.Determinant().ToString());
    Console.WriteLine("A Rank = " + A.Rank().ToString());

    Console.WriteLine("B Determinant = " + B.Determinant().ToString());
    Console.WriteLine("B Rank = " + B.Rank().ToString());


    Console.WriteLine(A);
    Console.WriteLine(B);
    Console.WriteLine(x);
    Console.WriteLine(z);

#3

Thanks so much!

As usual, it was my understanding of the problem that was flawed so thanks for helping to clarify it.

So yes, I check the rank and if it’s not equal to the number of columns then I try to remove columns and recalculate the rank. Then I can calculate the transformation matrix since I know the matrix will be invertible.

Thanks again!