# Iterative “solver” for $x^t \Sigma^{-1} x$

I can't imagine I'm the first to think about the following problem, so I'll be satisfied with a reference (but a complete, detailed answer is always appreciated):

Say you have a symmetric positive definite $\Sigma \in \mathbb{R}^{n \times n}$. $n$ is thought of as very large, so holding $\Sigma$ in memory is impossible. You can, however, evaluate $\Sigma x$, for any $x \in \mathbb{R}^{n}$. Given some $x \in \mathbb{R}^{n}$, you'd like to find $x^t\Sigma^{-1}x$.

The first solution that comes to mind is to find $\Sigma^{-1}x$ using (say) conjugate gradients. However, this seems somewhat wasteful - you seek a scalar and in the process you find a gigantic vector in $\mathbb{R}^{n}$. It seems to make more sense to come up with a method to calculate the scalar directly (i.e. without passing through $\Sigma^{-1}x$). I am looking for this kind of method.

• Does your matrix arise from $\Sigma = A^TA$ for some "short & wide" rectangular $A$? – GeoMatt22 Oct 13 '15 at 0:50
• @GeoMatt22 unfortunately not. But let's say it does - what would you suggest in that case? – Yair Daon Oct 13 '15 at 12:03
• Yair, I was just thinking if there is some smaller matrix to work with ... not sure it would help anyway. Have you tried googling "matrix free mahalanobis distance" or similar? Sorry to not be of more help! – GeoMatt22 Oct 13 '15 at 13:25
• Thanks @GeoMatt22, I wasn't able to find anything online. – Yair Daon Oct 13 '15 at 17:04

I don't think I've heard of any method that does what you want without actually solving $y=\Sigma^{-1}x$.
The only alternative I can offer is if you knew something about the eigenvectors and -values of $\Sigma$. Say you knew that they are $\lambda_i,v_i$, then you can represent $\Sigma=V^T L V$ where the columns of $V$ are the $v_i$, and $L$ is a diagonal matrix with the eigenvalues on the diagonal. Consequently, you have that $\Sigma^{-1}=V^T L^{-1} V$ and you get that $$x^T \Sigma^{-1} x = x^T V^T L^{-1} V x = \sum_i \lambda_i^{-1} (v_i^T x)^2.$$ This would of course require you to store all eigenvalues, i.e., a full matrix $V$. But, if you happened to know that only a few of the eigenvalues of $\Sigma$ are small, say the first $m$, and the rest is so large that you can neglect all terms with $\lambda^{-1}_i$ for $i>m$, then you can approximate $$x^T \Sigma^{-1} x = \sum_{i=1}^n \lambda_i^{-1} (v_i^T x)^2 \approx \sum_{i=1}^m \lambda_i^{-1} (v_i^T x)^2.$$ This then only requires you to store $m$ vectors, instead of all $n$ eigenvectors.
Of course, in practice it is often equally or more difficult to compute the eigenvalues and eigenvectors, compared to simply solving $y=\Sigma^{-1}x$ multiple times.
• But then you're left with the task of finding the $m$ smallest eigenvalues of a matrix, which is not an easy task... – Yair Daon Oct 13 '15 at 20:22
• Correct. But might be worth it if you need to evaluate $x^T\Sigma^{-1}x$ many times. – Wolfgang Bangerth Oct 14 '15 at 11:17
• Bangreth Thank you for the reference. I checked SLEPc (though not extremely thoroughly) and it seems that the way eigenpairs are found is by (perhaps modified) Lanczos iteration. If one seeks the smallest $m$ eigenpairs, one has to find and store the all of the eigenpairs in memory. This typically is not possible for matrix free problems - it takes as much memory as storing an $n \times n$ matrix. If one wishes to use an inverse iteration - no guarantee is made for the order of the eigenpair(s) found. Am I missing something? – Yair Daon Oct 15 '15 at 0:40