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.
 A: I was overthinking this.
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.)
