Eigen入门教程（4）线性方程组的求解

myleader

线性方程的求解是个复杂的命题，理论上非常复杂，不过这里讨论的是软件包的用法，所以让我们忽略掉那些艰深晦涩的理论，直接看看Eigen怎么使用好了

行列式

古典的线性代数求解都是用行列式相除来求的，或者说在中国的教育中，通常都是从行列式入手开始讲解线性代数，尤其是二元、三元线性方程组的求解都是用行列式相除入门的。所以本教程先从行列式入手。Eigen提供了一个简单的函数determinant()

1. double x=A.determinant();

复制代码

就可以求得A的行列式，只要针对对应的列进行代换就可以求得Ax=B的解

1. double a=A.determinant();
   
2. MatrixXd C=A;
   
3. C.col(0)=B;
   
4. double c=C.determinant();
   
5. x(0)=c/a;

复制代码

当然，通常我们都不这么干，这种古典方法除了入门以外在实际应用中几乎没有任何用途。这里只是介绍一下行列式的求法，大家千万不要用这种方法来解方程组。

矩阵求逆与除法

在上一篇Eigen入门教程（3）中介绍了矩阵的四则运算。矩阵的除法是非常常用的线性方程组求解方法。方程组Ax=B的解就是x=A\B，而且矩阵的除法其实就是求逆后再相乘。Eigen提供了inverse()函数用于求矩阵的逆，就像Eigen入门教程（3）中介绍的那样，我们很容易得到方程组的解

1. VectorXd x=A.inverse()*B;

复制代码

相对于行列式除法的解法，矩阵求逆和除法在矩阵求解中的用处就要大一些。按照Eigen官方手册所述，这种方法对于小尺寸线性方程组速度特别快，尤其是针对4元及以下方程组，比后面介绍的方法还要快很多，这是因为SSE指令集能一次操作4个数的原因，这一点对于涉及到空间点阵等3元方程组或者OpenGL数据的4元方程组时优势特别明显，笔者在这里也着力推荐：4元及以下方程组推荐采用矩阵求逆的方法求解。

矩阵的分解

对于大型矩阵而言，最常用的求解方法其实是分解法求解，针对不同的方程组和不同的应用环境，有多种不同的分解方法。本文只是讲述Eigen软件包的使用，对于其内在的分解方法这里不会赘述，有兴趣的童鞋可以自己去翻文献。

Eigen提供了多种矩阵分解的算法，这里列出7种常用的，它们的优缺点如下：

Decomposition 分解方法	Module 所在模块	Method 函数	Requirements on the matrix 对矩阵的要求	Speed 计算速度	Accuracy 计算精度
PartialPivLU	#include <Eigen/LU>	partialPivLu()	Invertible可逆	++	+
FullPivLU	#include <Eigen/LU>	fullPivLu()	None任意矩阵	-	+++
HouseholderQR	#include <Eigen/QR>	householderQr()	None任意矩阵	++	+
ColPivHouseholderQR	#include <Eigen/QR>	colPivHouseholderQr()	None任意矩阵	+	++
FullPivHouseholderQR	#include <Eigen/QR>	fullPivHouseholderQr()	None任意矩阵	-	+++
LLT	#include <Eigen/Cholesky>	llt()	Positive definite正定矩阵	+++	+
LDLT	#include <Eigen/Cholesky>	ldlt()	Positive or negative semidefinite 半正定或半负定矩阵	+++	++

由于计算机采用二进制表示数字，所以会有浮点损失，而这种损失会在不同的计算操作中产生累加，有的累加多，有的累加少，并由此带来计算精度不同程度的下降，这取决于你的需求。

首先必须要根据A矩阵的性质选取可用的分解方法，然后再根据计算能力和计算精度的需求来寻找合适的最合适的分解方法。

按照Eigen入门教程（1）中所述，你可以直接包含Dense头文件，也可以按需包含各模块，考虑到C++编译的模式，推荐按需包含各模块，这会使编译速度加快。在官方手册中讲解到矩阵的求解部分时，都是直接包含Dense的，楼主这里特意调整一下，在所提供的样例代码中只包含必须的模块。

