簡介

這里介紹線性系統(tǒng)的解析，如何進行各種分解計算，如LU，QR，SVD，特征值分解等。

簡單線性求解

在一個線性系統(tǒng)，常如下表示，其中A，b分別是一個矩陣，需要求x：
$$Ax=b$$

在Eigen中，我們可以根據(jù)需要的精確性要求，性能要求，從而從好幾種方法中選擇一種來求解這個線性系統(tǒng)。

下面的示例，可以看到一種求解方法。

//matrxi_decomp1.cpp #include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {Matrix3f A;Vector3f b;A << 1,2,3, 4,5,6, 7,8,10;b << 3, 3, 4;cout << "Here is the matrix A:\n" << A << endl;cout << "Here is the vector b:\n" << b << endl;Vector3f x = A.colPivHouseholderQr().solve(b);cout << "The solution is:\n" << x << endl; }

執(zhí)行：

$ g++ -I /usr/local/include/eigen3 matrxi_decomp1.cpp -o matrxi_decomp1 $ $ ./matrxi_decomp1 Here is the matrix A:1 2 34 5 67 8 10 Here is the vector b: 3 3 4 The solution is:-2 0.9999971

這個示例中，使用矩陣的colPivHouseholderQr()方法，其返回一個ColPivHouseholderQR實例對象。因為調(diào)用對象是一個Matrix3f類型，所以也可以如此設計

ColPivHouseholderQR<Matrix3f> dec(A); Vector3f x = dec.solve(b);

正如名字所示，這里ColPivHouseholderQR是一個采用列旋轉(zhuǎn)方法的QR分解。下面Eigen提供的表列出了你可以使用分解類型：

DecompositionMethodRequirements on the matrixSpeed(small-to-medium)Speed(large)Accuracy

PartialPivLU	partialPivLu()	Invertible(可逆)	++	++	+
FullPivLU	fullPivLu()	None	-	- -	+++
HouseholderQR	householderQr()	None	++	++	+
ColPivHouseholderQR	colPivHouseholderQr()	None	+	-	+++
FullPivHouseholderQR	fullPivHouseholderQr()	None	-	- -	+++
CompleteOrthogonalDecomposition	completeOrthogonalDecomposition()	None	+	-	+++
LLT	llt()	Positive definite(正定)	+++	+++	+
LDLT	ldlt()	Positive or negative semidefinite	+++	+	++
BDCSVD	bdcSvd()	None	-	-	+++
JacobiSVD	jacobiSvd()	None	-	- - -	+++

其中的部分分解公式：

• FullPivLU:
  $$A=P^{?1}LU{Q}^{?1}$$

• ColPivHouseholderQR:
  $$\mathbf{A}\mathbf{P}=\mathbf{Q}\mathbf{R}$$

• FullPivHouseholderQR:
  $$\mathbf{P}\mathbf{A}{\mathbf{P}}^{\prime }=\mathbf{Q}\mathbf{R}$$

• LLT:
  $$L{L}^T$$ Cholesky decomposition of a symmetric, positive definite matrix A such that
  $$A=L{L}^{?}=U^{?}U$$ , where L is lower triangular.

• LDLT
  $$LD{L}^T$$
  Cholesky decomposition:
  $$A=P^{T}LD{L}^{?}P$$
  where P is a permutation matrix, L is lower triangular with a unit diagonal and D is a diagonal matrix.

Note: Eigen的測試benchmark

---

如示例所示，所有的分解計算提供的solve()方法，此方法真正進行分解運算。

如果你的運算的矩陣是正定的，如上表中所展示的，更好的選擇是LLT和LDLT分解方法。下面的示例演示了使用方法，其中右邊使用一個普通的矩陣，而不是向量vector。

示例：

#include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() {Matrix2f A, b;A << 2, -1, -1, 3;b << 1, 2, 3, 1;cout << "Here is the matrix A:\n" << A << endl;cout << "Here is the right hand side b:\n" << b << endl;Matrix2f x = A.ldlt().solve(b);cout << "The solution is:\n" << x << endl; }

這里使用了A.ldlt().solve(b)來進行計算，得到結果x。

執(zhí)行結果：

Here is the matrix A:2 -1 -1 3 Here is the right hand side b: 1 2 3 1 The solution is: 1.2 1.4 1.4 0.8

檢查是否有解

計算存在誤差，只有我們知道了一個計算的解的可允許的誤差容限error margin，我們才能判斷一個解是否滿足要求。在Eigen內(nèi)，很容易做這種計算。

