Efficient sum along block diagonals

2 次查看（过去 30 天）
Alex Batts
Alex Batts 2023-5-4
移动：Matt J 2023-5-5
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!
  2 个评论
chicken vector
chicken vector 2023-5-4
编辑：chicken vector 2023-5-4
Ran in:
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
Alex Batts 2023-5-4
编辑：Alex Batts 2023-5-4
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?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-5-4
编辑：Matt J 2023-5-4
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
  2 个评论
Matt J
Matt J 2023-5-4
移动：Matt J 2023-5-5
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.

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2023-5-4
编辑：Matt J 2023-5-4
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.
  8 个评论
显示 6更早的评论
Alex Batts
Alex Batts 2023-5-4
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
Alex Batts 2023-5-4
Actually, I think there should be a permute(result,[2 1 3]). But other than that it looks good!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by