开源Math.NET基础数学类库使用(06)数值分析之线性方程组直接求解

　　在前几篇关于Math.NET的博客中(见上面链接)，主要是介绍了Math.NET中主要的数值功能，并进行了简单的矩阵向量计算例子，接着使用Math.NET的矩阵等对象，对3种常用的矩阵数据交换格式的读写。一方面可以了解Math.NET的使用，另一方面以后也可以直接读取和保存数据为这两种格式，给大家的科研或者工作带来便利。接下来的文章将继续对Math.NET的功能进行讲解和演示，并附带一些数学方面的基础知识。毕竟很多人没有精力去研究Math.NET，那我就把我的研究心得一一写出来，方便后来人。

1.数值分析与线性方程

　　数值分析的基本含义与特点：引用

　　数值分析(numerical analysis)是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科，是数学的一个分支，它以数字计算机求解数学问题的理论和方法为研究对象。为计算数学的主体部分。数值分析有如下特点：
1·面向计算机
2·有可靠的理论分析
3·要有好的计算复杂性
4·要有数值实验
5.要对算法进行误差分析　　

　　数值分析的主要内容：插值法，函数逼近，曲线拟和，数值积分，数值微分，解线性方程组的直接方法，解线性方程组的迭代法，非线性方程求根，常微分方程的数值解法。

　　所以我们今天要解决的就是数值分析的一个很小但又最常接触到的部分，方程(组)的求解。方程(组)的求解有2种方法，一种是直接求解，一种是迭代求解。直接方法比较好理解，相当与使用公式进行直接计算，结果也比较精确。而另外一种是迭代方法。从最初的猜测,迭代方法逐次近似形式收敛只在极限的精确解。

　　正如上面所说，方程（组）的形式很多，不同的形式以及实际要求的精度都可以使用不同的方法，而本文主要介绍最简单，也是最常见的线性方程的直接求解方法。

2.Math.NET解线性方程源码分析

　　本文首先对Math.NET中求解线性方程的相关源码进行分析，这样大家碰到了问题，也可以更好的查找源码解决，或者进行扩展，实现自己的一些特殊需求。Math.NET中MathNet.Numerics.LinearAlgebra.Factorization命名空间下有一个泛型接口：ISolver<T>，就是解AX = B类型的线性方程的接口类型。该接口功能很多，看看下面的接口源代码和注释(本人进行了简单的翻译)，就很清楚了。

 1 using System;
 2 
 3 namespace MathNet.Numerics.LinearAlgebra.Factorization
 4 {
 5     /// <summary>形如AX = B的线性方程组求解的接口类型</summary>
 6     /// <typeparam name="T">泛型参数，支持类型有：double, single, <see cref="Complex"/>, 
 7     /// 以及 <see cref="Complex32"/>.</typeparam>
 8     public interface ISolver<T> where T : struct, IEquatable<T>, IFormattable
 9     {
10         /// <summary>求解AX=B的线性方程组</summary>
11         /// <param name="input">右矩阵<c>B</c>.</param>
12         /// <returns>返回求解结果矩阵 <c>X</c>.</returns>
13         Matrix<T> Solve(Matrix<T> input);
14 
15         /// <summary>求解AX=B的线性方程组</summary>
16         /// <param name="input">右矩阵 <c>B</c>.</param>
17         /// <param name="result">求解结果矩阵, <c>X</c>.</param>
18         void Solve(Matrix<T> input, Matrix<T> result);
19 
20         /// <summary>求解AX=b的线性方程组</summary>
21         /// <param name="input">等式的右边向量 <c>b</c>.</param>
22         /// <returns>返回求解结果的左边向量 , <c>x</c>.</returns>
23         Vector<T> Solve(Vector<T> input);
24 
25         /// <summary>求解AX=b的线性方程组</summary>
26         /// <param name="input">等式的右边向量 <c>b</c>.</param>
27         /// <param name="result">求解结果矩阵, <c>x</c>.</param>
28         void Solve(Vector<T> input, Vector<T> result);
29     }
30 }