#include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {MatrixXd A = MatrixXd::Random(100,100);MatrixXd b = MatrixXd::Random(100,50);MatrixXd x = A.fullPivLu().solve(b);double relative_error = (A*x - b).norm() / b.norm(); // norm() is L2 normcout << "The relative error is:\n" << relative_error << endl; }

執(zhí)行結果：

The relative error is: 3.75098e-13

計算特征值和特征向量

在線性代數(shù)里計算特征值Eigen value和特征向量 EigenVector是非常常見的計算，計算前，要檢查確保你的矩陣是自伴隨self-adjoint矩陣。Eigen里提供了SelfAdjointEigenSolver，很容易地在矩陣對象上使用EigenSolver , ComplexEigenSolver。

特征值和特征向量的計算未必收斂，這種非收斂的情況很罕見，但確實存在。Eigen中的 SelfAdjointEigenSolver.info() 函數(shù)可以用來檢查。

#include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {Matrix2f A;A << 1, 2, 2, 3;cout << "Here is the matrix A:\n" << A << endl;SelfAdjointEigenSolver<Matrix2f> eigensolver(A);if (eigensolver.info() != Success) {cout << "converge error!"<<endl;abort(); }cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;cout << "Here's a matrix whose columns are eigenvectors of A \n"<< "corresponding to these eigenvalues:\n"<< eigensolver.eigenvectors() << endl; }

執(zhí)行結果：

$ g++ -I /usr/local/include/eigen3 matrxi_decomp3.cpp -o matrxi_decomp3 $ ./matrxi_decomp3 Here is the matrix A: 1 2 2 3 The eigenvalues of A are: -0.2360684.23607 Here's a matrix whose columns are eigenvectors of A corresponding to these eigenvalues: -0.850651 -0.5257310.525731 -0.850651

計算逆矩陣、行列式

首先，逆矩陣和行列式在數(shù)學上是一個基本的概念，但計算逆矩陣、計算行列式并不是一個好主意，因為這在處理問題時，并不是必須的。如上面介紹的，各種solve()更能滿足各種要求。而且，也可以不必使用計算行列式來判斷一個矩陣是否可逆。

但對一些小矩陣，計算逆矩陣或計算行列式并不是什么大問題。同時，一些分解算法也提供了inverse(), determinant()方法，來進行計算。

#include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {Matrix3f A;A << 1, 2, 1,2, 1, 0,-1, 1, 2;cout << "Here is the matrix A:\n" << A << endl;cout << "The determinant of A is " << A.determinant() << endl;cout << "The inverse of A is:\n" << A.inverse() << endl; }

執(zhí)行結果:

Here is the matrix A:1 2 12 1 0 -1 1 2 The determinant of A is -3 The inverse of A is: -0.667 1 0.3331.33 -1 -0.667-1 1 1

---

最小二乘法求解 Least squares solving

精度最高的最小二乘法求解是SVD分解法，而且Eigen提供了2中實現(xiàn)。推薦給大家的是 BDCSVD，它對大規(guī)模的求解問題能很好匹配，同時在處理小規(guī)模的問題時，自動回落到JacobiSVD處理模式。在設計上是一致性的，它們的solve()方法用于求解。

#include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {MatrixXf A = MatrixXf::Random(3, 2);cout << "Here is the matrix A:\n" << A << endl;VectorXf b = VectorXf::Random(3);cout << "Here is the right hand side b:\n" << b << endl;cout << "The least-squares solution is:\n"<< A.bdcSvd(ComputeThinU | ComputeThinV).solve(b) << endl; }

執(zhí)行結果：

$ g++ -I /usr/local/include/eigen3 matrxi_decomp4.cpp -o matrxi_decomp4 $ $ ./matrxi_decomp4 Here is the matrix A:-0.999984 -0.0826997-0.736924 0.06553450.511211 -0.562082 Here is the right hand side b: -0.9059110.3577290.358593 The least-squares solution is:0.46358 0.0429898

還有其他的求解方法，比如Cholesky分解法，QR分解法等，可以查看Eigen相關參考說明和文檔。

---

分離分解器（求解計算對象）的計算和構造

在上面的示例中，矩陣分解/線性求解的源碼，看上去在分解器對象的構造和計算同時在一個語句內(nèi)。其實我們可以把它們分開，這樣構造時，我們知道了需要分解的矩陣，而且可以對這個分解器對象重復使用。這樣的好處是：

• 所有的分解器有缺省的構造器。
• 所有的分解器有做計算的方法，可以重復調(diào)用，而且可以重新初始化。
  這樣會為程序帶來額外的好處。

下面的示例，演示了’LLT’的線性求解，矩陣分解方法。重復使用了LLT<Matrix2f> llt分解器對象。

