有效求解 Ax = b，其中 A 是 4x4 对称矩阵，b 是 4 x 1 向量 [关闭]答案

【问题标题】：Efficiently solve Ax = b where A is a 4x4 symmetric metrix and b is 4x1 vector [closed]有效求解 Ax = b，其中 A 是 4x4 对称矩阵，b 是 4 x 1 向量 [关闭]
【发布时间】：2014-05-03 14:26:33
【问题描述】：

我将求解一个小型线性系统 Ax = b，其中 A 是一个 4×4 对称矩阵，存储了 16 个 double 数字（实际上其中 10 个足以表示它） , b 是 4×1 向量。问题是，我必须运行这种系统数百万次。所以我正在寻找最有效的库来解决它。我在OpenCV 中尝试了cv::solve() 方法，但我仍然觉得它很慢。

由于矩阵A 是对称的，我记得Conjugate Gradient 算法由于其效率可能是一个很好的候选者。但是，我还没有找到一个库（Intel MKL 似乎有一个，但它是为稀疏矩阵设计的，不适合我的问题）。

谁能帮帮我？

【问题讨论】：

A 是否总是（或至少一次）相同？
不，每次我解决它，A 和 b 都会是新的。
如果你正在寻找一个图书馆，这是题外话，如果不是，这可能是题外话（和math.stackexchange.com 的话题）。
由于尺寸很小，而且您只有 4 个值可以猜测，我怀疑您最好对解决方案进行硬编码，而不是依赖可能无法获取所有细节的库（如对称) 考虑在内。
您是否对solve() 可以使用的所有可能算法进行了实验？对于不同大小的矩阵，它们的性能可能会有很大差异。您可以在此处查看所有选项：docs.opencv.org/modules/core/doc/…

标签： c++ opencv linear-equation

【解决方案1】：

由于矩阵维度是固定的，我认为你最好直接实现逆。这个任务有一个现成的公式。你有：

B 的条目由以下给出：

两个公式均取自this site。

您应该能够利用矩阵是对称的这一事实进一步简化这些条目的计算。如果你这样做，我认为你会比任何一般的矩阵逆实现更快。

那么你仍然需要将A^-1 应用到你的b，这是一个简单的矩阵向量乘法，你还应该硬编码，以获得最佳性能。

所有这些都假设您确实需要针对这个特定问题的最佳性能。如果矩阵维度发生变化，或者这不是算法中最关键的部分，请查看Eigen、Lapack/Blas 或任何其他线性代数库。求解稠密线性系统是基本任务，任何此类库都应包含该任务。

【讨论】：

记得在随机矩阵上针对其他解决方案（例如 openCV 逆矩阵解决方案）测试您的代码，以避免硬编码时出现索引错误；）
这是克莱默的规则。对于 $3\times3$ 矩阵可能没问题，但对于 $4\times 4$ 矩阵则不行。在这种情况下，LU 分解可能更快。 LU 分解通常在数值上也更稳定。 en.wikipedia.org/wiki/LU_decomposition#Solving_linear_equations 在早期阶段向学生教授矩阵行列式和矩阵逆，但通常可以避免两者。

【解决方案2】：

看在上帝的份上，请不要自己写。如果我理解正确，您正在寻求有效地解决密集线性系统。这正是 LAPACK 的用途。 netlib.org 上的版本（请参阅 this page 以获取有关您应该使用哪个例程的指导）非常快，但如果您需要一些真正令人尖叫的东西，请查看 MKL、ATLAS 或 GotoBLAS。

编辑：由于这是一个 C++ 论坛，我应该补充一点，您可以使用 Eigen 包来解决问题。它将使用 LAPACK 例程之一的一些实现。

【讨论】：

