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.

  • $\begingroup$ You might want to spell out what "MVM" means. $\endgroup$ Feb 5, 2017 at 7:31
  • $\begingroup$ @WolfgangBangerth done $\endgroup$
    – VF1
    Feb 5, 2017 at 13:13
  • 1
    $\begingroup$ 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. $\endgroup$ Mar 7, 2017 at 18:15
  • $\begingroup$ 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.... $\endgroup$ Mar 7, 2017 at 18:19
  • $\begingroup$ @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. $\endgroup$
    – VF1
    Mar 7, 2017 at 18:50

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

  • $\begingroup$ 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? $\endgroup$ Feb 7, 2017 at 7:22
  • $\begingroup$ 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$. $\endgroup$
    – VF1
    Feb 7, 2017 at 13:10
  • $\begingroup$ 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/… $\endgroup$ Nov 7, 2022 at 19:34
  • $\begingroup$ Relevant theorem from Vershynin probability book i.stack.imgur.com/coQrg.png , $m$ is number of samples $\endgroup$ Nov 7, 2022 at 19:36
  • $\begingroup$ @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. $\endgroup$
    – VF1
    Nov 8, 2022 at 3:09

Your Answer

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

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