线性代数与稠密矩阵分解

本章介绍如何利用不同的矩阵分解方法（LU，QR，SVD，特征分解等）求解线性系统

矩阵分解

矩阵分解（decomposition/factorization）是指将矩阵拆分为若干矩阵的乘积，包括三角分解（LU）、满秩分解、QR分解、Jordan分解和奇异值分解（SVD）等。Eigen提供了常用的三种矩阵分解：LU、QR、SVD

三角分解 LU

$A=LU$

将原方阵（square matrix）分解成一个上三角矩阵和一个下三角形矩阵。主要用于简化一个大矩阵的行列式的计算、求逆矩阵、求解方程组。这种分解法所得到的上下三角形矩阵并非唯一。

LLT分解（Cholesky分解）

若A为正定矩阵，则A有如下唯一的分解形式：

$A=LL^*$

其中 $L$ 为下三角矩阵，$L^$ 为其共轭转置矩阵，当A为实矩阵时，$L^=L^T$

LDLT分解

若A为一对称矩阵且其任意k阶主子阵均不为零，即半正定或半负定，则A有如下惟一的分解形式：

$A=LDL^T$

其中 $L$ 为下三角形单位矩阵（即主对角线元素皆为1），$D$ 为对角矩阵，$L^T$ 为 $L$ 的转置矩阵。LDLT分解法是Cholesky分解法的改进，Cholesky分解法虽然不需要选主元，但其运算过程中涉及到开方问题，而LDLT分解法则避免了这一问题，可用于求解线性方程组。

QR分解

对列向量线性无关的矩阵 A ：

$A=QR$

其中，$Q$ 为 $m\times m$ 的酉矩阵（实数域内称为正交矩阵），$R$ 为上三角矩阵

奇异值分解 SVD

对任意矩阵 A：

$A=U\Sigma V^*$

其中 $U$ 和 $V$ 均为酉矩阵（实数域内称为正交矩阵，$V^*=V^T$）且不唯一， $\Sigma$ 为对角矩阵，对角线上的元素称为A的奇异值。用于解最小二乘和数据压缩。

用矩阵分解求解线性方程组

Eigen为矩阵提供了如图所示的分解方法，调用第二列的方法后会返回一个第一列对应的类对象，再通过调用该对象的 solve() 方法实现对方程的求解，具体见后文代码示例。

LU 分解

partialPivLu()

基于部分消元的LU分解（必须是可逆方阵）。$A=PLU$ ，其中 $L$ 为单位下三角矩阵，$U$ 为上三角矩阵，$P$ 为置换矩阵，即只进行行变换
fullPivLu()

基于全消元的LU分解（对任意矩阵）。$A=P^{-1}LUQ^{-1}$ ，$L,U,P$ 同上，$Q$ 为置换矩阵，即同时进行行变换和列变换。速度比部分消元慢
llt()

标准Cholesky分解（必须是对称正定矩阵）。$A=LL^*$
ldlt()

有主元的鲁棒Cholesky分解（必须是半正定或半负定矩阵）。$A=P^TLDL^*P$

QR 分解

householderQr()

使用 Househoder 变换实现 $A=QR$ 的分解，其中 $Q$ 为酉矩阵，$R$ 为上三角矩阵。无主元，快，但不稳定
colPivHouseholderQr()

基于行变换（列主元消元法）的 Householder 变换，$AP=QR$ 。较慢，但更精确
fullPivHouseholderQr()

同时进行行变换和列变换的 Householder 变换，$PAP’=QR$ 。最慢，但最稳定
completeOrthogonalDecomposition()

完全正交分解。可以看作是QR分解的推广

SVD 分解

bdcSvd()

最快的SVD算法（推荐使用）
jacobiSvd()

对小矩阵速度快，大矩阵慢

表格给出了上述所有分解方法的效率对比

代码示例：

#include <Eigen/Eigen>
#include <ctime>
#include <iostream>

using namespace std;
using namespace Eigen;

#define MATRIX_SIZE 10

