[Eigen中文文档] 就地矩阵分解

文档总目录

英文原文(Inplace matrix decompositions)

从 Eigen 3.3 开始，LU、Cholesky 和 QR 分解可以就地操作，即直接在给定的输入矩阵内操作。当处理大矩阵时，或者当可用内存非常有限（嵌入式系统）时，此功能特别有用。

为此，必须使用 Ref<> 矩阵类型实例化相应的分解类，并且必须使用输入矩阵作为参数构造分解对象。作为一个例子，让我们考虑一个带有部分旋转的就地 LU 分解。

声明一个 2x2 矩阵 A：

 Eigen::MatrixXd A(2,2); 
 A << 2, -1, 
      1, 3;
 std::cout << "Here is the input matrix A before decomposition:\n" << A << "\n";

输出：

Here is the input matrix A before decomposition:
 2 -1
 1  3

然后，声明就地 LU 分解对象 lu，并检查矩阵 A 的内容：

Eigen::PartialPivLU<Eigen::Ref<Eigen::MatrixXd> > lu(A);
std::cout << "Here is the input matrix A after decomposition:\n" << A << "\n";

输出：

Here is the input matrix A after decomposition:
  2  -1
0.5 3.5

在这里，lu 对象计算 L 和 U 因子并将其存储在矩阵 A 所持有的内存中。A 的系数在分解过程中被破坏，并由 L 和 U 因子代替，可以验证：

std::cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << "\n";

输出：

Here is the matrix storing the L and U factors:
  2  -1
0.5 3.5

然后，可以使用 lu 对象，例如解决 $Ax=b$ 问题：

Eigen::MatrixXd A0(2,2); A0 << 2, -1, 1, 3;
Eigen::VectorXd b(2);    b << 1, 2;
Eigen::VectorXd x = lu.solve(b);
std::cout << "Residual: " << (A0 * x - b).norm() << "\n";

输出：

Residual: 0

由于A和lu共享内存，修改矩阵A将使lu无效。可以通过修改内容A并再次尝试解决初始问题来轻松验证：

A << 3, 4, -2, 1;
x = lu.solve(b);
std::cout << "Residual: " << (A0 * x - b).norm() << "\n";

输出：

Residual: 15.8114

请注意，这里没有共享指针，用户需要使输入矩阵 A 和 lu 保持相同的生命周期。

如果想用修改后的 A 更新因式分解，则必须像往常一样调用 compute 方法：

A0 = A; // save A
lu.compute(A);
x = lu.solve(b);
std::cout << "Residual: " << (A0 * x - b).norm() << "\n";

输出：

Residual: 0

请注意，调用 compute 不会更改 lu 对象引用的内存。因此，如果使用不同于 A 的另一个矩阵 A1 调用 compute 方法，则不会修改 A1 的内容。这仍然是 A 的内容，它将用于存储矩阵 A1 的 L 和 U 因子。验证如下：

Eigen::MatrixXd A1(2,2);
A1 << 5,-2,3,4;
lu.compute(A1);
std::cout << "Here is the input matrix A1 after decomposition:\n" << A1 << "\n";

输出：

Here is the input matrix A1 after decomposition:
 5 -2
 3  4

矩阵 A1 没有改变，并且可以求解 $A1\ast x=b$，也可以不复制 A1 直接求得残差：

x = lu.solve(b);
std::cout << "Residual: " << (A1 * x - b).norm() << "\n";

输出：

Residual: 2.48253e-16

以下是支持这种就地机制的矩阵分解列表：

• class LLT
• class LDLT
• class PartialPivLU
• class FullPivLU
• class HouseholderQR
• class ColPivHouseholderQR
• class FullPivHouseholderQR
• class CompleteOrthogonalDecomposition