原文：Eigen官网-Linear algebra and decompositions

本篇文章介绍了线性方程求解、矩阵分解，包括LU分解法，QR分解法，SVD（奇异值分解）、特征值分解等。

The problem：对于一般形式如下的线性系统：
[ Ax : = : b ]

The solution：解决上述方程的方式一般是将矩阵A进行分解，你可以根据A的形式选择不同的分解方法。当然最基本的方法是高斯消元法。
示例如下：

#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;
}

结果如下：

Here is the matrix A:
 1  2  3
 4  5  6
 7  8 10
Here is the vector b:
3
3
4
The solution is:
-2
 1
 1

Eigen内置的解线性方程组的算法如下表所示：

如上所示的所有分解方法都提供了一个solve()方法，如上例中所示。

比如，假设有一个正定的矩阵，如上表中所示，最好的选择是使用LLT或LDLT分解方法。

#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;
}

结果如下：

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

2. 确认解是否真的存在

你可以计算误差并确认得到的解是否是有效的，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 norm
   cout << "The relative error is:\\n" << relative_error << endl;
}

结果如下：

The relative error is:
2.31495e-14

3. 计算特征值和特征向量

首先要检查矩阵是否是self-adjoint。
SelfAdjointEigenSolver 可以通过EigenSolver或者ComplexEigenSolver方法适用于一般矩阵。

特征值和特征向量的计算不一定收敛，但这种不收敛的情况非常少见。调用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) 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;
}

结果如下：

Here is the matrix A:
1 2
2 3
The eigenvalues of A are:
-0.236
  4.24
Here's a matrix whose columns are eigenvectors of A 
corresponding to these eigenvalues:
-0.851 -0.526
 0.526 -0.851

4. 计算逆矩阵和矩阵行列式

Eigen 也提供了求逆矩阵和求矩阵行列式的算法，但是这两种算法对于大型矩阵来说都是非常不经济的算法，当需要对大型矩阵做这种的操作时，需要自己判断到底需不需这样做。但是对于小型矩阵 则可以没有顾虑地使用。

PartialPivLU和FullPivLU提供了inverse()和determinant()方法，你也可以直接对矩阵使用inverse()和determinant()方法。如果矩阵是一个很小的固定尺寸矩阵（至少是4*4），允许Eigen阻止使用LU分解，在这种小矩阵上使用公式将会更有效。

#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  1
 2  1  0
-1  1  2
The determinant of A is -3
The inverse of A is:
-0.667      1  0.333
  1.33     -1 -0.667
    -1      1      1

5. 最小二乘法求解

最小二乘问题最准确的方法是SVD分解。
Eigen也提供了解最小二乘问题的解法，并给出两种实现，分别是BDCSVD和JacobiSVD，BDCSVD可以很好地处理大问题，而对于较小的问题则自动返回到JacobiSVD类。因此推荐使用的一种是BDCSVD。这两个类的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;
}

结果如下：

Here is the matrix A:
  0.68  0.597
-0.211  0.823
 0.566 -0.605
Here is the right hand side b:
 -0.33
 0.536
-0.444
The least-squares solution is:
-0.67
0.314

6. 将计算与构造分离

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

使这成为可能的是：

• 所有分解都有一个默认构造函数，
• 所有分解都有一个compute（matrix）方法来执行计算，并且可以在已计算的分解上再次调用该方法，重新初始化它。

#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;
}

结果如下：

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.29
    1 0.571

最后，您可以告诉分解构造函数为分解给定大小的矩阵预先分配存储空间，以便在随后分解此类矩阵时，不执行动态内存分配（当然，如果使用的是固定大小的矩阵，则根本不进行动态内存分配）。这是通过将大小传递给分解构造函数来完成的，如本例所示：

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

7. 揭示秩的分解

一些分解是可以揭示秩的，也就是说可以计算矩阵的秩。
** Rank-revealing decompositions** 提供了rank()方法，isInvertible()，有一些还提供了计算kernel（null-space)和image(column-space)的方法，示例如下：

#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
}

结果如下：

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.5
   1
-0.5
Here is a matrix whose columns form a basis of the column-space of A:
5 1
4 2
3 3

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

#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 2
With threshold 1e-5, the rank of A is found to be 1

参考：

“Eigen学习之简单线性方程与矩阵分解”