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

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..

• This looks like a stretch goal. Let's assume we use the fastest Skylake-X Xeon Platinum 8180 with a theoretical peak throughput of 4 single-precision TFLOPs, and that one 10x10 system requires about 700 (roughly 2n**3/3) floating-point operations to be solved. Then a batch of 1M such systems could theoretically be solved in 175 microseconds. That is a cannot-exceed speed-of-light number. Can you share what performance you are currently achieving with your fastest existing code? BTW, is the data single precision or double precision? – njuffa Mar 2 '19 at 9:07
• @njuffa yes I aimed to achieve close to 1ms but micro is another story. For micro I considered exploiting incremental inverse structure in the batch by detecting similar matrices, which occur often. Perf is curremtly at 10-500ms range depending on processor. Precision is double or even complex double. Single precision does slower. – rfabbri Mar 2 '19 at 10:32
• @njuffa I can reduce or up precision for speed – rfabbri Mar 2 '19 at 10:46
• It seems like precision/accuracy isn't your priority. For your goal, maybe an iterative method truncated at a relatively small number of evaluations is useful? Especially if you have a reasonable initial guess. – Spencer Bryngelson Mar 3 '19 at 5:37
• Do you pivot? Could you do a QR factorization instead of Gaussian elimination. Do you interleave you systems so that you can use SIMD instructions and do several systems at once? Do you write straight-line programs with no loops and no indirect addressing? What accuracy do you want and how I'll conditioned are your system? Do they have any structure which could be exploited. – Carl Christian Mar 3 '19 at 12:18

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.

• tx. I already use Eigen like Map<const Matrix<complex, 10, 10> > AA(A) successfully. will check into the other stuff. – rfabbri Mar 2 '19 at 1:06
• Eigen also has AVX and even a complex.h header for it. Why PETSc for this? It is hard to compete with Eigen in this case. I specialized Eigen even more for my problem and with an approximate pivot strategy that instead of taking max over a column, swaps a pivot immediately when it finds another that is 3 orders of magnitude bigger. – rfabbri Mar 8 '19 at 3:25
• @rfabbri I wasn't suggesting that you use PETSc for this, only that what they do in that particular instance could be instructive. I've edited the answer to make that clearer. – Daniel Shapero Mar 8 '19 at 16:52

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).

• since the matrices are so small, perhaps pivoting can be done away if they are pre-scaled. Not even pre-pivoting the matrices. All we need is that entries are within 2-3 orders of magnitude of each other. – rfabbri Mar 14 '19 at 3:45

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.

• So far Eigen already optimizes heavily, uses SEE, AVX, etc and I tried iterative methods in a preliminary test and they did not help. I tried Intel MKL but no better than Eigen with optimized GCC flags. I am currently trying to handcraft something better and simpler than Eigen and to do more detailed tests with iterative methods. – rfabbri Mar 7 '19 at 13:51

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).