我看到的主要问题是，LAPACK 没有针对具有 4x4 矩阵的小型系统进行优化，而是针对具有 1000x1000 矩阵的系统进行了优化。因此，就性能而言，我认为它们不会接近硬编码实现。如果这值得努力实现硬编码，则取决于计算的频率。另一方面，开始使用 LAPACK/BLAS 并不一定比直接实现更容易。
如果您查看Eigen main page 的底部，他们包含的优化的 4x4 矩阵反转代码归功于英特尔。我敢打赌，这已经足够快了。正如@Micka 指出的那样，测试是一个问题，您可以打赌 Eigen 代码经过充分测试。
很好，我不知道 Eigen（或者在这种情况下是英特尔）实际上对这种小情况有高度优化的代码。我想它通常用于出现这种小问题的一些分治算法。对于该特定问题，您可能不会击败 SSE 优化代码，这可能是对的。对于测试，我假设 OP 有一个可以测试的工作求解器，并且只需要一个高性能变体。否则，Eigen 肯定是一个很好的测试对象。

【解决方案3】：

这个主题可能已经过时了，但是我解决了非常相似的问题，这个主题可能对其他人有用。

4 维非常小，因此直接计算矩阵求逆的算法是有效的。但是，如果您只需要一个解决方案的 A 矩阵，则无需计算整个反演。可以看出除以行列式（矩阵求逆中的一种运算）对于逆矩阵的所有项都不是必需的。您只能划分解的 4 个项，而不是所有 16 个矩阵项。还有一些其他的优化。这是我的 4D-solver 实现（c++ 代码）

void getsub(double* sub, double* mat, double* vec)
{
*(sub++) = *(mat  ) * *(mat+5) - *(mat+1) * *(mat+4);
*(sub++) = *(mat  ) * *(mat+6) - *(mat+2) * *(mat+4);
*(sub++) = *(mat  ) * *(mat+7) - *(mat+3) * *(mat+4);
*(sub++) = *(mat  ) * *(vec+1) - *(vec  ) * *(mat+4);
*(sub++) = *(mat+1) * *(mat+6) - *(mat+2) * *(mat+5);
*(sub++) = *(mat+1) * *(mat+7) - *(mat+3) * *(mat+5);
*(sub++) = *(mat+1) * *(vec+1) - *(vec  ) * *(mat+5);
*(sub++) = *(mat+2) * *(mat+7) - *(mat+3) * *(mat+6);
*(sub++) = *(mat+2) * *(vec+1) - *(vec  ) * *(mat+6);
*(sub  ) = *(mat+3) * *(vec+1) - *(vec  ) * *(mat+7);
}

void solver_4D(double* mat, double* vec)
{
    double a[10], b[10]; // values of 20 specific 2D subdeterminants

    getsub(&a[0], mat, vec);
    getsub(&b[0], mat+8, vec+2);

    *(vec++) = a[5]*b[8] + a[8]*b[5] - a[6]*b[7] - a[7]*b[6] - a[4]*b[9] - a[9]*b[4];
    *(vec++) = a[1]*b[9] + a[9]*b[1] + a[3]*b[7] + a[7]*b[3] - a[2]*b[8] - a[8]*b[2];
    *(vec++) = a[2]*b[6] + a[6]*b[2] - a[0]*b[9] - a[9]*b[0] - a[3]*b[5] - a[5]*b[3];
    *(vec  ) = a[0]*b[8] + a[8]*b[0] + a[3]*b[4] + a[4]*b[3] - a[6]*b[1] - a[1]*b[6];

    register double idet = 1./(a[0]*b[7] + a[7]*b[0] + a[2]*b[4] + a[4]*b[2] - a[5]*b[1] - a[1]*b[5]);

    *(vec--) *= idet;
    *(vec--) *= idet;
    *(vec--) *= idet;
    *(vec  ) *= idet;
}

*mat表示指向A矩阵的指针

0 1 2 3
4 5 6 7
8 9 10 11
12 13 14 15

和*vec 以0 1 2 3 的形式表示惯性矢量。方法丢弃惯性向量并用新的替换它们。

您只需要 74 次乘法和 1 次除法。如果您在 SSE 优化方面经验丰富，此方法可能适用于 SSE 编码器（大多数操作是双重的）。

一般来说，这对所有可逆 A 矩阵都很有效。您可以将身份 a[7] = b[0] 用于对称 A 矩阵。不多，但是一般4维可逆A矩阵的原码还是挺快的。

【讨论】：