如果你要计算结构动力学，通常都是正定矩阵，那么采用LLT分解速度快，内存占用小，如果是进行精确计算，那就用LDLT分解方法，多消耗一点内存，可以提高一些精度。对于漂浮在空中的飞行器可能是半正定矩阵，那就最好直接用LDLT方法

如果你要计算流体力学问题，通常都是非正定的，不过一般都是可逆的，那么如果你需要较快速度就可以采用PartialPivLU

至于QR类分解，除了FullPivHouseholderQR收敛性稍好以外，其他的都不太好，而FullPivHouseholderQR精度高，速度慢，这就决定了其适用范围比较狭小。除非你确定你的方程组更适合使用QR分解，楼主这里不建议使用

http://eigen.tuxfamily.org/dox/TopicLinearAlgebraDecompositions.html
提供了一些建议，针对不同的问题，最好采用什么求解器。楼主本人是搞流体的，所以基本上都在使用PartialPivLU，对于其他的求解器仅仅知道其用法，对于特性不甚了了。

分解器求解的方法使用起来也非常简单，来看一段样例代码

1. #include <iostream>
   
2. #include <Eigen/Cholesky>
   
3. using namespace std;
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. Matrix2f A, b;
   
8. A << 2, -1, -1, 3;
   
9. b << 1, 2, 3, 1;
   
10. cout << "Here is the matrix A:\n" << A << endl;
    
11. cout << "Here is the right hand side b:\n" << b << endl;
    
12. Matrix2f x = A.ldlt().solve(b);
    
13. cout << "The solution is:\n" << x << endl;
    
14. }

复制代码

A.ldlt().solve(b)就可以得到方程的解了。例子中是ldlt分解，如果你需要其他的分解方法，使用方法一样。
这段代码的输出结果是：

1. Here is the matrix A:
   
2. 2 -1
   
3. -1  3
   
4. Here is the right hand side b:
   
5. 1 2
   
6. 3 1
   
7. The solution is:
   
8. 1.2 1.4
   
9. 1.4 0.8

复制代码

诸位看官可以用MATLAB验证一下看看结果怎么样。当然这只是一个例子，对于这种尺寸的方程组，还是推荐使用求逆的方法来解，这里只是演示一下分解器的用法。

下面这些代码演示了一些其他的求解器的用法，这种用法可以返回一个bool值，以判断方程是否有解，不仅如此，还可以输出分解后的矩阵，这个矩阵也许你在后面的运算中还用得着，这样可以节省一些运算

1. // Solve Ax = b. Result stored in x. Matlab: x = A \ b.
   
2. bool solved;
   
3. solved = A.ldlt().solve(b, &x));  // A sym. p.s.d.    #include <Eigen/Cholesky>
   
4. solved = A.llt() .solve(b, &x));  // A sym. p.d.      #include <Eigen/Cholesky>
   
5. solved = A.lu()  .solve(b, &x));  // Stable and fast. #include <Eigen/LU>
   
6. solved = A.qr()  .solve(b, &x));  // No pivoting.     #include <Eigen/QR>
   
7. solved = A.svd() .solve(b, &x));  // Stable, slowest. #include <Eigen/SVD>
   
8. // .ldlt() -> .matrixL() and .matrixD()
   
9. // .llt()  -> .matrixL()
   
10. // .lu()   -> .matrixL() and .matrixU()
    
11. // .qr()   -> .matrixQ() and .matrixR()
    
12. // .svd()  -> .matrixU(), .singularValues(), and .matrixV()

复制代码

这里的lu()和qr()其实就是收敛性较好的fullPivLu()和fullPivHouseholderQr()

构造分解器

其实你也可以把求解器单独构造，这样也可以

1. #include <iostream>
   
2. #include <Eigen/Dense>
   
3. using namespace std;
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. Matrix2f A, b;
   
8. LLT<Matrix2f> llt;
   
9. A << 2, -1, -1, 3;
   
10. b << 1, 2, 3, 1;
    
11. cout << "Here is the matrix A:\n" << A << endl;
    
12. cout << "Here is the right hand side b:\n" << b << endl;
    
13. cout << "Computing LLT decomposition..." << endl;
    
14. llt.compute(A);
    
