「C++: Eigen」第二章 第一节 Linear algebra and decompositions

参考链接

Eigen支持分解矩阵及其对应方法表：Catalogue of decompositions offered by Eigen
Eigen各矩阵分解方法基准参考：矩阵分解基准

2 密集线性问题和分解

第二章主要解决线性相关问题，以及一些矩阵分解问题。第二章主要包括五个小节

2.1 线性代数和分解

本节主要介绍如何使用Eigen应对线性系统中的问题，以及计算各种分解LU、QR、SVD等

2.1.1 基础线性求解

对于线性方程组Ax = b，我们可以根据A的属性选择不同的分解方法，并且决定分解追求速度还是精度。下面是一个线性方程求解的例子：

#include 
#include 
 
int main()
{
   Eigen::Matrix3f A;
   Eigen::Vector3f b;
   A << 1,2,3,  4,5,6,  7,8,10;
   b << 3, 3, 4;
   std::cout << "Here is the matrix A:\n" << A << std::endl;
   std::cout << "Here is the vector b:\n" << b << std::endl;
   Eigen::Vector3f x = A.colPivHouseholderQr().solve(b);
   std::cout << "The solution is:\n" << x << std::endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14

其中colPivHouseholderQr()方法返回一个colPivHouseholderQR对象，因此上面的求解还可以写成：

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

其中，ColPivHouseholderQR是一个带列旋转的QR分解（QR分解- wiki百科）。

所有分解方法均带有一个solve()方法，用于分解调用。

如果在了解矩阵性质的前提下，可以根据性质选择最佳分解方法。例如，满秩非对称矩阵可以使用 PartialPivLU 方法；正定矩阵（positive-definite matrix）则LLT或LDLT更合适。

2.1.2 最小二乘求解

最小二乘：通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小

在最小二乘场景中，解决欠定（方程数小于变量个数，有无穷多解）方程组和超定（方程数大于变量个数方程组）方程组最通用、精确的方法是SVD分解。

Eigen提供两个方法，BDCSVD适用于比较大的问题，JacobiSVD适用于比较小的问题。除了SVD分解，还有一个与SVD有相同精度但速度更快的替代分解方法：完全正交分解，CompleteOrthogonalDecomposition

2.1.3 奇异矩阵检查

奇异矩阵：非满秩矩阵

奇异矩阵检查一般需要确定允许误差范围，例如下面这个例子，它求解出x后，计算了Ax - b的误差：

#include 
#include 
 
using Eigen::MatrixXd;
 
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 norm
   std::cout << "The relative error is:\n" << relative_error << std::endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13

2.1.4 计算特征值和特征向量

伴随矩阵：与矩阵的逆之间只相差一个系数

特征值和特征向量的计算不一定收敛，但这种不收敛的情况很少见。调用 info() 是为了检查这种可能性。下面是一个求解特征值和特征向量的例子：

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

2.1.5 计算逆和行列式

对于确定分解的方法 PartialPivLU 和 FullPivLU，它们提供了inverse()和determinant()方法来求解矩阵的逆和行列式。当然也可以直接对矩阵使用inverse()和determinant()方法。在eigen中，如果矩阵大小小于等于4，那么可以通过避免使用LU分解，转而使用inverse()和determinant()方法来求解：

#include 
#include 
 
int main()
{
   Eigen::Matrix3f A;
   A << 1, 2, 1,
        2, 1, 0,
        -1, 1, 2;
   std::cout << "Here is the matrix A:\n" << A << std::endl;
   std::cout << "The determinant of A is " << A.determinant() << std::endl;
   std::cout << "The inverse of A is:\n" << A.inverse() << std::endl;
}
1
2
3
4
5
6
7
8
9
10
11
12
13

2.1.6 将计算与构造分开

在上面的示例中，分解是在构造分解对象的同时计算的。然而，在某些情况下，我们可能希望将这两件事分开，例如在构造时不知道要分解的矩阵；或者如果想重用现有的分解对象。

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

最后，我们可以告诉分解构造函数为给定大小的矩阵分解预先分配存储，这样当随后分解这些矩阵时，就不会进行动态内存分配（当然，如果使用的是固定大小的矩阵，则不会进行动态内存分配）分配发生）。这是通过将大小传递给分解构造函数来完成的，例如：

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

2.1.7 等级揭示分解

某些分解是根据秩来计算的，也就是说，在分解时计算了矩阵的秩。这张表展示了各方法中是否计算矩阵秩。

等级揭示分解方法中，有提供rank()方法，除此之外还提供了isInvertible()方法，有些还提供计算矩阵的内核（零空间）和图像（列空间）的方法，例如FullPivLU。

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

当然，任何秩计算都取决于任意阈值的选择，因为实际上没有任何浮点矩阵是秩亏的。Eigen选择一个合理的默认阈值，这取决于分解，但通常是对角线大小乘以机器 epsilon。虽然这是我们可以选择的最佳默认值，但只有您知道您的应用程序的正确阈值是多少。您可以通过在调用 rank() 或任何其他需要使用此类阈值的方法之前对分解对象调用 setThreshold() 来设置它。分解本身，即 compute() 方法，与阈值无关。更改阈值后无需重新计算分解