CODEWITHSUNDEEP
matlabmultidimensional-arrayvectorization
shashashamti2008
shashashamti2008

Reputation: 2327

Need help in using bsxfun

I have two arrays in MATLAB:

A; % size(A) = [NX NY NZ 3 3]
b; % size(b) = [NX NY NZ 3 1]

In fact, in the three dimensional domain, I have two arrays defined for each (i, j, k) which are obtained from above-mentioned arrays A and b, respectively and their sizes are [3 3] and [3 1], respectively. Let's for the sake of example, call these arrays m and n.

m; % size(m) = [3 3]
n; % size(n) = [3 1]

How can I solve m\n for each point of the domain in a vectorize fashion? I used bsxfun but I am not successful.

solution = bsxfun( @(A,b) A\b, A, b );

I think the problem is with the expansion of the singleton elements and I don't know how to fix it.

Upvotes: 3

Views: 128

Answers (1)

Daniel
Daniel

Reputation: 36710

I tried some solutions, it seems that a for loop is acutally the fastest possibility in this case.

A naive approach looks like this:

%iterate
C=zeros(size(B));
for a=1:size(A,1)
    for b=1:size(A,2)
        for c=1:size(A,3)
            C(a,b,c,:)=squeeze(A(a,b,c,:,:))\squeeze(B(a,b,c,:));
        end
    end
end

The squeeze is expensive in computation time, because it needs some advanced indexing. Swapping the dimensions instead is faster.

A=permute(A,[4,5,1,2,3]);
B=permute(B,[4,1,2,3]);
C2=zeros(size(B));
for a=1:size(A,3)
    for b=1:size(A,4)
        for c=1:size(A,5)
            C2(:,a,b,c)=(A(:,:,a,b,c))\(B(:,a,b,c));
        end
    end
end
C2=permute(C2,[2,3,4,1]);

The second solution is about 5 times faster.

/Update: I found an improved version. Reshaping and using only one large loop increases the speed again. This version is also suitable to be used with the parallel computing toolbox, in case you own it replace the for with a parfor and start the workers.

A=permute(A,[4,5,1,2,3]);
B=permute(B,[4,1,2,3]);
%linearize A and B to get a better performance
linA=reshape(A,[size(A,1),size(A,2),size(A,3)*size(A,4)*size(A,5)]);
linB=reshape(B,[size(B,1),size(B,2)*size(B,3)*size(B,4)]);
C3=zeros(size(linB));
for a=1:size(linA,3)
    C3(:,a)=(linA(:,:,a))\(linB(:,a));
end
%undo linearization
C3=reshape(C3,size(B));
%undo dimension swap
C3=permute(C3,[2,3,4,1]);

Upvotes: 2

Related Questions