Can Math.NET solve any matrix?

Question

I am trying to use Math.NET to solve the following system:

Coefficient Matrix A:

var matrixA = DenseMatrix.OfArray(new[,] {
    { 20000, 0, 0, -20000, 0, 0, 0, 0, 0 },
    { 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 }, 
    { 0, 2000,  8000,  0,  -2000, 4000, 0, 0, 0 },
    { -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 },
    { 0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 },
    { 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 },
    { 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 },
    { 0, 0, 0, 0, -20000, 0, 0, 20000, 0 },
    { 0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 }});

Result Vector b:

double[] loadVector = { 0, 0, 0, 5, 0, 0, 0, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);

I pulled these numbers from a Finite Element Analysis example problem, so the answer I expect based on that example is:

[0.01316, 0, 0.0009199, 0.01316, -0.00009355, -0.00188, 0, 0, 0]

However, Math.NET and an online Matrix Calculator I found mostly give me zeros (from iterative solvers), NaN, or large incorrect numbers (from direct solvers) as the solution.

In Math.NET, I have tried plugging my matrices in to the examples provided including:

Iterative Example:

namespace Examples.LinearAlgebra.IterativeSolversExamples
{
/// <summary>
/// Composite matrix solver
/// </summary>
public class CompositeSolverExample : IExample
{
    public void Run()
    {
        // Format matrix output to console
        var formatProvider = (CultureInfo)CultureInfo.InvariantCulture.Clone();
        formatProvider.TextInfo.ListSeparator = " ";

        // Solve next system of linear equations (Ax=b):
        // 5*x + 2*y - 4*z = -7
        // 3*x - 7*y + 6*z = 38
        // 4*x + 1*y + 5*z = 43

        // Create matrix "A" with coefficients
        var matrixA = DenseMatrix.OfArray(new[,] { { 20000, 0, 0, -20000, 0, 0, 0, 0, 0 }, { 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 }, 
                                                { 0, 2000,  8000,  0,  -2000, 4000, 0, 0, 0 }, { -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 },
                                                {0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 }, { 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 },
                                                { 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 }, { 0, 0, 0, 0, -20000, 0, 0, 20000, 0 },
                                                 {0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 }});


        Console.WriteLine(@"Matrix 'A' with coefficients");
        Console.WriteLine(matrixA.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // Create vector "b" with the constant terms.
        double[] loadVector = {0,0,0,5,0,0,0,0,0};
        var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
        Console.WriteLine(@"Vector 'b' with the constant terms");
        Console.WriteLine(vectorB.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // Create stop criteria to monitor an iterative calculation. There are next available stop criteria:
        // - DivergenceStopCriterion: monitors an iterative calculation for signs of divergence;
        // - FailureStopCriterion: monitors residuals for NaN's;
        // - IterationCountStopCriterion: monitors the numbers of iteration steps;
        // - ResidualStopCriterion: monitors residuals if calculation is considered converged;

        // Stop calculation if 1000 iterations reached during calculation
        var iterationCountStopCriterion = new IterationCountStopCriterion<double>(500000);

        // Stop calculation if residuals are below 1E-10 --> the calculation is considered converged
        var residualStopCriterion = new ResidualStopCriterion<double>(1e-10);

        // Create monitor with defined stop criteria
        var monitor = new Iterator<double>(iterationCountStopCriterion, residualStopCriterion);

        // Load all suitable solvers from current assembly. Below in this example, there is user-defined solver
        // "class UserBiCgStab : IIterativeSolverSetup<double>" which uses regular BiCgStab solver. But user may create any other solver
        // and solver setup classes which implement IIterativeSolverSetup<T> and pass assembly to next function:
        var solver = new CompositeSolver(SolverSetup<double>.LoadFromAssembly(Assembly.GetExecutingAssembly()));

        // 1. Solve the matrix equation
        var resultX = matrixA.SolveIterative(vectorB, solver, monitor);
        Console.WriteLine(@"1. Solve the matrix equation");
        Console.WriteLine();

        // 2. Check solver status of the iterations.
        // Solver has property IterationResult which contains the status of the iteration once the calculation is finished.
        // Possible values are:
        // - CalculationCancelled: calculation was cancelled by the user;
        // - CalculationConverged: calculation has converged to the desired convergence levels;
        // - CalculationDiverged: calculation diverged;
        // - CalculationFailure: calculation has failed for some reason;
        // - CalculationIndetermined: calculation is indetermined, not started or stopped;
        // - CalculationRunning: calculation is running and no results are yet known;
        // - CalculationStoppedWithoutConvergence: calculation has been stopped due to reaching the stopping limits, but that convergence was not achieved;
        Console.WriteLine(@"2. Solver status of the iterations");
        Console.WriteLine(monitor.Status);
        Console.WriteLine();

        // 3. Solution result vector of the matrix equation
        Console.WriteLine(@"3. Solution result vector of the matrix equation");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 4. Verify result. Multiply coefficient matrix "A" by result vector "x"
        var reconstructVecorB = matrixA*resultX;
        Console.WriteLine(@"4. Multiply coefficient matrix 'A' by result vector 'x'");
        Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
        Console.WriteLine();
        Console.Read();
    }
}
}

Direct Example:

namespace Examples.LinearAlgebraExamples
{
/// <summary>
/// Direct solvers (using matrix decompositions)
/// </summary>
/// <seealso cref="http://en.wikipedia.org/wiki/Numerical_analysis#Direct_and_iterative_methods"/>
public class DirectSolvers : IExample
{
    /// <summary>
    /// Gets the name of this example
    /// </summary>
    public string Name
    {
        get
        {
            return "Direct solvers";
        }
    }

