Efficient sum along block diagonals
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
chicken vector
le 4 Mai 2023
Modifié(e) : chicken vector
le 4 Mai 2023
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?
Réponse acceptée
Matt J
le 4 Mai 2023
Modifié(e) : Matt J
le 4 Mai 2023
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 commentaires
Plus de réponses (1)
Voir également
Catégories
En savoir plus sur Operating on Diagonal Matrices dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!