Efficient/feasible sparse matrix inversion in R

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.

• having protected data doesn't/shouldn't mean you can't provide a reproducible example ... make up some data that looks like yours. If it doesn't produce the same error, you have a clue that your data has special structure. If it does, you've got a reproducible example. – Ben Bolker Oct 3 '16 at 20:27
• Stop! Whenever you see a matrix inverse $A^{-1}$ times something, that DOES NOT mean you should form the inverse. Instead, solve the linear system. – Matthew Gunn Oct 3 '16 at 20:27
• @MatthewGunn Right, I'm trying to avoid direct inversion, but doing a two-argument solve doesn't work either. Unless you had something else in mind? – Ray Oct 3 '16 at 20:30
• Whenever you see $A^{-1}b$ you instead want to solve the linear system $Ax = b$. In some sense, you really have a programming problem here. One overall question is whether you want to use "sparse" matrices. Using a sparse matrix, you can form G. You have to be incredibly careful though to only call functions that preserve sparsity and utilize the sparsity. Alternatively, you can break down the algorithm yourself, represent G as a list of matrices or something. I don't know what's better, but this is a programming problem... – Matthew Gunn Oct 3 '16 at 20:45
• I am sympathetic. I solved a problem like this recently and it was almost this large, too. After a day of working at it, exploring various algorithms, I finally found a solution that fit into available RAM and was fast. In the end it required all of eight lines of R code (which mostly manipulated arrays with four indexes each: the trick to vectorization often is getting the right numbers into exactly the right places). – whuber Oct 3 '16 at 21:35