matlab中欠定方程組超定方程組_一篇文章入門大規模線性方程組求解

前面介紹過主要的線性方程組求解庫，參考附錄。求解大規模線性方程組是仿真軟體求解器的底層技術，求解器時間基本都消耗在方程組求解上。線性方程組的解法比較成熟，方法也有很多，而且不同的方法對應不同類型方程組，是以在方法選擇上實際很講究。

商業軟體通常将方法封裝起來，使用者包括開發人員都接觸不到線性方程組求解方法。商業軟體内部一般會根據求解規模，求解類型等選擇合适的線性方程組求解方法。

有些商業軟體開放了部分接口供使用者選擇；開源軟體比如OpenFOAM，以及使用開源軟體的平台Simscale等提供了許多選項供使用者選擇。

本文簡單介紹下線性方程組的常用解法。通常将線性方程組表示為：

A*x=b

A為已知N*N的矩陣，通常稱為剛度矩陣（剛度是力學中的概念，電磁，熱等也習慣性這麼稱呼），b為已知向量，x為待求向量。解線性方程組的操作基本圍繞矩陣A展開。

首先介紹一些相關術語：

1.矩陣條件數

條件數是一個表征矩陣穩定特性的标志，條件數越大，說明矩陣越不穩定，即當矩陣中資料出現微小變化時，x結果變化非常大。Matlab中可直接使用指令cond(A)檢視矩陣條件數。

2.滿秩矩陣

用初等行變換将矩陣A化為階梯形矩陣, 矩陣中非零行的個數就定義為這個矩陣的秩，記做r(A)

即如果矩陣A為N*N，r(A)=N，則矩陣A為滿秩矩陣

在邊界元，矩量法等計算方法中，最終形成的A為滿秩矩陣。

3.稀疏矩陣

矩陣中絕大部分元素都是0

4.對稱矩陣

矩陣的上三角和下三角關于對角線對稱

5.求解線性方程組直接法

先求出矩陣A的逆矩陣，再乘以向量b

6.求解線性方程組疊代法

簡單講，就是給出一個初始解x'，帶入原方程中，每次評價差距逐漸修正，直到最終A*x'-b 接近0為止。

7. OOC（out of core）

對于大型方程組，記憶體通常無法裝下，在求解過程中需要将部分資料寫到硬碟上，再次使用的時候再讀回來。

常用直接法：

直接法的本質是要計算出A的逆矩陣，通常在求解小規模和特征值問題是可以考慮使用直接法。

1. 高斯消去法/Doolittle /三角分解法/追趕法

這是最基本的方法，時間複雜度和空間複雜度都是N的三次方，軟體一般都不會使用。

2. 矩陣分解相關

2. 1.LU分解法

LU分解就是将矩陣分解成機關下三角矩陣L和上三角矩陣U，本質上仍然屬于高斯消去法

2.2. Cholesky分解

Cholesky 分解是把一個對稱正定的矩陣表示成一個下三角矩陣L和其轉置的乘積的分解。如果矩陣是正定的，使用 Cholesky分解會比LU分解更加高效。

2.3. LDLT分解法

Cholesky 分解法的改進

2.4.QR分解

QR分解是把矩陣分解成一個正交矩陣與一個上三角矩陣的積。

2.5. Schur分解

2.6. SVD/GSVD

奇異值分解/一般奇異值分解

在求解大規模線性方程組中，一般不會使用直接法求解，但在使用疊代法過程中需要使用直接法中的方法加工資料。

間接法的本質是疊代，不同方法的差別在于如何選取初始值以及疊代方法。

1. 牛頓疊代，Jacobi(雅克比)疊代，Gauss-Seidel疊代

線性疊代方法，一般情況下收斂太慢，大規模方程組不推薦使用。

2. JOR -- Jacobi Over Relaxation

Jacobi 方法加入松弛因子

3. Lanczos方法

适用于稀疏矩陣特征值問題

4. 共轭梯度方法以及改進方法

共轭梯度法

Conjugate gradient method--CG

雙共轭梯度法

Bi-conjugate gradient method--BCG

穩定雙共轭梯度法

Bi-conjugate gradient stabilized method--GCGS

預條件共轭梯度法

Preconditioned gradient stabilized method--PCG

相比牛頓和Jacobi，通過優化使其共轭的求解向量和方向，加速了求解性能。

共轭梯度法的思想就是找到N個兩兩共轭的方向，每次沿着一個方向優化得到該方向上的極小值，後面再沿其它方向求極小值的時候，不會影響前面已經得到的沿哪些方向上的極小值，是以理論上對n個方向都求出極小值就得到了N維問題的極小值。

5. GMRAS廣義最小殘量

(Generalized Minimal Residual Algorithm)

