我正在使用 Eigen一个求解极小方阵(4X4)线性方程的C++程序。
我的测试代码是这样的
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
}
我测试了一些 EigenSolver其中包括:
- FullPixLU
- PartialPivLU
- HouseholderQR
- ColPivHouseholderQR
- ColPivHouseholderQR
- 完全正交分解
- 低密度脂蛋白
- 正反
正逆是:
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连接起来,结果也没有改变。
这是测试结果
- FullPixLU:477 毫秒
- PartialPivLU:468 毫秒
- HouseholderQR:849 毫秒
- ColPivHouseholderQR:766 毫秒
- ColPivHouseholderQR:857 毫秒
- 完全正交分解:832 毫秒
- LDLT:477 毫秒
- 正反:88 毫秒
它全部使用 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
的自伴 View 显式复制到临时文件中来暗示这一点,并希望编译器足够聪明以找到公共(public)子表达式:
Eigen::Matrix4d tmp = A.selfadjointView<Eigen::Upper>();
Eigen::Matrix4d inv = tmp.inverse();
为了减少一些除法,您还可以使用 -freciprocal-math
(在 gcc 或 clang 上)进行编译,这当然会稍微降低准确性。
如果这确实对性能至关重要,请尝试实现手动调整的 inverse_4x4_symmetric
方法。
利用 inv * b
的对称性不太可能对如此小的矩阵有益。
关于c++ - 非常小方阵的特征线性求解器,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/50909385/