Compute covariance matrix via rank-1 update to $M^\top M$

I have a large, sparse matrix $M\in\mathbb{R}^{n\times p}$. Centering $M$ to compute the covariance matrix $\Sigma$ would, in general, destroy the "zeros aren't stored" property of sparse matrices. Instead, we can center $\Sigma$ after computing the second moment. There are any number of ways to do this, but I want one that is cheap and numerically stable.

The goal of estimating $\Sigma$ is to compute many evaluations of $f(x) = x^\top \Sigma^{-1} x$.

My first thought was to compute $M^\top M$ and apply the Bunch-Nielsen-Sorensen formula for a rank-1 update to a spectral decomposition, which allows one to compute the eigenvalues and eigenvectors to $A$ for $D$ diagonal and $\rho\in\mathbb{R}$ $$A = D + \rho zz^\top$$ (Some algebra extends this to updating $\Sigma$.) However, BNS requires that $z_i>0 ~\forall ~ i$, and that does not hold for the mean vector $\mu$ in general.

As a work-around, I could apply two successive BNS updates, one for a $\tilde{\mu}$ which shifts the nonpositive elements of $\mu$ to be positive, and a second which offsets that shift. This is (probably) more efficient than two spectral decompositions.

Another way would form a pseudoinverse and apply the Sherman-Morrison formula, but we all know that one should never explicitly compute a matrix inverse, hence my interest in factorizations.

• How many of those quadratics are you planning to compute? If it's not many (relative to the size of your $S$), then solving $\Sigma^{-1}x$ each time may be cheaper than computing a pseudoinverse and storing it for repeated use. For the few-quadratics case, there are solvers that rely on repeated multiplications of the form $\frac{1}{n-1}(M^TM-n\mu\mu^T)v$. These multiplications can be done cheaply without forming $\frac{1}{n-1}(M^TM-n\mu\mu^T)$ explicitly. Where to find a software package that does this, though, I'm not sure. – eric_kernfeld Oct 6 '17 at 16:36
• $S$ is approximately $10^5 \times 10^5$ and I'll be evaluating $f$ on the order of $10^7$ times. I'm not afraid of coding stuff up myself, as long as I have a good strategy to solving the problem. :-) – Sycorax Oct 6 '17 at 16:42
• I'm still sorting through this mentally, and I realized that "using versus not using decompositions" is a slightly different decision than "explicitly or not explicitly forming an inverse". I will focus on fast, stable solutions by any means. Or any covariances. – eric_kernfeld Oct 6 '17 at 17:16
• @eric_kernfeld It might be that the best approach is the simplest -- compute the Cholesky factorization $\Sigma = LL^\top$ and forwardsolve the linear system to evaluate $f$. This can be done directly, or indirectly if there exists a cheaper way to apply a rank-1 update to $L$. – Sycorax Oct 6 '17 at 18:38
• That's what I'm thinking of. There may be a bonus in there -- if $A^TA=B$, the $R$ from the QR decomposition of $A$ is the same as the $R$ from the Cholesky factorization of $B$. If you use Givens rotators to compute the QR of $A$, you can take advantage of the sparsity by only rotating the nonzero elements. I know that my A-B scenario is not identical to your M-mu-sigma problem, but maybe some tweak would allow this trick to apply. – eric_kernfeld Oct 6 '17 at 18:43

The canonical way to compute $f(x)$ is to compute the Cholesky factor of $\Sigma=\tilde{R}^\top \tilde{R}$ where $\tilde{R}$ is and upper-triangular. Then we can write $$f(x)=x \tilde{R}^{\top}\tilde{R}^{-\top} x = y^\top y$$ and solve the triangular system for $y$. So we should focus on how to get to a Cholesky factor of $\Sigma$.
Using the equation $$\Sigma=\frac{1}{n-1}M^\top M - \frac{n}{n-1}\mu\mu^\top$$ I observe that this is just a rank-one down-date to $M^\top M$. Or, in terms of Cholesky factors, we can use the factorization $M=QR$ for orthonormal $Q$ and upper triangular $R$ to write $$M^\top M = R^\top Q^\top QR = R^\top R$$
We can update $R$ to be Cholesky factors of $\Sigma$ using hyperbolic rotations, as outlined in Golub & Van Loan, Matrix Computations 4th edition, p. 339-341. The procedure covers three pages of text, so it would be impractical to reproduce here. The point is that it is a down-dating framework similar to Givens rotations. (To efficiently compute $M=QR$ for sparse $M$ is what Golub & Van Loan call the "The Sparse QR Challenge" and it is discussed in pp. 606-608.)