    /// <summary>
    /// Gets the description of this example
    /// </summary>
    public string Description
    {
        get
        {
            return "Solve linear equations using matrix decompositions";
        }
    }

    /// <summary>
    /// Run example
    /// </summary>
    public void Run()
    {
        // Format matrix output to console
        var formatProvider = (CultureInfo) CultureInfo.InvariantCulture.Clone();
        formatProvider.TextInfo.ListSeparator = " ";

        // Solve next system of linear equations (Ax=b):
        // 5*x + 2*y - 4*z = -7
        // 3*x - 7*y + 6*z = 38
        // 4*x + 1*y + 5*z = 43

         matrixA = DenseMatrix.OfArray(new[,] { { 20000, 0, 0, -20000, 0, 0, 0, 0, 0 }, { 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 }, 
                                                { 0, 2000,  8000,  0,  -2000, 4000, 0, 0, 0 }, { -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 },
                                                {0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 }, { 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 },
                                                { 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 }, { 0, 0, 0, 0, -20000, 0, 0, 20000, 0 },
                                                 {0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 }});

        Console.WriteLine(@"Matrix 'A' with coefficients");
        Console.WriteLine(matrixA.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // Create vector "b" with the constant terms.
        double[] loadVector = { 0, 0, 0, 5000, 0, 0, 0, 0, 0 };
        var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);
        Console.WriteLine(@"Vector 'b' with the constant terms");
        Console.WriteLine(vectorB.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 1. Solve linear equations using LU decomposition
        var resultX = matrixA.LU().Solve(vectorB);
        Console.WriteLine(@"1. Solution using LU decomposition");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 2. Solve linear equations using QR decomposition
        resultX = matrixA.QR().Solve(vectorB);
        Console.WriteLine(@"2. Solution using QR decomposition");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 3. Solve linear equations using SVD decomposition
        matrixA.Svd().Solve(vectorB, resultX);
        Console.WriteLine(@"3. Solution using SVD decomposition");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 4. Solve linear equations using Gram-Shmidt decomposition
        matrixA.GramSchmidt().Solve(vectorB, resultX);
        Console.WriteLine(@"4. Solution using Gram-Shmidt decomposition");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 5. Verify result. Multiply coefficient matrix "A" by result vector "x"
        var reconstructVecorB = matrixA*resultX;
        Console.WriteLine(@"5. Multiply coefficient matrix 'A' by result vector 'x'");
        Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // To use Cholesky or Eigenvalue decomposition coefficient matrix must be 
        // symmetric (for Evd and Cholesky) and positive definite (for Cholesky)
        // Multipy matrix "A" by its transpose - the result will be symmetric and positive definite matrix
        var newMatrixA = matrixA.TransposeAndMultiply(matrixA);
        Console.WriteLine(@"Symmetric positive definite matrix");
        Console.WriteLine(newMatrixA.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 6. Solve linear equations using Cholesky decomposition
        newMatrixA.Cholesky().Solve(vectorB, resultX);
        Console.WriteLine(@"6. Solution using Cholesky decomposition");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 7. Solve linear equations using eigen value decomposition
        newMatrixA.Evd().Solve(vectorB, resultX);
        Console.WriteLine(@"7. Solution using eigen value decomposition");
        Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
        Console.WriteLine();

        // 8. Verify result. Multiply new coefficient matrix "A" by result vector "x"
        reconstructVecorB = newMatrixA*resultX;
        Console.WriteLine(@"8. Multiply new coefficient matrix 'A' by result vector 'x'");
        Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));
        Console.WriteLine();
        Console.Read();
    }
}
}

The numbers from the example problem may very well be wrong, but I need to be sure that I am using Math.NET correctly before proceeding. Am I using these solver examples the way they were meant to be used? Is there anything else I can try that these examples do not cover?

The Finite Element Analysis Example Problem (p.8, Example 1):

They seem to have messed up on the units somewhere, so in order to get my matrix to match there's, we had to use the following input:

Member  A (mm^2)    E (N/mm^2)  I (mm^4)    L (mm)
AB     600000000    0.0002      60000000      6
BC     600000000    0.0002      60000000      6

Also note that they have eliminated some rows and columns that should naturally go away in the course of the calculation. These rows and columns are still present in the matrix I am using

Nabala

The answer I get with x = inv(A)*b is x={-757.5757,0,0,-757.5757,0,0,-757.5757,0,0} and A is very ill conditioned. Please check your math first. – Carillon

Answer 1

The last two sentences in your question are the source of your problem:

Also note that they have eliminated some rows and columns that should naturally go away in the course of the calculation. These rows and columns are still present in the matrix I am using.

In your example problem, you have joints that are fixed against movement in certain directions (called boundary conditions). Sometimes when doing finite element analysis, if you don't remove the appropriate rows and columns from your stiffness matrix and load matrix according to these boundary conditions, you will end up with a system that cannot be solved, which is the case here.

Try the DirectSolver again with:

var matrixA = DenseMatrix.OfArray(new[,] { {20000,  0,  -20000, 0,  0}, {0, 8000    ,0, -2000   ,4000},
                                                   {-20000, 0,  20666.667   ,0, 2000}, {0,  -2000   ,0, 20666.67,   -2000},
                                                    {0, 4000    ,2000   ,-2000, 16000}});

and

double[] loadVector = { 0, 0, 5, 0, 0 };
var vectorB = MathNet.Numerics.LinearAlgebra.Vector<double>.Build.Dense(loadVector);

To answer your question, yes you are using the methods correctly, but you are solving the wrong system. Fix your input and you should get the output you are looking for.

I should also point out that the reason I suggest using the Direct Solver Example is because it seems like you are looking for an exact solution. The Iterative Solvers save computation time by merely approximating a solution.

Sugarcoat

Answer 2

Can Math.NET solve any matrix?

No it can't. Specifically, it can't solve a system of equations that has no solution, and nor can any other solver.

In this case your matrix A is singular, i.e. it doesn't have an inverse. This means that your system of equations either has no solution, i.e. it is inconsistent, or it has infinite solutions (see section 6.5 in Introduction to Numerical Methods for examples). A singular matrix has a zero determinant. You can see this in mathnet as follows using the Determinant method:

Console.WriteLine("Determinant {0}", matrixA.Determinant());

This gives

Determinant 0

A condition for A being singular is when a linear combination of its rows (or columns) is zero. For example here the sum of the 2nd, 5th and 8th rows is zero. These aren't the only rows that add together to give zero. (You'll see another example later. In fact there are three different ways of doing it, which technically means that this 9x9 matrix is "rank 6" rather than "rank 9".).

Remember that all you're doing when you're trying to solve Ax=b is to solve a set of simultaneous equations. In two dimensions you might have a system such as

A = [1 1   b = [1 
     2 2],      2]

and solving this is equivalent to finding x0 and x1 such that