15. cout << "The solution is:\n" << llt.solve(b) << endl;
    
16. A(1,1)++;
    
17. cout << "The matrix A is now:\n" << A << endl;
    
18. cout << "Computing LLT decomposition..." << endl;
    
19. llt.compute(A);
    
20. cout << "The solution is now:\n" << llt.solve(b) << endl;
    
21. }

复制代码

先用LLT构造一个求解器，然后用compute()函数计算求解矩阵，然后用solve()函数求解。如果你有多组方程组要解，这样可以节省一些内存。对于其他的分解器，也有类似的用法。

这种方法可以给定分解器的尺寸，对于指定尺寸的求解器，不能重新分配内存，这样可以节省一些内存分配的操作，程序运行更快。

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

复制代码

总之，构造分解器的方法会在一些情况下节省内存和内存操作，速度更快；但是它需要先在内存中开辟一块空间，这又需要额外的内存操作；如何使用才能加快运算速度取决于具体的问题，要由程序员自己把握。

求解的验算

如果你想了解一下计算的精度，可以这样

1. #include <iostream>
   
2. #include <Eigen/LU>
   
3. using namespace std;
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. MatrixXd A = MatrixXd::Random(100,100);
   
8. MatrixXd b = MatrixXd::Random(100,50);
   
9. MatrixXd x = A.fullPivLu().solve(b);
   
10. double relative_error = (A*x - b).norm() / b.norm(); // norm() is L2 norm
    
11. cout << "The relative error is:\n" << relative_error << endl;
    
12. }

复制代码

其实就是把求得的结果带回原方程组求方程右侧的值，和原始值对比一下，输出其差异，上面这段代码的输出结果是

1. The relative error is:
   
2. 2.29835e-14

复制代码

这个例子中的精度也还够用。但是童鞋们要注意了，如果这个精度差别特别大你就要考虑是不是方程求解出错。你可以验算一下行列式，如果行列式结果为0，那么很可能本来这个方程组就是无解的，那无论如何都会出错。如果你要写的求解器非常重要，推荐你在求解之后再验算一下这个差值，以便确定计算确实没错，这个验算过程只是一个简单的乘法后求模，计算量和分解、求解比差得多，对于有品质要求的计算这种验算的代价其实是值得的。

楼主曾经自己编译LAPACK，在最新的3.4.2版中增加了求解精度的测试项目。在Intel扣肉M520的CPU上，如果采用官方blas后端，这个精度蛮高的，如果采用openblas做后端，复数计算的精度下降就很厉害。

最小二乘法

楼主个人印象中，这种方法使用的非常少，不过既然官方手册当要点特意说明了，那就拿来一叙吧。看一下样例代码

1. #include <iostream>
   
2. #include <Eigen/SVD>
   
3. using namespace std;
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. MatrixXf A = MatrixXf::Random(3, 2);
   
8. cout << "Here is the matrix A:\n" << A << endl;
   
9. VectorXf b = VectorXf::Random(3);
   
10. cout << "Here is the right hand side b:\n" << b << endl;
    
11. cout << "The least-squares solution is:\n"
    
12. << A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << endl;
    
13. }

复制代码

使用起来和其他的分解器一样res=A.jacobiSvd().solve(b)就可以了。要特意说明的是：jacobiSvd()函数支持一个bitset参数，如果你不给定，那么它会有默认值（不过楼主没认真看过代码，所以不知道默认值是什么），这个bitset参数可以用位运算（通常都是用或），实际上Eigen提供了一些枚举参数，由这些枚举参数来进行位运算比较安全，如果未来Eigen升级改变了bitset参数的含义，你可能就需要修改你的代码，而且修改量可能会很大，采用枚举参数则不必修改你自己的代码，Eigen会做好枚举与bitset的对应。
http://eigen.tuxfamily.org/dox/group__enums.html
提供了Eigen所有的枚举参数，适用于jacobiSvd的参数也在其中，如果需要请自行查阅，楼主不用最小二乘法，所以对这个也没有什么研究。

