# Iterative single variable solutions in large linear systems

I have a system where $A$ is a large $n\times n$ marix with fast MVMs. It may have many nonzero entries (albeit in a structured way so as to allow fast MVMs), and is not necessarily diagonally dominant.

However, $A$ is weakly positive definite.

$Y$ is a $n\times m$ dense matrix.

I understand iterative Krylov-subspace-based methods exist for finding $X=A^{-1}Y$, and they perform well. Are there any optimizations that can be implemented if I am only interested in finding: $$\textbf{x}=\text{diag }Y^\top A^{-1}Y$$ In other words, I only want to solve for the $i$-th entry of the solution $Y^\top\textbf{x}_i$ where $A\textbf{x}_i=\textbf{y}_i$ for $Y=[\textbf{y}_i]_i$.

MVM stands for matrix-vector multiplication. In my application this takes only ~$n$ runtime compared to $n^2$ when dense.

• You might want to spell out what "MVM" means. Feb 5, 2017 at 7:31
• @WolfgangBangerth done
– VF1
Feb 5, 2017 at 13:13
• I believe the general case is still a research question at this point. I can think of several applications in statistics, ML, and physics that will immediately benefit from an efficient solution. Mar 7, 2017 at 18:15
• Do you know anything specific about the eigenvalues of the matrix $A$? If they are clustered, then we may use Krylov methods to construct a polynomial approximation $B = p(A)$ minimizing $\|BA - I \|_2$ that may be used as a cheap surrogate for the matrix inverse. Otherwise, Krylov methods may not be of much use.... Mar 7, 2017 at 18:19
• @RichardZhang That's a shame; you're probably right regarding it being an open question (but I suppose that's an answer to my question - would you consider posting it?) In my case there is no guarantee about clustering of eigenvalues.
– VF1
Mar 7, 2017 at 18:50

An approach by Papandreou and Yuille for diagonal $$A$$ relates variance estimation to the expectation of a quadratic form. The logic follows more generally: since $$A$$ is PD so is its inverse. Then $$Z\sim N(0, A^{-1})$$ is well-formed, and $$\textbf{x}=\mathbb{E}[(Y^\top Z)^2]$$.

If we can sample $$Z$$ then an approximation is possible (and converges). The linked paper has success when we can find the square root of $$A$$: solving $$AZ=A^{1/2}G$$ for samples $$G\sim N(0,I)$$ offers a way of sampling $$Z$$ since $$A^{-1/2}G\sim N(0,A^{-1})$$.

When an operator $$A^{1/2}$$ is available, this provides an estimation opportunity for variance.

Edit (5 years later):

Given only fast MVMs for $$A$$ we can actually achieve this. Using a matrix countour integral we can approximate MVMs of $$A^{-1/2}$$ using essentially a single iterative solve of a linear systems in $$A$$ (it's actuall $$O(1)$$ concurrent solves, where this $$O(1)$$ depends on quadrature error from the integral, and here concurrent solves means that Krylov iterates are shifts of one another, so we don't require additional applications of the operator $$A$$). As a result, we can directly sample $$Z=A^{-1/2}G$$ via this iterative routine.

The question then becomes, how many such solves do we need for our particular variance terms of interest $$\mathbb{E}[(Y^\top Z)^2]$$? At the end of the day, this is essentially asking what's the efficiency of estimating the covariance matrix for some $$W$$ where $$W\sim N(0,Y^\top A^{-1} Y)$$, which admits Wishart statistics.

Edit edit:

Thanks to @YaroslavBulatov we can add some meat to the last statement, which is that worst-case relative $$\epsilon$$ error in expectation in estimating $$\mathrm{cov}\ W$$ demands $$O(\epsilon^{-2} r \log m)$$ samples, where $$r=\mathrm{tr}\ \Sigma/\lambda_{\max}(\Sigma)$$ with $$\Sigma = Y^\top A^{-1} Y$$.

While $$r$$ is unknown, it may be practically feasible to estimate using Hutchinson's trick. Sampling new $$Z$$ and independent Rademacher vectors $$\rho$$, $$\mathbb{E}[\rho^\top Y^\top Z Z^\top Y \rho]=\mathrm{tr}\ \Sigma$$ and $$\lambda_{\max}$$ we can estimate with inverse power iteration. The composition $$Y\rho$$ is a single vector, so this second expectation is (conditionally on $$\rho$$) a chi-squared statistic with degrees of freedom equal to the number of samples $$Z$$ we use for this error estimation, which is sub-exponential and concentrates.

• How do you plan to solve $AZ=A^{1/2}G$, when the entire procedure has to be cheaper than solving a linear system with $A$, otherwise there would be no point in using it? Feb 7, 2017 at 7:22
• Right, I'm saying this is an approach only when you have cheap access to $A^{1/2}$ (say, $A$ is a large Kronecker product). If so, then you're only solving one linear system, $AZ=A^{1/2}G$, where $Z,A^{1/2}G$ are vectors. This is definitely better than the multiple linear systems required for $A^{-1}Y$, where $Y$ is a matrix. That's why I didn't accept my answer. Of course, this is only cheaper if the number of samples of $G$ that you make is less than the width of $Y$.
– VF1
Feb 7, 2017 at 13:10
• Regarding the edit part, Vershynin has some dimension-free way to estimate sample complexity of covariance estimation. If "stable rank" is small, the problem is easy -- mathoverflow.net/questions/395098/… Nov 7, 2022 at 19:34
• Relevant theorem from Vershynin probability book i.stack.imgur.com/coQrg.png , $m$ is number of samples Nov 7, 2022 at 19:36
• @YaroslavBulatov That's a good link, thank you. One can imagine by choice of $Y$ stable rank can be as large as $\min(m, n)$, which is somewhat disappointing since it feels like the fact that we care only of the diagonal and not the full covariance, we should have an easier time. Indeed, across $k$ samples each mean sum-of-squares estimate is $\chi_k^2$ which is sub-exponential. However, the maximum among each such sub-exponential isn't clearly bounded, because the samples aren't independent. I still feel like there's wiggle room in the case $n > m$ beyond the high-dimensional result.
– VF1
Nov 8, 2022 at 3:09