Using matrix structure to speed up matlab

Question

Suppose that I have an N-by-K matrix A, N-by-P matrix B. I want to do the following calculations to get my final N-by-P matrix X.

X(n,p) = B(n,p) - dot(gamma(p,:),A(n,:))

where

gamma(p,k) = dot(A(:,k),B(:,p))/sum( A(:,k).^2 )

In MATLAB, I have my code like

for p = 1:P
    for n = 1:N
        for k = 1:K
            gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
        end
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

which are highly inefficient since it uses three for loops! Is there a good way to speed up this code?

Answer 1

Use bsxfun for the division and matrix multiplication for the loops:

gamma = bsxfun(@rdivide, B.'*A, sum(A.^2));
x = B - A*gamma.';

And here is a test script

N = 3;
K = 4;
P = 5;

A = rand(N, K);
B = rand(N, P);

for p = 1:P
    for n = 1:N
        for k = 1:K
            gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
        end
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

gamma2 = bsxfun(@rdivide, B.'*A, sum(A.^2));
X2 = B - A*gamma2.';

isequal(x, X2)
isequal(gamma, gamma2)

which returns

ans =
     1
ans =
     1

Answer 2

It looks to me like you can hoist the gamma calculations out of the loop; at least, I don't see any dependencies on N in the gamma calculations.

So something like this:

for p = 1:P
    for k = 1:K
        gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
    end
end
for p = 1:P
    for n = 1:N
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

I'm not familiar enough with your code (or matlab) to really know if you can merge the two loops, but if you can:

for p = 1:P
    for k = 1:K
        gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
    end
    for n = 1:N
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end