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.