How to vectorize the following operations?

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?

Answer 1

[Voss]

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]). ]

Answer 2

[Sean de Wolski]

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]