#include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {Matrix2f A, b;LLT<Matrix2f> llt;A << 2, -1, -1, 3;b << 1, 2, 3, 1;cout << "Here is the matrix A:\n" << A << endl;cout << "Here is the right hand side b:\n" << b << endl;cout << "Computing LLT decomposition..." << endl;llt.compute(A);cout << "The solution is:\n" << llt.solve(b) << endl;A(1,1)++;cout << "The matrix A is now:\n" << A << endl;cout << "Computing LLT decomposition..." << endl;llt.compute(A);cout << "The solution is now:\n" << llt.solve(b) << endl; }

執(zhí)行結果：

Here is the matrix A:2 -1 -1 3 Here is the right hand side b: 1 2 3 1 Computing LLT decomposition... The solution is: 1.2 1.4 1.4 0.8 The matrix A is now:2 -1 -1 4 Computing LLT decomposition... The solution is now:1 1.291 0.571

而且，我們可以指定分解器的處理矩陣尺寸，讓編譯器編譯時，預分配存儲空間，這樣在后面執(zhí)行時，無需動態(tài)分配內(nèi)存空間，從而提高執(zhí)行速度和執(zhí)行效率。
比如：

HouseholderQR<MatrixXf> qr(50,50); MatrixXf A = MatrixXf::Random(50,50); qr.compute(A); // no dynamic memory allocation

---

秩分解 Rank-revealing decompositions

某些分解可以用來對矩陣求秩，可以稱為秩分解。代表性的是，奇異矩陣的非全秩矩陣的分解。在Catalogue of dense decompositions上可以看到Eigen的分解器是否支持Rank-Revealing。

具有矩陣秩分解的分解器，至少會提供rank()函數(shù)。它們可能還會提供isInvertible();還有一些提供了計算空域的核kernel的方法，也有提供列空間的像image。比如FullPivLU,如下示例：

#include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() {Matrix3f A;A << 1, 2, 5,2, 1, 4,3, 0, 3;cout << "Here is the matrix A:\n" << A << endl;FullPivLU<Matrix3f> lu_decomp(A);cout << "The rank of A is " << lu_decomp.rank() << endl;cout << "Here is a matrix whose columns form a basis of the null-space of A:\n"<< lu_decomp.kernel() << endl;cout << "Here is a matrix whose columns form a basis of the column-space of A:\n"<< lu_decomp.image(A) << endl; // yes, have to pass the original A }

執(zhí)行：

$ g++ -I /usr/local/include/eigen3 matrxi_decomp5.cpp -o matrxi_decomp5 $ $ $ ./matrxi_decomp5 Here is the matrix A: 1 2 5 2 1 4 3 0 3 The rank of A is 2 Here is a matrix whose columns form a basis of the null-space of A:0.51 -0.5 Here is a matrix whose columns form a basis of the column-space of A: 5 1 4 2 3 3

當然,任何矩陣的秩rank的計算依賴于一個選擇的任意閾值, 實際上，非浮點數(shù)的矩陣是有秩缺陷的rank-deficient.。

計算特征Eigen時，可以選擇合理的默認閾值,這取決于分解算法，但通常的選擇是對角線尺寸乘以ε。雖然這是我們可以選擇的最好的缺省值。但針對您的應用程序，只有你才知道什么是正確的閾值。在分解計算對象調(diào)用rank()或任何其他方法之前，你可以設置通過調(diào)用setThreshold(),來設置一個合適的閾值。而分解計算方法自身，即compute()方法,獨立于閾值 – 你無需在修改了閾值后，再重新計算分解對象。

當然，任何矩陣秩的計算，有賴于一個閾值的選擇。特別是因為計算機計算的問題，浮點類型的矩陣是無秩的。Eigen中的分解器會選擇一個合理的閾值，通常是于機器精度的數(shù)倍于對角陣中的尺寸 – 這個說法很繞。如果你明確知道閾值，你可以使用setThreshold()函數(shù)來設置閾值。

注意: 矩陣分解器的computer()計算方法，并不依賴于這個閾值。

看一個示例：

#include <iostream> #include <Eigen/Dense>using namespace std; using namespace Eigen;int main() {Matrix2d A;A << 2, 1,2, 0.9999999999;FullPivLU<Matrix2d> lu(A);cout << "By default, the rank of A is found to be " << lu.rank() << endl;lu.setThreshold(1e-5);cout << "With threshold 1e-5, the rank of A is found to be " << lu.rank() << endl; }

執(zhí)行結果：

By default, the rank of A is found to be 2 With threshold 1e-5, the rank of A is found to be 1

總結

以上是生活随笔為你收集整理的Eigen密集矩阵求解 1 - 线性代数及矩阵分解的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。