Efficient sum along block diagonals

Hello,

I have a matrix that has two levels of block-diagonal structure. The matrix is of size M*N*L-by-M*N*L, and I restructure the matrix using reshape() into a 6-D array of size

([M M N N L L])

In the original matrix, there are L^2 L-by-L blocks each containing an M*N-by-M*N block diagonal matrix.

The idea is to perform the outer block diagonal sum so that, in terms of the 6-D array, the dimensions would be

([M M N N 1 2*L-1])

Then, the inner block diagonal sum would give dimensions

([M M 1 2*N-1 1 2*L-1])

And finally, a sum along the diagonals to give dimensions

([1 2*M-1 1 2*N-1 1 2*L-1])

From here, I would use squeeze() to return a

([2*M-1 2*N-1 2*L-1])

cube.

I would like something equivalent to

sum(spdiags(A))

for each of these extra dimensions (or block diagonals, whichever way you prefer to view it). Additionally, because of the high-dimensionality, I would prefer to not use for loops, but would rather have something scalable and vectorized. Does anyone know of a way to do this, or have any insights for how to get started?

Thank you,

Alex

EDIT:

Here is some code that tries to implement what I am looking for in a vectorized but non-scalable fashion.

M = 2;
N = 2;
L = 2;
A = randn(M*N*L);

% Outer block diagonal sum
Out1 = A(5:8,1:4);
Out2 = A(1:4,1:4) + A(5:8,5:8);
Out3 = A(1:4,5:8);

% Inner block diagonal sum
In1_1 = Out1(3:4,1:2);
In1_2 = Out1(1:2,1:2) + Out1(3:4,3:4);
In1_3 = Out1(1:2,3:4);

In2_1 = Out2(3:4,1:2);
In2_2 = Out2(1:2,1:2) + Out1(3:4,3:4);
In2_3 = Out2(1:2,3:4);

In3_1 = Out3(3:4,1:2);
In3_2 = Out3(1:2,1:2) + Out1(3:4,3:4);
In3_3 = Out3(1:2,3:4);

% Diagonal sum
B1 = [sum(spdiags(In1_1)); sum(spdiags(In1_2)); sum(spdiags(In1_3))];
B2 = [sum(spdiags(In2_1)); sum(spdiags(In2_2)); sum(spdiags(In2_3))];
B3 = [sum(spdiags(In3_1)); sum(spdiags(In3_2)); sum(spdiags(In3_3))];

% Final result
C = cat(3,cat(3,B1,B2),B3);

EDIT 2:

Here is a latex equation that describes what I am looking to achieve:

I know this can be achieved by using the "find()" function to make this scalable, but that is super slow.

Any help is greatly appreciated!

[chicken vector]

Could you create a MRE because, at least for me, the problem is a bit hard to understand.

For instance:

"In the original matrix, there are L^2 L-by-L blocks each containing an M*N-by-M*N block diagonal matrix."

1) I don't understand how is it possible to have L^2 (instead of L) M-by-N submatrices?

Your original structure is something like this?

M = 2; N = 2; L = 3;
mnBlocks = repmat({ones(M,M)},N,1);
mnDiag = blkdiag(mnBlocks{:})

mnDiag = 4×4

1 1 0 0 1 1 0 0 0 0 1 1 0 0 1 1

mnlBlocks = repmat({mnDiag},L,1);
mnlDiag = blkdiag(mnlBlocks{:})

mnlDiag = 12×12

1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0

"The idea is to perform the outer block diagonal sum"

"Then, the inner block diagonal sum would give dimensions"

"And finally, a sum along the diagonals to give dimensions"

2) What do you mean by outer, inner and along block diagonal sum?

[Alex Batts]

My apologies, here's some code to help reproduce what I'm trying to say.

M = 2;
N = 2;
L = 2;
A = randn(M*N*L);

% Outer block diagonal sum
Out1 = A(5:8,1:4);
Out2 = A(1:4,1:4) + A(5:8,5:8);
Out3 = A(1:4,5:8);

% Inner block diagonal sum
In1_1 = Out1(3:4,1:2);
In1_2 = Out1(1:2,1:2) + Out1(3:4,3:4);
In1_3 = Out1(1:2,3:4);

In2_1 = Out2(3:4,1:2);
In2_2 = Out2(1:2,1:2) + Out1(3:4,3:4);
In2_3 = Out2(1:2,3:4);

