Efficient matrix multiplication with weights

Let A and B be two matrices, say square NxN matrices. Ordinary matrix multiplication A*B implements (A*B)_{ij} = Sum_k A_{ik} B_{kj}. Is there an efficient way in Matlab to implement a weighted version of this product, where we have a matrix of weights W and we want to do :

Weighted(A*B)_{ij} = Sum_k A_{ik} B_{kj} W_{i-j,k}

(let's say here that A and B are triangular so that only i>=j need be considered).

How can I efficiently express Weighted(A*B), avoiding, if possible, for loops and the like ? I would like to keep everything vectorialized / use only matrix products and elements wise products etc.

[Matt J]

I would like to keep everything vectorialized / use only matrix products and elements wise products etc.

Even if for-loops are faster?

[Matt J]

Also, are A and B Toeplitz as well as triangular?

[Pierre-Louis Giscard]

Actually the fastest option is the best so you are right if the for loops are faster I would use them. In general A and B are not Toeplitz. In applications A and B are rather large (say 1000x1000) so memory usage could also be an issue.

Answer 1

[Matt J]

A more memory efficient solution is as follows. It has a loop, but is still highly vectorized.

Wt=W.';
At=A.';
T=toeplitz(1:N,[1,zeros(1,N-1)]);

result=zeros(N);
for i=1:N

    result(T==i)=sum(  At(:,1:end+1-i).*Wt(:,i).*B(:,i:end) ,1);
    
end

[Pierre-Louis Giscard]

Thank you for your codes ! I think this second proposition will be more suited to my applications as I am worried about the memory usage. I will try to see how fast this is but it seems it will be much faster than with all the nested for loops of the naive approach.

[Matt J]

You're welcome. If it works as you need it to, though, please Accept-click the answer.

Answer 2

[Matt J]

Using sepblockfun() from,

https://www.mathworks.com/matlabcentral/fileexchange/48089-separable-block-wise-operations

T=toeplitz(1:N);

WW=W.';
WW=reshape(WW(:,T), N^2,N);

BB=repmat(B,N^2,1);
AA=repmat( reshape(A.',[],1)  ,1,N^2);

result=sepblockfun(AA.*WW.*BB, [N,1] , 'sum' ); %

[Matt J]

For N=1000, you would need a lot of RAM for this to work. You might be able to mitigate RAM requiements by using single floats inputs. The result could still be obtained in doubles with,

result=sepblockfun(AA.*WW.*BB, [N,1] , @(x,d)sum(x,d,'double') ); %