在前幾篇關于Math.NET的部落格中(見上面連結),主要是介紹了Math.NET中主要的數值功能,并進行了簡單的矩陣向量計算例子,接着使用Math.NET的矩陣等對象,對3種常用的矩陣資料交換格式的讀寫。一方面可以了解Math.NET的使用,另一方面以後也可以直接讀取和儲存資料為這兩種格式,給大家的科研或者工作帶來便利。接下來的文章将繼續對Math.NET的功能進行講解和示範,并附帶一些數學方面的基礎知識。畢竟很多人沒有精力去研究Math.NET,那我就把我的研究心得一一寫出來,友善後來人。
本部落格所有文章分類的總目錄:【總目錄】本部落格博文總目錄-實時更新
開源Math.NET基礎數學類庫使用總目錄:【目錄】開源Math.NET基礎數學類庫使用總目錄
前言
在前幾篇關于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
今天的用法就到此為止,請關注部落格,後續還有更多的分析和使用文章。
4.資源
如果本文資源或者顯示有問題,請參考 本文原文位址:http://www.cnblogs.com/asxinyu/p/4285245.html
.NET資料挖掘與機器學習,作者部落格:
http://www.cnblogs.com/asxinyu
E-mail:[email protected]