Fundamental solution of homogeneous linear equations: Ax=0 with Det(A)=0 with MathNet

Question

I am trying to solve a system of homogeneous linear equations like Ax=0. Here is an example matrix which has already been reduced for simplicity:

1 2 | 0
3 6 | 0

The solution I am hoping to get is at least [ 2, -1 ]. But the fundamental solution is [2C; -1C]. You can see that Det(A) = 0 and Rank(A) = 1. Of course you know that such systems have trivial solution [0,0].

I am trying:

Matrix<double> A = Matrix<double>.Build.DenseOfArray(new double[,]
{
    { 1, 2 },
    { 3, 6 }
});    
Vector<double> B = Vector<double>.Build.Dense(new double[] { 0, 0 });
var result = A.Solve(B); //result = Nan, Nan.

This solution doesn't work for my situation (B = 0, Det(A) = 0).

Answer 1

To solve linear equations you can try Matrix<T>.SolveIterative(Vector<T> input, IIterativeSolver<T> solver, Iterator<T> iterator = null, IPreconditioner<T> preconditioner = null).

To get available solvers:

var solvers = Assembly.GetAssembly(typeof(Matrix<double>))
    .GetTypes()
    .Where(t => t.GetInterfaces().Contains(typeof(IIterativeSolver<double>)) &&
                t.GetConstructors().Any(ctor => ctor.GetParameters().Count() == 0))
    .Select(t => Activator.CreateInstance(t))
    .Cast<IIterativeSolver<double>>();

It gives:

MathNet.Numerics.LinearAlgebra.Double.Solvers.BiCgStab
MathNet.Numerics.LinearAlgebra.Double.Solvers.GpBiCg
MathNet.Numerics.LinearAlgebra.Double.Solvers.MlkBiCgStab
MathNet.Numerics.LinearAlgebra.Double.Solvers.TFQMR

Try them all with you data:

Matrix<double> A = Matrix<double>.Build.DenseOfArray(new double[,]
{
    { 1, 2 },
    { 3, 6 }
});
Vector<double> B = Vector<double>.Build.Dense(new double[] { 0, 0 });

Like this:

foreach (var solver in solvers)
{
    try
    {
        Console.WriteLine(solver);
        Console.WriteLine(A.SolveIterative(B, solver));
    }
    catch (Exception e)
    {
        Console.WriteLine(e.Message);
    }
}

In your case it's no luck. The result is:

MathNet.Numerics.LinearAlgebra.Double.Solvers.BiCgStab
Algorithm experience a numerical break down

MathNet.Numerics.LinearAlgebra.Double.Solvers.GpBiCg
DenseVector 2-Double
NaN
NaN

MathNet.Numerics.LinearAlgebra.Double.Solvers.MlkBiCgStab
Algorithm experience a numerical break down

MathNet.Numerics.LinearAlgebra.Double.Solvers.TFQMR
DenseVector 2-Double
0
0

Though it could work like here.

Anyway MathNet is so good that one can easily build the solution.

Using linear algebra from college and the library functionality:

static class MatrixExtension
{
    public static Vector<double>[] SolveDegenerate(this Matrix<double> matrix, Vector<double> input)
    {
        var augmentedMatrix =
            Matrix<double>.Build.DenseOfColumnVectors(matrix.EnumerateColumns().Append(input));

        if (augmentedMatrix.Rank() != matrix.Rank())
            throw new ArgumentException("Augmented matrix rank does not match coefficient matrix rank.");

        return augmentedMatrix.SolveAugmented();
    }

    private static Vector<double>[] SolveAugmented(this Matrix<double> matrix)
    {
        var rank = matrix.Rank();
        var cut = matrix.CutByRank(rank);
        // [A|R]x[X] = [B]            
        var A = Matrix<double>.Build.DenseOfColumnVectors(cut.EnumerateColumns().Take(rank));
        var R = Matrix<double>.Build.DenseOfColumnVectors(cut.EnumerateColumns().Skip(rank).Take(cut.ColumnCount - rank - 1));
        var B = cut.EnumerateColumns().Last();

        var vectors = Matrix<double>.Build.DenseDiagonal(R.ColumnCount, 1)
            .EnumerateColumns().ToArray();

        return vectors.Select(v => A.Solve(B - R * v))
            .Zip(vectors, (x, v) => Vector<double>.Build.DenseOfEnumerable(x.Concat(v)))
            .ToArray();
    }

    private static Matrix<double> CutByRank(this Matrix<double> matrix, int rank)
    {
        var result = Matrix<double>.Build.DenseOfMatrix(matrix);
        while (result.Rank() < result.RowCount)
            result = result.EnumerateRows()
                           .Select((r, index) => result.RemoveRow(index))
                           .FirstOrDefault(m => m.Rank() == rank);
        return result;
    }
}

And now:

Console.WriteLine(A.SolveDegenerate(B).First());

Gives:

DenseVector 2-Double
-2
 1

Another example:

Matrix<double> A = Matrix<double>.Build.DenseOfArray(new double[,]
{
    { 1, 2, 1 },
    { 3, 6, 3 },
    { 4, 8, 4 }
});
Vector<double> B = Vector<double>.Build.Dense(new double[3]);
foreach (var x in A.SolveDegenerate(B))
    Console.WriteLine(x);

Gives:

DenseVector 3-Double
-2
 1
 0

DenseVector 3-Double
-1
 0
 1