How to vectorise code involving matrix multiplications with vectors

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:

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.

Answer 1

[James Tursa]

See this link:

https://www.mathworks.com/matlabcentral/answers/776812-how-to-vectorize-a-b-operation-on-slices-of-3d-matrices?s_tid=srchtitle

[Bruno Luong]

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]