由于求解线性方程组主要用到了矩阵的分解，Math.NET实现了5种矩阵分解的算法：LU，QR，Svd，Evd，Cholesky。而GramSchmidt是继承QR的，每一个都是实现了ISolver<T>接口，因此就可以直接使用矩阵的分解功能，直接进行线性方程组的求解。为了方便，我们举一个LU的源码例子，简单的看看源码的基本情况：

 1 public abstract class LU<T> : ISolver<T> where T : struct, IEquatable<T>, IFormattable
 2 {
 3     static readonly T One = BuilderInstance<T>.Matrix.One;
 4 
 5     readonly Lazy<Matrix<T>> _lazyL;
 6     readonly Lazy<Matrix<T>> _lazyU;
 7     readonly Lazy<Permutation> _lazyP;
 8 
 9     protected readonly Matrix<T> Factors;
10     protected readonly int[] Pivots;
11 
12     protected LU(Matrix<T> factors, int[] pivots)
13     {
14         Factors = factors;
15         Pivots = pivots;
16 
17         _lazyL = new Lazy<Matrix<T>>(ComputeL);
18         _lazyU = new Lazy<Matrix<T>>(Factors.UpperTriangle);
19         _lazyP = new Lazy<Permutation>(() => Permutation.FromInversions(Pivots));
20     }
21 
22     Matrix<T> ComputeL()
23     {
24         var result = Factors.LowerTriangle();
25         for (var i = 0; i < result.RowCount; i++)
26         {
27             result.At(i, i, One);
28         }
29         return result;
30     }
31 
32     /// <summary>Gets the lower triangular factor.</summary>
33     public Matrix<T> L
34     {
35         get { return _lazyL.Value; }
36     }
37     /// <summary>Gets the upper triangular factor.</summary>
38     public Matrix<T> U
39     {
40         get { return _lazyU.Value; }
41     }
42     /// <summary>Gets the permutation applied to LU factorization.</summary>
43     public Permutation P
44     {
45         get { return _lazyP.Value; }
46     }
47     /// <summary>Gets the determinant of the matrix for which the LU factorization was computed.</summary>
48     public abstract T Determinant { get; }
49     public virtual Matrix<T> Solve(Matrix<T> input)
50     {
51         var x = Matrix<T>.Build.SameAs(input, input.RowCount, input.ColumnCount);
52         Solve(input, x);
53         return x;
54     }    
55     public abstract void Solve(Matrix<T> input, Matrix<T> result);
56     
57     public virtual Vector<T> Solve(Vector<T> input)
58     {
59         var x = Vector<T>.Build.SameAs(input, input.Count);
60         Solve(input, x);
61         return x;
62     }    
63     public abstract void Solve(Vector<T> input, Vector<T> result);    
64     public abstract Matrix<T> Inverse();
65 }

　　大家可能会注意到，上面是抽象类，这和Math.NET的实现是有关的。最终都是实现相应版本的Matrix矩阵，然后实现对应版本的类型的分解方法。下面例子会介绍具体使用，大家有疑问，可以拿着源码和例子，调试一番，知道上面的2个实现过程，就比较简单了

3.Math.NET求解线性方程的实例

　　假设下面是要求解的线性方程组(Ax=b):

5*x + 2*y - 4*z = -7
3*x - 7*y + 6*z = 38
4*x + 1*y + 5*z = 43

测试代码，由于求解的方法很多，只列举了几种，其他的以此类推：

 1 var formatProvider = (CultureInfo) CultureInfo.InvariantCulture.Clone();
 2 formatProvider.TextInfo.ListSeparator = " ";
 3 
 4 //先创建系数矩阵A 
 5 var matrixA = DenseMatrix.OfArray(new[,] {{5.00, 2.00, -4.00}, {3.00, -7.00, 6.00}, {4.00, 1.00, 5.00}});
 6 
 7 //创建向量b
 8 var vectorB = new DenseVector(new[] {-7.0, 38.0, 43.0});
 9 
10 // 1.使用LU分解方法求解
11 var resultX = matrixA.LU().Solve(vectorB);
12 Console.WriteLine(@"1. Solution using LU decomposition");
13 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
14 
15 // 2.使用QR分解方法求解
16 resultX = matrixA.QR().Solve(vectorB);
17 Console.WriteLine(@"2. Solution using QR decomposition");
18 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
19 
20 // 3. 使用SVD分解方法求解
21 matrixA.Svd().Solve(vectorB, resultX);
22 Console.WriteLine(@"3. Solution using SVD decomposition");
23 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
24 
25 // 4.使用Gram-Shmidt分解方法求解
26 matrixA.GramSchmidt().Solve(vectorB, resultX);
27 Console.WriteLine(@"4. Solution using Gram-Shmidt decomposition");
28 Console.WriteLine(resultX.ToString("#0.00\t", formatProvider));
29 
30 // 5.验证结果，就是把结果和A相乘，看看和b是否相等
31 var reconstructVecorB = matrixA*resultX;
32 Console.WriteLine(@"5. Multiply coefficient matrix ‘A‘ by result vector ‘x‘");
33 Console.WriteLine(reconstructVecorB.ToString("#0.00\t", formatProvider));

结果如下：

1. Solution using LU decomposition
DenseVector 3-Double
3.00
1.00
6.00

2. Solution using QR decomposition
DenseVector 3-Double
3.00
1.00
6.00

3. Solution using SVD decomposition
DenseVector 3-Double
3.00
1.00
6.00

4. Solution using Gram-Shmidt decomposition
DenseVector 3-Double
3.00
1.00
6.00

5. Multiply coefficient matrix ‘A‘ by result vector ‘x‘
DenseVector 3-Double
-7.00
38.00
43.00