Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: Efficient, foolproof calculation of matrix quadratic form with block-diagonal middle matrix? (Mata)


From   Venable <venablito@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   st: Efficient, foolproof calculation of matrix quadratic form with block-diagonal middle matrix? (Mata)
Date   Fri, 25 Sep 2009 10:24:30 -0400

Dear Statalisters,

I am using Mata to calcuate a quadratic form, the middle matrix of
which is block-diagonal, with all square and identical blocks. This is
most compactly expressed in Mata as:

A=X'*(I(N)#Om_Inv_Block)*X

X is (NxT) x K and is made up of N cross-sectional units, each of
which has T observations (so, a balanced panel of N with T
observations for each unit. You can think of X as X1 \ X2 \ ... \ XN,
with each Xn being K by T. Om_Inv_Block is TxT.

This is extremely simple to code, but Mata needs to create the NTxNT
matrix I(N)#Om_Inv_Block and with N and / or T large, this becomes a
very, very large matrix which occasionally causes Mata to crash.
(Interestingly, when Mata does not crash it doesn't take long at all
to calculate.)

I was wondering if there were a way for Mata to recognize the special
structure of the block-diagonal matrix and get around the need to
create this very large middle matrix explicitly.

Alternatively, I know that I could do the sum of the N terms
Xn'*Om_Inv_Block*Xn but I am worried that I would somehow mess this
up. So, another solution, I suppose, would be some idiot-proof way to
do this sum.

I apologize in advance if this is obvious, has been discussed before,
or is in the Mata manuals somewhere. I searched the Statalist using
terms "block diagonal" and "quadatic form" (individually and together)
and did not find a solution.

Many thanks.
*
*   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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index