Eigen没有提供高斯消元法，这种算法不适合计算机做数值运算，不仅除法过多带来的浮点损失大，而且速度慢，还不方便编程中使用矢量指令集加速。楼主绝没有瞧不起先贤的意思，只是这种3个世纪以前的算法和行列式一样确实仅适用于入门，在计算机时代已经不合适了。

接下来稍微介绍一下跟求解有关的一些高级内容

收敛阀值

分解和求秩运算的收敛性判断与收敛准则密切相关，纯理论上适用于任何精度，但在计算机中精度都是有限的，如果是整数矩阵，而且这个阀值也取整数的话，那收敛性就不好说了。Eigen软件包对不同类型的矩阵设定了默认的阀值，按照官方手册所述：通常这个阀值都是对角线尺寸与机器所支持的最小值的乘积。如果使用构造分解器的方法，你可以使用setThreshold()函数设定这个阀值用以调整收敛性和计算精度。

根据官方手册所述，这个阀值在方程求解中作用不是特别大，但是在求秩时更有用。楼主实际上不是十分理解，因为对于实际应用而言，我们其实不太关心秩的精度，我们更关心解的精度。没关系，吐啊吐啊就习惯了。

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

复制代码

看一下输出结果

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这个英文单词的意思就是特征值，这个软件包自然是能求特征值的。线性代数理论中对特征值的解释非常简单易懂，但是真正要计算起来才发现，教科书中的算法计算量不知道有多大，对于计算程序而言绝不是教科书上那么简单。

首先这个矩阵必须有伴随矩阵，否则不存在特征值。即使如此，求特征值的计算也不一定收敛，不信的话你自己拿MATLAB试试，不过Eigen软件包还好，提供了一个info()函数用以判断收敛性，如果收敛再计算特征值，这样程序的稳定性就会好一些。

1. #include <iostream>
   
2. #include <Eigen/Eigenvalues>
   
3. using namespace std;
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. Matrix2f A;
   
8. A << 1, 2, 2, 3;
   
9. cout << "Here is the matrix A:\n" << A << endl;
   
10. SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
    
11. if (eigensolver.info() != Success) abort();
    
12. cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
    
13. cout << "Here's a matrix whose columns are eigenvectors of A \n"
    
14. << "corresponding to these eigenvalues:\n"
    
15. << eigensolver.eigenvectors() << endl;
    
16. }

复制代码

先用已知矩阵定义一个SelfAdjointEigenSolver，在用info()函数检测一下是否有解，然后用eigenvalues()函数求特征值，如果需要还可以用eigenvectors()函数输出特征向量

这段代码的输出结果是

1. Here is the matrix A:
   
2. 1 2
   
3. 2 3
   
4. The eigenvalues of A are:
   
5. -0.236
   
6. 4.24
   
7. Here's a matrix whose columns are eigenvectors of A
   
8. corresponding to these eigenvalues:
   
9. -0.851 -0.526
   
10. 0.526 -0.851

复制代码

除了使用SelfAdjointEigenSolver构造求解器以外，还可以使用EigenSolver或ComplexEigenSolver构造求解器，下面是一段样例代码

1. MatrixXd A = MatrixXd::Random(6,6);
   
2. cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
   
3. EigenSolver<MatrixXd> es(A);
   
4. cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl;
   
5. cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
   
6. complex<double> lambda = es.eigenvalues()[0];
   
7. cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
   
8. VectorXcd v = es.eigenvectors().col(0);
   
9. cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
   
10. cout << "... and A * v = " << endl << A.cast<complex<double> >() * v << endl << endl;
    
11. MatrixXcd D = es.eigenvalues().asDiagonal();
    
12. MatrixXcd V = es.eigenvectors();
    
13. cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;

复制代码

后面求特征值和特征向量的用法和SelfAdjointEigenSolver一样，只不过这两个求解器的收敛性差一些。

在上一篇Eigen入门教程（3）中求迹的运算中其实没有使用特征值来计算，而是直接使用了四则运算，所以速度和收敛性都要更好一些，童鞋们求迹的时候千万不要用特征值求和的方法呀，不仅慢，而且还可能不收敛。