  x0 +   x1 = 1
2*x0 + 2*x1 = 2

Here there are infinite solutions satisfying x1 = 1 - x0, i.e. along the line x0 + x1 = 1. Alternatively for

A = [1 1   b = [1 
     1 1],      2]

which is equivalent to

  x0 +  x1 = 1
  x0 +  x1 = 2

there is clearly no solution because we can subtract the first equation from the second to get 0 = 1!

In your case the 1st, 4th, and 7th equations are

 20000*x0 -20000               *x3                                          = 0
-20000*x0 +20666.66666666666663*x3 +2000*x5 -666.66666666666663*x6 +2000*x8 = 5
            -666.66666666666663*x3 -2000*x5 +666.66666666666663*x6 -2000*x8 = 0

Adding these together gives 0=5, and hence your system has no solution.

It's easiest to explore matrices in an interactive environment like Matlab or R. As Python is available in Visual Studio and it provides a Matlab-like environment via numpy I've demonstrated the above with some code in Python. I would recommend Python tools for visual studio, which I've used successfully in both Visual Studio 2012 and 2013.

# numpy is a Matlab-like environment for linear algebra in Python
import numpy as np

# matrix A
A = np.matrix ([
    [ 20000, 0, 0, -20000, 0, 0, 0, 0, 0 ],
    [ 0, 666.66666666666663, 2000, 0, -666.66666666666663, 2000, 0, 0, 0 ], 
    [ 0, 2000,  8000,  0,  -2000, 4000, 0, 0, 0 ],
    [ -20000, 0, 0, 20666.66666666666663, 0, 2000, -666.66666666666663, 0, 2000 ],
    [ 0, -666.66666666666663, -2000, 0, 20666.66666666666663, -2000, 0, -20000, 0 ],
    [ 0, 2000, 4000, 2000, -2000, 16000, -2000, 0, 4000 ],
    [ 0, 0, 0, -666.66666666666663, 0, -2000, 666.66666666666663, 0, -2000 ],
    [ 0, 0, 0, 0, -20000, 0, 0, 20000, 0 ],
    [ 0, 0, 0, 2000, 0, 4000, -2000, 0, 7999.9999999999991 ]])

# vector b
b = np.array([0, 0, 0, 5, 0, 0, 0, 0, 0])
b.shape = (9,1)

# attempt to solve Ax=b
np.linalg.solve(A,b)

This fails with an informative error message: LinAlgError: Singular matrix. You can see that A is singular by, e.g., showing that the sum of the 2nd, 5th and 8th rows is zero

A[1,]+A[4,]+A[7,]

noting that rows are zero-indexed.

To demonstrate that the 1st, 4th, and 7th equations lead to 0=5 form the augmented matrix by appending the columnar vector b onto A, and then adding the corresponding (0-indexed) rows together

Aaug = np.append(A,b,1)

Aaug[0,] + Aaug[3,] + Aaug[6,]

Finally even if your matrix isn't singular you can still have a numerically unstable problem: in this case the problem is known as ill-conditioned. Check the condition number of the matrix to see how to do this (wikipedia, np.linalg.cond(A), matrixA.ConditionNumber()).

Preterite

This answer is a very, very patient way of saying "you should actually learn linear algebra and probably some numerical analysis" as well as teaching it. Bravo. – Typecase

Why would you give an answer using python and numpy to a question clearly tagged .net mathnet? – Bounteous

@Bounteous thanks for commenting. Reason is that this is more general than any one programming language -- it can't be solved in any of them. I'm more experienced in C# than Python, but Python/numpy was the easiest way to explore the problem. I've been thinking today that I'd like to add to add one or two mathnet specific points. In particular I'm interested in calculating the determinant (0 for a singular matrix) and also, ideally, the condition number. But it took me a while to address the math last night and I ran out of time-- please check back after the weekend. – Preterite

Answer 3

No it cannot solve singular matrices. But neither does any other code out there as there is no solution here.

For your particular case the matrix A posted is singular. The size is 9×9 but the rank is 6. MATLAB reports a condition of 1.9e17. So please check your stiffness matrix composition before you can expect a reasonable answer. Maybe you need to normalize your matrix, i.e. extract the E I coefficients to bring the numbers down from 1e5 to near 1 which is more numerically palatable.

FYI

If you don't like Math.NET, or you want to validate code use pure c#. Read this MSDN Magazine article by James McCaffrey and use the code listed.

var A = new [,] { ... };
var b = new [] { ... };
var x = LU.SustemSolve(A,b);

Carillon