非對稱系統的線性方程組的數值解疊代法，該方法與最小殘量的 Krylov 子空間中向量來逼近解，是求解大規模線性方程組的常用算法之一,具有收斂速度快、穩定性好等優點。

改進的GMRAS：

SGMRAS Simpler GMRAS

PGMRAS Preconditioned GMRAS

6. 快速多級（Fast Multi pole Method）

當矩陣為滿秩矩陣時，傳統計算方法資源需求和求解規模呈指數級上升，大規模的系統求解需要使用快速多級方法。本技術部落格和公衆号中有詳細介紹。

改進方法：多層快速多級

7. Successive Over Relaxation(SOR)連續松弛

連續過度松弛（SOR）方法是高斯-賽德爾（Gauss-Seidel）方法的一種變體，用于求解線性方程組，進而可以更快地收斂。任何緩慢收斂的疊代過程都可以使用類似的方法。

改進的SOR：

Accelerated Over relaxation (AOR)過松弛

Preconditioned AOR -- PAOR 預條件過松弛

Quasi AOR -- QAOR準過松弛

8. Krylov子空間疊代法

将矩陣A分解成多個子矩陣，并用一系列線性表達式組合表示。将整個系統降維，利于并行計算，是大規模線性方程組有效的一種解法。其中涉及到了Lanczos方法和Arnoldi方法。

按照矩陣A的特點，我們可以做如下分類：

1. 是否是滿秩矩陣；

2. 是否是稀疏矩陣；

3. 是否是對稱矩陣；

4. 是否是病态矩陣；

5. 是否是正定矩陣

6. 矩陣規模（矩陣中N的大小）

在選擇求解方法的時候，需要考慮到方程組矩陣以上特點。根據作者的經驗，方程組的求解性能高度依賴矩陣特征，矩陣規模和硬體。

問題：

求解線性方程組 Ax=b.

Matlab 中提供了方程的求解方法：

x=Ab （ 反除法符号）

最簡單的問題，一群兔子和鴨子，一共有40條腿，20張嘴，問兔子鴨子各多少隻，兩個變量，兩個方程組成一個簡單的線性方程組：

4*x+2*y =40

x+y =20

在Matlab中輸入:

>>A =

4 2

1 1

>> b=[40, 15]'

b =

40

15

>> x=Ab

x =

5

10

當方程求解規模達到1000*1000的時候，Matlab還能應付，速度也不錯。當規模達到10w*10w的時候，Matlab會運作相當長時間，然後彈出一個錯誤。這是因為對于10w*10w的矩陣，用原始的資料結構，性能會下降很快。由有限元計算的矩陣特征可知，矩陣是資料大部分為0的稀疏矩陣，Matlab提供了稀疏矩陣功能，在使用x=Ab 之前 使用如下指令：

x=sparse(x);

10萬*10萬的方程組1秒之内求出結果。

對于大型線性方程組的求解，使用Matlab做驗證是個不錯的選擇

回到主題:

1. BLAS(Basic Linear Algebra Subprograms)

blas 是許多數值計算軟體庫的核心，用Fortran語言寫成，庫一共有三個level，第一level包含了vector計算的一些function；第二個level則是矩陣與vector計算；第三個level是矩陣與矩陣的計算。BLAS接口稍微複雜，一般很少直接使用。

開源

2. Intel MKL

商用是不錯的選擇。

Intel數學核心函數庫（MKL）是一套高度優化、線程安全的數學例程、函數，面向高性能的工程、科學與财務應用。英特爾 MKL 的叢集版本包括 ScaLAPACK 與分布式記憶體快速傅立葉轉換，并提供了線性代數 (BLAS、LAPACK 和Sparse Solver)、快速傅立葉轉換、矢量數學 (Vector Math) 與随機号碼生成器支援。

主要包括：

LAPACK (線形代數工具linear algebra package)

DFTs (離散傅立葉變換 Discrete Fourier transforms)

VML (矢量數學庫Vector Math Library)

VSL (矢量統計庫Vector Statistical Library)

MKL的主要功能

2.1）BLAS 和 LAPACK

在英特爾處理器中部署經過高度優化的基本線性代數例程BLAS（Basic Linear Algebra Subroutines）和 線性代數包LAPACK（Linear Algebra Package） 例程，它們提供的性能改善十分顯著。

2.2）ScaLAPACK

ScaLAPACK是一個并行計算軟體包，适用于分布存儲的MIMD并行機。ScaLAPACK提供若幹線性代數求解功能，具有高效、可移植、可伸縮、高可靠性的特點，利用它的求解庫可以開發出基于線性代數運算的并行應用程式。ScaLAPACK 的英特爾MKL 實施可提供顯著的性能改進，遠遠超出标準 NETLIB 實施所能達到的程度。

2.3）PARDISO稀疏矩陣解算器

