4
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ 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. $\endgroup$
    – eric_kernfeld
    Commented Oct 6, 2017 at 16:36
  • $\begingroup$ $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. :-) $\endgroup$
    – Sycorax ♦
    Commented Oct 6, 2017 at 16:42
  • 2
    $\begingroup$ 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. $\endgroup$
    – eric_kernfeld
    Commented Oct 6, 2017 at 17:16
  • $\begingroup$ @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$. $\endgroup$
    – Sycorax ♦
    Commented Oct 6, 2017 at 18:38
  • 1
    $\begingroup$ 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. $\endgroup$
    – eric_kernfeld
    Commented Oct 6, 2017 at 18:43

1 Answer 1

Reset to default
2
$\begingroup$

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.)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.