Solving a large number (thousands) of different "Ax=b" equations when the matrices A are all the same small size (3x3 to 6x6), but different numbers.

Hi all,

Wondering about the fastest way, within the Matlab programming environment, to solve a large number (in the thousands) of linear equation sets "Ax=b"; when all the sets are the same size, but the "A" matrices are different for each set.

Any insight will be appreciated.

Thanks,

Jahir Pabon

Answer 1

[Matt J]

I suspect that a for loop will be fastest, but you could also try combining them into a single set of equations using,

A=blkdiag(sparse(A1),sparse(A2),...) 
    b=[b1;b2;...]

and then

    x=A\b;

[Walter Roberson]

I do not know about algorithms for sparse matrices, but if I recall correctly, this kind of situation for "full" matrices should be detected as a band-limited matrix to trigger solvers that would be a bit more efficient than the general case.

[Jahir Pabon]

Thank you Matt. That seems like a very simple alternative to consider. I am just concerned that assembling the A matrix will require something like a for loop going for a few thousand iterations; and for loops typically result on considerable slowdown.

I will try it anyway.

Thanks,

Jahir

[Jahir Pabon]

To add to my question, what I have right now is a vectorized version of the gauss-elimination and back propagation that works on all the sets at the same time. That is, each element of the 6x6 matrix is really a vector. The same with "b". I notice there is significant speedup doing it this way than going through a for loop. But I don't know if there are faster ways to do this in Matlab.

[Walter Roberson]

If you build a cell array of the individual matrices, then you can use cell expansion with blkdiag() to build the A in one step.

In some cases it makes sense to arrayfun or cellfun to create the cell of A* values. However it is common that a for loop can be optimized more than arrayfun or cellfun.

[Jahir Pabon]

Thanks Walter. I need to become more familiar with cell arrays. Will look into it.

Jahir

Answer 2

[Walter Roberson]

Solutions up to 6 x 6 are feasible to calculate analytically.

a = sym('a', [6 6]);
B = sym('B', [6 1]);
tic; sol6 = a\B; time_for_symbolic_solution = toc

tic; solfun6no = matlabFunction(sol6, 'vars', {a, B}, 'file', 'sol6.m', 'optimize', false); time_to_create_unoptimized_function = toc

tic; solfun6op = matlabFunction(sol6, 'vars', {a, B}, 'file', 'sol6.m', 'optimize', true); time_to_create_optimized_function = toc

A = randn(6,6);
b = randn(6,1);
tic; xn = A\b; time_for_plain_numeric_calculation = toc

tic; xf = solfun6no(A, b); time_for_using_unoptimized_function = toc
tic; xf = solfun6op(A, b); time_for_using_optimized_function = toc

On my system, I got

time_for_symbolic_solution =

0.116915761

time_to_create_unoptimized_function =

1.470677184

time_to_create_optimized_function =

24.179678466

time_for_plain_numeric_calculation =

0.000126835

time_for_using_unoptimized_function =

0.167501003

time_for_using_optimized_function =

0.069751491

You might notice that this is significantly slower than using the high-performance numeric libraries!

However -- if I recall correctly, for the 3 x 3 case you can get higher performance using these kinds of techniques.

Also, the code generated by matlabFunction tossed in a lot of unnecessary overhead starting in R2021a (your release), so earlier releases might have higher performance.

[Jahir Pabon]

Thank you walter. I am looking for fast alternatives. I am trying to avoid using for loops to go over each set.

Jahir

[Jahir Pabon]

I meant "Walter"

Answer 3

[Matt J]

If you have the Parallel Computing Toolbox, you can use

pagefun(@mldivide,A,b)

https://www.mathworks.com/help/parallel-computing/pagefun.html

[Jahir Pabon]

Thanks Matt. I don't have that toolbox. But good to know the tip.

Jahir

Answer 4

[James Tursa]

See this FEX submission by Bruno Luong:

https://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver?s_tid=srchtitle

[Matt J]

Only for 2x2 and 3x3 systems, unfortuantely, but with Walter's answer you could probably derive a version for 6x6.

[Jahir Pabon]

Thanks for the link James. Will give Bruno's approach a try and compare speed with the vectorized version I have.

Matt, not sure why you write that it only works for 2x2 and 3x3 systems. It looks general enough to me. It is very much what you had suggested at the beginning. Just shows a way to build the sparse matrix I was not familiar with.

Thank you all,

Jahir

[Matt J]

Matt, not sure why you write that it only works for 2x2 and 3x3 systems.

I was thinking of this one, which does not use sparse algebra,

https://www.mathworks.com/matlabcentral/fileexchange/27762-small-size-linear-solver?s_tid=prof_contriblnk

Because it uses closed form solutions, I suspect it might be faster than what James posted, but it would need to be generalized to larger than 3x3.

A few additional nots:

(1) The small size solver no longer requires James' mtimesx package if you have R2020b or higher. You can use pagemtimes() instead.

(2) In the 4x4, 5x5 and 6x6 cases, if you happen to know that the upper-left or lower right 3x3 submatrix will always be invertible, you can use inv3 and Schur complements to build the higher order matrix inverses.