5

我一直在寻找 C++ 中的库来解决像这样的未确定系统

q 是一个向量，w, x, y, z 变量和 a,b,c,d 常数。

argmin_q MAX(q) - MIN(q)

英石

q[1] = a - w - y
q[2] = b - w - z
q[3] = c - x - y
q[4] = d - x - z

找到求解器、算法等会非常有用。我发现了几个库能够解决未确定的系统，但另外我需要最小化系数之间的距离。

先感谢您

阿尔托伯

4

1 回答 1

3

尽管您可能没有意识到，但您的问题可能更像是线性代数问题而不是编程问题，尽管它确实具有编程维度。

在这个答案中，我假设您的实际问题可能比您在插图中提出的 N = 8 草图大得多。首先是数学讨论。示例代码如下。

线性代数问题的本质似乎是它们要么实际上可以通过一般线性代数技术处理，要么不能。LU分解和奇异值分解等一般线性代数技术自1960年代以来得到了彻底的发展，并具有坚实的理论基础。问题在于，这种通用技术在计算处理中倾向于 O(N*N*N)，也就是说，使你的向量长 10 倍需要 1000 倍的计算机时间来处理。超出一些最大实际大小，您的问题在计算上变得难以处理 - 至少通过一般线性代数技术。

如果您的问题非常大

幸运的是，有许多令人眼花缭乱的专用技术可以将 O(N*N*N) 问题简化为 O(N*N) 问题。这些技术大多出现在 1960 年代，是一个活跃的研究领域。不幸的是，要知道应用哪种技术以及如何应用，需要对问题有深入的了解，并且对矩阵数学有相当好的把握。在 O(N*N*N) 时，一般技术是已知的。在 O(N*N) 时，没有通用技术，只有许多适用于受限类型系统的专用技术。如果您的工作时间足够长，您通常可以找到一个合适的限制，但它绝不像将矩阵转储到某个通用求解器中那么简单。不在 O(N*N) 处。

我读过的关于这个主题的最好的书是 Henk A. van der Vorst 的Iterative Krylov Methods for Large Linear Systems。 令人高兴的是，这本书的价格是合理的。当然，在处理 van der Voorst 之前，您需要对一般线性代数技术有坚实的基础。我建议乔尔·N·富兰克林 1968 年的短书《矩阵理论》，该书有 1993 年多佛的廉价平装本。（Van der Vorst 没有碰巧提到 Kapp 和 Brown 的聪明、简单且经常有效的有序多重交互方法，我反复发现它很有用，所以我在这里提一下，以便您在需要时查找。披露：我个人认识布朗，因此可能会偏袒；但是，仍然是他和卡普

由于我不完全理解的原因，大多数最近的技术（尽管不是 Kapp 和 Brown 的）似乎更喜欢在实数中工作，将复数视为特殊情况或完全忽略它们。您没有提到您的数字是否复杂，但如果是，那么这可能会在一定程度上限制您的选择。

作为富兰克林和范德沃斯特之间的中间水平，如果你没有时间或兴趣去解决后者，你至少应该查一下并熟悉一下共轭梯度的方法。

如果您的问题只是比较大

如果你的问题只是中等大——比如，N < 20,000 左右——那么你的任务就容易多了。在这种情况下忘记范德沃斯特：你不需要他。根据您是否希望以毫秒、秒、分钟或小时为单位（您没有提到哪个，但它只影响 N 的实际限制），您可以容忍 O(N*N*N) 性能，并且在这个案例富兰克林从 1960 年代开始的一般技术非常强大。

快速、高效、彻底和正确的 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) 欠定系统的解决方案。这超出了富兰克林书的范围，但是，如果你已经吸收了富兰克林的材料，那么我自己关于伪逆的笔记可能会对你有所帮助。

显然，如果您正在寻找快速解决方案，我没有。可能有一个预装的软件包可以完全满足您的需求，但我的经验表明，这对您不利。在实践中出现了数百种不同类型的矩阵问题，预装软件无法破解。但是，LAPACK 和 BLAS 以及 Franklin、van der Vorst 和您现在正在阅读的答案提供了解决特定问题所需的所有工具。

祝你好运。

于 2012-03-17T15:20:44.043 回答