Re: st: Mata question: putting matrices along the diagonal

[Stas Kolenikov <skolenik@gmail.com>]

On Thu, May 27, 2010 at 1:02 PM,  <mjbaker@hunter.cuny.edu> wrote:
> Dear Listers:
>
> I am interested in making one large "diagonal" matrix out of a bunch of matrices in Mata. For example, suppose I have 3 2x2 matrices labeled m1, m2, and m3. I want to construct the matrix M, where M is something like:
> M=
> (m1[1,1],m1[1,2],     0,        0,           0, 0               \
>  m1[2,1],m1[2,2],     0,        0,           0, 0               \
>      0    ,   0,       m2[1,1],m2[1,2],     0, 0               \
>      0    ,   0,       m2[2,1],m2[2,2],     0,  0              \
>      0    ,   0,             0,        0,      m3[1,1],m3[1,2] \
>      0    ,   0,             0,        0,      m3[2,1],m3[2,2] )
>
> That is, I seek to place each of the individual matrices along the diagonal, zeros elsewhere. In my problem, both the number of matrices and the dimensions of each matrix may be rather large.

For a fixed number of matrices:

M = m1
M = (M, J(rows(M),cols(m2),0) \ J(rows(m2),cols(M),0), m2 )
M = (M, J(rows(M),cols(m3),0) \ J(rows(m3),cols(M),0), m3 )

There might be a more elegant solution, but this is a quick and dirty fix.

For a varying number of matrices, I am not sure what a good solution
might be. In Stata, a natural solution would be to use locals, but
Mata does not have them. So you can switch back and forth between
Stata and Mata; or you can set up some pointers to matrices and cycle
over them.

If you have an upper bound on the number of matrices, say 5, then you
can write a couple of functions to work this through:

real matrix BlockDiag( real matrix M1, real matrix M2, | real matrix
M3, real matrix M4, real matrix M5 ) {
   real matrix M;

   M = DiagAppend(M1, M2)
   if ( ! missing(M3) ) M = DiagAppend( M, M3 )
   if ( ! missing(M4) ) M = DiagAppend( M, M4 )
   if ( ! missing(M5) ) M = DiagAppend( M, M5 )
   return(M)
}

real matrix DiagAppend(real matrix M1, real matrix M2) {
   return( (M1, J(rows(M1),cols(M2),0) \ J(rows(M2),cols(M1),0), M2) )
}

-- 
Stas Kolenikov, also found at http://stas.kolenikov.name
Small print: I use this email account for mailing lists only.

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/