c# - 齐次线性方程组的基本解 : Ax=0 with Det(A)=0 with MathNet

Question

我正在尝试求解类似 Ax=0 的齐次线性方程组。 这是一个示例矩阵，为简单起见，它已经被简化:

1 2 | 0
3 6 | 0

我希望得到的解决方案至少是[ 2, -1 ]。但根本的解决方案是 [2C; -1C]。您可以看到 Det(A) = 0 和 Rank(A) = 1。你当然知道这样的系统有简单的解决方案 [0,0]。

我正在尝试:

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不适用于我的情况(B = 0，Det(A) = 0)。

Answer 1

要求解线性方程，您可以尝试 Matrix<T>.SolveIterative(Vector<T> input, IIterativeSolver<T> solver, Iterator<T> iterator = null, IPreconditioner<T> preconditioner = null) .

获取可用的求解器:

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>>();

它给出:

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

用你的数据试试它们:

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

像这样:

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

在你的情况下，这不是运气。结果是:

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

虽然它可以像 here 那样工作.

无论如何MathNet非常好，可以轻松构建解决方案。

使用大学的线性代数和图书馆的功能:

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;
    }
}

现在:

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

给予:

DenseVector 2-Double
-2
 1

另一个例子:

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);

给予:

DenseVector 3-Double
-2
 1
 0

DenseVector 3-Double
-1
 0
 1