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

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) 的复杂度(其中 ^ 表示“的幂”)。

计算行列式有多种替代方法，例如评估特征值集(方阵的行列式等于其特征值的乘积)。对于一般矩阵，计算完整的特征值集实际上也是 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) 更大的矩阵(我不知道，从来没有理由调查这样的事情)可以解释你的观察。