MathNet矩阵分解

矩阵分解是解线性方程组的一个方法，比较容易。

有LU,QR,SVD等等方法。在c#中用MathNet操作起来也是比较简单。

https://blog.csdn.net/c914620529/article/details/50393223/

https://www.cnblogs.com/asxinyu/p/4285245.html

https://www.cnblogs.com/asxinyu/p/Dotnet_Opensource_MathNet_Basic_Statistics_10.html

程序如下：

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using MathNet.Numerics.LinearAlgebra.Double;

namespace SimpleCalcModule
{
    static public class MatrixCalc
    {
        //LU
        public static bool calc1(int n, double[,] A, double[] b, ref double[] X)
        {
            if (n < 3)
                return false;
            if (b.Length != n)
                return false;
            if (A.GetLength(0) != n && A.GetLength(1) != n)
                return false;
            try
            {
                var matrixA = new DenseMatrix(n, n);
                var vectorB = new DenseVector(n);
                for (int i = 0; i < n; i++)
                {
                    vectorB[i] = b[i];
                    X[i] = .0;
                    for (int j = 0; j < n; j++)
                    {
                        matrixA[i, j] = A[i, j];
                    }
                }
                //1.LU
                var resultX = (DenseVector)matrixA.LU().Solve(vectorB);
                X = resultX.ToArray();
                if (X.Length == n)
                    return true;
                else
                    return false;
            }
            catch
            {
                return false;
            }
            
        }
        //QR
        public static bool calc2(int n, double[,] A, double[] b, ref double[] X)
        {
            if (n < 3)
                return false;
            if (b.Length != n)
                return false;
            if (A.GetLength(0) != n && A.GetLength(1) != n)
                return false;
            try
            {
                var matrixA = new DenseMatrix(n, n);
                var vectorB = new DenseVector(n);
                for (int i = 0; i < n; i++)
                {
                    vectorB[i] = b[i];
                    X[i] = .0;
                    for (int j = 0; j < n; j++)
                    {
                        matrixA[i, j] = A[i, j];
                    }
                }
                //1.LU
                var resultX = (DenseVector)matrixA.QR().Solve(vectorB);
                X = resultX.ToArray();
                if (X.Length == n)
                    return true;
                else
                    return false;
            }
            catch
            {
                return false;
            }
        }
        //SVD
        public static bool calc3(int n, double[,] A, double[] b, ref double[] X)
        {
            if (n < 3)
                return false;
            if (b.Length != n)
                return false;
            if (A.GetLength(0) != n && A.GetLength(1) != n)
                return false;
            try
            {
                var matrixA = new DenseMatrix(n, n);
                var vectorB = new DenseVector(n);
                for (int i = 0; i < n; i++)
                {
                    vectorB[i] = b[i];
                    X[i] = .0;
                    for (int j = 0; j < n; j++)
                    {
                        matrixA[i, j] = A[i, j];
                    }
                }
                //1.LU
                var resultX = (DenseVector)matrixA.Svd().Solve(vectorB);
                X = resultX.ToArray();
                if (X.Length == n)
                    return true;
                else
                    return false;
            }
            catch
            {
                return false;
            }
        }
        //Gram-Shmidt
        public static bool calc4(int n, double[,] A, double[] b, ref double[] X)
        {
            if (n < 3)
                return false;
            if (b.Length != n)
                return false;
            if (A.GetLength(0) != n && A.GetLength(1) != n)
                return false;
            try
            {
                var matrixA = new DenseMatrix(n, n);
                var vectorB = new DenseVector(n);
                for (int i = 0; i < n; i++)
                {
                    vectorB[i] = b[i];
                    X[i] = .0;
                    for (int j = 0; j < n; j++)
                    {
                        matrixA[i, j] = A[i, j];
                    }
                }
                //1.LU
                var resultX = (DenseVector)matrixA.GramSchmidt().Solve(vectorB);
                X = resultX.ToArray();
                if (X.Length == n)
                    return true;
                else
                    return false;
            }
            catch
            {
                return false;
            }
        }
    }
}

测试函数
 

if (commOperate == 99)
            {
                //double[,] a = { { 8.1, 2.3, -1.5, 6.1 }, { 0.5, -6.23, 0.87, 2.3 }, { 2.5, 1.5, 10.2, 1.8 } };
                //double[,] a = { { 2, -1, 3, 1 }, { 4, 2, 5, 4 }, { 1, 2, 0, 7 } };

                double[,] A = { { .0, .0, 1.0 }, { 50.0, .0, 1.0 }, { .0, 50.0, 1.0 } };
                int n = A.GetLength(0);//数组a的第一维长度，即行数,3
                double[] b1 = new double[n];//存放解的数组，初始值为0
                double[] b2 = new double[n];//存放解的数组，初始值为0
                b1[0] = 37.5;
                b1[1] = 62.5;
                b1[2] = 37.5;
                double[] X1 = new double[n];
                double[] X2 = new double[n];
                double[] X3 = new double[n];
                double[] X4 = new double[n];
                b2[0] = 38.75;
                b2[1] = 38.75;
                b2[2] = 63.75;
                double[] X5 = new double[n];
                double[] X6 = new double[n];
                double[] X7 = new double[n];
                double[] X8 = new double[n];
                bool rst1 = SimpleCalcModule.MatrixCalc.calc1(n, A, b1, ref X1);
                bool rst2 = SimpleCalcModule.MatrixCalc.calc2(n, A, b1, ref X2);
                bool rst3 = SimpleCalcModule.MatrixCalc.calc3(n, A, b1, ref X3);
                bool rst4 = SimpleCalcModule.MatrixCalc.calc4(n, A, b1, ref X4);
                bool rst5 = SimpleCalcModule.MatrixCalc.calc1(n, A, b2, ref X5);
                bool rst6 = SimpleCalcModule.MatrixCalc.calc2(n, A, b2, ref X6);
                bool rst7 = SimpleCalcModule.MatrixCalc.calc3(n, A, b2, ref X7);
                bool rst8 = SimpleCalcModule.MatrixCalc.calc4(n, A, b2, ref X8);
                double x4 = 50;
                double y4 = 50;
                double Tx1 = X1[0] * x4 + X1[1] * y4 + X1[2];
                double Tx2 = X2[0] * x4 + X2[1] * y4 + X2[2];
                double Tx3 = X3[0] * x4 + X3[1] * y4 + X3[2];
                double Tx4 = X4[0] * x4 + X4[1] * y4 + X4[2];

                double Ty1 = X5[0] * x4 + X5[1] * y4 + X5[2];
                double Ty2 = X6[0] * x4 + X6[1] * y4 + X6[2];
                double Ty3 = X7[0] * x4 + X7[1] * y4 + X7[2];
                double Ty4 = X8[0] * x4 + X8[1] * y4 + X8[2];

            }

多谢，亲爱的美美。