Efficient sum along block diagonals
4 次查看(过去 30 天)
显示 更早的评论
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1375239/image.png)
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
2023-5-4
编辑:chicken vector
2023-5-4
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{:})
mnlBlocks = repmat({mnDiag},L,1);
mnlDiag = blkdiag(mnlBlocks{:})
"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?
采纳的回答
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 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!