How to vectorize the following operations?

1 次查看（过去 30 天）
Wei HU
Wei HU 2019-3-10
编辑： Sean de Wolski 2022-6-6
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?

请先登录，再回答此问题。

回答（2 个）

Voss
Voss 2022-6-4
Ran in:
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 2022-6-4
编辑：Sean de Wolski 2022-6-6
Ran in:
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]

类别

在 Help Center 和 File Exchange 中查找有关 Matrices and Arrays 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by