.net - Math.NET 可以解决任何矩阵吗？

Question

我正在尝试使用 Math.NET 来解决以下系统:

系数矩阵 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 }});

结果向量 b:

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

我从有限元分析示例问题中提取了这些数字，因此我期望基于该示例的答案是:

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

但是，Math.NET 和 online Matrix Calculator我发现大多给我零(来自迭代求解器)、NaN 或大的不正确数字(来自直接求解器)作为解决方案。

在 Math.NET 中，我尝试将矩阵插入提供的示例中，包括:

迭代示例:

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

直接示例:

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

示例问题中的数字很可能是错误的，但在继续之前，我需要确保我正确使用了 Math.NET。我是否按照预期的方式使用这些求解器示例？还有什么我可以尝试但这些例子没有涵盖的吗？

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

他们似乎在某个地方搞砸了单位，所以为了让我的矩阵匹配那里，我们不得不使用以下输入:

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

还要注意，他们已经消除了一些在计算过程中应该自然消失的行和列。这些行和列仍然存在于我使用的矩阵中

Answer 1

Can Math.NET solve any matrix?

不，不能。具体来说，它无法求解无解的方程组，任何其他求解器也无法求解。

在这种情况下，您的矩阵 A是单数，即它没有逆。这意味着您的方程组要么没有解，即不一致，要么有无限解(参见 Introduction to Numerical Methods 中的第 6.5 节)。奇异矩阵的行列式为零。您可以使用 Determinant 在 mathnet 中看到这一点，如下所示方法:

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

这给

Determinant 0

A 是奇异的一个条件是其行(或列)的线性组合为零。例如这里第 2、第 5 和第 8 行的总和为零。这些不是唯一相加得到零的行。 (稍后您将看到另一个示例。实际上有三种不同的方法，这在技术上意味着这个 9x9 矩阵是“等级 6”而不是“等级 9”。)。

请记住，当您试图解决问题时所做的一切 Ax=b是求解一组联立方程组。在二维中，您可能有一个系统，例如

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

解决这个问题相当于找到 x0和 x1以至于

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

这里有无穷解满足 x1 = 1 - x0 ，即沿线 x0 + x1 = 1 .或者为

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

这相当于

  x0 +  x1 = 1
  x0 +  x1 = 2

显然没有解，因为我们可以从第二个方程中减去第一个方程得到 0 = 1 !

在您的情况下，第 1、第 4 和第 7 个方程是

 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

将这些加在一起得到 0=5 ，因此您的系统没有解决方案。

在像 Matlab 或 R 这样的交互式环境中探索矩阵是最简单的。由于 Python 在 Visual Studio 中可用，并且它通过 numpy 提供了一个类似 Matlab 的环境，我已经用 Python 中的一些代码演示了上面的内容。我会推荐 Python tools for visual studio ，我在 Visual Studio 2012 和 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)

这失败并显示一条信息性错误消息: LinAlgError: Singular matrix .你可以看到 A是单数，例如，显示第 2、第 5 和第 8 行的总和为零

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

注意行是零索引的。

为了证明第 1、第 4 和第 7 个方程导致 0=5通过附加列向量形成增广矩阵 b到 A ，然后将相应的(0 索引)行相加

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

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

最后，即使您的矩阵不是奇异的，您仍然会遇到数值不稳定的问题:在这种情况下，该问题被称为病态。检查矩阵的条件编号以了解如何执行此操作( wikipedia 、 np.linalg.cond(A) 、 matrixA.ConditionNumber() )。