I am looking to perform a 2-stage least-squares estimation with sparse matrices in R, in the style of Bramoulle et al (J. Econometrics 2009). Specifically, let:
G be a very sparse block-diagonal matrix, roughly 63,000 x 63,000
X be a design matrix, not sparse but full of dummy variables (with the relevant excluded variables to avoid the obvious collinearity problem), 63,000 x 16
y be a vector of outcomes, 63,000 x 1
The analysis has me calculating a number of the standard X(X'X)^-1 X' regression operations. For example, define S to be 3 horizontally appended sparse matrices (so S is 63,000 x 48):
S=[(I-G)X, (I-G)GX, (I-G)G^2 X]
P=S (S'S)^-1 S'
For simplicity, I want to perform the following calculation:
estimate = (X'PX)^-1 X' Py
I've been trying to use the 'Matrix' package to perform the relevant inversions directly (particularly the S'S inversion in calculating P), and have been running into memory issues:
P=S %*% solve(t(S) %*% S) %*% t(S) Error in LU.dgC(a) : cs_lu(A) failed: near-singular A (or out of memory)
I can obviously increase the memory allowed, but that only helps to a point.
My question is: are there best practices to exploit the structure of the problem? For example
1) I tried running
solve(t(S) %*% S, t(S)) to try to 'avoid' inverting, but got the same error message as above (out of memory)
2) S'S is symmetric, so doesn't this mean we can use the Cholesky decomposition? I tried and got an error:
chol2inv(chol(t(S)%*%S)) Warning message: In .local(x, ...) : Cholmod warning 'not positive definite' at file ../Cholesky/t_cholmod_rowfac.c, line 431 Error in chol2inv(chol(t(S) %*% S)) : error in evaluating the argument 'x' in selecting a method for function 'chol2inv': Error in .local(x, ...) : internal_chm_factor: Cholesky factorization failed
This is a bit strange since by construction S'S is symmetric.
I must admit I'm not the most well-versed in best-practices on inverting matrices in R, so any suggestions are welcome (also any other packages that might be better suited). Unfortunately I cannot provide a reproducible example as the data I'm using is protected.