计算行列式的最快方法是什么？

Question

我正在编写一个库，我希望在其中具有一些没有任何依赖关系的基本 NxN 矩阵功能，这是一个学习项目。我正在将我的表现与 Eigen 进行比较。我已经能够相当平等，甚至在 SSE2 和 AVX2 的几个方面击败了它的性能（它只使用 SSE2，所以并不令人惊讶）。

我的问题是我正在使用高斯消元法创建上对角化矩阵，然后乘以对角线以获得 N < 300 的 determinant.I 击败本征，但之后本征让我大吃一惊，情况变得更糟随着矩阵变大。考虑到所有内存都是按顺序访问的，而且编译器反汇编看起来并不糟糕，我不认为这是一个优化问题。

可以进行更多优化，但计时看起来更像是算法计时复杂性问题，或者我没有看到 SSE 的主要优势。简单地展开循环对我来说并没有太大帮助。

是否有更好的计算行列式的算法？

标量码

/*
    Warning: Creates Temporaries!
*/
template<typename T, int ROW, int COLUMN> MML_INLINE T matrix<T, ROW, COLUMN>::determinant(void) const
{
    /*
    This method assumes square matrix
    */
    assert(row() == col());
    /*
    We need to create a temporary
    */
    matrix<T, ROW, COLUMN> temp(*this);
    /*We convert the temporary to upper triangular form*/
    uint N = row();
    T det = T(1);
    for (uint c = 0; c < N; ++c)
    {
         det = det*temp(c,c);
        for (uint r = c + 1; r < N; ++r)
        {
            T ratio = temp(r, c) / temp(c, c);
            for (uint k = c; k < N; k++)
            {
                temp(r, k) = temp(r, k) - ratio * temp(c, k);
            }
        }
    }

    return det;
}

AVX2

template<> float matrix<float>::determinant(void) const
{
    /*
    This method assumes square matrix
    */
    assert(row() == col());
    /*
    We need to create a temporary
    */
    matrix<float> temp(*this);
    /*We convert the temporary to upper triangular form*/
    float det = 1.0f;

    const uint N = row();
    const uint Nm8 = N - 8;
    const uint Nm4 = N - 4;

    uint c = 0;
    for (; c < Nm8; ++c)
    {
        det *= temp(c, c);
        float8 Diagonal = _mm256_set1_ps(temp(c, c));

        for (uint r = c + 1; r < N;++r)
        {
            float8 ratio1 = _mm256_div_ps(_mm256_set1_ps(temp(r,c)), Diagonal);

            uint k = c + 1;
            for (; k < Nm8; k += 8)
            {
                float8 ref = _mm256_loadu_ps(temp._v + c*N + k);
                float8 r0 = _mm256_loadu_ps(temp._v + r*N + k);

                _mm256_storeu_ps(temp._v + r*N + k, _mm256_fmsub_ps(ratio1, ref, r0));
            }

            /*We go Scalar for the last few elements to handle non-multiples of 8*/
            for (; k < N; ++k)
            {
                _mm_store_ss(temp._v + index(r, k), _mm_sub_ss(_mm_set_ss(temp(r, k)), _mm_mul_ss(_mm256_castps256_ps128(ratio1),_mm_set_ss(temp(c, k)))));
            }
        }
    }

    for (; c < Nm4; ++c)
    {
        det *= temp(c, c);
        float4 Diagonal = _mm_set1_ps(temp(c, c));

        for (uint r = c + 1; r < N; ++r)
        {
            float4 ratio = _mm_div_ps(_mm_set1_ps(temp[r*N + c]), Diagonal);
            uint k = c + 1;
            for (; k < Nm4; k += 4)
            {
                float4 ref = _mm_loadu_ps(temp._v + c*N + k);
                float4 r0 = _mm_loadu_ps(temp._v + r*N + k);

                _mm_storeu_ps(temp._v + r*N + k, _mm_sub_ps(r0, _mm_mul_ps(ref, ratio)));
            }

            float fratio = _mm_cvtss_f32(ratio);
            for (; k < N; ++k)
            {
                temp(r, k) = temp(r, k) - fratio*temp(c, k);
            }
        }
    }

    for (; c < N; ++c)
    {
        det *= temp(c, c);
        float Diagonal = temp(c, c);
        for (uint r = c + 1; r < N; ++r)
        {
            float ratio = temp[r*N + c] / Diagonal;
            for (uint k = c+1; k < N;++k)
            {
                temp(r, k) = temp(r, k) - ratio*temp(c, k);
            }
        }
    }

    return det;
}

Answer 1

通过高斯消元将 n×n 矩阵简化为上（或下）三角形式的算法通常具有 O(n^3) 的复杂度（其中 ^ 表示 "to power of"）。

计算行列式有多种替代方法，例如评估特征值集（方阵的行列式等于其特征值的乘积）。对于一般矩阵，完整特征值集的计算实际上也是 O(n^3).

然而，理论上，特征值集的计算具有 n^w 的复杂性，其中 w 介于 2 和 2.376 之间 - 这意味着对于（大得多）更大的矩阵，它将比使用高斯消元法更快。查看 James Demmel、Ioana Dumitriu 和 Olga Holtz 在 Numerische Mathematik 中发表的文章 "Fast linear algebra is stable"，第 108 卷，第 1 期，第 59-91 页，2007 年 11 月。如果 Eigen 使用复杂度小于 O 的方法(n^3) 对于更大的矩阵（我不知道，从来没有理由调查这些事情）可以解释你的观察结果。

Answer 2

大多数地方的答案似乎使用 Block LU Factorization 在同一内存中创建下三角和上三角矩阵 space。它是 ~O(n^2.5) 取决于您使用的块的大小。

这里是莱斯大学的一个 power point，它解释了该算法。

www.caam.rice.edu/~timwar/MA471F03/Lecture24.ppt

除以矩阵意味着乘以它的逆矩阵。

这个想法似乎是显着增加 n^2 操作的数量，但减少 m^3 的数量，这实际上降低了算法的复杂性，因为 m 是一个固定的小尺寸。

需要花点时间以高效的方式编写它，因为要高效地完成它需要 'in place' 个我尚未编写的算法。

计算行列式的最快方法是什么？

What is the fastest way to calculate determinant?

c++

algorithm

algebra

matrix

linear-algebra