Efficient Summation of two matrices multiplied

Hello,

I have a summation that I have been trying to vectorize using block processing, but have been unable to find an underlying structure to compute this. The equation is:

I understand that this can be implemented using for loops and conditional statements, however, this is very slow and doesn't work efficiently as the matrices get larger.

It's also fine to loop through α and β as the B matrix will typically not get so big to where it would bog down a double for loop. However, if it can be fully vectorized including α and β that would be even better.

Thank you in advance!

Alex

EDIT:

As was pointed out, the equation can likely instead be solved for:

which allows for the appropriate indexing of the "4D" matrix.

Here is some code that works:

A = randn(8,8)+1i*randn(8,8);
data = randn(16,16) + 1i*randn(16,16);
presums = size(A);
A = repmat(A,8,8);
out = foo(data,A,idx);

function out = foo(data,A,idxd)

k = 0:idxs(1)-1;
l = k;
p = 0:idxs(2)-2;
q = p;
t = (idxs(1)):size(data,1)-1;
s = (idxs(2)):size(data,2)-1;
a = 1:size(data,1);
b = 1:size(data,2);
out = zeros(size(data));

for aa = a
    for bb = b
        for tt = t
            for kk = k
                if tt - kk == aa 
                    for ss = s
                        for pp = p
                            if ss - pp == bb
                                for qq = q
                                    for ll = l
                                        out(aa,bb) = out(aa,bb) + A(kk+idxs(1)*(pp)+1,ll+idxs(1)*(qq)+1)*...
                                            data(tt-ll+1,ss-qq+1);
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

end

[Walter Roberson]

are you sure k=0..m and not k=0..m-1?

If it were m-1 then the equations would correspond to simulating a 4d matrix in a 2d matrix but with the way it is now, it is overlapping on the boundaries

[Alex Batts]

It’s very possible that is the case. The equation as is is the way I derived it. However, A is supposed to be a block diagonal matrix which could be represented as a 4d matrix. Would that then also make t=[m … M-1]?

Answer 1

[Sugandhi]

Hi,

I understand that you are trying to perform a sum over a 4D array, where each element of the array is a 2x2 matrix. You want to sum up all the elements of the 4D array, but only when certain conditions are met.

The equation you provided is correct, but there might be a mistake in the bounds of the loop. If the dimensions of the 4D array are M x N x P x Q, then the loops should go from 0 to M-1, 0 to N-1, 0 to P-1, and 0 to Q-1.

To avoid looping over the entire 4D array, you can use advanced indexing to extract the desired subarray and then sum it up.