Is there any method to accelerate many small matrix and vector's "mldivide" (4*4)?

I am trying to solve many(4000000) mldivide evaluation. All of them are in the form of "x = A\b". A is a 4*4 matrix, and b is 4*1 vector. e.g.

num = 4000000;
A = rand(4,4,num);
b = rand(4,num);
x = zeros(4,num);
for i=1:num
    x(:,i) = A(:,:,i)\b(:,i);
end

Should I use the parfor? Or do it in the vectorized way? Any suggestion would be appreciated.

----------------------------------------------------------------------------------------------------------------------------

The following is my time profile of the answers. My computer is somehow old. only 2GZ cpu.

num = 4000000;
A = rand(4,4,num);
b = rand(4,num);
%% method1 : for loop (slowest)
tic
x1 = zeros(4,num);
for i=1:num
    x1(:,i) = A(:,:,i)\b(:,i);
end
toc % Elapsed time is 61.896299 seconds.
%% method2 : parfor with 10 workers
% parpool(10);
tic
x2 = zeros(4,num);
parfor i=1:num
    x2(:,i) = A(:,:,i)\b(:,i);
end
toc % Elapsed time is 6.148422 seconds.
%% method 3 (@Walter Roberson)
tic
syms tA [4 4]
syms tb [4 1];
tx = tA\tb;
X = matlabFunction(tx, 'vars', {[tA(:); tb(:)]});
A3 = reshape(A,16,[]);
x3 = X(reshape([A3;b], 20, []));
toc % Elapsed time is 2.351828 seconds.
%% method 4 (@ Bruno Luong)
tic
x4 = MultiSolver(A, b);
toc % Elapsed time is 9.759959 seconds.
%% the difference ratio between answers
max(max(abs((x1-x2)./x1))) % 0
max(max(abs((x1-x3)./x1))) % 5.9513e-09
max(max(abs((x1-x4)./x1))) % 2.9376e-10

If any problem or suggestions, please leave a message.

[Walter Roberson]

Note that once the function X had been created, it could be re-used; there is definitely set-up cost that could be amortized over repeated use.

Answer 1

[Walter Roberson]

syms tA [4 4]
syms tb [4 1];
tx = tA\tb;
X = matlabFunction(tx, 'vars', {[tA(:); tb(:)]});

x = X(reshape([a,b], 20, []));

This creates a vectorized function X, that accepts a 20 x N matrix and returns a 4 x N matrix.

[wei zhang]

I have a further problem when using your code, please help me.

I create a function like below.

function X = test1
syms tA [4 4]
syms tb [4 1]
tx = tA\tb ;
X = matlabFunction(tx, 'vars', {[tA(:); tb(:)]});
end

I always got a red wavy line below " tA\tb" in the function editor, which says "the variable "tA" ("tb") might be used before it is defined". Is it necessary to correct it? If yes, how to correct it?

[Walter Roberson]

function X = test1
tA = sym('tA', [4 4]);
tb = sym('tB', [4 1]);
tx = tA\tb ;
X = matlabFunction(tx, 'vars', {[tA(:); tb(:)]});
end

will shut up the warning.

Answer 2

[Bruno Luong]

Use my FEX file of MultiSolver

num = 4000000;
A = rand(4,4,num);
b = rand(4,num);
x = zeros(4,num);

tic
for i=1:num
    x(:,i) = A(:,:,i)\b(:,i);
end
toc % Elapsed time is 19.004386 seconds.

% https://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver
tic
X = MultiSolver(A, b);
toc % Elapsed time is 4.601745 seconds.

[wei zhang]

Thank you for you toolbox. But sorry I accepted the ealier answer. I collected your toolbox, which is smart to transfor small linear equations in to a huge sparse solver. I would like to use later. Thank you very much for your answer.

[Bruno Luong]

Just wonder how much tic/toc run on the example with Walter symbolic solution? I don't have the tbx to test.

[wei zhang]

I had submitted a time profile with answers. Please check it in the edited question part. If you have any suggestions, please let me know.

[Bruno Luong]

Thanks. that is useful.

Answer 3

[Bruno Luong]

I have submitted fast method FEX file MultipleQRSolver (C-compiler required) that can carried out the job in 1.5 second, 15 times faster than MATLAB for-loop (23 secs)

>> TestMultipleQRSolve
size(A) = [4 4 4000000]
size(y) = [4 1 4000000]
MultipleQRSolve time   = 1.47891 [s]
Matlab loop time       = 23.1515 [s]

The test cript is:

nA = 4;
mA = 4;
nY = 1;
nP = 4000000;

szA = [nA,mA,nP];
szY = [nA,nY,nP];
A = rand(szA);
y = rand(szY);

tic
x = MultipleQRSolve(A,y);
t1=toc;

tic
xx = zeros(size(x));
for k=1:nP
    xx(:,:,k) = A(:,:,k)\y(:,:,k);
end
t2=toc;

fprintf('size(A) = %s\n', mat2str(size(A)));
fprintf('size(y) = %s\n', mat2str(size(y))); 
fprintf('MultipleQRSolve time   = %g [s]\n', t1) 
fprintf('Matlab loop time       = %g [s]\n', t2)