Re: st: Mata - generalised cross products of the form X'SX

[wgould@stata.com (William Gould, StataCorp LP)]

Gordon Hughes <G.A.Hughes@ed.ac.uk> writes, 

> [In Mata] I have identified that most of the time is consumed in loops
> in my likelihood evaluation that calculate matrices of the form Z=X'SX
> where S is symmetric but not a diagonal matrix.  [...]
> 
> This is a standard form in any GLS type calculation.  They are not 
> directly amenable to use with the Mata function --cross()-- but they 
> could be got into that form by forming by factorising S=QQ', then 
> forming P=Q'X, and finally using P=cross(P,P).

Gorden asks, 

> A.  Is there any Mata function that does this more directly?  Maybe 
>      Stata Corp might consider extending -cross()- to handle such 
>      cases.  The command --matrix glsaccum-- does this in Stata.

There is currently no command to perform directly what you ask.


> B.  Roughly, under what conditions might this sequence of steps 
>     reduce execution time and/or memory use on the assumption 
>     that the dimension of S is large relative to cols(X)?  

I'm sorry, I don't have an answer, except to say that the answer will
depend on whether you are using Stata/MP because MP is parallelizes
matrix multiplication, and thus performs it faster than does
Stata/SE.  In cases were rows>>cols or cols>>rows, multiplication
time approaches theoretical limits, which is to say, execution time 
becomes proportional to 1/cores.

> In a related context I have to calculate X'B'BX where B is square
> but not symmetric and found that use of --cross()-- does not save
> much (if any) time, though it is probably more efficient on memory.

Sometimes the Mata optimizer is not as smart as we wish that it were.
To calculate X'B'B'X as quickly as possible, code 

		T = B'X
		R = T'T

That will execute more quickly than R = X'B'B'X.  When you code
X'B'B'X, Mata does not notice that X'B'B'x = (B'X)'(B'X).  Calculation
of X'B'B'X results in Mata proforming three transposes and four matrix
multiplications.

My general recomendation is that you code long matrix expressions
carefully, performing the optmization yourself. The Mata optimizer is
good at combining separate calculations efficiently; it is poor at
splitting out compound calculations into seprate pieces.

In addition, use parenthesis where necessary to force the order of the
operations.  Concerning order of operations, Mata does not know what
you know, namely that rows >> cols.  Mata makes all of its
optimization decisions at compile time, and therefore can can choose
inefficient order of operations.  In this case, for memory
effinciency reasons, you do not want to calculate B'B to get to
X'B'B'X, you want to calculate X'B and B'X.  That right there will
make a differece in execution time, even without exploiting X'B
= (B'X)', which can be used additionally to save an entire multiplication.

-- Bill
wgould@stata.com
*
*   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/