YNWA
YNWA

Reputation: 700

Is there any short way to solve a system Ax=B, when A is [3x3xN], i.e more then one matrix?

I want to try and implement it without a loop...

I have A, as A[3x3xN] or [3,3,N] and for each N it is a different matrix.

and B as [3x1xN] ofcourse...

How can I solve it without doing a loop and make A^-1 * B every time ?

Upvotes: 4

Views: 123

Answers (1)

petrichor
petrichor

Reputation: 6569

Z = cellfun(@(a,b) a\b, ... %# Solve for each pair
            num2cell(A,[1 2]), ... %# Make a cell array containing each slice
            num2cell(B,[1 2]), ... %# Make a cell array containing each slice
            'UniformOutput',false);
Z = cat(3,Z{:}); %# Merge the results to a 3x1xN array

Please see num2cell and cellfun documentation for further details on the functions used.

Let us compare its speed to a for loop:

clc, clear

N = 100000;
D = 10;
A = rand(D,D,N);
B = rand(D,1,N);

tic
Z = cellfun(@(a,b) a\b, ...
            num2cell(A,[1 2]),num2cell(B,[1 2]),'UniformOutput',false);
Z = cat(3,Z{:});
toc

tic
Z2 = zeros(D,1,N);
for i = 1:N
    Z2(:,:,i) = A(:,:,i) \ B(:,:,i);
end
toc

all(isequal(Z,Z2))

My results are as follows:

Elapsed time is 2.130507 seconds.
Elapsed time is 1.306873 seconds.

I tried it with different D values and got always similar ratios. Jonas' bet is correct!

Upvotes: 2

Related Questions