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

Posted 2023-02-22

【中文标题】计算行列式的最快方法是啥？

【英文标题】：What is the fastest way to calculate determinant?计算行列式的最快方法是什么？

【发布时间】：2016-03-28 00:18:29

【问题描述】：

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

我的问题是我正在使用高斯消元法创建一个上对角矩阵，然后将对角线相乘以获得行列式。我在 N

可以进行更多优化，但时序看起来更像是算法时序复杂性问题，或者我没有看到主要的 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;

【问题讨论】：

鉴于您显然已准备好深入研究细节，并且 Eigen 是开源的...why not look at what it does 或逐步完成...？

他们的方法对我来说没有多大意义，因此很难适应我的图书馆的运作方式。如果我理解它背后的数学推理，我就能很容易地适应它。我认为这与我开始研究的部分旋转有关。其他方法对我来说很有意义，但这是我第一次无法理解其背后的方法。一般来说，只是好奇另一个大脑是否有“最好的方法”的想法。即使在发布了这个问题之后，我仍然非常关注它，当我找到更好的方法时会发布我的代码。

this 你可能会感兴趣。

我认为您需要进行旋转以处理 (0 1; 1 0) 等具有行列式 -1 的矩阵，但我认为您的方法将在该矩阵上失败。

是的，它会的，我正在努力弄清楚我想在这样做之前使用什么算法。

【参考方案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) 更大的矩阵（我不知道，从来没有理由调查这些事情）可以解释你的观察。

【讨论】：

感谢这篇论文很有趣。我认为“块 LU 分解”是 Eigen 用于 8x8 子块的方法。那篇论文说时间复杂度约为 O(n^2.5)。我也会更多地研究特征值选项。

【参考方案2】：

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

这是莱斯大学的一个简报，用于解释该算法。

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

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

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

要花一点时间以一种有效的方式写出来，因为要有效地写出来需要我还没有写过的“就地”算法。