Iterative single variable solutions in large linear systems

Question

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.

Answer 1

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.