How to vectorize the following operations?

1 view (last 30 days)
Wei HU
Wei HU on 10 Mar 2019
Edited: Sean de Wolski on 6 Jun 2022
For example, suppose we are going to calculate serveral quadratic forms:
A is a 3-d array with shape (d,d,n). x y are 2d vectors with shape (d,n). I need to calculate,
b(i) = x(:,i)'*A(:,:,i)*y(:,i)
But can we do this without such a loop? can we vectorize it?

Answers (2)

Voss
Voss on 4 Jun 2022
Yes, this can be vectorized.
d = 4;
n = 5;
A = rand(d,d,n);
x = rand(d,n);
y = rand(d,n);
% for-loop b calculation
b_for = zeros(1,n);
for i = 1:n
b_for(i) = x(:,i)'*A(:,:,i)*y(:,i);
end
disp(b_for)
2.5627 2.7269 3.9463 1.5213 1.8338
% vectorized b calculation
b_vec = permute(sum( ...
repmat(permute(x,[1 3 2]),[1 d 1]).*A.*repmat(permute(y,[3 1 2]),[d 1 1]), ...
[1 2]),[1 3 2]);
disp(b_vec)
2.5627 2.7269 3.9463 1.5213 1.8338
% check the difference
max(abs(b_for-b_vec))
ans = 4.4409e-16
[ Here I assumed x is real (so that x(:,i)' could've been written x(:,i).'). If x is complex, use permute(conj(x),[1 3 2]) instead of permute(x,[1 3 2]). ]

Sean de Wolski
Sean de Wolski on 4 Jun 2022
Edited: Sean de Wolski on 6 Jun 2022
This ought to do it for you in one shot. pagemtimes
d = 4;
n = 5;
A = rand(d,d,n);
x = rand(d,n);
y = rand(d,n);
% for-loop b calculation
b_for = zeros(1,n);
for i = 1:n
b_for(i) = x(:,i)'*A(:,:,i)*y(:,i);
end
disp(b_for)
0.2415 1.1436 1.9688 2.1402 2.2604
For comparison
q = squeeze(pagemtimes(pagemtimes(reshape(x,1,d,n),A),reshape(y,d,1,n)))
q = 5×1
0.2415 1.1436 1.9688 2.1402 2.2604
Note, this may not be faster than the for-loop, but it's "vectorized".
[Edit Monday morning for better implementation]

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!