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

2 views (last 30 days)
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.
  1 Comment
Walter Roberson
Walter Roberson on 23 Sep 2020
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.

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 23 Sep 2020
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.
  2 Comments
wei zhang
wei zhang on 23 Sep 2020
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
Walter Roberson on 23 Sep 2020
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.

Sign in to comment.

More Answers (2)

Bruno Luong
Bruno Luong on 23 Sep 2020
Edited: Bruno Luong on 23 Sep 2020
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.
  4 Comments
wei zhang
wei zhang on 23 Sep 2020
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.

Sign in to comment.


Bruno Luong
Bruno Luong on 19 Dec 2020
Edited: Bruno Luong on 19 Dec 2020
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)

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!