What is the fastest way to calculate determinant?

Question

I'm writing a library where I want to have some basic NxN matrix functionality that doesn't have any dependencies and it is a bit of a learning project. I'm comparing my performance to Eigen. I've been able to be pretty equal and even beat its performance on a couple front with SSE2 and with AVX2 beat it on quite a few fronts (it only uses SSE2 so not super surprising).

My issue is I'm using Gaussian Elimination to create an Upper Diagonalized matrix then multiplying the diagonal to get the determinant.I beat Eigen for N < 300 but after that Eigen blows me away and it just gets worse as the matrices get bigger. Given all the memory is accessed sequentially and the compiler dissassembly doesn't look terrible I don't think it is an optimization issue.

There is more optimization that can be done but the timings look much more like an algorithmic timing complexity issue or there is a major SSE advantage I'm not seeing. Simply unrolling the loops a bit hasn't done much for me when trying that.

Is there a better algorithm for calculating determinants?

Scalar code

/*
    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

Algorithms to reduce an n by n matrix to upper (or lower) triangular form by Gaussian elimination generally have complexity of O(n^3) (where ^ represents "to power of").

There are alternative approaches for computing determinate, such as evaluating the set of eigenvalues (the determinant of a square matrix is equal to the product of its eigenvalues). For general matrices, computation of the complete set of eigenvalues is also - practically - O(n^3).

In theory, however, calculation of the set of eigenvalues has complexity of n^w where w is between 2 and 2.376 - which means for (much) larger matrices it will be faster than using Gaussian elimination. Have a look at an article "Fast linear algebra is stable" by James Demmel, Ioana Dumitriu, and Olga Holtz in Numerische Mathematik, Volume 108, Issue 1, pp. 59-91, November 2007. If Eigen uses an approach with complexity less than O(n^3) for larger matrices (I don't know, never having had reason to investigate such things) that would explain your observations.

Answer 2

The answer most places seem to use Block LU Factorization to create an Lower triangle and Upper triangle matrix in the same memory space. It is ~O(n^2.5) depending on the size of block you use.

Here is a power point from Rice University that explains the algorithm.

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

Division by a matrix means multiplication by its inverse.

The idea seems to be to increase the number of n^2 operations significantly but reduce the number m^3 which in effect lowers the complexity of the algorithm since m is of a fixed small size.

Going to take a little bit to write this up in an efficient manner since to do it efficiently requires 'in place' algorithms I don't have written yet.