简介
这里介绍线性系统的解析，如何进行各种分解计算，如LU，QR，SVD，特征值分解等。
简单线性求解
在一个线性系统，常如下表示，其中A，b分别是一个矩阵，需要求x：
Ax=bAx \:= \: b Ax=b
在Eigen中，我们可以根据需要的精确性要求，性能要求，从而从好几种方法中选择一种来求解这个线性系统。
下面的示例，可以看到一种求解方法。
//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;}
执行：
$ 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 10Here is the vector b:334The solution is:-20.9999971
这个示例中，使用矩阵的colPivHouseholderQr()方法，其返回一个ColPivHouseholderQR实例对象。因为调用对象是一个Matrix3f类型，所以也可以如此设计
ColPivHouseholderQR<Matrix3f> dec(A);Vector3f x = dec.solve(b);
正如名字所示，这里ColPivHouseholderQR是一个采用列旋转方法的QR分解。下面Eigen提供的表列出了你可以使用分解类型：
其中的部分分解公式：
FullPivLU:
A=P−1LUQ−1A = P^{-1} L U Q^{-1} A=P−1LUQ−1
ColPivHouseholderQR:
AP=QR\mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} AP=QR
FullPivHouseholderQR:
PAP′=QR\mathbf{P} \, \mathbf{A} \, \mathbf{P}' = \mathbf{Q} \, \mathbf{R} PAP′=QR
LLT:
LLTLL^TLLT Cholesky decomposition of a symmetric, positive definite matrix A such that
A=LL∗=U∗UA = LL^* = U^*U A=LL∗=U∗U , where L is lower triangular.
LDLT
LDLTLDL^T LDLT
Cholesky decomposition:
A=PTLDL∗PA = P^TLDL^*P A=PTLDL∗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。
执行结果：
Here is the matrix A:2 -1-1 3Here is the right hand side b:1 23 1The solution is:1.2 1.41.4 0.8
检查是否有解
计算存在误差，只有我们知道了一个计算的解的可允许的误差容限error margin，我们才能判断一个解是否满足要求。在Eigen内，很容易做这种计算。
#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;}
执行结果：
The relative error is:3.75098e-13
计算特征值和特征向量
在线性代数里计算特征值Eigen value和特征向量 EigenVector是非常常见的计算，计算前，要检查确保你的矩阵是自伴随self-adjoint矩阵。Eigen里提供了SelfAdjointEigenSolver，很容易地在矩阵对象上使用EigenSolver , ComplexEigenSolver。
特征值和特征向量的计算未必收敛，这种非收敛的情况很罕见，但确实存在。Eigen中的SelfAdjointEigenSolver.info()函数可以用来检查。
#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;}
执行结果：
$ g++ -I /usr/local/include/eigen3 matrxi_decomp3.cpp -o matrxi_decomp3$ ./matrxi_decomp3Here is the matrix A:1 22 3The eigenvalues of A are:-0.2360684.23607Here's a matrix whose columns are eigenvectors of A corresponding to these eigenvalues:-0.850651 -0.5257310.525731 -0.850651
计算逆矩阵、行列式
首先，逆矩阵和行列式在数学上是一个基本的概念，但计算逆矩阵、计算行列式并不是一个好主意，因为这在处理问题时，并不是必须的。如上面介绍的，各种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;}
执行结果:
Here is the matrix A:1 2 12 1 0-1 1 2The determinant of A is -3The inverse of A is:-0.6671 0.3331.33-1 -0.667-111
最小二乘法求解Least squares solving
精度最高的最小二乘法求解是SVD分解法，而且Eigen提供了2中实现。推荐给大家的是BDCSVD，它对大规模的求解问题能很好匹配，同时在处理小规模的问题时，自动回落到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;}
执行结果：
$ g++ -I /usr/local/include/eigen3 matrxi_decomp4.cpp -o matrxi_decomp4$ $ ./matrxi_decomp4Here is the matrix A:-0.999984 -0.0826997-0.736924 0.06553450.511211 -0.562082Here is the right hand side b:-0.9059110.3577290.358593The least-squares solution is:0.463580.0429898
还有其他的求解方法，比如Cholesky分解法，QR分解法等，可以查看Eigen相关参考说明和文档。
分离分解器（求解计算对象）的计算和构造
在上面的示例中，矩阵分解/线性求解的源码，看上去在分解器对象的构造和计算同时在一个语句内。其实我们可以把它们分开，这样构造时，我们知道了需要分解的矩阵，而且可以对这个分解器对象重复使用。这样的好处是：
所有的分解器有缺省的构造器。所有的分解器有做计算的方法，可以重复调用，而且可以重新初始化。
这样会为程序带来额外的好处。
下面的示例，演示了’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;pute(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;pute(A);cout << "The solution is now:\n" << llt.solve(b) << endl;}
执行结果：
Here is the matrix A:2 -1-1 3Here is the right hand side b:1 23 1Computing LLT decomposition...The solution is:1.2 1.41.4 0.8The matrix A is now:2 -1-1 4Computing LLT decomposition...The solution is now:1 1.291 0.571
而且，我们可以指定分解器的处理矩阵尺寸，让编译器编译时，预分配存储空间，这样在后面执行时，无需动态分配内存空间，从而提高执行速度和执行效率。
比如：
HouseholderQR<MatrixXf> qr(50,50);MatrixXf A = MatrixXf::Random(50,50);pute(A); // no dynamic memory allocation
秩分解 Rank-revealing decompositions
某些分解可以用来对矩阵求秩，可以称为秩分解。代表性的是，奇异矩阵的非全秩矩阵的分解。在Catalogue of dense decompositions上可以看到Eigen的分解器是否支持Rank-Revealing。
具有矩阵秩分解的分解器，至少会提供rank()函数。它们可能还会提供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}
执行：
$ g++ -I /usr/local/include/eigen3 matrxi_decomp5.cpp -o matrxi_decomp5$ $ $ ./matrxi_decomp5Here is the matrix A:1 2 52 1 43 0 3The rank of A is 2Here is a matrix whose columns form a basis of the null-space of A:0.51-0.5Here is a matrix whose columns form a basis of the column-space of A:5 14 23 3
当然,任何矩阵的秩rank的计算依赖于一个选择的任意阈值, 实际上，非浮点数的矩阵是有秩缺陷的rank-deficient.。
计算特征Eigen时，可以选择合理的默认阈值,这取决于分解算法，但通常的选择是对角线尺寸乘以ε。虽然这是我们可以选择的最好的缺省值。但针对您的应用程序，只有你才知道什么是正确的阈值。在分解计算对象调用rank()或任何其他方法之前，你可以设置通过调用setThreshold(),来设置一个合适的阈值。而分解计算方法自身，即compute()方法,独立于阈值 – 你无需在修改了阈值后，再重新计算分解对象。
当然，任何矩阵秩的计算，有赖于一个阈值的选择。特别是因为计算机计算的问题，浮点类型的矩阵是无秩的。Eigen中的分解器会选择一个合理的阈值，通常是于机器精度的数倍于对角阵中的尺寸 – 这个说法很绕。如果你明确知道阈值，你可以使用setThreshold()函数来设置阈值。
注意:矩阵分解器的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;}
执行结果：
By default, the rank of A is found to be 2With threshold 1e-5, the rank of A is found to be 1
