cs.LGSep 16, 2026

Accurate Trace Estimation with Fewer Random Bits via Recursive TensorSketch

Authors: Mohammad Azhar KhanRameshwar PratapAmit Sharma

Abstract

We consider the problem of estimating the trace of an implicit matrix ARdp×dp\mathbf{A} \in \mathbb{R}^{d^p\times d^p} that can only be accessed through matrix-vector products queries. The \textit{Hutchinson trace estimator}% \cite{Girard1987algorithme, article-hutchinson} is a classical sketching method for this problem. Their estimator, Hm(A)=1mi=1mz(i)TAz(i),where  z(i)RdpH_{m}(\mathbf{A}) = \frac{1}{m} \sum_{i=1}^{m} {\mathbf{z}^{(i)}}^T \mathbf{A} \mathbf{z}^{(i)}, \quad \text{where } \ {\mathbf{z}^{(i)}}\in \mathbb{R}^{d^p}, and zj(i)N(0,1),j[dp]z^{(i)}_j \in {N}(0, 1), j\in [d^p], satisfies the following guarantees: (i) E[Hm(A)]=tr(A)\mathbb{E}[H_{m}(\mathbf{A})]=\operatorname{tr}(\mathbf{A}), and (ii) Var[Hm(A)]=2mAF2\mathrm{Var}[H_{m}(\mathbf{A})]=\frac{2}{m}||\mathbf{A}||_F^2. Generating one query vector z(i)\mathbf{z}^{(i)} requires O(dp)O(d^p) random bits; thus, mm queries require O(mdp)O(md^p) random bits, which can be prohibitive in large-scale applications. Recent work by Meyer et al.\cite{meyer2025hutchinsonsestimatorbadkroneckertraceestimation} proposes a variant of the Hutchinson trace estimator in which each query vector in Rdp\mathbb{R}^{d^p} is constructed as the Kronecker product of pp random vectors in Rd\mathbb{R}^d, requiring O(mpd)O(mpd) random bits for mm query vectors. The estimator of~\cite{meyer2025hutchinsonsestimatorbadkroneckertraceestimation} is unbiased; however, its variance grows exponentially with pp. In this work, we address this limitation by proposing a sketching-based estimator that requires O ⁣(p(d+m)logm)O\!\big(p (d + m)\log m\big) random bits, yields an unbiased estimate of the trace, and simultaneously achieves a variance bound that grows polynomially with pp.

Explore similar work

Jun 14, 2026stat.ML

Stochastic trace estimation with tensor train random vectors

Stochastic trace estimation is a standard tool for approximating the trace of a large-scale matrix available only through matrix-vector products. However, in tensor-structured settings, unstructured Gaussian or Rademacher test vectors may be prohibitively expensive to store and compute with, while cheaper rank-one tensor-product vectors can require sample complexities that grow exponentially with the tensor order. This work studies Gaussian random tensor train vectors as a structured alternative for stochastic trace estimation. We show that, with a suitable choice of the tensor train rank, random tensor train vectors recover dimension-independent guarantees for the Girard--Hutchinson estimator. In particular, a median-of-means variant with tensor train rank rd1r \geq d-1 achieves the same dependence on the accuracy ε\varepsilon and failure probability δδ as the classical estimator based on unstructured Gaussian vectors. We further prove an oblivious subspace injection result for sketches formed from independent Gaussian random tensor train vectors: tensor train rank rd1r\geq d-1 and O(ε2(k+log(1/δ)))\mathcal{O}(\varepsilon^{-2}(k+\log(1/δ))) samples suffice for a kk-dimensional target subspace. Finally, we investigate the use of such sketches within the Nyström++ framework. We show that the resulting estimator can achieve the desired O(ε1)\mathcal{O}(\varepsilon^{-1}) sample complexity under an additional spectral-tail condition. These results provide clarififcation on both the potential and the limitations of random tensor train vectors in stochastic trace estimation.
Zvonimir Bujanović, Daniel Kressner, Hrvoje Olić
Aug 11, 2026cs.DS

Improving TensorSketch Using Complex Random Variables

\texttt{TensorSketch} by~\cite{pham2013fast,kar2012random} provides efficient sketching algorithms for high-dimensional polynomial kernels xpRdp\vec{x}^{\otimes p} \in \R^{d^p}. \cite{kar2012random} uses dense Johnson-Lindenstrauss (JL)-type projections with computational cost O(pDd)O(pDd), where DD denotes the sketch dimension, whereas~\cite{pham2013fast} extends the sparse \texttt{CountSketch}\citep{count_sketch} algorithm, yielding a faster algorithm for high-dimensional sparse inputs with running time O(p(\nnzx+DlogD))O\big(p(\nnz{\vec{x}} + D \log D)\big). However, the variance of both estimators grows exponentially with the polynomial degree pp, scaling as 3p/D3^{p}/D. Recent work by\cite{pmlr-v206-wacker23a} showed that using complex-valued distribution reduces this dependence to 2p/D2^{p}/D for the approach of~\cite{kar2012random}. However, their method relies on dense JL-type projections with computational cost O(pDd)O(pDd) and does not extend to the algorithm of~\cite{pham2013fast}. In this work, we introduce a simple variant of \texttt{TensorSketch}\citep{pham2013fast} that achieves the same variance bound as\cite{pmlr-v206-wacker23a}, while retaining its advantage of the input-sparsity running time. We validate our results with supporting experiments on synthetic and real-world datasets.
Amit Sharma, Mohammad Azhar Khan, Rameshwar Pratap +1
Jun 22, 2026cs.LG

Dynamic estimation of slowly varying sequences

We consider the problem of sequentially approximating functions of each element in a slowly-varying sequence, i.e. one where the magnitude αiα_i of the difference between the elements at positions ii and i1i-1 is small. Recent work on implicit trace estimation shows that when αtα_t is small, reusing queries to past sequence elements can reduce the overall cost [Dharangutte & Musco, NeurIPS2021; Woodruff et al., NeurIPS2022]. We introduce a framework generalizing this to a variety of linear and nonlinear functions on diverse vector spaces, obtaining novel sequential estimation results for matrix powers, spectral densities, Monte Carlo integration, and a boundary value problem from partial differential equations~(PDEs). Furthermore, we develop a novel algorithm for use with this framework that locally scales the estimation budget with αtα_t, obtaining sharper path-length-style variation bounds of form O(i=1mαi)\mathcal O(\sum_{i=1}^mα_i) on the cost of estimating a sequence of length mm. This improves upon the previous implicit trace estimation bound of O(mmaxiαi)\mathcal O(m\cdot\max_iα_i) [Dharangutte & Musco, NeurIPS~2021], which is achieved by fixing the query budget using the worst-case αiα_i and is thus inefficient for stable sequences with rare bursts. Lastly, while all past work assumes a known bound on αiα_i, we show in certain cases how the changes can be estimated on-the-fly with (nearly) no added cost. In summary, our framework makes the sequential approximation toolkit general-purpose and adaptive while improving upon state-of-the-art-guarantees for dynamic trace estimation.
Prashant Gokhale, Mikhail Khodak, Sandeep Silwal