How to vectorise code involving matrix multiplications with vectors
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
KostasK
el 7 de Abr. de 2021
Comentada: Bruno Luong
el 7 de Abr. de 2021
Hi all,
I have a large number of
matrices and an equally large number of corresponding
vectors. For the
matrix and vector I have to solve the following linear sysem: ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/576582/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/576567/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/576572/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/576577/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/576582/image.png)
With that, I have been trying to find a way to vectorise this operation such that I can preferably avoid the for loop, and here's what I have come up with:
clear ; clc
n = 5 ; % Size of matrices and vectors
Nm = 1e3 ; % Number of matrices and vectors
for i = 1:Nm
A(:,:,i) = rand(n) ; % Generate matrices
b(:,i) = rand(5, 1) ; % Generate vectors
end
tic
% Option 1
for i = 1:Nm
v1(:,i) = A(:,:,i)\b(:,i) ;
end
toc1 = toc ;
tic
% Option 2
Ac = mat2cell(A, n, n, ones(Nm, 1)) ;
v2 = blkdiag(Ac{:})\b(:) ;
toc2 = toc ;
% Performance
fprintf('Option 1: %.3fs\nOption 2: %.3fs\n', toc1, toc2)
The performance of the for loop seems to be extremely fast, however I can't help but think that there must be an efficient way to vectorise the \ operation without using blkdiag, that is not really vectorising as it includes a for loop, and is incredibly slow.
With that said, is there a way to vectorise the above operation?
KMT.
PS. This is by no means a 'life or death' problem, if anything the for loop is clearly better and that should be the end of that. However its more out of stubborness that I think there must be way to vectorise this operation that I am not able to see, and thus I would like to ask for you help to see if this is true.
0 comentarios
Respuesta aceptada
James Tursa
el 7 de Abr. de 2021
1 comentario
Bruno Luong
el 7 de Abr. de 2021
Benchmark with the some methods in James's link
nA = 5;
mA = 5;
nY = 1;
nP = 1e3;
checkresult = false;
szA = [nA,mA,nP];
szY = [nA,nY,nP];
A = randn(szA);
y = randn(szY);
% https://www.mathworks.com/matlabcentral/fileexchange/68976-multipleqr
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)
try %#ok
% https://www.mathworks.com/matlabcentral/fileexchange/24260-multiple-same-size-linear-solver?s_tid=srchtitle
tic
xxx = SliceMultiSolver(A,y);
t3=toc;
fprintf('SliceMultiSolver time = %g [s]\n', t3)
end
tic
% Option 2
n = nA;
Nm = nP;
Ac = mat2cell(A, n, n, ones(Nm, 1)) ;
v2 = blkdiag(Ac{:})\y(:) ;
t4 = toc ;
fprintf('BlkDiag time = %g [s]\n', t4)
if checkresult
for k=[1,nP]
x(:,:,k)
xx(:,:,k)
try %#ok
xxx(:,:,k)
end
end
err = norm(x(:)-xx(:),Inf)/norm(xx(:),Inf)
if err >= 1e-6
fprintf('Bug detected\n');
keyboard
end
end
The tic/toc results on my system is:
>> TestMultipleQRSolve
size(A) = [5 5 1000]
size(y) = [5 1 1000]
MultipleQRSolve time = 0.0014875 [s]
Matlab loop time = 0.0112222 [s]
SliceMultiSolver time = 0.002807 [s]
BlkDiag time = 1.09904 [s]
Más respuestas (0)
Ver también
Categorías
Más información sobre Logical en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!