c – 非常小的方阵的特征线性解算器

2024-02-05 18:26:28

我在C程序上使用Eigen来求解非常小的方阵(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
}

我测试了一些EigenSolver,其中包括：

> FullPixLU
> PartialPivLU
> HouseholderQR
> ColPivHouseholderQR
> ColPivHouseholderQR
> CompleteOrthogonalDecomposition
> LDLT
>直接逆

直接反转是：

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毫秒
> CompleteOrthogonalDecomposition：832 ms
> LDLT：477毫秒
>直接反转：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的selfadjoint视图复制到临时中来暗示它,并希望编译器足够聪明以找到公共子表达式：

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

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

如果这对性能至关重要,请尝试实现手动调整的inverse_4x4_symmetric方法.
利用inv * b的对称性对于这样的小矩阵来说不太可能是有益的.

码农公寓

相关文章