Fastest way to implement a GLS estimator with a known covariance matrix Omega.

[Malcolm Wardlaw]

I'm trying to write a routine in Mata to calculate a GLS estimator where the covariance matrix Omega is estimated from some set of historical pre-period data.
So I have Omega and I'm trying to calculate (X'Omega^-1X)^-1(X'Omega^-1y).

I'm running into speed issues because the coefficients have to be calculated around 200 times using different X an y matricies with the same Omega matrix, and it has become computationally expensive no matter what I have done.

In Matlab the code

coeff = (X'*OmegaInv*X)\(X'*OmegaInv*Y)

where OmegaInv is the pre-calculated inverse of Omega, performs incredibly quickly, with the 200 or so independent calculations finishing in a few seconds. For the life of me, I can't seem to figure out what Matlab is doing under the hood to optimize this.

My initial attempt was to do a Cholesky decomposition on Omega followed by a triagular solve:

L = cholesky(Omega)
X = solvelower(L,X)
y = solvelower(L,y)
XX = quadcross(X,X)
Xy = quadcross(X,y)
beta = cholsolve(XX,Xy)'

This was quite fast in a single pass, but performing the 200 independent calculations required calling solvelower(L,X) and solvelover(L,y) 200 times each which seems to be the primary bottleneck. On an equivalent computer this clocks in at well over 20-25 seconds, which is not terrible but it's not quite fast enough to do the simulations I want.

The second thing I tried was to just brute force it with a similar code to the Matlab approach:

OmegaInv = cholinv(Omega)
beta    = cholsolve(X'*OmegaInv*X, X'*OmegaInv*y)

This approach doesn't appear to be any better.(As a side note, I'm aware that it's generally inadvisable to actually invert a large matrix, but I was trying to match what had been done in Matlab.)

I know that the Matlab interpreter is somehow optimizing the (X'Omega^-1X)^-1(X'Omega^-1y) calcualtion in order to avoid costly inversions and using solvers where possible to optimize the computation, but I can't really tell what it's doing that is making the Matlab code so fast and my Mata code so slow.

Is there something I'm missing? Is there an obvious way to code the calculation of (X'Omega^-1X)^-1(X'Omega^-1y), with or without pre-inverting Omega, that is more performant?

[Joerg Luedicke (StataCorp)]

Try this:

beta = cholsolve(X'*cross(X,OmegaInv)', (y'*cross(X,OmegaInv)')')

This should speed it up substantially. However, orders of magnitude faster still would be to do the linear transform desribed here followed by simple OLS. For example, something like:

L    = cholesky(Omega)
Li   = cholinv(L)
ys   = Li*y
Xs   = Li*X
beta = cross(invsym(Xs'*Xs),Xs')*ys