# Solving for auto-regression coefficients when covariance is singular?

I need coefficients of autoregressive model of $$x$$: matrix $$B$$ with zero diagonal which minimizes squared error:

$$\text{argmin}_B E[||x - Bx||^2]$$

For non-singular covariance, I can get solution directly from $$E[xx']^{-1}$$ using following procedure:

1. Let $$D2$$ be the inverse second moment with off-diagonal terms set to zero $$D2=E[xx']^{-1}_d$$

2. Read coefficients $$\{a_1\ldots a_n\}$$ for the model $$x_i=a_0 x_0 + \ldots a_n x_n$$ from $$i$$th row of following matrix

$$B=I-(X'X)^{-1} D2^{-1}$$

What should I do for singular $$E[xx']$$?

For instance, suppose

$$X=\left( \begin{array}{cccc} 2 & -1 & -2 & 1 \\ -2 & 1 & -2 & 1 \\ 2 & 1 & 1 & 3 \\ \end{array} \right)$$

I can turn this into 4 regression problems, solve each problem using $$\beta=(X'X)^{-1}(Y'X)$$, and get following coefficients of auto-regressive model with $$x=Bx$$ $$B=\left( \begin{array}{cccc} 0 & -\frac{1}{2} & \frac{7}{4} & \frac{7}{8} \\ -2 & 0 & \frac{7}{2} & \frac{7}{4} \\ \frac{4}{7} & \frac{2}{7} & 0 & -\frac{1}{2} \\ \frac{8}{7} & \frac{4}{7} & -2 & 0 \\ \end{array} \right)$$

1. This requires 4 matrix inversions, whereas for non-singular covariance matrix I had to only do 1

2. This is a well-posed problem, but in general it may be ill-posed, so one would use some linear solver instead of $$(X'X)^{-1}(Y'X)$$ formula. However, calling linear solver for each column is too expensive.

Edit

Aksakal suggested to use regularized estimator for inverting the covariance, trying this out here

Using shrunk covariance estimator with very small lambda gives good results. Numerically inverting singular covariance without regularization and sticking it into the formula gives slightly better error.

Error direct:  1.4051584874249273e-30
Error inverse:  3.2047474274603605e-30
Error Ledoit-Wolf:  9.351834226497054
Error shrunk covariance:  0.7203818230423954
Error very-shrunk covariance:  0.00952058588215488

• Your comparisons are problematic. The better way to compare would be to start with a known lag structure, then generate the set of rank deficient samples from it. Then study properties of different estimates of lag coefficients. Currently you are looking at a single sample so any conclusions you make are going to be highly uncertain Jun 21 at 22:14
• I think it's an issue of Ledoit-Wolf optimizing for different objective. To get coefficients, I divide entries of precision matrix by entries of the diagonal of the precision matrix. Small bias in the precision matrix estimation produces large error in the coefficients. Jun 22 at 22:40

your problem is rank rank-deficient, not necessarily singular when you have more variables than observations.

You could regularize the covariance matrix before operating it, e.g. using Ledoit-Wolf method:

Sig = np.cov(Y.T)
iSig = np.linalg.inv(Sig)

from sklearn.covariance import LedoitWolf

RSig = LedoitWolf().fit(Y)
iRSig = np.linalg.inv(RSig.covariance_)
print('ordinary:',iSig, '\nregilarized: ', iRSig)


Output

ordinary: [[ 1.40737488e+14  2.81474977e+14  7.50599938e+14 -1.40737488e+15]
[ 2.81474977e+14  5.62949953e+14  1.50119988e+15 -2.81474977e+15]
[-1.68884986e+15 -3.37769972e+15  3.00239975e+15 -1.12589991e+15]
[ 2.25179981e+15  4.50359963e+15 -6.00479950e+15  4.50359963e+15]]
regilarized:  [[ 0.45698883  0.17402552 -0.14171221 -0.0944748 ]
[ 0.17402552  0.80670551 -0.15040682 -0.10027121]
[-0.14171221 -0.15040682  0.6620951  -0.22126292]
[-0.0944748  -0.10027121 -0.22126292  0.84648087]]

• Thanks for the idea! I actually get much better results by replacing LedoitWolf with ShrunkCovariance, updated question with link to code. Both of them are worse than just numerically inverting the singular covariance matrix and sticking the result into formula -- formula divides by diagonal, so the resulting coefs become reasonable. This is a general problem with regularization, optimal lambda is problem specific, hard to check for correctness Jun 20 at 20:09
• Ledoit wolf shrinks to a certain kind of a matrix. You can use the idea to shrink to something else, if you have a better target matrix. For instance, you could use a diagonal matrix or scaled identity matrix Jun 20 at 21:01