int main()
{
     /* 假设矩阵A为方阵，待求解方程为 Ax=b */
     // 首先定义一个随机矩阵A，随机列向量b和待求解向量x
     MatrixXd A;
     Matrix<double, MATRIX_SIZE, 1> b;
     A = MatrixXd::Random(MATRIX_SIZE, MATRIX_SIZE);
     b = MatrixXd::Random(MATRIX_SIZE, 1);
     // 定义时间变量，方便比较效率
     clock_t start_t;
     double time;

     /* 方法一：直接右乘A的逆矩阵 x = A^{-1}*b */
     Matrix<double, MATRIX_SIZE, 1> x_inv;
     start_t = clock();
     x_inv = A.inverse()*b;
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_inv time = " << time << " ms" << endl;
     cout << "x_inv = " << x_inv.transpose() << endl;

     /* 方法二：普通LU分解 
      * Eigen 提供了四种LU分解的方法
      * - partialPivLu(): A = PLU -> 要求A必须可逆
      * - fullPivLu(): A = P^{-1}LUQ^{-1}
      * - llt(): A = LL^T -> 要求A必须对称正定
      * - ldlt(): A = LDL^T -> 要求A必须半正定或半负定
      */
     // 2.1 partialPivLu()
     Matrix<double, MATRIX_SIZE, 1> x_plu;
     start_t = clock();
     x_plu = A.partialPivLu().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_plu time = " << time << " ms" << endl;
     cout << "x_plu = " << x_plu.transpose() << endl; 
     // 2.2 fullPivLu()
     Matrix<double, MATRIX_SIZE, 1> x_flu;
     start_t = clock();
     x_flu = A.fullPivLu().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_flu time = " << time << " ms" << endl;
     cout << "x_flu = " << x_flu.transpose() << endl; 

     /* 方法三：LLT分解 */
     Matrix<double, MATRIX_SIZE, 1> x_llt;
     start_t = clock();
     // 因为llt分解要求 A 必须为正定矩阵，所以构造新的方程：(A^T*A)*x = A^T*b 来求解x
     x_llt = (A.transpose()*A).llt().solve(A.transpose()*b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_llt time = " << time << " ms" << endl;
     cout << "x_llt = " << x_llt.transpose() << endl; 

     /* 方法四：LDLT分解 */
     Matrix<double, MATRIX_SIZE, 1> x_ldlt;
     start_t = clock();
     x_ldlt = (A.transpose()*A).llt().solve(A.transpose()*b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_ldlt time = " << time << " ms" << endl;
     cout << "x_ldlt = " << x_ldlt.transpose() << endl; 

     /* 方法五：QR分解
      * Eigen 提供了四种QR分解的方法
      * - householderQr()
      * - colPivHouseholderQr()
      * - fullPivHouseholderQr()
      * - completeOrthogonalDecomposition()
      * 常用2和3
      */
     // 5.1 colPivHhouseholderQr()
     Matrix<double, MATRIX_SIZE, 1> x_cqr;
     start_t = clock();
     x_cqr = A.colPivHouseholderQr().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_cqr time = " << time << " ms" << endl;
     cout << "x_cqr = " << x_cqr.transpose() << endl; 
     // 5.2 fullPivHhouseholderQr()
     Matrix<double, MATRIX_SIZE, 1> x_fqr;
     start_t = clock();
     x_fqr = A.fullPivHouseholderQr().solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_fqr time = " << time << " ms" << endl;
     cout << "x_fqr = " << x_plu.transpose() << endl; 

     /* 方法六：SVD分解
      * Eigen 提供了两种SVD分解方法
      * - bdcSvd()
      * - jacobiSvd()
      */
     // 6.1 bdcSvd()
     Matrix<double, MATRIX_SIZE, 1> x_bsvd;
     start_t = clock();
     x_bsvd = A.bdcSvd(ComputeThinU|ComputeThinV).solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_bsvd time = " << time << " ms" << endl;
     cout << "x_bsvd = " << x_bsvd.transpose() << endl; 
     // 6.2 jacobiSvd()
     Matrix<double, MATRIX_SIZE, 1> x_jsvd;
     start_t = clock();
     x_jsvd = A.jacobiSvd(ComputeThinU|ComputeThinV).solve(b);
     time = 1000 * (clock() - start_t) / (double)CLOCKS_PER_SEC;
     cout << "x_jsvd time = " << time << " ms" << endl;
     cout << "x_jsvd = " << x_jsvd.transpose() << endl; 
}

计算结果：

计算特征值与特征向量

如图，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) 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

计算逆矩阵与行列式

尽管逆和行列式是基本的数学概念，但它们在数值线性代数中不像在纯数学中那么常用。求逆通常由 solve() 代替，而行列式一般并不是检测矩阵是否可逆的好方法。不过对于比较小的矩阵，逆和行列式还是比较有用的。

尽管Eigen提供了上述矩阵分解的方法，我们仍然可以直接调用 inverse() 方法和 determinant() 方法。若矩阵尺寸较小（不超过4x4），那么Eigen可以避免使用LU分解，而是使用数学公式，这样更高效。

求解最小二乘

求解最小二乘最精确的方法是使用SVD分解。Eigen提供了两种方法（见上文），推荐使用 bdcSvd() ，对大规模问题友好，同时对小规模问题能自动降级到 jacobiSvd() 方法；

另一种方法是使用 Cholesky 分解或 QR 分解，速度可能会快点，但准确率会下降；

第三种方法是将方程转化为 $A^TAx=A^Tb$ 并使用 ldlt() 方法求解。不过，如果 A 是一个病态矩阵的话，这种方法就不会太好，因为 $A^TA$ 的条件数是 $A$ 的平方，这意味着这种方法会比其他方法损失两倍数值。

将求解和构造分开

上面的例子将矩阵分解对象的构造和计算写在同一句中，当然也有方法将它们分开写，每种矩阵分解类都有一个默认的构造函数和一个 compute(matrix) 方法，该方法可以对不同的 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

也可以通过向构造函数中传入矩阵大小来预分配分解矩阵所需要的内存空间大小：

1
2
3

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

秩显式分解

某些分解方法可以同时计算出矩阵的秩，这些通常也是在面对非满秩矩阵（方阵时称为奇异矩阵）时表现最好的方法。包括 FullPivLU 、ColPivHouseholderQR 、FullPivHouseholderQR 、BDCSVD 、JacobiSVD 、SelfAdjointEigenSolver 、ComplexEigenSolver 和 EigenSolver

这些分解类至少提供了一个 rank() 方法用来求矩阵的秩。此外，还提供了一些方便的方法，例如 isInvertible() ；有些还提供了计算矩阵的内核(零空间)和图像(列空间)方法，例如 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
}

输出：

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

当然，任何秩计算都取决于对随机阈值的选择，因为实际上没有哪个浮点矩阵是完全秩亏的。Eigen会选择一个合理的默认阈值。我们也可以自己设定想要的正确阈值：在调用 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;
}

输出：

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