尽管您可能没有意识到,您的问题可能更多是线性代数问题,而不是编程问题,尽管它确实具有编程维度。
在这个答案中,我假设您的实际问题可能比您在插图中提出的 N = 8 草图大得多。数学讨论是第一位的。示例代码如下。
线性代数问题的本质似乎是它们要么实际上可以用一般的线性代数技术来处理,要么不能。 L-U 分解和奇异值分解等通用线性代数技术自 20 世纪 60 年代以来已得到彻底发展,并具有坚实的理论基础。问题在于,这种通用技术在计算处理上往往是 O(N*N*N),也就是说,让你的向量长 10 倍需要 1000 倍的计算机处理时间。超出某个最大实际尺寸,您的问题在计算上变得难以解决——至少通过一般线性代数技术是这样。
如果您的问题非常严重
幸运的是,存在多种令人眼花缭乱的专用技术可以将 O(N*N*N) 问题简化为 O(N*N) 问题。这些技术大多自 20 世纪 60 年代以来出现,是一个活跃的研究领域。不幸的是,要知道应用哪种技术以及如何应用,需要对问题有深入的了解并且对矩阵数学有相当好的掌握。在 O(N*N*N) 下,通用技术是已知的。在 O(N*N) 中,没有通用技术,只有许多适用于有限类型系统的专用技术。如果您工作足够长的时间,您通常可以找到合适的限制,但它永远不会像将矩阵转储到某些通用求解器中那么简单。不是 O(N*N)。
我读过的关于这个主题的最好的书是 Henk A. van der Vorst 的大型线性系统的迭代 Krylov 方法。令人高兴的是,这本书的价格很合理。当然,在解决范德沃斯特问题之前,您需要在一般线性代数技术方面有坚实的基础。我推荐乔尔·N·富兰克林 (Joel N. Franklin) 1968 年出版的短书矩阵理论,多佛 1993 年推出廉价平装本。 (Van der Vorst 并没有提到 Kapp 和 Brown 的巧妙、简单且通常有效的有序多重交互方法,我多次发现该方法很有用,所以让我在这里提一下,以便您需要时可以查找它。披露:我个人认识布朗,因此可能会偏袒;但是,他和卡普的技术仍然是一项很好的技术,范德沃斯特确实忽略了这一点。)
由于我不完全理解的原因,大多数最近的技术(尽管不是卡普和布朗的)似乎更喜欢在实数中工作,将复数视为特殊情况或完全忽略它们。您没有提到您的数字是否复杂,但如果是,那么这可能会在一定程度上限制您的选择。
作为富兰克林和范德沃斯特之间的中间水平,如果您缺乏时间或兴趣来解决后者,您至少应该查找并熟悉共轭梯度方法。
如果您的问题只是中等程度
如果您的问题只是中等大小(例如,N
LAPACK 库及其低级助手 BLAS 快速、高效、彻底且正确,是该领域的标准。你应该使用它们。
LAPACK、BLAS 和 C++
我和我的同事都尝试过 LAPACK 和 BLAS 的各种 C 和 C++ 接口。由于各种原因,我们发现他们中没有一个是完全令人满意的,尽管他们确实试图提供帮助。我们每个人最终决定完全跳过接口并直接使用 LAPACK 和 BLAS。这就是我倾向于向您推荐的。
LAPACK 和 BLAS 应该从 FORTRAN 77 中调用,而不是从 C++ 中调用。尽管如此,从 C++ 调用它们并不是那么困难,而且您不需要 C++ 接口来完成它。事实上,您可能不想要 C++ 接口,至少不想要其他人准备的通用接口。如果你不把 LAPACK 和 BLAS 合二为一的话,他们会更高兴。 (请记住,FORTRAN 77 的链接功能虽然有效但有限,但缺乏 C 风格的头文件。)
直接使用LAPACK和BLAS首先要知道的是这个,不明白的话你会很惨:矩阵按列优先存储。也就是说,矩阵的每一列(而不是每一行)在内存中顺序表示为一个单元。因此,矩阵元素直接存储在其上方和下方的元素之间,而不是直接存储在其左侧和右侧的元素之间。事实上,LAPACK和BLAS以这种方式存储矩阵是谨慎的,因为发明C的Kernighan和Ritchie似乎对线性代数从来没有很感兴趣。 C 是一种很好的语言,但 C 的默认矩阵存储约定可能从一开始就是一个错误。从数学家的角度来看,LAPACK 和 BLAS 以应有的存储方式存储矩阵,即列优先;而且,如果您正在编写线性代数代码,那么您也应该以这种方式存储矩阵。忽略 C 的默认行优先约定:它与矩阵数学的约定无关,并且您不需要它。
这是一个完整的例子:
#include <iostream>
extern "C" {
void dgesv_(
const int *N,
const int *NRHS,
double *A,
const int *LDA,
int *IPIV,
double *B,
const int *LDB,
int *INFO
);
}
int main()
{
const int N = 2;
// Let A = | 3.0 6.1 |
// |-1.2 1.7 |
double A[N*N] = {3.0, -1.2, 6.1, 1.7};
// Let b = | 5.5 |
// | 0.4 |
double x[N] = {5.5, 0.4};
// (Why is b named x? Answer: because b and x share
// the same storage, because Lapack will write x over b.)
// Solve the linear system A*x = b.
const int NRHS = 1;
const int LDA = N;
int IPIV [N];
const int LDB = N;
int INFO;
dgesv_(&N, &NRHS, A, &LDA, IPIV, x, &LDB, &INFO);
// Uncomment the next line if you wish to see
// the error code Lapack returns.
//std::cout << "INFO == " << INFO << "\n";
// Output x.
for (int i = 0; i < N; ++i) std::cout << x[i] << "\n";
return 0;
}
【为什么函数名为dgesv_()?答案:不是,至少在 FORTRAN 77 中不是,它被命名为 DGESV。然而,FORTRAN 77 不区分大小写,当我们将其转换为 C 的链接约定时,我们得到 dgesv_()。至少,这是我和我的同事在我们尝试过的每台 Linux、BSD 或 OSX 机器上得到的结果。在 Debian 或 Ubuntu 上,shell 命令“readelf --symbols /usr/lib/liblapack.so | grep -i DGESV”会发现这一点。在 Microsoft 平台上,C 连接符号可能有其他名称:您必须调查这是否与您相关。]
LAPACK 和 BLAS 非常擅长他们的工作。你应该使用它们。
将数据存储在 C 样式数组(如示例所示)中还是 C++ 样式 std::valarrays 中取决于您。 LAPACK 和 BLAS 并不关心,只要存储是列优先的即可。
系数最小化
您没有提供足够的信息让我准确地告诉我您希望如何处理系数,但有人怀疑您可能需要的是 N×(2*N) 欠定系统的解决方案。这超出了富兰克林书的范围,但是,如果您已经吸收了富兰克林的材料,那么我自己关于伪逆的笔记 http://derivations.org/可能会帮助你。
显然,如果您正在寻找快速解决方案,我没有。可能有一个预装软件包可以完全满足您的要求,但我的经验表明,这样做的可能性对您不利。实践中会出现数百种不同类型的矩阵问题,预装软件无法破解。然而,LAPACK 和 BLAS 以及 Franklin、van der Vorst 以及您现在正在阅读的答案,提供了解决您的特定问题所需的所有工具。
祝你好运。