fastest linear system solve for small square matrices (10x10)

Question

I am very interested in optimizing the hell out of linear system solving for small matrices (10x10), sometimes called tiny matrices. Is there a ready solution for this? The matrix can be assumed nonsingular.

This solver is to be executed in excess of 1 000 000 times in microseconds on an Intel CPU. I am talking to the level of optimization used in computer games. No matter if I code it in assembly and architecture-specific, or study precision or reliability tradeoffs reductions and use floating point hacks (I use the -ffast-math compile flag, no problem). The solve can even fail for about 20% of the time!

Eigen's partialPivLu is the fastest in my current benchmark, outperforming LAPACK when optimized with -O3 and a good compiler. But now I am at the point of handcrafting a custom linear solver. Any advice would be greatly appreciated. I will make my solution open source and I'll akcnowledge key insights in publications, etc..

Related: Speed of solving linear system with block diagonal matrix What is the fastest method to invert millions of matrices? https://stackoverflow.com/q/50909385/1489510

Answer 1

Using an Eigen matrix type where the number of rows and columns is encoded into the type at compile time gives you an edge over LAPACK, where the matrix size is known only at runtime. This extra information allows the compiler to do full or partial loop unrolling, eliminating lots of branch instructions. If you're looking at using an existing library rather than writing your own kernels, having a data type where the matrix size can be included as C++ template parameters will probably be essential. The only other library I know of that does this is blaze, so that might be worth benchmarking against Eigen.

If you decide to roll your own implementation, you might find what PETSc does for its block CSR format to be a useful example, although PETSc itself probably won't be the right tool for what you have in mind. Rather than write a loop, they write out every single operation for small matrix-vector multiplies explicitly (see this file in their repository). This guarantees that there are no branch instructions like you might get with a loop. The versions of the code with AVX instructions are a good example of how to actually use vector extensions. For example, this function uses the __m256d data type to simultaneously operate on four doubles at the same time. You could get an appreciable performance boost by explicitly writing out all the operations using vector extensions, only for LU factorization instead of matrix-vector multiply. Rather than actually write the C code by hand, you'd be better off using a script to generate it. It might also be fun to see if there's an appreciable performance difference when you reorder some of the operations to better take advantage of instruction pipelining.

You might also get some mileage out of the tool STOKE, which will randomly explore the space of possible program transformations to find a faster version.

Answer 2

Another idea could be to use a generative approach (a program writing a program). Author a (meta)program that spits out the sequence of C/C++ instructions to perform unpivoted** LU on a 10x10 system.. basically taking the k/i/j loop nest and flattening it into O(1000) or so lines of scalar arithmetic. Then feed that generated program into whichever optimizing compiler. What I think is sort of interesting here, is removing the loops exposes every data dependency and redundant subexpression, and gives the compiler maximum opportunity to reorder instructions so that they map well to actual hardware (eg number of execution units, hazards/stalls, so on).

If you happen to know all of the matrices (or even just a few of them), you can improve throughput by calling SIMD intrinsics/functions (SSE/AVX) instead of scalar code. Here you'd be exploiting the embarrassing parallelism across the instances, instead of chasing any parallelism within a single instance. For instance, you could perform 4 double precision LU's simultaneously using AVX256 intrinsics, by packing 4 matrices "across" the register and doing the same operations** on all of them.

** Hence the focus on unpivoted LU. Pivoting spoils this approach in two ways. First, it introduces branches due to pivot selection, meaning your data dependencies are not as perfectly known. Second, it means that different SIMD "slots" would have to do different things, because instance A might pivot differently than instance B. So if you pursue any of this, I'd suggest statically pivoting your matrices prior to computation (permute largest entry of each column to diagonal).

Answer 3

Your question leads to two different considerations.

First, you need to pick the right algorithm. Hence, the question if the matrices have any structure, should be considered. E.g., when the matrices are symmetric, a Cholesky decomposition is more efficient than LU. When you only need a limited amount of accuracy an iterative method can be faster.

Second, you need to implement the algorithm efficiently. To do so, you need to know the bottleneck of your algorithm. Is your implementation bound by the speed of the memory transfer or by the speed of the computation. Since you consider only $10 \times 10$ matrices, your matrix should fit into the CPU cache completely. Thus, you should make use of the SIMD units (SSE, AVX, etc.) and cores of your processor, to do as many computations per cycle as possible.

In all, the answer to your question heavily depends on the hardware and matrices that you consider. There probably is no definite answer and you have to try out a few things to find an optimal method.

Answer 4

I would try blockwise inversion.

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen uses an optimized routine to calculate the inverse of a 4x4 matrix, which is probably the best you're going to get. Try using that as much as possible.

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

Top left: 8x8. Top right: 8x2. Bottom left: 2x8. Bottom right: 2x2. Invert the 8x8 using the optimized 4x4 inversion code. The rest is matrix products.

EDIT: Using 6x6, 6x4, 4x6, and 4x4 blocks has shown to be a bit faster than what I described above.

using namespace Eigen;

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
    Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;

    Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
    Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<br_size, br_size>() = DCAB_inv;

    return result;
}

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
    const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
    const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
    const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
    const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();

    return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
    const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
    const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
    const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();

    Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
    Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<2, 2>() = DCAB_inv;

    return result;
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
    const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
    const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
    const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();

    Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
    Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();

    result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<4, 4>() = DCAB_inv;

    return result;
}

Here are the results of one bench mark run using one million Eigen::Matrix<double,10,10>::Random() matrices and Eigen::Matrix<double,10,1>::Random() vectors. In all my tests, my inverse is always faster. My solve routine involves computing the inverse and then multiplying it by a vector. Sometimes its faster than Eigen, sometimes its not. My bench marking method may be flawed (didn't disable turbo boost, etc). Also, Eigen's random functions may not representative of real data.

• Eigen partial pivot inverse: 3036 milliseconds
• My inverse with 8x8 upper block: 1638 milliseconds
• My inverse with 6x6 upper block: 1234 milliseconds
• Eigen partial pivot solve: 1791 milliseconds
• My solve with 8x8 upper block: 1739 milliseconds
• My solve with 6x6 upper block: 1286 milliseconds

I am very interested to see if anyone can optimize this further, as I have a finite element application which inverts a gazillion 10x10 matrices (and yes, I do need individual coefficients of the inverse so directly solving a linear system is not always an option).