看过Eigen软件包对于其求特征值的介绍之后，楼主实际上有些吐血，这个自称为“特征值”的软件包在计算“特征值”这项功能上居然如此不自信，实在是盛名之下其实难副。当然本来特征值就是比较难算的，LAPACK其实处理的也不是那么好。所以楼主也就不再苛求了。越深入使用越了解她，越了解她的缺点，不过也更愿意包容她。真是爱之深恨之切。

还好他的LU分解、求逆、矩阵乘法的速度都蛮快，还支持自动多线程，所以楼主就忍了。

求矩阵的秩

有些分解算法如LU分解的计算过程中就可以生成矩阵的秩，如此来看只要直接输出就可以了。

1. #include <iostream>
   
2. #include <Eigen/LU>
   
3. using namespace std;
   
4. using namespace Eigen;
   
5. int main()
   
6. {
   
7. Matrix3f A;
   
8. A << 1, 2, 5,
   
9. 2, 1, 4,
   
10. 3, 0, 3;
    
11. cout << "Here is the matrix A:\n" << A << endl;
    
12. FullPivLU<Matrix3f> lu_decomp(A);
    
13. cout << "The rank of A is " << lu_decomp.rank() << endl;
    
14. cout << "Here is a matrix whose columns form a basis of the null-space of A:\n"
    
15. << lu_decomp.kernel() << endl;
    
16. cout << "Here is a matrix whose columns form a basis of the column-space of A:\n"
    
17. << lu_decomp.image(A) << endl; // yes, have to pass the original A
    
18. }

复制代码

实际上只有LU和QR能够在分解时顺便算秩，但是考虑到QR的适应性不好，推荐LU

要注意的是，凡是能够在分解过程中能够得到矩阵秩的算法，收敛性都比较好，尤其是对于非满秩矩阵。本来非满秩矩阵就是有无穷多解的，所以求解时容易发散，这种分解方法会在计算过程中把对角线为0的元素算出来并避免更进一步的计算错误。

如果你实在放心不下，官方手册中建议：可以用isInvertible()先检查一下矩阵，然后再计算矩阵。楼主并没有看过那段代码，以楼主个人所学，实在是难以理解，这个函数是用来检查矩阵是否可逆的，可逆矩阵都是满秩的，没必要再算秩了，理解不了。另外楼主高度怀疑那个函数其实是判断行列式是否为0

另外样例中的分解器还提供了kernel()函数返回矩阵中的null-space，以及image()函数——注意和复数矩阵取虚部的函数只差一个e——返回矩阵中的null-column

实话实说，楼主的线性代数水平还没到可以理解这两个术语的程度，如果有谁能指点一下，请在下面跟帖，楼主不胜感激。

看一下上面这段代码的计算结果

1. Here is the matrix A:
   
2. 1 2 5
   
3. 2 1 4
   
4. 3 0 3
   
5. The rank of A is 2
   
6. Here is a matrix whose columns form a basis of the null-space of A:
   
7. 0.5
   
8. 1
   
9. -0.5
   
10. Here is a matrix whose columns form a basis of the column-space of A:
    
11. 5 1
    
12. 4 2
    
13. 3 3

复制代码

到这里方程求解的常用功能应该都介绍到了。各位童鞋应该足以解决一般的方程组了。

和前面讲解的矩阵的常规运算不同，常规运算都是有确定结果的，而方程的求解很可能是无解或者是无穷多解的，这就带来了额外的变数——收敛性、稳定性等等问题，而且实际上考虑到计算机本身的精度极限和数学原理，这些问题很可能根本就是无法解决的。楼主在使用中吐过几次血，但是当我使用其他的软件包来尝试的时候，发现其实别人也解决不了这个问题，终于才明白是楼主误解了Eigen。所谓“爱之深恨之切”，深入了解她不仅仅意味着了解她的优点，更有缺点。工具在我们手里能发挥出多大的功效来，不仅仅跟工具有关，更与我们使用者的水平有关，喷工具的缺点或者把工具的优点捧上天都是不对的，了解工具的优缺点并有意的扬长避短才是我们应该做的。楼主在几篇文章中都特意强调了Eigen的缺点就是希望大家在使用中能够有效的规避风险，发挥出更大的优势来。

常规矩阵的部分基本上就这些内容了，接下来的一篇准备介绍稀疏矩阵，敬请期待