【问题标题】:C++ library for undetermined equation systems [closed]用于未定方程系统的 C++ 库
【发布时间】:2011-10-17 20:18:38
【问题描述】:

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

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

argmin_q MAX(q) - MIN(q)

s.t.

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

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

提前谢谢你

阿尔托伯

【问题讨论】:

    标签: c++ solver equation-solving


    【解决方案1】:

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

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

    线性代数问题的本质似乎是它们要么可以通过一般线性代数技术实际处理,要么不能。 L-U因式分解和奇异值分解等一般线性代数技术自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 之前,您需要对一般线性代数技术有坚实的基础。我推荐 Joel N. Franklin 1968 年的短书Matrix Theory,在 1993 年 Dover 的廉价平装本中有售。 (Van der Vorst 没有碰巧提到 Kapp 和 Brown 的聪明、简单且经常有效的有序多重交互方法,我反复发现它很有用,所以我在这里提一下,以便您在需要时查找。披露:我个人认识布朗,因此可能会偏袒;但是,他和卡普的技术仍然很好,范德沃斯特确实省略了。)

    由于我不完全理解的原因,大多数最近的技术(尽管不是 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) 欠定系统的解决方案。这超出了富兰克林的书的范围,但是,如果你已经吸收了富兰克林的材料,那么my own notes on the pseudoinverse 可能会帮助你。

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

    祝你好运。

    【讨论】:

      猜你喜欢
      • 2011-12-30
      • 2013-07-26
      • 1970-01-01
      • 2023-03-09
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多