Speed of solving linear system with block diagonal matrix

Question

I have a bunch of 3x3 linear systems of the form $Ax=b$. In general, would it be faster to solve each individual system, or to formulate it as a giant block diagonal system and solve that?

I expect to have a million or so systems to solve. I also was planning on using the Eigen C++ library, if that makes any differences, speed-wise.

Answer 1

The numerical costs of the problem itself are linearly dependent on the number of rows. You reach these minimal costs if you solve the blocks separately. The usual matrix solvers that do not exploit the internal structure of the system cause numerical costs that grow with the third power of the number of rows. So, they are more expensive.

But, there are general solvers that detect the matrix structure and then solve the system with linear costs. The detection of the matrix structure has costs that are quadratic in the number of rows but have a very low impact since the required time for testing one matrix element is very low. We are using such a solver for SimulationX. With this solver the time for solving the linear system is rarely the limiting factor in the simulation.

Answer 2

If there's only a million of them, you should be able to easily solve all of them in about 1 second on a decent laptop or desktop, so performance shouldn't even be a concern (unless 1 second is still too slow for your intended application).

To illustrate this, consider the following Mathematica code, which solves 1 million random $3\times 3$ systems:

First@Timing@Table[LinearSolve[RandomReal[{0, 1}, {3, 3}]], {1000000}]

On my ancient 7-year-old budget laptop using a single CPU core, this takes 30 seconds. With 4-6 parallel CPU cores on a modern desktop computer, this should take about 1 second. In CUDA with a GPU exploiting the parallelism, you could probably do it in milliseconds.

Answer 3

If you know something about a matrix, it almost always pays off to use it. In your case, use the fact that you know that it is composed of many $3\times 3$ blocks and solve them individually -- which you can do in parallel, for example. If you just give a large block diagonal matrix to an iterative or sparse direct solver, the best they could do is to analyze the matrix and do exactly what you could do yourself, but this analysis step is (i) not usually done, at least in iterative solvers, (ii) expensive.

Answer 4

I think it's best to solve the individual 3 x 3 systems. The problems of doing such a large matrix, containing millions of 3 x 3 systems:

1. computational cost of assembly: assembling such a large matrix could be take a substantial amount of time (even with a sparse solver).
2. parallelization overheads: should there be a need to parallelize the code, the large matrix would have to be distributed over many processors and there would be communication overheads for passing information between processors. On the other hand, for individual 3 x 3 systems, the coefficients could easily be distributed over the various processors, after which it would be an "embarrassingly parallel" operation (i.e. no communication between processors would be necessary). Furthermore, parallelization using OpenMP or MPI would not be too difficult. On the other hand, for the large matrix, it would probably be easier to make use of a library (such as PETSc: e.g. the example ex2.c is a nice start), rather than to write your own code.

Lastly, I would like to point out that for 3 x 3 systems, it's easy to work out an analytical formula using computer algebra systems (or using Wolfram Alpha). As an illustration, here's how it can be done with the Python package, sympy:

from sympy import *
x1,x2,x3=symbols('x1 x2 x3')
b1,b2,b3=symbols('b1 b2 b3')
a11,a12,a13=symbols('a11 a12 a13')
a21,a22,a23=symbols('a21 a22 a23')
a31,a32,a33=symbols('a31 a32 a33')
eq1=Eq(a11*x1+a12*x2+a13*x3,b1)
eq2=Eq(a21*x1+a22*x2+a23*x3,b2)
eq3=Eq(a31*x1+a32*x2+a33*x3,b3)
x=solve([eq1,eq2,eq3], [x1,x2,x3])
x.get(x1)
x.get(x2)
x.get(x3)

Just replace x1, x2, x3, b1, b2, b3, etc. with the actual variables in your existing code. If generation of C code is required:

[(c_name, c_code), (h_name, c_header)] = codegen(("x1",x.get(x1)), "C", "test", header=False, empty=False)
print(c_code)
[(c_name, c_code), (h_name, c_header)] = codegen(("x2",x.get(x2)), "C", "test", header=False, empty=False)
print(c_code)
[(c_name, c_code), (h_name, c_header)] = codegen(("x3",x.get(x3)), "C", "test", header=False, empty=False)
print(c_code)

Sample output (for just x1):

>>> print(c_code)
#include "test.h"
#include <math.h>
double x1(double a11, double a12, double a13, double a21, double a22, double a23, double a31, double a32, double a33, double b1, double b2, double b3) {
   return (a12*a23*b3 - a12*a33*b2 - a13*a22*b3 + a13*a32*b2 + a22*a33*b1 - a23*a32*b1)/(a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);
}

And here it may make sense to pre-calculate the determinant to further reduce the computational cost, i.e.:

determinant = (a11*a22*a33 - a11*a23*a32 - a12*a21*a33 + a12*a23*a31 + a13*a21*a32 - a13*a22*a31);

Answer 5

I agree with Mr. DumpsterDoofus's comments.

On a fairly modest 5-year-old CPU (AMD Phenom II X4 830), this Eigen code:

Eigen::PartialPivLU< Eigen::Matrix3d > lu(a);
Eigen::Vector3d u = lu.solve(b);

requires about .37 seconds to execute one million times.