Solving for auto-regression coefficients when covariance is singular?

Question

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

Answer 1

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]]