用 C++ 求解正规方程组答案

【问题标题】：Solving normal equation system in C++用 C++ 求解正规方程组
【发布时间】：2011-06-02 20:43:08
【问题描述】：

我想求解线性方程组：

 Ax = b

A 是 n x m 矩阵（不是正方形），b 和 x 都是 n x 1 向量。在已知 A 和 b 的情况下，n 大约为 50-100，m 约为 2（换句话说，A 可以是最大值 [100x2]）。

我知道x的解决方案： $x = \inv(A^T A) A^T b$

我找到了几种解决方法：uBLAS (Boost)、Lapack、Eigen 等，但我不知道使用这些包的“x”的 CPU 计算时间有多快。我也不知道这在数值上是否快速为什么要解决'x'

对我来说重要的是，由于我是新手，所以 CPU 计算时间将尽可能短，并且有良好的文档。

解出正规方程Ax = b 后，我想使用回归和可能稍后应用卡尔曼滤波器来改进我的近似值。

我的问题是哪个 C++ 库更健壮、更快速地满足我上面描述的需求？

【问题讨论】：

如何将 n x m 矩阵乘以 n 维列向量？大概 x 实际上是 m 维的。
另外，您是否有一些规定最低限度符合流行语的要求？
@Eagle 我不认为 Boost uBLAS 库实现了这一点，但如果我错了，请纠正我。 uBLAS 似乎为您提供了向量、矩阵和基本运算（乘法、加法），但没有像 LU、QR、SVD 或矩阵求逆，更不用说 OLS 实现了。然而，它可能是一个很好的库来实现这样的算法。同样，如果我错了，或者您是否找到了一个好的 Boost uBLAS OLS 实现，请告诉我...
我错了，lu.hpp中有LU分解。除了包含的三角形solver，它还可以让您做一些事情

【解决方案1】：

uBlas 未经优化，除非您将其与优化的 BLAS 绑定一起使用。

以下针对多线程和 SIMD 进行了优化：

顺便说一句，我不知道你在做什么，但通常正态方程不是进行线性回归的正确方法。除非您的矩阵条件良好，否则应该首选 QR 或 SVD。

【讨论】：

还有用于 AMD 芯片的 ACML。我相信这个是免费的。
我不确定优化的多线程版本对于像这样微不足道的矩阵是否会有这么大的好处。
@David Heffernan。 100x100 并不是那么小。
@watson @Alexandre 我很高兴我的手机可以实时执行此操作，但我认为 OP 使用的是真正的计算机！
@David Heffernan。 100x100 对于 QR 来说足够大了。 SVD 更贵。多线程的大量机会。 100x100 双打不适合 L1 缓存。使用优化实现的更多理由。

【解决方案2】：

如果你真的需要专门化，你可以使用 Skilling 方法来近似矩阵求逆（到任意精度）。它仅使用阶 (N^2) 运算（而不是通常矩阵求逆的阶 N^3 - LU 分解等）。

它在 Gibbs 的论文中有所描述，链接到这里（大约第 27 页）：

【讨论】：

【解决方案3】：

这是一个最小二乘解，因为你有比方程更多的未知数。如果 m 确实等于 2，这告诉我一个简单的线性最小二乘对你来说就足够了。公式可以写成封闭形式。你不需要图书馆。

如果 m 是个位数，我仍然会说您可以使用 A(transpose)*A*X = A(transpose)*b 轻松解决此问题。求解系数的简单 LU 分解就足够了。这应该是一个比你想象的要简单得多的问题。

【讨论】：

他谈到了卡尔曼滤波器。我认为他对线性代数和 OLS 尤其熟悉。他想要优化的库。
@duffymo，你是对的，现在我正在寻找 $x = \inv(A^T A) A^T b$ 的解决方案。卡尔曼滤波器可能是为了未来的发展。对我来说重要的是我应该使用哪个库（支持逆、转置、矩阵乘法等）（Boost、Eigen、Lapack 等）
@duffymo，我想我需要 Eigen 或 Boost 来进行矩阵转置、逆运算和乘法运算。如果我现在不需要这些东西，那我需要什么？
@Eagle 你不需要计算逆。
特征值问题不同于求解方程组。我将不得不回过头来回顾一下我对卡尔曼滤波的了解，以回忆特征值是否对其实施至关重要。（对不起，我的书在家里。）

【解决方案4】：

如果授权不是问题，你可以试试 gnu 科学库

它带有一个 blas 库，如果您以后需要，您可以换成一个优化的库（例如 intel、ATLAS 或 ACML（AMD 芯片）库。

【讨论】：

GSL 线性代数例程未优化。
@watson 那又怎样，你不需要花哨的 100x2 优化吗？
@watson 但是它提供了一个到底层 BLAS 库的接口？如果你发现你真的需要优化，你可以用链接而不是代码来交换你最喜欢的 BLAS 库
@David Heffernan。缓存没问题，因为 100x2 适合 L1，但使用 SSE，它可以比 uBLAS 快 4-8 倍（如果以单精度完成）。下载 Eigen 并亲自查看。

【解决方案5】：

如果您可以使用 MATLAB，我建议您使用它的 C 库。

【讨论】：

嗯，一个小问题的残酷解决方案！
AFAIK，Matlab C 库（至少是线性代数例程）使用/基于一些众所周知的公开可用的库 (LAPACK)。
@watson 所有的线性代数库本质上都是一样的，都来源于手册
@David Heffernan。希望你是在开玩笑。 MKL、Eigen、ATLAS 经过优化，可以更有效地利用高速缓存。库存 LAPACK（你指的是 LAPACK 手册？）仅提供（未优化的）参考实现和高级例程。
@David Heffernan。数学公式相同，但实现方式不同。调整内存层次结构的方式不同。 LAPACK 和 uBLAS 不做任何调整。