Why does Hutchinson's trace estimator reduce computation complexity? Given a matrix $A$, we want to compute its trace, in which we can use a trick name Hutchinson's trace estimator
\begin{align}
tr(A) = tr(A\mathbb{E}[\epsilon \epsilon^T])=\mathbb{E}[tr(A \epsilon \epsilon^T)]=\mathbb{E}[tr(\epsilon^T A \epsilon)]=\mathbb{E}[\epsilon^T A \epsilon],
\end{align}
where $\epsilon \sim \mathcal{N}(0,1)$ is a standard normal distribution. Then the trace of $A$ can be estimated using Monte Carlo sampling.
Some literature said, by transforming the computation into quadratic form, the computation complexity of the trace calculators can be reduced. I can't understand it. Why can such a stochastic quadratic way reduce the complexity of trace computation?

It depends a lot on what accuracy you need as well as the form of the matrix.  In the application I had, there was a matrix $A=GP$ that was the product of a large (say $N\times M$) sparse matrix $G$ of genotypes and a projection matrix $P$ on to the residual space of a linear regression. This means $A$ wasn't itself sparse, but computing matrix-vector products $Ax$ or $A^TAx$ was still fast.
I needed the trace of $B=A^TA$ and $B^TB$ as part of an approximation to the distribution of a quadratic form, and they didn't need to be all that accurate. The trace of $B$ can be computed directly in $MN$ operations, but $tr(B^TB)$ requires explicitly computing all of $B$, taking $M^2N$ (assuming $M<N$, with both in the thousands).  Hutchinson's trace estimator with sample size $k$ could be computed in the time needed for $2k$ matrix multiplications by $A$ (small compared to $MNk$). The error in the estimator was $O_p(k^{-1/2})$, but that was ok (with $k\sim 500$)
The rest of the approximation involved the first few (~100)leading eigenvalues of $B$ and stochastic SVD gives an algorithm whose time is dominated by doing ~100 matrix-vector multiplications by $B$, so it was desirable for the trace part of the whole computation not to be hugely slower than this.
