Sparse Cholesky decomposition of factorized matrix

Question

I want the diagonal of a matrix $Y^TA^{-1}Y$ where $A=X^TX$ and $X$ is very sparse with dimensions ~1e6 x ~1e5 (so $A$ is 1e5 by 1e5). $Y$ is something like 1e5 by 1e4 (also sparse). Currently I'm explicitly calculating $A$ then getting its Cholesky decomposition $LL^T$ via sksparse.cholmod which exploits its sparsity in $A$ nicely. Then I can solve(L,Y) and get my diagonal terms.

However, I feel like I should be able to do something with $X$ directly, e.g. if I had a QR decomposition $X=QR$ then I would have $A = R^T Q^T Q R = R^T R$. This would presumably save me some computation by avoiding the explicit calculation of $A$, but there don't seem to sparse QR methods widely available (https://github.com/yig/PySPQR fails to compile for me): I'm limited to Python since this is part of a package.

Is there some more straightforward approach I'm missing (e.g. via least squares)?

(also please let me know if another stackexchange like https://scicomp.stackexchange.com/ would be more appropriate)

Answer 1

You'll find the LSQR as an option for solving large sparse linear systems that max be over or under determined. This is included in scipy.sparse.linalg.lsqr. However, this solver doesn't seem applicable at for your case (unless your $Y$ has a deeper story, like $Y=X^T b$?). Iterative solvers wouldn't be my pick for solving 1e4 RHS' regardless.

Note that your $A$ has the major advantage of being 10 times smaller than $X$, and the factorization of sparse systems aren't necessarily themselves (very) sparse. Computing $A$ is basically free compared to everything else we need to solve here. There is possibly some deeper cleverness beyond me that might be found here but I suspect the boring answer that you aren't going to do a lot better than what you already do.