2

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

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

有人可以帮我吗？

4

3 回答 3

5

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

逆矩阵公式

的条目由B以下给出： B_ij 的计算

这两个公式都取自这个网站。

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

然后你仍然需要申请A^-1你的b，这是一个简单的矩阵向量乘法，你也应该硬编码，以获得最佳性能。

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

于 2014-05-03T17:04:27.263 回答
1

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

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

于 2014-05-07T14:32:04.777 回答
0

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

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

于 2016-02-16T21:00:50.013 回答