利用 PARDISO 直接稀疏矩陣解算器解算大型的稀疏線性方程組，該解算器獲得了巴塞爾大學的授權，是一款易于使用、具備線程安全性、高性能的記憶體高效型軟體庫。英特爾? MKL 還包含共轭梯度解算器和 FGMRES 疊代稀疏矩陣解算器。

2.4）快速傅立葉變換 (FFT)

充分利用帶有易于使用的新型 C/Fortran 接口的多元 FFT 子程式（從 1 維至 7 維）。英特爾? MKL 支援采用相同 API 的分布式記憶體叢集，支援将工作負載輕松地分布到大量處理器上，進而實作大幅的性能提升。此外，英特爾 MKL 還提供了一系列 C 語言例程（“wrapper”），這些例程可模拟 FFTW 2.x 和 3.0 接口，進而支援目前的 FFTW 使用者将英特爾 MKL 內建到現有應用中。

2.5）矢量數學庫（VML）

矢量數學庫（Vector Math Library）借助計算密集型核心數學函數（幂函數、三角函數、指數函數、雙曲函數、對數函數等）的矢量實施顯著提升應用速度。

2.6）矢量統計庫—随機數生成器（VSL）

利用矢量統計庫（Vector Statistical Library）随機數生成器加速模拟，進而實作遠遠高于标量随機數生成器的系統性能提升

商用

3. Eigen

首頁：

http://eigen.tuxfamily.org/index.php?title=Main_Page

Eigen 是一個輕量級的庫，使用隻需頭檔案，接口比較齊全，文檔也很全面，一般的矩陣操作運算都能滿足。提供了DenseMatrix 和 SparseMatrix，對于稀疏矩陣的效率也不錯。

開源

4. Spoolse

早期很多人使用，性能不錯，使用簡單。維護的不是很好，文檔使用ps格式。

http://www.netlib.org/linalg/spooles/spooles.2.2.html

開源

5. Lapack

使用最廣泛的線性代數包，底層調用 BLAS，能求解 AX=b, 矩陣分解、求逆，求矩陣特征值、奇異值等。文檔齊全，使用靈活，穩定性好，用作研究和商業開發都是不錯的選擇。提供SVN源碼下載下傳

http://www.netlib.org/lapack/

開源

6. Pardiso

PARDISO是為解決大 稀疏對稱和非對稱線性系統的方程的軟體包，如下特點：

線程安全，,高性能的, 和記憶體使用率低，支援共享記憶體和記憶體多處理器。

http://www.pardiso-project.org/

商用，學術用免費

7. Mumps

使用多波前法的稀疏矩陣求解庫，支援并行計算。

http://mumps.enseeiht.fr/

開源

8. PETSc

PETSc是一個高大上的科學計算庫，求解線性方程是其中一個功能，支援記憶體共享并行計算機，支援多線程，GPU加速等，支援求解稀疏矩陣

http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html

開源

9. SuperLU

SuperLU是一個求解大規模，稀疏，非對稱系統對線性方程組的通用庫,c語言編寫，可供Fortran和C調用。最新版本為4.3 （2014/10/1）支援并行計算和分布式計算。

http://www.cs.berkeley.edu/~demmel/SuperLU.html

開源

10. Umfpack

UMFPACK是 用來求解不對稱稀疏線性系統軟體包, Ax = b,使用非對稱多波方法，SuitSparse有此算法

使用比較簡單，Windows下需要自己編譯。

http://www.cise.ufl.edu/research/sparse/umfpack/

開源

11. TAUCS

求解稀疏矩陣線性方程軟體包，最新版本支援多線程

http://www.tau.ac.il/~stoledo/taucs/

開源

12. SuitSparse

SuiteSparse是一組C、Fortran和MATLAB函數集，用來生成空間稀疏矩陣資料。在SuiteSparse中幾何多種稀疏矩陣的處理方法，包括矩陣的LU分解，QR分解，Cholesky分解，提供了解非線性方程組、實作最小二乘法等多種函數代碼

http://faculty.cse.tamu.edu/davis/suitesparse.html

開源

13. Sparselib++

用C++寫的 稀疏矩陣庫，可以和 IML++一起使用做線性方程組的疊代求解。

http://math.nist.gov/sparselib++/

開源

14. Trilinos

Trilinos也是個高大上的東西，求解線性方程隻是其中一部分功能

http://trilinos.org/capability-areas/#ScalableLinearAlgebra

開源

15.IML++

使用C++模闆庫 和疊代方法 求解對稱，非對稱矩陣庫

http://math.nist.gov/iml++/

開源

使用推薦：

1. 簡單快速上手 Eigen

2. 商業使用 Intel MKL

3. 研究使用 Lapack，PETSc，SuitSparse