Triple loop into tensor matlab

Question

I am aware of similar questions, but couldn't find a useful one for my problem

Is there a way to speed up my code by transforming this triple loop into matrix/tensor operations?

% preallocate
X_tensor = zeros(N, N, N);

% loop
for i = 1: N
    for k = 1: N
        for j = 1:N
            X_tensor(i,k,j) = X(i,j) * B(i,k) * Btilde(k,j)  / B(i,j);
        end
    end
end

EDIT I am sorry, I forgot one important information: X, B and Btilde are all NxN matrices

The first index is the one that changes most rapidly. Therefore, you should switch the order of the loops (i should be the inner loop, and j the outer loop). You can probably also vectorize the operation, which might or might not improve the speed. — Cris Luengo
– Cris Luengo, Commented Sep 2, 2022 at 14:51

beaker · Accepted Answer · 2022-09-02 15:15:04Z

2

You can vectorize the operation by permuting the matrices and using element-wise multiplication/division.

i, k, and j are dimensions [1,2,3]. Looking at the first term in the equation, that means we want to take the [1,2,3] order of X (where the third dimension is a singleton) and rearrange it to [1,3,2]. This corresponds to an index of [i,j,k]. B in the second term is already in the required [1,2,3] order.

So the loops can be replaced by:

X_tensor2 = permute(X,[1,3,2]) .* B .* permute(Btilde,[3,1,2]) ./ permute(B,[1,3,2])

Here's a test program to verify correctness:

N = 5;
P = primes(400);
X = reshape(P(1:25),N,N);
B = reshape(P(26:50),N,N);
Btilde = reshape(P(51:75),N,N);

% preallocate
X_tensor = zeros(N, N, N);
X_tensor2 = zeros(N, N, N);

% loop
for i = 1: N
    for k = 1: N
        for j = 1:N
            X_tensor(i,k,j) = X(i,j) * B(i,k) * Btilde(k,j)  / B(i,j);
        end
    end
end

X_tensor2 = permute(X,[1,3,2]) .* B .* permute(Btilde,[3,1,2]) ./ permute(B,[1,3,2]);

isequal(X_tensor, X_tensor2)

answered Sep 2, 2022 at 15:15

beaker

16.9k3 gold badges36 silver badges53 bronze badges

Sign up to request clarification or add additional context in comments.

7 Comments

Ander Biguri Over a year ago

Great answer! Is it actually faster?

beaker Over a year ago

@AnderBiguri In Octave with N=100, it's about 13.0s vs. 0.005s. But I would expect the difference to be much less pronounced with MATLAB and JIT.

jack Over a year ago

@AnderBiguri thank you! I am sorry I forgot the dimension of the right hand side, does that change the approach? You created NxN matrices so it shouldn't, right? Just to be sure

beaker Over a year ago

@jack I assumed that they were NxN matrices because that's how you were indexing them in your loops. You should be able to verify that both approaches produce the same output (within floating point rounding errors) with your data.

jack Over a year ago

@beaker Sorry, one last follow up question: why isn't the Btilde's permuntation [2,3,1] intead of [3,1,2]? I see that the test code would give a different result, but I can't see how that relates to the intuition. thanks again!

|

Collectives™ on Stack Overflow

Triple loop into tensor matlab

1 Answer 1

7 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

7 Comments

Your Answer

Sign up or log in

Post as a guest

Related