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

Question

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

我的问题是我正在使用高斯消元法创建一个上对角化矩阵，然后乘以对角线以获得行列式。我在 N < 300 上击败了 Eigen，但在那之后 Eigen 让我大吃一惊，随着矩阵变大它变得更糟。鉴于所有内存都是按顺序访问的，并且编译器反汇编看起来并不可怕，我认为这不是优化问题。

可以进行更多优化，但时间看起来更像是算法时间复杂性问题，或者我没有看到主要的 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，第 108 卷，第 1 期，第 59-91 页，2007 年 11 月发表的一篇文章“快速线性代数是稳定的”。如果 Eigen 使用复杂度较低的方法比 O(n^3) 更大的矩阵（我不知道，从来没有理由调查这样的事情）可以解释你的观察。