用于非常小方阵的特征线性求解器

Eigen Linear Solver for very small square matrix

本文关键字：线性 特征 非常 方阵 用于      更新时间：2023-10-16

我正在C++程序上使用本征求解非常小的方阵(4X4(的线性方程。

我的测试代码是这样的

template<template <typename MatrixType> typename EigenSolver>
Vertor3d solve(){
//Solve Ax = b and A is a real symmetric matrix and positive semidefinite
... // Construct 4X4 square matrix A and 4X1 vector b
EigenSolver<Matrix4d> solver(A);
auto x = solver.solve(b);
... // Compute relative error for validating
}

我测试了一些特征求解器，其中包括：

  1. 全像素LU
  2. 部分皮夫LU
  3. 房主二维码
  4. 科尔皮夫户主QR
  5. 科尔皮夫户主QR
  6. 完全正交分解
  7. LDLT
  8. 直接逆

直接逆是：

template<typename MatrixType>
struct InverseSolve
{
private:
MatrixType  inv;
public:
InverseSolve(const MatrixType &matrix) :inv(matrix.inverse()) {
}
template<typename VectorType>
auto solve(const VectorType & b) {
return inv * b;
}
};

我发现快速方法是DirectInverse，即使我将 Eigen 与 MKL 链接，结果也没有改变。

这是测试结果

  1. 全像素LU ： 477 ms
  2. 部分皮夫LU ： 468 ms
  3. 户主QR ： 849 ms
  4. ColPivHouseholderQR ： 766 ms
  5. ColPivHouseholderQR ： 857 ms
  6. 完全正交分解 ： 832 ms
  7. LDLT ： 477 ms
  8. 直接反转 ： 88 ms

它们都使用 1000000 个矩阵，从均匀分布 [0,100] 随机双精度。我小心翼翼地构造上三角形，然后复制到下三角形。

DirectInverse的唯一问题是它的相对误差略大于其他求解器，但可以接受。

我的程序是否有更快或更优雅的解决方案？DirectInverse 是我的程序的快速解决方案吗？

DirectInverse不使用对称信息，那么为什么DirectInverse比LDLT快得多呢？

尽管许多人建议当你只想求解线性系统时，永远不要显式计算逆，但对于非常小的矩阵，这实际上是有益的，因为有使用协因子的闭式解。

您测试的所有其他替代方案都会变慢，因为它们将进行透视(这意味着分支(，即使对于小型固定大小的矩阵也是如此。此外，它们中的大多数会导致更多的除法，并且不能像直接计算那样向量化。

为了提高精度(如果需要，这种技术实际上可以独立于求解器使用(，您可以通过使用残差再次求解系统来优化初始解：

Eigen::Vector4d solveDirect(const Eigen::Matrix4d& A, const Eigen::Vector4d& b)
{
Eigen::Matrix4d inv = A.inverse();
Eigen::Vector4d x = inv * b;
x += inv*(b-A*x);
return x;
}

我不认为 Eigen 直接提供了一种利用这里A对称性的方法(对于直接计算的逆(。您可以尝试通过将A的自伴随视图显式复制到临时视图中来暗示这一点，并希望编译器足够聪明以找到常见的子表达式：

Eigen::Matrix4d tmp = A.selfadjointView<Eigen::Upper>();
Eigen::Matrix4d inv = tmp.inverse();

为了减少一些划分，你也可以用-freciprocal-math(在 gcc 或 clang 上(进行编译，这当然会略微降低准确性。

如果这确实对性能至关重要，请尝试实现手动调整的inverse_4x4_symmetric方法。 利用inv * b的对称性不太可能对如此小的矩阵有益。

相关文章：