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

2 次查看（过去 30 天）
wei zhang
wei zhang 2020-9-23
编辑： Bruno Luong 2020-12-19
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 个评论
Walter Roberson
Walter Roberson 2020-9-23
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2020-9-23
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 个评论
wei zhang
wei zhang 2020-9-23
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 2020-9-23
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.

请先登录，再进行评论。

更多回答（2 个）

Bruno Luong
Bruno Luong 2020-9-23
编辑：Bruno Luong 2020-9-23
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 个评论
显示 2更早的评论
wei zhang
wei zhang 2020-9-23
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
Bruno Luong 2020-12-19
编辑：Bruno Luong 2020-12-19
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)

类别

在 Help Center 和 File Exchange 中查找有关 Code Performance 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by