使拉普拉斯展开更高效

Question

我创建了一个能够在 C++ 中计算矩阵行列式的小程序。我使用了 laplace-expansion，虽然我知道有更有效的方法来做到这一点：

double getDeterminantLaplace(const std::vector<std::vector<double>> vect) {

    int dimension = vect.size();

    if(dimension == 0) {
        return 1;
    }

    if(dimension == 1) {
        return vect[0][0];
    }

    //Formula for 2x2-matrix
    if(dimension == 2) {
        return vect[0][0] * vect[1][1] - vect[0][1] * vect[1][0];
    }

    double result = 0;
    int sign = 1;
    for(int i = 0; i < dimension; i++) {

        //Submatrix
        std::vector<std::vector<double>> subVect(dimension - 1, std::vector<double> (dimension - 1));
        for(int m = 1; m < dimension; m++) {
            int z = 0;
            for(int n = 0; n < dimension; n++) {
                if(n != i) {
                    subVect[m-1][z] = vect[m][n];
                    z++;
                }
            }
        }

        //recursive call
        result = result + sign * vect[0][i] * getDeterminantLaplace(subVect);
        sign = -sign;
    }

    return result;
}

我现在的问题是：如何让这个算法更高效？

我的一个想法是不创建“子矩阵”而只使用原始矩阵，但我真的不知道该怎么做。你怎么看这个想法？如何在 C++ 中做到这一点？

你还有什么想法吗？

Answer 1

首先，微不足道的优化是在当前元素为零时不递归。这将使您在稀疏矩阵上立即加速。

下一个优化是您已经建议的：不要创建所有子矩阵。您可以通过创建索引向量来做到这一点。例如，如果您的原始矩阵有 4×4 元素，则使用以下索引向量进行递归：

0: {1, 2, 3}
1: {0, 2, 3}
2: {0, 1, 3}
3: {0, 1, 2}

您不需要每次都从头开始创建索引向量：从作为当前向量的子向量开始，然后用当前子向量的第 i 个条目覆盖第 i 个位置。

当您访问子矩阵的元素s[r][c] 时，访问原始矩阵的元素a[r + top][col[c]]。可以根据当前列向量的维度和原始矩阵的维度来确定顶行的索引。

您永远不会创建子矩阵，只会创建子列向量。将您的函数一分为二：一个作为前端的公共函数，它调用递归工作函数。

这将在一定程度上加快计算速度，但不幸的是，当您的矩阵增长时，这种改进不会给您带来太多收益。让我们再看一下 4×4 矩阵。在第一个递归步骤中，您将考虑这些 3×3 子矩阵：

1, 2, 3     0, 2, 3     0, 1, 3     0, 1, 2

从那里，您将计算这些 2×2 子矩阵：

   2, 3        2, 3        1, 3        1, 2
1,    3     0,    3     0,    3     0,    2
1, 2        0, 2,       0, 1,       0, 1,  

请注意，这 12 个索引实际上只是 6 个不同的对。您将对它们中的每一个进行两次计算。你的原始矩阵越大，这会变得更糟。一个解决方案是记忆：一旦你计算了某个子矩阵的行列式，将值存储在一个关联的数组中。在计算子矩阵之前，请检查您是否已经这样做了，如果是，则返回您之前计算的值。

这会加速你的函数，但它是有代价的：它会在相关的数组中创建许多条目。

无论如何，这是实现我所描述的所有优化的代码：

#include <vector>
#include <map>
#include <iostream>

double subdet(const std::vector<std::vector<double> > &a,
    const std::vector<int> &col,
    std::map<std::vector<int>, double> &memo)
{
    int dim = col.size();
    int top = a.size() - dim;

    if (memo.find(col) != memo.end()) {
        return memo[col];
    }   

    if (dim == 2) return a[top + 0][col[0]] * a[top + 1][col[1]]
                       - a[top + 0][col[1]] * a[top + 1][col[0]];

    double result = 0.0;
    int sign = 1;

    std::vector<int> ncol(&col[1], &col[dim]);

    for (int i = 0; i < dim; i++) {
        if (a[top][col[i]]) {
            double d = subdet(a, ncol, memo);

            result = result + sign * a[top][col[i]] * d;
        }

        sign = -sign;

        if (i + 1 < dim) ncol[i] = col[i];
    }

    memo[col] = result;
    return result;
}

double det(const std::vector<std::vector<double> > a)
{
    int dim = a.size();

    if (dim == 0) return 1.0;
    if (dim == 1) return a[0][0];

    std::vector<int> col(dim);
    std::map<std::vector<int>, double> memo;

    for (unsigned i = 0; i < a.size(); i++) col[i] = i;

    return subdet(a, col, memo);
}

注意：映射（具有 O(log n) 查找的二叉树）实际上应该是未排序的映射（具有 O(1) 查找的哈希表），但我无法得到它可以工作，因为我不擅长 C++。对此感到抱歉。

可能还有优化查找键的空间：可以枚举可能的索引向量或使用位掩码，也许，从而节省哈希映射中的内存。它不是对列索引向量的良好字符串引用，因为它是短暂的，我们经常在其中交换，所以它不是恒定的。

当然，其他算法更适合寻找大矩阵的行列式。我的回答侧重于改进现有方法。