In3_1 = Out3(3:4,1:2);
In3_2 = Out3(1:2,1:2) + Out1(3:4,3:4);
In3_3 = Out3(1:2,3:4);

% Diagonal sum
B1 = [sum(spdiags(In1_1)); sum(spdiags(In1_2)); sum(spdiags(In1_3))];
B2 = [sum(spdiags(In2_1)); sum(spdiags(In2_2)); sum(spdiags(In2_3))];
B3 = [sum(spdiags(In3_1)); sum(spdiags(In3_2)); sum(spdiags(In3_3))];

C = cat(3,cat(3,B1,B2),B3);

Here, C is a 3x3x3 3-D array.

Does this help?

Answer 1

[Matt J]

I believe this works! If you want to put this into a separate answer I will vote for this!

I'm glad, but if it does work, please Accept-click it as well.

temp=oneStep(X,M*N,L);
temp=oneStep(temp,M,N);
temp=oneStep(temp,1,M);

result=reshape(temp,2*M-1,2*N-1,2*L-1);

function Y=oneStep(X,p,q)
%p - blocklength
%q - length of matrix, counted in blocks

  Y=blkReshape(X,[p,p],1,q^2,[]);% from https://www.mathworks.com/matlabcentral/fileexchange/115340-block-transposition-permutation-and-reshaping-of-arrays 

  Y=pagemtimes( reshape(Y,p^2,q^2,[]) , summationMatrix(q) );
  
  Y=reshape(Y,p,p,[]);
  
end

function S=summationMatrix(q)
%q - length of matrix, counted in blocks

  T=toeplitz(0:-1:1-q,0:q-1);

   S=(T(:)==(1-q:q-1));

end

[Matt J]

If you're going to be repeating this for different array, but with the same M,N,L parameters, it will save computation if you pre-compute the summationMatrix results one time only.

[Alex Batts]

Thank you for your help!

Answer 2

[Matt J]

Using this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/115340-block-transposition-permutation-and-reshaping-of-arrays

You can do the outer sum as follows, where X is your matrix in its original form,

Xr=blkReshape(X,[M*N,M*N],1,1,[]);
outersum=sum(Xr(:,:,1:L+1:end),3);

and then proceed likewise for the inner sums.

[Alex Batts]

I made an edit to the original question and added code to clarify what I am looking for. In this case, "outersum" gives the same answer as "Out2".

I am aware that this method will give me what I am looking for if I loop over the third index of "Xr", but that is what I am trying to avoid.

[Matt J]

You can vectorize the different sums over Xr(:,:,i) using,

outersums = permute( tensorprod(diagSum(L),Xr,2,3), [2,3,1])

function S=diagSums(L)

  T=toeplitz(0:-1:1-L,0:L-1);

   S=(T(:)'==(1-L:L-1)');

end

[Alex Batts]

What is "T" in that function?

[Matt J]

 T=toeplitz(0:-1:1-L,0:L-1);

[Alex Batts]

The code worked, but I am now trying to adapt it to the inner sums. For some reason, I am getting a 3x2x2x3 array, with incorrect values, when running the code:

temp = blkReshape(outersums,[M M],1,1,N^2,[]); % Reshapes outersums to a 2x2x4x3 or MxMx(N^2)x(2*L-1) array
innersums = tensorprod(diagSums(N),temp,2,3);

Would I have to reshape each MNxMN matrix within outersums individually prior to doing the tensorprod? Or is there something I am missing?

[Matt J]

I think this might be the way to go,

temp=oneStep(X,M*N,L);
temp=oneStep(temp,M,N);
temp=oneStep(temp,1,M);

result=reshape(temp,2*M-1,2*N-1,2*L-1);

function Y=oneStep(X,p,q)
%p - blocklength
%q - length of matrix, counted in blocks

  Y=blkReshape(X,[p,p],1,q^2,[]);

  Y=pagemtimes( reshape(Y,p^2,q^2,[]) , summationMatrix(q) );
  
  Y=reshape(Y,p,p,[]);
  
end

function S=summationMatrix(q)
%q - length of matrix, counted in blocks

  T=toeplitz(0:-1:1-q,0:q-1);

   S=(T(:)==(1-q:q-1));

end

[Alex Batts]

Ah, interesting with the pagemtimes. I believe this works! If you want to put this into a separate answer I will vote for this!

[Alex Batts]

Actually, I think there should be a permute(result,[2 1 3]). But other than that it looks good!