topic Reshaping block matrix in block diagonal form in SAS/IML Software and Matrix Computations

Reshaping block matrix in block diagonal form

Daniel_Paul — Wed, 27 Jun 2018 14:24:43 GMT

Hello,

I have a block matrix of the form A = {A1, A2, ..., AN} with symmetric matrices A1, .., AN each of the same order 9*9 (N could be as large as 400). I want to shape the matrix in the block-diagonal form B = {C 0 0, 0 D 0, ..., 0 0 E} with matrices of the form C = {2A1 A1 A1, A1 2A1 A1, A1 A1 2A1} and D and E of the same form with A2 vs. AN instead of A1. If I allow loops, I can figure out a solution, but I try to avoid loops. Any suggestions?

Bye, Daniel

Re: Reshaping block matrix in block diagonal form

Rick_SAS — Wed, 27 Jun 2018 14:52:42 GMT

Use Kronecker products to construct the C, D, and E matrices from A1, A2, and A3.

Use the BLOCK function to get the block diagonal matrix of C, D, and E.

proc iml;
N = 3;
N2 = N*(N+1)/2;
A1 = sqrsym(1:N2);
A2 = sqrsym(N2:1);

M = j(N, N, 1) + I(N);
C = M @ A1;
D = M @ A2;
print C, D;

B = block(C, D);
print B;

Re: Reshaping block matrix in block diagonal form

Daniel_Paul — Wed, 27 Jun 2018 15:01:25 GMT

I see the principal idee, but N can be large, so I have to loop in the block part. I try to avoid that loop.

proc iml;
p = 3;
N = 2;

start TBlock(AN, p); 
M = j(p,p,1) + I(p);
C = M @ AN;
return (C);
finish TBlock;

start AD(N, A, p);
t = 1;
do rs=1 to N;
	B = A[t:(t+2),];
	BB = TBlock(B, p);
	AF = block(AF, BB);
	t = t + p;
end;
return (AF);
finish AD;

A1 = {1 2 3, 
	  2 4 5,
      3 5 7};
A2 = {1 7 9,
	  7 8 1,
	  9 1 1};

AC = A1//A2;

AF = AD(N, AC, p);

print AF;
quit;

Re: Reshaping block matrix in block diagonal form

IanWakeling — Thu, 28 Jun 2018 10:31:13 GMT

It looks too complicated to be able to avoid the loop. What might be slowing things down, is the fact that you are 'growing' the matrix AF by steps within the loop. After each iteration a new version of AF is declared, and ever more data needs to be moved around in memory. It would be more efficient to declare AF as a zero matrix, at it's final size, before the loop, and then write the block diagonal sub-matrices using row and column indexing.