Фундаментальное решение однородных линейных уравнений: Ax = 0 с Det (A) = 0 с MathNet - PullRequest
Новая
       55

Фундаментальное решение однородных линейных уравнений: Ax = 0 с Det (A) = 0 с MathNet

0 голосов
/ 15 ноября 2018

Я пытаюсь решить систему однородных линейных уравнений типа 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.

Это решение не работает для моей ситуации (B= 0, Det (A) = 0).

1 Ответ

0 голосов
/ 31 июля 2019

Для решения линейных уравнений вы можете попробовать 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

Хотя это может работать как здесь .

В любом случае 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
...