arXiv CSopen access

Accurate Trace Estimation with Fewer Random Bits via Recursive TensorSketch

· arxiv_cs

arXiv CS · Papers · License: Open Access

machine learning, deep learning, neural networks

Accurate Trace Estimation with Fewer Random Bits via Recursive TensorSketch Mohammad Azhar Khan1*, Rameshwar Pratap1 and Amit Sharma1

arXiv:2609.18577v1 [cs.LG] 16 Sep 2026

Indian Institute of Technology Hyderabad, Kandi, Sangareddy, 502284, Telangana, India. *Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected]; [email protected]; Abstract We consider the problem of estimating the trace of an implicit matrix A ∈ p p Rd ×d that can only be accessed through matrix-vector products queries. The Hutchinson trace estimator [1, 2] is a classical sketching method for this problem. p 1 Pm (i) T Their estimator, Hm (A) = m Az(i) , where z(i) ∈ Rd , and i=1 z (i)

zj ∈ N (0, 1), j ∈ [dp ], satisfies the following guarantees: (i) E[Hm (A)] = 2 ||A||2F . Generating one query vector z(i) tr(A), and (ii) Var[Hm (A)] = m p requires O(d ) random bits; thus, m queries require O(mdp ) random bits, which can be prohibitive in large-scale applications. Recent work by Meyer et al. [3] proposes a variant of the Hutchinson trace estimator in which each query vector p in Rd is constructed as the Kronecker product of p random vectors in Rd , requiring O(mpd) random bits for m query vectors. The estimator of [3] is unbiased; however, its variance grows exponentially with p. In this work, we address this limitation by proposing a sketching-based estimator that requires O p(d + m) log m random bits, yields an unbiased estimate of the trace, and simultaneously achieves a variance bound that grows polynomially with p. Keywords: Trace estimation, Randomized Algorithms, Numerical Linear Algebra, Sketching Algorithms, Implicit linear operators

1 Introduction A central problem in scientific computing is the estimation of the trace of a large matrix p p A ∈ Rd ×d when explicit access to its entries is restricted. Instead, the matrix is accessible only through an oracle that returns matrix–vector products Ax for arbitrary p vectors x ∈ Rd . Under this restricted access model, the goal is to develop efficient algorithms that approximate the trace of matrix A while minimizing the number of oracle queries. The matrix–vector oracle model, also referred to as the implicit matrix model, is a widely adopted computational framework in the numerical linear algebra community [4–9]. The trace estimation problem in implicit matrix model can be solved exactly using D = dp oracle queries by using standard basis vectors PD the T e1 , e2 , . . . , eD via the following estimator tr(A) = i=1 ei Aei . Each term in the summation corresponds to a single diagonal entry of A, leading to a total of O(dp ) matrix–vector queries. However, the computational cost associated with such a large number of oracle queries is prohibitive. The seminal algorithm due to Girard and Hutchinson [1, 2], known as the Hutchinson trace estimator, provides an efficient approximation for trace estimation. Given p p an implicit matrix A ∈ Rd ×d , the Hutchinson trace estimator is defined as, p H (A) = zT Az, where z ∈ Rd with i.i.d. entries zi ∼ N (0, 1), for i ∈ [dp ]. The estimator satisfies the following guarantee E[H (A)] = tr(A), Var[H (A)] = 2||A||2F . Furthermore, to reduce the variance, the above procedure is repeated independently m times, and the final estimator is defined as the mean of these m estimators, that is, m

1 X (i) T Hm (A) = z Az(i) , m i=1

(i)

where z(i) ∈ Rd , zj ∈ N (0, 1), and j ∈ [dp ].

(1)

The estimator satisfies the following guarantee

E[Hm (A)] = tr(A),

Var[Hm (A)] =

2 ||A||2F . m

(2)

Subsequent work further improved the sample-complexity analysis of classical trace estimators. [10] derived sharper bounds for Gaussian, Rademacher, and unit-vector estimators, including a Hutchinson bound without the rank-dependent term appearing in the earlier analysis. Under the quadratic-form query model, [11] characterized optimal linear nonadaptive estimators and established lower bounds for multiplicative trace approximation. More recently, [12] studied nearly optimal high-probability traceestimation sketches under matrix-vector access. These works primarily seek to reduce the number of oracle queries, whereas our work studies the complementary objective of reducing the randomness required to construct queries in kronecker-structured spaces. The Hutchinson trace estimator, as stated in Equation (1), requires m random p vectors z(i) ∈ Rd . Consequently, the total number of random bits required by the estimator is O(mdp ). Structured random queries based on Kronecker-structured random vectors for trace estimation were proposed by [13]. Building on this idea, [3] addresses the challenge of random bits and suggests an estimator that requires sigp nificantly fewer random bits. Their estimator construct a query vector x ∈ Rd as a 2

Kronecker product of p independent random vectors in Rd , that is, x = x1 ⊗ · · · ⊗ xp , where xi ∈ Rd for i ∈ [p]. Therefore, generating a single random vector x requires O(dp) random bits, and the final estimator - formed by averaging m such estimators requires O(mdp) random bits, in contrast to the O(mdp ) random bits required by the Hutchinson trace estimator [2]. However, the main limitation of their approach is that p

the variance of their estimator grows exponentially with p, that is, O 3m (tr (A)) - making the estimator less accurate. This motivates the problem considered in this paper, which we state as follows: p p Problem Statement: Given an implicit matrix A ∈ Rd ×d , the goal is to design a trace estimation algorithm that requires asymptotically fewer random bits and simultaneously provides an accurate trace estimation. We draw inspiration from recent advances in sketching techniques to address this problem. In particular, [14] introduced Recursive TensorSketch, an efficient sketching method for approximating high-degree polynomial kernels. Their approach enables effective compression of polynomial kernels using a sketching dimension that scales only polynomially with the degree of the kernel function. We leverage Recursive p p TensorSketch to design a trace estimator for an implicit matrix A ∈ Rd ×d . We prove that the proposed estimator is unbiased, requires asymptotically fewer random bits than the classical Hutchinson trace estimator [2], and admits variance bounds with only polynomial dependence on the p. This constitutes an exponential improvement in the dependence on p over the variance bounds established for Kronecker-Hutchinson estimators in [3]. We summarize our key contributions as follows:

Our Contribution: • We propose a novel trace estimator for an implicit matrix A ∈ Rd ×d . Our esti p p ⊤ mator, tr Π A(Π ) (see Definition 5), leverages the Recursive TensorSketch p m×dp matrix Π ∈ R proposed by [14]. • We show that the proposed estimator is unbiased and derive a variance bound that is a polynomial of degree p. Further, when the input matrix A is a Positive Semi-Definite (PSD) matrix then the variance of our estimator achieves exponential improvement over the estimator proposed in [3]. Furthermore, the number of random bits required by our estimator is O p(d + m) log m , which is asymptotically better to that of required in [3]. Also, it is exponentially smaller than the Hutchinson trace estimator [2], which requires O(mdp ) random bits. • We further propose a complex-valued analogue of our estimator (Definition 6), in p which the entries of the Recursive TensorSketch matrix Πp ∈ Cm×d are sampled from complex random variables. We show that this variant achieves variance that is exponentially smaller than that of [3] conditioned that the input matrix is a PSD matrix, while simultaneously requiring asymptotically fewer random bits. p

There are two complementary approaches to reducing the amount of randomness required by randomized sketching algorithms. One approach is to redesign the sketching construction so that its random choices are shared through an underlying structure. This is the approach pursued in this paper through Recursive TensorSketch. A different approach is to retain an existing sketching construction while reducing its 3

randomness by implementing the underlying hash functions using tabulation-based hashing [15]. Prior work has shown that simple and double tabulation hashing can provide strong concentration guarantees for a variety of randomized algorithms and data structures, including Minwise Independent Permutations and Cuckoo Hashing [16, 17]. However, hash functions generated by tabulation hashing are generally not 4-wise independent and, therefore, cannot be directly used in standard trace estimation algorithms [1, 2], where such independence is required by the analysis. Nevertheless, tabulation hashing may offer an alternative approach for reducing the random seed required by existing trace-estimation sketches. An interesting direction for future work is to investigate whether such implementations can preserve the required JL moment properties and variance guarantees. Trace estimation is a fundamental primitive with numerous large-scale applications. Hutchinson trace estimator [2] has been used extensively as a key subroutine in various applications such as sublinear-time spectral density estimation [18], faster eigenvalue approximation [19], counting triangles in large graphs [20, 21], approximating spectral sums [22], and estimating ||A||F (using the well-known identity ||A||2F = tr AT A ) to name a few. Our proposed estimator can be plugged into these applications in place of [2] to yield a randomness-efficient algorithm with almost the same accuracy.

2 Related Work Trace estimation has a long history in randomized numerical linear algebra. The seminal algorithm by Girard and Hutchinson [1, 2], known as the Hutchinson trace estimator, provides an efficient method for approximating the trace. Given an p p implicit matrix A ∈ Rd ×d , the Hutchinson trace estimator is defined as H (A) = p z⊤ Az, where z ∈ Rd is a random vector with i.i.d. entries, typically drawn from N (0, 1) or a Rademacher distribution. The estimator satisfies E[H (A)] = tr(A) and Var(H (A)) = 2∥A∥2F . To reduce the variance, the estimator is repeated independently p mP times. Let z(1) , . . . , z(m) ∈ Rd be independent copies of z, and define Hm (A) = m 2 1 (i) ⊤ (i) 2 i=1 (z ) Az . Then, E[Hm (A)] = tr(A) and Var(Hm (A)) = m ∥A∥F . Each m dp p query requires generating a random vector z ∈ R , which uses O(d ) random bits, leading to a total randomness of O(mdp ) for m samples. In the classical setting, trace estimators are based on the form X := z⊤ Az, where z is a random query vector. Avron and Toledo [23] study several such estimators that differ in the choice of the distribution of z. In particular, they analyze the Hutchinson estimator under different choices of the query vector z, including the case where its entries are i.i.d. N (0, 1), the variant with i.i.d. Rademacher entries, and unit-vectorbased estimators in which z is sampled uniformly from the standard basis. They also study a mixed unit-vector estimator of the form XM := e⊤ FAF⊤ e, where each e is sampled uniformly from the standard basis vectors, and F is a fixed orthogonal transform (e.g. Hadamard matrix). Their work provides high-probability guarantees and highlights the trade-off between variance and the number of random bits used in these estimators. A recent work by Meyer et al. [3] proposed a Kronecker-structured trace estimator to reduce the number of random bits required for trace estimation. Their query

Estimator Hutchinson [23]

Variance (Gaussian) =

Hutchinson (Rademacher) [23] Normalized quotient [23] Unit [23]

vector

2 ∥A∥2F m 2 m p

Rayleigh = dm estimator = dm

O(mdp )

∥A∥2F −

n X

! A2ii

O(mdp )

i=1

P p d i

P p p

Mixed unit vector esti- – mator [23]

Randomness

d i

A2ii − (tr(A))2

O(mdp )

A2ii − (tr(A))2

O(mp log d)

2 3p Kronecker-Hutchinson ≤ tr(A) m (real) [3] 2 2p Kronecker-Hutchinson ≤ tr(A) m (complex) [3] 2 10p 100p2 Recursive TensorSketch ≤ + tr(A) m m2 (real) [this paper] 2 4p 16p2 Recursive TensorSketch ≤ + tr(A) 2 m m (complex) [this paper]

O(mp log d)

Bound on # samples for an (ε, δ)-approx. 2 20ε−2 ln δ 2r −2 6ε ln δ n2 κ2f (A) 2 ln δ 2r2 ε2 2 rD (A) 2 ln δ 2ε2 2 4n 4 −2 8ε ln ln δ δ

O(mpd)

3p 1 ln ε2 δ

O(mpd)

2p 1 ln ε2 δ

O p(d + m) log m

20p ε2 δ

O p(d + m) log m

8p ε2 δ

Table 1 Comparison of trace estimators for a fixed nonzero symmetric positive semidefinite p p matrix A ∈ Rd ×d , where n = dp , r = rank(A), κf (A) = λmax (A)/λ+ min (A), and rD (A) = n maxi Aii / tr(A). An (ε, δ)-approximation b t satisfies Pr[|b t − tr(A)| ≤ ε tr(A)] ≥ 1 − δ. The table highlights the trade-off among variance, randomness, and the number of samples required for an (ε, δ)-approximation. Classical Hutchinson estimators require O(mdp ) randomness, whereas Kronecker-Hutchinson estimators reduce this requirement to O(mpd) but incur an exponential dependence on p in both variance and sample complexity. In contrast, the proposed Recursive TensorSketch estimators require O p(d + m) log m randomness and have polynomial, rather than exponential, dependence on p.

vector x ∈ Rd is constructed as the Kronecker product of p independent random vectors in Rd , namely, x = x1 ⊗ · · · ⊗ xp , where xi ∈ Rd for each i ∈ [p]. Consequently, generating a single query vector x requires only O(dp) random bits, and an estimator obtained by averaging m independent samples requires O(mdp) random bits. This is substantially smaller than the O(mdp ) random bits required by the classical Hutchinson trace estimator [2]. However, this reduction in randomness comes at the p cost of increased variance. In particular, the variance of the estimator scales as O 3m , p whereas a complex-valued variant improves this dependence to O 2m . Thus, although the Kronecker-structured approach significantly reduces the randomness requirement by constructing each query vector from p independent vectors in Rd , the exponential dependence of the variance on p can make the estimator increasingly inaccurate as p grows. In this work, we address this challenge by proposing an estimator that requires asymptotically fewer random bits while ensuring that its variance grows only polynomially with p. Our work is inspired by the work of [14], which proposed a recursive sketching algorithm Recursive TensorSketch for compressing polynomial kernels. We show that Recursive TensorSketch can also be leveraged to design a trace estimator that significantly reduces the number of random bits p p required while maintaining low variance. For an implicit matrix A ∈ Rd ×d , our 5

p p ⊤ estimator T ( A ) := tr Π A (Π ) is unbiased and admits a variance bound of 2 10p 100p2 O tr(A) , which is polynomial of p. Our estimator yields a expom + m2 nential over Kronecker-Hutchinson estimators [3], whose variance scales pimprovement 3 (tr(A))2 2p (tr(A))2 as O and O in the real and complex settings, respectively. m m Moreover, the randomness complexity of our estimator is O p(d + m) log m , which is asymptotically smaller than the O(mpd) randomness required by KroneckerHutchinson estimator [3] and exponentially smaller than the O(mdp ) randomness required by classical Hutchinson estimator [1, 2]. We further extend our framework to a complex-valued setting, obtaining improved variance bounds with fewer random bits than the corresponding estimator of [3]. A standard way to evaluate a randomized trace estimator is through an (ε, δ )approximation guarantee [23]. For a fixed nonzero positive semidefinite matrix A, an estimator b t is called an (ε, δ )-approximation of tr(A) if

Pr b t − tr(A) ≤ ε tr(A) ≥ 1 − δ.

(3)

Here, ε specifies the allowed relative error and δ specifies the failure probability. This guarantee is important because it translates variance or concentration bounds into a required sample or sketch size, thereby allowing different trace estimators to be compared in terms of accuracy. We summarize our comparison with the baseline methods in Table 1, which highlights the trade-offs among variance, randomness, and the number of samples required to obtain an (ε, δ )-approximation.

3 Background Notation. We denote vectors by lowercase bold letters (e.g., x) and matrices by uppercase bold letters (e.g., M). For a matrix M ∈ Rn×n , tr(M) denotes its trace and ∥M∥F its Frobenius norm. We write M ⪰ 0 to indicate that M is symmetric positive semi-definite. For a positive integer d, we denote [d] := {1, 2, . . . , d}. Kronecker product p is denoted by ⊗, and for vectors x1 , . . . , xp ∈ Rd , we write x = x1 ⊗ · · · ⊗ xp ∈ Rd . We use E[·] and Var(·) to denote expectation and variance, respectively. Throughout the paper, i ∈ [m] indexes sketch dimensions. For a complex vector or matrix in the field C , we denote by (·)∗ its conjugate transpose. Finally, 1[·] denotes the indicator function. We first state the classical Hutchinson Trace Estimator and its concentration guarantee. p

Theorem 1 (Hutchinson Trace Estimator [1, 2]) Let A ∈ Rd ×d be any implicit matrix. Then, the trace estimator is defined as H(A) := p z⊤ Az, where z ∈ Rd such that zi ∼ N (0, 1). p Let z(1) , . . . , z(m) ∈ Rd be i.i.d. copies of z, then the final estimator is defined as follows m 1 X (i) ⊤ (i) Hm (A) = z Az . (4) m i=1

Then,

E[Hm (A)] = tr(A),

Var(Hm (A)) =

2 ∥A∥2F . m

(5)

Theorem 2 (High-Probability Error Bound [21]) Let A ⪰ 0 and let Hm (A) be the estimator defined in Equation (4) using Rademacher or Gaussian vectors. Then for any ε, δ ∈ (0, 1), it suffices to choose m = O

log(1/δ) ε2

samples to guarantee

Pr(|Hm (A) − tr(A)| ≤ ε tr(A)) ≥ 1 − δ.

(6)

3.1 Trace Estimation via Kronecker-Matrix vector product [3] considered a variant of Hutchinson trace estimator where the problem is estimatp p ing the trace of an implicit matrix A ∈ Rd ×d that can only be accessed through Kronecker-matrix-vector products. That is, for any Kronecker-structured vector that is, x = x1 ⊗· · ·⊗xp , where random vector xi ∈ Rd for i ∈ [p], Kronecker-matrix-vector product Ax can be computed. Their estimator is termed as Kronecker-Hutchinson estimator and defined as follows: T := x⊤ Ax. They propose several estimators, each corresponding to different choices of distributions from which the random vectors xi are sampled. Theorem 3 (Variance for real-valued Kronecker random vectors [3, Theorem 5.4]) Let A ∈ p p Rd ×d be a PSD matrix. Let x = x1 ⊗ · · · ⊗ xp , where x1 , . . . , xp ∈ Rd are independent and identically distributed random vectors. Then, all the following estimators are unbiased, and satisfy the following variance bounds Gaussian:

if xi ∼ N (0, Id ), Var[x⊤ Ax] ≤ 3p (tr(A))2 ,

Rademacher:

Uniform sphere:

if entries of each xi are i.i.d. in {−1, +1}, p 2 Var[x⊤ Ax] ≤ 3 − (tr(A))2 , d if each xi is drawn uniformly from Sd−1 , p 6 Var[x⊤ Ax] ≤ 3 − (tr(A))2 . d+2

The bounds in Theorem 3 exhibit an exponential dependence on the parameter p, rendering the Kronecker–Hutchinson estimator inefficient for large values of p. They further demonstrate that using complex-valued random vectors leads to improved bounds with a smaller exponential factor. The following theorem summarizes those guarantees. Theorem 4 (Variance for Complex-Valued Structures (Theorem 6.2 and Lemma 6.3 of [3])) p p Let A ∈ Rd ×d be a PSD matrix, and x = x1 ⊗ · · · ⊗ xp , where x1 , . . . , xp ∈ Cd are i.i.d. random vectors. Then, all the following estimators are unbiased, and satisfy the following variance bounds Complex Gaussian:

if xi = √1 (ri + imi ), with ri , mi ∼ N (0, Id ), 2

∗

Var[x Ax] ≤ 2p (tr(A))2 ,

Π8 = S21 (w1 ⊗ w2 )

Sℓi

→ Tensor Sketch (p=2)

Tj x

→ Count Sketch

S21

w1 := S41 (v1 ⊗ v2 )

w2 := S42 (v3 ⊗ v4 )

S41

v1 := S81 (T1 x ⊗ T2 x)

S42

v2 := S82 (T3 x ⊗ T4 x)

S81

T1 x

v3 := S83 (T5 x ⊗ T6 x)

v4 := S84 (T7 x ⊗ T8 x)

S83

S84

S82

T2 x

T3 x

T4 x

T5 x

T6 x

T7 x

T8 x

Fig. 1 Recursive TensorSketch construction for p = 8. Each Tj denotes a CountSketch, while Siℓ denotes a TensorSketch operator. Intermediate vectors are combined recursively.

Complex Rademacher:

Complex sphere:

if each entry of xi is drawn i.i.d. from p 1 (tr(A))2 , Var[x∗ Ax] ≤ 2 − d

± √1 , ± √i 2

if each xi is uniformly distributed on the complex sphere , p 2 Var[x∗ Ax] ≤ 2 − (tr(A))2 . d+1

We address the limitations of the Kronecker–Hutchinson estimator by designing an alternative estimator that leverages Recursive TensorSketch proposed by [14]. In their work, the authors develop this technique in the context of sketching high-degree polynomial kernels, demonstrating that tensor product structures can be efficiently compressed via recursive linear-mappings while preserving inner-product similarity to a high degree of accuracy. We state their sketching algorithm in the following subsection.

3.2 Introduction to Recursive TensorSketch We begin by presenting the CountSketch [24] algorithm, which enables fast dimensionality reduction for high-dimensional vectors. We then describe TensorSketch [25, 26] of degree 2, which extends this idea to efficiently compress vectors formed via the Kronecker product of two vectors. Definition 1 (CountSketch [24]) Given an input vector y ∈ Rd , the CountSketch is a randomized linear map T ∈ Rm×d that maps y to a lower-dimensional vector z = T y ∈ Rm . The CountSketch matrix T is constructed by two hash functions: (a) h : [d] → [m] a 3-wise independent hash function, and (b) s : [d] → {1, −1} a 4-wise independent random sign function. The j th entry of vector z ∈ Rm is computed as, X zj = s(i) yi , ∀j ∈ {1, . . . , m}. h(i)=j

The time complexity of computing the CountSketch is O(nnz(y)), which in the worst case can be O(d). Furthermore, CountSketch provides an unbiased estimator and variance of this estimator is Var[∥Ty∥22 ] ≤

2∥y∥42 m .

TensorSketch extends the idea of CountSketch to tensor products and allows them to be sketched efficiently. Definition 2 (TensorSketch of Degree Two [25, 26]) Let h1 , h2 : [d] → [m] be 3-wise independent hash functions, and σ1 , σ2 : [d] → {−1, +1} be 4-wise independent random sign 2 functions. Then the TensorSketch of degree two S ∈ Rm×d is defined ∀ r ∈ [m], i1 , i2 ∈ [d], as follows Sr,(i1 ,i2 ) = σ1 (i1 ) · σ2 (i2 ) · 1 [h1 (i1 ) + h2 (i2 ) ≡ r (mod m)] . (7) TensorSketch provides an unbiased estimator of the squared ℓ2 -norm, whose variance is h i bounded by Var ∥S(x ⊗ x)∥22 ≤

8∥x∥42 d m . Furthermore, for any x ∈ R , the sketch S(x ⊗ x)

can be computed in O(m log m + nnz(x)) time using the Fast Fourier Transform (FFT). p

Given a vector x ∈ Rd of the form x = x1 ⊗ · · · ⊗ xp , Recursive TensorSketch provides an efficient sketching procedure that avoids the explicit construction of x. The method proceeds by first applying independent CountSketch transformations to each component vector xi , for i ∈ [p], and then recursively combining the resulting sketches using degree-two TensorSketch operations, producing a hierarchical tree-structured representation refer to Figure 1. p

Definition 3 (Recursive TensorSketch [14]) Given a vector x ∈ Rd , where p is a power of two, the Recursive TensorSketch is a randomized linear map p

Πp : Rd → Rm ,

defined as

Πp := Qp · Tp , where

• Tp = T1 ⊗ T2 ⊗ · · · ⊗ Tp , with each Ti ∈ Rm×d for i ∈ [p] a CountSketch matrix (Definition 1), ℓ/2 ℓ • Qp = S2 · S4 · · · Sp/2 · Sp , with each Sℓ ∈ Rm ×m a Kronecker product of matrices 2 Sjℓ ∈ Rm×m , • each Sℓj is a TensorSketch matrix of degree 2 (Definition 2), and Sℓ = Sℓ1 ⊗ Sℓ2 ⊗ · · · ⊗ Sℓℓ/2 . When x is given in the form of Kronecker product of p vectors, i.e., x = x1 ⊗ · · · ⊗ xp with xi ∈ Rd for all i ∈ [p], the Recursive TensorSketch can be computed efficiently in p time O(p m log m + pd). In contrast, when x is an arbitrary vector in Rd without explicit p p Kronecker structure, computing Π x requires O(md ) time. Definition 4 (Definition 18 of [14]: JL Moment Property) For every positive integer t and every δ, ε ≥ 0, a distribution over random matrices M ∈ Rm×d has the (ε, δ, t)-JL Moment Property if h i ∥Mx∥22 − 1 Lt ≤ ε δ 1/t and E ∥Mx∥22 = 1

for every unit vector x ∈ Rd .

Now, we state some useful results from [14] that will be used in our proofs. Lemma 5 (Lemma 9 of [14]: Two-vector JL Moment Property) For any x, y ∈ Rd , if M has the (ε, δ, t)-JL Moment Property, then ⟨Mx, My⟩ − ⟨x, y⟩ Lt ≤ ε δ 1/t ∥x∥2 ∥y∥2 . Lemma 6 (Lemma 12 of [14]: Factorisation of Πp ) For any integer p which is a power of p two, let Πp : Rd → Rm be Recursive TensorSketch defined in Definition 3, for sketches 2 (j) Sℓi : Rm → Rm and Tj : Rd → Rm . Then there exist matrices M(i) i∈[p−1] , M′ j∈[p]

and integers (ki )i∈[p−1] , (ki′ )i∈[p−1] , (lj )j∈[p] , (lj′ )j∈[p] , such that

Πp = M(p−1) · · · M(1) · M′(p) · · · M′(1) , (j)

where M(i) = Iki ⊗Sℓi ⊗Iki′ and M′ = Iℓj ⊗Tj ⊗Iℓ′j , with Sℓi and Tj independent instances of TensorSketch of Degree-2 and CountSketch, respectively, for every i ∈ [p − 1] and j ∈ [p]. Lemma 7 (Lemma 14 of [14]: JL Moment Property under tensor wraps) If the matrix S has the (ε, δ, t)-JL Moment Property, then for any positive integers k, k′ , the matrix M = Ik ⊗ S ⊗ Ik′ has the (ε, δ, t)-JL Moment Property. Lemma 8 (Lemma 15 of [14]: Composition lemma for the second moment) For any ε, δ ≥ 0 d2 ×d1 , . . . , M(k) ∈ Rdk+1 ×dk are independent random matrices and any integer k, if M(1) ∈R each with the √ε , δ, 2 -JL Moment Property, then the product matrix M = M(k) · · · M(1) 2k satisfies the (ε, δ, 2)-JL Moment Property. Corollary 9 (Corollary 16 of [14]: Second moment property for Πp ) For any power-ofp two integer p, let Πp : Rd → Rm be defined in Definition 3, where both base distributions 2 ε Siℓ : Rm → Rm and Tj : Rd → Rm satisfy the √4p+2 , δ, 2 -JL Moment Property. Then Πp satisfies the (ε, δ, 2)-JL Moment Property.

The exponential variance growth of the Kronecker-Hutchinson estimator motivates alternative approaches for tensor-structured trace estimation. Although complexvalued structures offer partial improvement, they do not remove this dependence. The Recursive TensorSketch provides a structured way to compress tensor products, suggesting a more efficient estimator. We introduced trace estimator using Recursive TensorSketch and analyse its variance in the following section.

4 Trace Estimator using Recursive TensorSketch In this section, Definition 5 introduces a Recursive TensorSketch-based trace estimator for positive semidefinite matrices given in implicit form. Our analysis uses the 10

independent-layer factorization established in Lemma 6. Lemma 10 first establishes expectation and variance bounds for a single sketch satisfying the second-moment JL property, and Lemma 12 extends these bounds to the Kronecker-wrapped sketching operators arising in the recursive construction. Theorem 11 then applies these bounds conditionally across the 2p − 1 independent layers to prove unbiasedness and derive the variance bound. Finally, Lemma 13 analyzes the number of random bits required to construct the estimator, and Theorem 14 presents the corresponding concentration guarantee. p

Definition 5 (Recursive TensorSketch (RTS) Trace Estimator) Let Πp ∈ Rm×d denote thep Recursive TensorSketch matrix stated in Definition 3. For an implicit PSD matrix A ∈ p Rd ×d , its trace estimator is defined as follows: T (A) := tr Πp A(Πp )⊤ . (8) Lemma 10 (Expectation and Variance bound for a single-layer of Sketch) Let S ∈ Rm×di be (i) (i) (i) 2 (i) ci . δ0 ≤ m matrix satisfying the (ε0 , δ0 , 2)-JL Moment Property (Definition 4), with ε0 di ×di e Then, for every positive semidefinite matrix B ∈ R , h i e ⊤ ) = tr(B), e and Var tr(SBS e ⊤ ) ≤ ci tr(B) e 2. E tr(SBS (9) m

e The unbiasedness follows from linearity Proof The proof uses the eigendecomposition of B. of the trace, while the variance bound follows by applying the second-moment JL property to each eigenvector and then using Minkowski’s inequality. Let R X e = λr ur u⊤ B r r=1

e where λr ≥ 0 and {ur }R is an orthonormal set. By linearity be an eigendecomposition of B, r=1 of the trace, R X e ⊤ = tr SBS λr ∥Sur ∥22 . r=1

h i The unbiasedness property of S gives E ∥Sur ∥22 = 1. Therefore, R h i h i X e ⊤ = E tr SBS λr E ∥Sur ∥22 r=1

R X

e λr = tr(B).

r=1

Define Zr := ∥Sur ∥22 − 1. By unbiasedness property, E[Zr ] = 0, and the second-moment JL property gives (i) (i) 1/2 ∥Zr ∥L2 ≤ ε0 δ0 .

1/2 Here, for any random variable Y , ∥Y ∥L2 := E[|Y |2 ] denotes its L2 norm. Since R X e ⊤ − tr(B) e = tr SBS λr Zr , r=1

Minkowski’s inequality, which is the triangle inequality for the L2 norm, gives R X r=1

≤

λr Zr L2

R X

λr ∥Zr ∥L2

r=1 (i)

≤ ε0

(i) 1/2

δ0

R X

λr

r=1 (i)

= ε0

(i) 1/2

δ0

e tr(B).

This argument does not require the random variables Zr to be independent. Since the sum P R r=1 λr Zr has mean zero, its squared L2 norm equals its variance. Squaring the preceding inequality therefore yields e 2 e ⊤ ≤ ε(i) 2 δ (i) tr(B) Var tr SBS 0 0 c e 2. ≤ i tr(B) m □

Having established the expectation and variance bounds for a single sketching layer in Lemma 10, we now state the main guarantee and subsequently prove it for the complete Recursive TensorSketch trace estimator. Theorem 11 (Unbiasedness and Variance of RTS Trace Estimator) Let T (A) be the estimator defined in Definition 5. Then, 2 10p 100p2 tr(A) . (10) E[T (A)] = tr(A), and Var(T (A)) ≤ + m m2

To prove Theorem 11, recall from Lemma 6 that the Recursive TensorSketch matrix Πp can be expressed as a product of 2p − 1 independent sketching matrices. We first establish expectation and variance bounds for one Kronecker-wrapped sketching layer in Lemma 10. We then apply these bounds successively to all 2p − 1 layers through a composition argument to obtain the variance bound for Πp . Lemma 12 (Unbiasedness and Variance Guarantees for a Layer) Let S ∈ Rm×di satisfy the assumptions of Lemma 10. Let da and db be positive integers, and let Ida ∈ Rda ×da and Idb ∈ Rdb ×db denote the identity matrices acting on the tensor components before and after the component sketched by S, respectively. Define M(i) = Ida ⊗ S ⊗ Idb . Then, for every positive semidefinite matrix B ∈ Rda di db ×da di db , we have h i E tr M(i) B(M(i) )⊤ = tr(B), 2 c Var tr M(i) B(M(i) )⊤ ≤ i tr(B) . m

(11) (12)

Proof Let matrix of suitable dimension be B ∈ Rda di db ×da di db where da , di , and db are suitable dimensions. Let the full space be indexed by the tuple (a, j, b) corresponding to the dimensions da , ′ ′ di , and db respectively. We partition the matrix B into blocks B(a,a ,b,b ) ∈ Rdi ×di by fixing the outer dimensions at indices (a, a′ ) ∈ [da ]2 and (b, b′ ) ∈ [db ]2 . The sketching operator at layer i is defined as the Kronecker product: M(i) = Ida ⊗ S ⊗ Idb . Given that the base sketch S ∈ Rm×di and the identity matrices are Ida ∈ Rda ×da and Idb ∈ Rdb ×db , the dimensions of M(i) multiply across the tensor product. Thus, M(i) has dimensions: M(i) ∈ R(da mdb )×(da di db ) . When we sketch B using M(i) , the matrix multiplication aligns as follows:

• M(i) is of size (da mdb ) × (da di db ) • B is of size (da di db ) × (da di db ) • (M(i) )⊤ is of size (da di db ) × (da mdb )

′

We now partition B into di × di blocks denoted by B(a,a ,b,b ) , such that: B=

da X

db X

′

(a,a ,b,b ) ⊗ (eb e⊤ (ea e⊤ a′ ) ⊗ B b′ ).

a,a′ =1 b,b′ =1

We now apply the sketching operator M(i) to B. Using the mixed-product property of Kronecker products, (X ⊗ Y)(U ⊗ V) = (XU ⊗ YV), we obtain: M(i) B(M(i) )⊤ = Ida ⊗ S ⊗ Idb B Ida ⊗ S⊤ ⊗ Idb X (a,a′ ,b,b′ ) ⊤ S ⊗ Idb eb e⊤ = Ida ea e⊤ a′ Ida ⊗ SB b′ Idb a,a′ ,b,b′

(a,a′ ,b,b′ ) ⊤ (ea e⊤ S ⊗ (eb e⊤ a′ ) ⊗ SB b′ ).

X a,a′ ,b,b′

Finally, we apply the trace operator. The trace of a Kronecker product is the product of the traces, i.e., tr(X ⊗ Y) = tr(X) tr(Y). Applying this to our summation gives: da X tr M(i) B(M(i) )⊤ =

db X

(a,a′ ,b,b′ ) ⊤ tr(ea e⊤ S · tr(eb e⊤ a′ ) · tr SB b′ ).

(13)

a,a′ =1 b,b′ =1

Recall that the trace of an outer product of basis vectors is the inner product of the vectors: ′ ⊤ tr(ea e⊤ a′ ) = ⟨ea′ , ea ⟩. Thus, tr(ea ea′ ) is 1 if a = a and 0 otherwise. Therefore we can write: db da X X tr M(i) B(M(i) )⊤ = tr SB(a,a,b,b) S⊤ .

(14)

a=1 b=1

Using the linearity of matrix addition, we define a matrix B̃ ∈ Rdi ×di on the single subspace where the random sketch S operates: B̃ :=

db da X X a=1 b=1

B(a,a,b,b) .

(15)

Substituting B̃ back into Equation (14), the trace of the entire high-dimensional Kronecker layer simplifies to a standard matrix sketch trace on the lower-dimensional space: tr(M(i) B(M(i) )⊤ ) = tr(SB̃S⊤ ).

(16)

Applying Lemma 10 and evaluating it further gives h h i i e E tr M(i) B(M(i) )⊤ = E tr(SB̃S⊤ ) = tr(B)   db db da X da X X X (a,a,b,b) = = tr  B tr B(a,a,b,b) a=1 b=1 db X da X di X

(17) (18)

a=1 b=1

B(a,j,b),(a,j,b) = tr(B),

(19)

e ⊤) Var tr M(i) B(M(i) )⊤ = Var tr(SBS c e 2 ≤ i tr(B) m 2 c = i tr(B) . m

(20)

a=1 b=1 j=1

and similarly,

(21) (22) □

We now apply bounds of single layer established in Lemma 12 successively to all 2p− 1 layers through a composition argument to obtain the variance bound for Πp .

4.1 Proof of Theorem 11 via Composition The proof of Theorem 11 relies on the factorization of the Recursive TensorSketch matrix Πp into a sequence of mutually independent random sketching matrices. The argument applies the single-layer trace moment bounds conditionally at each layer and uses the PSD structure of the input matrix. Proof Let A ⪰ 0 and let k = 2p − 1. By Lemma 6, the Recursive TensorSketch matrix admits the independent-layer factorization Πp = M(k) M(k−1) · · · M(1) . Here, the first p factors correspond to the CountSketch maps at the leaf level, while the remaining p − 1 factors correspond to the degree-2 TensorSketch maps at the internal nodes of the recursive tree. Define A0 := A, Ai := M(i) Ai−1 (M(i) )⊤ ,

i = 1, . . . , k,

and let Xi := tr(Ai ). Since A0 ⪰ 0 and each matrix Ai is obtained from Ai−1 by multiplication with M(i) and its transpose, positive semidefiniteness is preserved at every step. Therefore, Ai ⪰ 0

for every i = 0, . . . , k.

(23)

Let Fi−1 represent all the random choices made in the first i−1 layers. Putting Condition on this, the matrix Ai−1 is fixed and positive semidefinite, while M(i) remains independent and random. Therefore, Lemma 12 gives E[Xi | Fi−1 ] = Xi−1 , (24) ci 2 Var(Xi | Fi−1 ) ≤ Xi−1 . (25) m where, by the variance of CountSketch and TensorSketch of degree-2 given in Definitions 1 and 2, respectively, ( 2, i = 1, . . . , p, ci = 8 = 32 − 1, i = p + 1, . . . , 2p − 1. Here, 2 and 8 are the corresponding second-moment JL constants. Expectation. Taking expectations in Equation (24) and applying the tower property yields E[Xi ] = E[Xi−1 ]. Iterating over all k layers gives E[Xk ] = E[X0 ] = tr(A). Since Xk = T (A), the trace estimator is unbiased: E[T (A)] = tr(A).

(26)

Variance. Using the conditional second-moment identity together with Equation (24) and (25), we obtain E[Xi2 | Fi−1 ] = Var(Xi | Fi−1 ) + (E[Xi | Fi−1 ])2 c 2 2 + Xi−1 ≤ i Xi−1 m c 2 = 1 + i Xi−1 . m Taking expectations and iterating from i = 1 to i = k gives E[Xk2 ] ≤

k Y c 1 + i X02 m

i=1

k Y 2 c 1 + i tr(A) . m

(27)

i=1

Combining Equation (26) and (27), we obtain 2 Var(T (A)) = E[Xk2 ] − E[Xk ] " k # Y 2 ci ≤ − 1 tr(A) . 1+ m i=1

The sum of the layer constants is k X

ci = 2p + 8(p − 1)

i=1

= 10p − 8.

(28)

Using 1 + x ≤ ex for x ≥ 0, we have k Y c 1 + i ≤ exp m

! k 1 X ci m i=1 i=1 10p − 8 = exp . m Consequently, the general variance bound is 2 10p − 8 − 1 tr(A) . Var(T (A)) ≤ exp m If m ≥ 10p − 8, then 10p − 8 0≤ ≤ 1. m x 2 The inequality e − 1 ≤ x + x , valid for 0 ≤ x ≤ 1, therefore gives 2 (10p − 8)2 10p − 8 Var(T (A)) ≤ + tr(A) 2 m m 2 10p 100p2 + ≤ tr(A) . m m2 This proves the claimed variance bound.

(29)

□

We now analyze our trace estimator’s randomness complexity. The following lemma separately counts the random bits required for the leaf-level CountSketch matrices in Tp and the internal degree-2 TensorSketch matrices in Qp . Lemma 13 (Randomness Complexity of Recursive TensorSketch) Let Πp = Qp Tp be the Recursive TensorSketch matrix as defined in Definition 3. Then, the total number of random bits required to construct Πp are O p(d + m) log m . Consequently, the estimator T (A) := tr Πp A(Πp )⊤ can be implemented using O p(d + m) log m random bits.

Proof We decompose the randomness required to construct Πp = Qp Tp into two parts. (1) Randomness for Tp . Recall that Tp = T1 ⊗ · · · ⊗ Tp , where each Ti ∈ Rm×d is a CountSketch matrix. Each Ti is specified by:

• a hash function hi : [d] → [m], requiring ⌈log2 m⌉ bits per coordinate (to store the index j ∈ [d] is mapped into which index j ′ ∈ [m]), • a sign function si : [d] → {±1}, requiring 1 bit per coordinate.

Thus, each Ti requires d ⌈log2 m⌉ + 1 random bits, and over all p matrices, Number of bits in Tp = pd ⌈log2 m⌉ + 1 .

(30)

(2) Randomness for Qp . By definition, Qp = S2 · S4 · · · Sp , 2

where each Sℓ is a Kronecker product of ℓ/2 matrices Sℓj ∈ Rm×m , each being a degree2 TensorSketch as defined in Definition 2. In particular, each Sℓj is constructed using two 3-wise independent hash functions and two 4-wise independent random sign functions, as specified in Definition 2. Thus, each Sℓj requires:

• two hash functions h1 , h2 : [m] → [m], requiring O(log m) bits per coordinate for each hash function, and • two sign functions σ1 , σ2 : [m] → {−1, +1}, requiring 1 bit per coordinate for each sign function. Since each Sℓj acts on m2 coordinates but is implemented implicitly via hash functions, its description requires O(m(log m + 1)) random bits. At level ℓ, there are ℓ/2 such matrices, hence Number of bits in Sℓ = O ℓ m(log m + 1) . Summing over levels ℓ = 2, 4, . . . , p, Number of bits in Qp =

O ℓ m(log m + 1)

ℓ

= O pm(log m + 1) .

(31)

(3) Total randomness. Combining both parts, we have O p(d + m) log m .

(32)

This proves the stated bound. The final asymptotic form follows immediately.

□

We conclude this section analysis by deriving a concentration guarantee for the estimator. The following theorem combines the unbiasedness and variance bound from Theorem 11 with Chebyshev’s inequality to obtain a relative failure-probability bound and a sufficient condition on the sketch dimension m for an (ε, δ )-approximation. p

Theorem 14 [Concentration Analysis of RTS Trace Estimator] Let A ∈ Rd ×d be a fixed nonzero symmetric positive semidefinite matrix and let T (A) be the trace estimator defined in Definition 5. Then, for every ε > 0, 1 10p 100p2 Pr[|T (A) − tr(A)| ≥ ε tr(A)] ≤ 2 . (33) + m ε m2 Moreover, for 0 < ε ≤ 1 and δ ∈ (0, 1), T (A) is an (ε, δ)-approximation whenever 20p m≥ 2 . ε δ

(34)

Proof From Theorem 11, we have E[T (A)] = tr(A), 2 100p2 10p + tr(A) . Var(T (A)) ≤ 2 m m

(35) (36)

Since A is nonzero and positive semidefinite, tr(A) > 0. Therefore, Chebyshev’s inequality gives Var(T (A)) 1 10p 100p2 Pr[|T (A) − tr(A)| ≥ ε tr(A)] ≤ ≤ + . 2 m ε2 m2 ε2 tr(A) To make the failure probability at most δ, it is sufficient that 1 10p 100p2 + ≤ δ. m ε2 m2

(37)

Multiplying both sides by m2 ε2 gives ε2 δ m2 − 10pm − 100p2 ≥ 0. Solving this quadratic inequality for m yields p 5p m ≥ 2 1 + 1 + 4ε2 δ . ε δ √ √ Since 0 < ε ≤ 1 and δ ∈ (0, 1), 1 + 1 + 4ε2 δ ≤ 1 + 5 < 4. Hence, p 5p 20p 1 + 1 + 4ε2 δ < 2 . 2 ε δ ε δ Therefore, the condition 20p m≥ 2 ε δ is sufficient to make the failure probability at most δ.

(38)

(39)

□

While the Recursive TensorSketch yields favourable variance bounds in the real-valued setting, further improvements can be obtained by considering complexvalued sketching constructions. As observed in prior work [3], complex random projections often exhibit improved concentration properties and reduced variance compared to their real-valued counterparts. Motivated by this, we extend the Recursive TensorSketch framework to the complex domain and analyze the resulting trace estimator in the section below.

5 Trace Estimator using Complex Recursive TensorSketch In this section, Definition 6 introduces the Complex Recursive TensorSketch trace estimator. Using the independent-layer factorization from Lemma 6, Theorem 15 establishes the unbiasedness and variance bound of the estimator. Lemma 16 then analyzes the number of random bits required to construct the complex sketch. Finally, Theorem 17 derives the corresponding concentration guarantee. p

Definition 6 (Complex Recursive TensorSketch (RTS) Trace Estimator) Let Πp ∈ Cm×d denote the Complex Recursive TensorSketch matrix constructed as in Definition 3, so that Πp = Qp Tp , except that the real-valued sign functions used in the CountSketch and degree-2 TensorSketch matrices are replaced by independent hash functions whose values arep unip formly distributed over the fourth roots of unity {1, i, −1, −i}. For a PSD matrix A ∈ Rd ×d , we define TC (A) := tr Πp A(Πp )∗ , (40) where (·)∗ denotes the conjugate transpose.

Theorem 15 [Unbiasedness and Variance of Complex Recursive TensorSketch Trace Estimator] Let TC (A) be the estimator defined in Definition 6. Then 2 4p 16p2 E[TC (A)] = tr(A), and Var(TC (A)) ≤ + tr(A) . (41) m m2

Proof Let k = 2p − 1. By Lemma 6, the Complex Recursive TensorSketch matrix admits the independent-layer factorization Πp = M(k) M(k−1) · · · M(1) . Each factor is a Kronecker-wrapped sketching matrix of the form M(i) = Id(i) ⊗ K(i) ⊗ Id(i) , a

where the identity matrices act on the tensor components that remain unchanged and  Ti ∈ Cm×d , i = 1, . . . , p, K(i) = Sℓi ∈ Cm×m2 , i = p + 1, . . . , 2p − 1. ji Here, Ti is a Complex CountSketch matrix given in Appendix A.1 acting at the ith leaf, whereas Sℓjii is a degree-2 Complex TensorSketch matrix given in Appendix A.2 acting at an internal node. In particular, M(k) contains the degree-2 Complex TensorSketch transformation at the root node. The random choices used in the 2p − 1 factors are mutually independent. The proof of Lemma 6 depends only on the recursive Kronecker structure and matrix multiplication. It therefore applies over C after replacing the real sign functions by random hash function drawn from fourth root of unity i.i.d. Define A0 := A, Ai := M(i) Ai−1 (M(i) )∗ ,

i = 1, . . . , k,

and let Xi := tr(Ai ). Since A0 ⪰ 0 and each Ai is obtained from Ai−1 by multiplication with M(i) and its conjugate transpose, positive semidefiniteness is preserved at every step. Therefore, Ai ⪰ 0

for every i = 0, . . . , k.

(42)

Consequently, each Xi is real and nonnegative. Let Fi−1 represent all the random choices made in the first i − 1 layers. Conditional on this information, Ai−1 is fixed and positive semidefinite, while M(i) remains independent and random. The proof of Lemma 12 applies over C after replacing the transpose by the conjugate transpose. Hence, E[Xi | Fi−1 ] = Xi−1 , c 2 Var(Xi | Fi−1 ) ≤ i Xi−1 , m

(43) (44)

where ( ci =

i = 1, . . . , p,

i = p + 1, . . . , 2p − 1.

The constant 1 is the second-moment JL constant for Complex CountSketch, as established in Theorem 18 of Appendix A.1. The constant 3 is the corresponding constant for degree-2 Complex TensorSketch, obtained from Theorem 19 of Appendix A.2 by setting the degree equal to 2. Expectation. Taking expectations in Equation (43) and applying the tower property gives E[Xi ] = E[Xi−1 ].

Iterating over all k layers yields E[Xk ] = E[X0 ] = tr(A). Since Xk = TC (A), it follows that E[TC (A)] = tr(A).

(45)

Variance. Using the conditional second-moment identity together with Equations (43) and (44), we obtain E[Xi2 | Fi−1 ] = Var(Xi | Fi−1 ) + (E[Xi | Fi−1 ])2 c 2 2 + Xi−1 ≤ i Xi−1 m c 2 = 1 + i Xi−1 . m Taking expectations and iterating gives E[Xk2 ] ≤

k Y

ci 2 X0 m

2 ci tr(A) . m

i=1

k Y i=1

(46)

Combining Equations (45) and (46), we obtain 2 Var(TC (A)) = E[Xk2 ] − E[Xk ] " k # Y 2 c ≤ 1 + i − 1 tr(A) . m

(47)

i=1

Using ex − 1 ≤ x + x2 for 0 ≤ x ≤ 1, we obtain 2 (4p − 3)2 4p − 3 Var(TC (A)) ≤ + tr(A) m m2 2 4p 16p2 ≤ + tr(A) . 2 m m This proves the stated unbiasedness and variance bounds.

□

Lemma 16 (Randomness Complexity of Complex Recursive TensorSketch) Let Πp = Qp Tp be the complex Recursive TensorSketch matrix. Then, the total number of random bits p required to construct Π is O p(d + m) log m . Consequently, the estimator TC (A) := tr Πp A(Πp )∗ can be implemented using O p(d + m) log m random bits.

Proof The proof follows along the structure as in Lemma 13. In the complex setting, each random variable can be expressed in the form a + ib, where a and b are real-valued random variables. Thus, compared to the real case, the construction involves at most a constant factor increase in the number of underlying random variables. Therefore, the total number of random bits required remains O p(d + m) log m . □

Theorem 17 [Concentration Analysis of Complex RTS Trace Estimator] Let A ∈ Rd ×d be a fixed nonzero symmetric positive semidefinite matrix and let TC (A) be the trace estimator defined in Definition 6. Then, for every ε > 0, 1 4p 16p2 . (48) Pr[|TC (A) − tr(A)| ≥ ε tr(A)] ≤ 2 + m ε m2 Moreover, for 0 < ε ≤ 1 and δ ∈ (0, 1), TC (A) is an (ε, δ)-approximation whenever m≥

8p . ε2 δ

(49)

Proof From Theorem 15, we have E[TC (A)] = tr(A), 2 4p 16p2 Var(TC (A)) ≤ + tr(A) . 2 m m

(50) (51)

Since A is nonzero and positive semidefinite, tr(A) > 0. Therefore, Chebyshev’s inequality gives Var(TC (A)) 2 ε2 tr(A) 1 4p 16p2 ≤ 2 + . m ε m2

Pr[|TC (A) − tr(A)| ≥ ε tr(A)] ≤

To make the failure probability at most δ, it is sufficient that 1 4p 16p2 + ≤ δ. ε2 m m2

(52)

(53)

Multiplying both sides by m2 ε2 gives ε2 δ m2 − 4pm − 16p2 ≥ 0. Solving this quadratic inequality for m yields p 2p m ≥ 2 1 + 1 + 4ε2 δ . ε δ √ √ Since 0 < ε ≤ 1 and δ ∈ (0, 1), 1 + 1 + 4ε2 δ ≤ 1 + 5 < 4. Hence, p 2p 2 δ < 8p . 1 + 1 + 4ε ε2 δ ε2 δ Therefore, the condition 8p m≥ 2 ε δ is sufficient to make the failure probability at most δ.

(54)

(55)

□

6 Conclusion In this paper, we introduce a trace estimation algorithm for an implicit matrix p p p A ∈ Rd ×d based on Recursive TensorSketch Πp ∈ Rm×d proposed in [14]. Our estimator leverages structured random projections and requires significantly fewer random bits than existing baselines, while maintaining strong theoretical guarantees. We show that the proposed estimator is unbiased and admits a variance bound of 21

2 100p2 O 10p + tr( A ) . We further introduce a complex-valued variant, in which 2 m m the entries of Πp are sampled from complex random variables, leading to improved variance bounds. In contrast to the Kronecker-Hutchinson estimator of [3], our estimator avoids exponential dependence on p in variance bounds of the respective estimators. These properties make the proposed approach well-suited for high-dimensional settings. Our work also suggests several directions for future investigation. Several improved variants of the Hutchinson trace estimator have been proposed that offer additional variance reduction, such as Hutch++ [8, 9] and Krylov-aware trace estimation [5]. It would be interesting to investigate whether our approach can be combined with these techniques to achieve further variance reduction. Extensions of the Hutchinson trace estimator have also been developed for related problems, including the estimation of tr (f (A)) [27], multivariate trace estimation [28], partial trace estimation [29], and trace estimation for tensor data [30]. It would be of interest to explore whether our technique can be integrated into these frameworks to yield randomness-efficient estimators for the corresponding problems. Finally, to obtain an (ε, δ )-approximation, our current analysis relies on Chebyshev’s inequality, resulting in a sample complexity with suboptimal dependence on δ . An important direction for future work is to improve this dependence by leveraging higher-moment analysis of the estimator.

Appendix A Analysis of Complex Sketches This section provides the theoretical analysis of the complex-valued sketching constructions used in this work. We first analyze Complex CountSketch, deriving its unbiasedness, variance, and sketching-time guarantees. We then analyse Complex TensorSketch for degree-p polynomial kernels and establish the corresponding expectation, variance, and computational bounds. Finally, we prove an auxiliary complex AMS moment result used in the analysis of Complex TensorSketch.

A.1 Theoretical Analysis of Complex CountSketch Theorem 18 (Unbiasedness, Variance, and Sketching Time of Complex CountSketch Inner-Product Estimator) Let x, y ∈ Rd , and let C ∈ CD×d be a Complex CountSketch matrix. Define the inner-product estimator by b kC (x, y) = ⟨Cx, Cy⟩C , where ⟨a, b⟩C = a∗ b denotes the Hermitian inner product. Then h i E b kC (x, y) = ⟨x, y⟩, (56) ! d h i X 1 ∥x∥22 ∥y∥22 − x2i yi2 . (57) Var b kC (x, y) = D i=1

Moreover, the sketches Cx and Cy can be computed in O(nnz(x)) and O(nnz(y)) time, respectively.

Proof We first outline the structure of the proof. To prove unbiasedness, we begin by expanding the inner product expression of the sketched vectors obtained from Complex CountSketch. We compute the expectation using the independence of the hash function and complex random hash funtions along with the moment properties E[s(i)s(r)] = 0 for i ̸= r and E[|s(i)|2 ] = 1 due to which all cross terms vanish and we obtain the unbiased estimation of actual inner product. For the variance, we expand the second moment of the estimator and 2 analyze the non-zero terms. Since E[s(i)2 ] = E[s(i) ] = 0 and E[s(i)s(i)] = E[|s(i)|2 ] = 1 for all i ∈ [d] imply that all terms vanish except those corresponding to index configurations with pairwise matchings. Combining these contributions provides a closed-form expression for the second moment, and subtracting the squared mean gives a variance of order 1/D, completing the proof. We now provide the detailed argument. By expanding the estimator, we obtain k̂C (x, y) = ΦC (x)∗ ΦC (y) = ⟨Cx, Cy⟩, =

D X

(58)

(Cx)j (Cy)j ,

j=1

D X d d X X ( s(i)1h(i)=j xi )( s(r)1h(r)=j yr ), r=1

j=1 i=1

D X d X d X

s(i)s(r) 1h(i)=j 1h(r)=j xi yr .

j=1 i=1 r=1

(59)

Computing Expectation: We compute the expected value of Equation (59). D X d h i X i h i h E k̂C (x, y) = xi yr E s(i)s(r) E 1h(i)=j 1h(r)=j , j=1 i,r=1

D X d X

h i h i xi yi E |s(i)|2 E 12h(i)=j + · · ·

j=1 i=1

··· +

D X d X

i h i h xi yr E s(i)s(r) E 1h(i)=j 1h(r)=j .

(60)

j=1 i,r=1 i̸=r

By independence and symmetry of the functions h(.) and s(.), we have E[|s(i)|2 ] = 1 and E[s(i)s(r)] = 0 for i ̸= r. Moreover, since 12h(i)=j = 1{h(i)=j} with h(i) uniform on [D], 1 . D Substituting these identities into (60) vanishes cross term and we get D d h i X 1 X E k̂C (x, y) = xi yi = ⟨x, y⟩. D E[12h(i)=j ] = E[1h(i)=j ] =

j=1

(61)

i=1

This completes the proof of unbiasedness. We next turn to the analysis of the variance of the estimator. Computing Variance: The variance of the complex estimator can be expressed as h i2 h i 2 (62) − E k̂C (x, y) . Var k̂C (x, y) = E k̂C (x, y) To evaluate the first term, we expand it using Equation (59) as follows

k̂C (x, y)

= ⟨Cx, Cy⟩

D X d X

s(i) s(r) 1h(i)=j 1h(r)=j xi yr

j=1 i,r=1

D X

d X

s(i) s(r) s(p) s(q) 1h(i)=j 1h(r)=j 1h(p)=j ′ 1h(q)=j ′ xi yr xp yq . (63)

j,j ′ =1 i,r,p,q=1

Taking expectations with respect to the randomness in h(.) and s(.), we obtain 2 E k̂C (x, y) =

D X

d X

h i h i xi yr xp yq E s(i)s(r)s(p)s(q) E 1h(i)=j 1h(r)=j 1h(p)=j ′ 1h(q)=j ′ .

(64)

j,j ′ =1 i,r,p,q=1

We begin with the case j = j ′ : Since the random variables {s(i)}di=1 are i.i.d. with E[s(i)] = 0, E[|s(i)|2 ] = 1, and E[s(i)2 ] = 0, the fourth-order moment h i E s(i) s(r) s(p) s(q) is nonzero only when each index appears an even number of times. Following terms which are non-zero:

1 . E[14h(i)=j ] = E[1h(i)=j ] = D 1 2 2 E[1h(i)=j 1h(p)=j ] = D2 . E[12h(i)=j 12h(r)=j ] = D12 .

E[|s(i)|4 ] = 1, E[|s(i)|2 |s(p)|2 ] = 1, E[|s(i)|2 |s(r)|2 ] = 1,

(a) i = r = p = q : (b) i = r ̸= p = q : (c) i = p ̸= r = q :

All other cases are zero. Adding the contributions from the above cases, we obtain  d X

x2i yi2 +

i=1



d d X X  1   xi yi xp yp + x2i yr2   . D i,p=1 i̸=p

(65)

i,r=1 i̸=r

We next consider the case j ̸= j ′ : Since the same index cannot hash to two different buckets, all terms vanish except the following case.

(a) i = r ̸= p = q : E[|s(i)|2 |s(p)|2 ] = 1,

E[12h(i)=j 12h(p)=j ′ ] = D12 .

Therefore we have, D−1 X xi yi xp yp . D

(66)

i̸=p

Combining Equation (65) and Equation (66), we get     X d X X X 2 1 D − 1  E k̂C (x, y) = x2i yi2 +  xi yi xp yp + x2i yr2  + xi yi xp yp  , (67) D D i

i̸=p

i̸=r

i̸=p

1 X 2 2 xi yr . = ⟨x, y⟩ + D 2

(68)

i̸=r

Substituting Equation (68) and Equation (61) into Equation (62), we get   i X 1  Var k̂C (x, y) = ⟨x, y⟩2 + x2i yr2  − ⟨x, y⟩2 , D i̸=r X 2 2 1 ∥x∥22 ∥y∥22 − xi yi . = D h

(69) (70)

□ Remark 1 (Sketching time for Complex CountSketch) For a vector x ∈ Rd , the Complex CountSketch sketch Cx can be computed in O(nnz(x)) time. This is because each nonzero entry xi contributes to exactly one bucket h(i) with a single multiplication by the corresponding complex random variable s(i) and a single addition, while zero entries require no computation.

A.2 Theoretical Analysis of Complex TensorSketch Theorem 19 (Unbiasedness and Variance of Complex TensorSketch for Degree-p Polynop p mial Kernel) Let x, y ∈ Rd and let x⊗p , y⊗p ∈ Rd . Let C ∈ CD×d denote a Complex

TensorSketch matrix. Define the degree-p polynomial kernel estimator by b kC (x, y) = Cx⊗p , Cy⊗p C , where ⟨a, b⟩C = a∗ b denotes the Hermitian inner product. Then h i D E E b kC (x, y) = x⊗p , y⊗p = ⟨x, y⟩p , (71)   !p d h i X 1 2 2 2 2 2p  ∥x∥2 ∥y∥2 − xi yi − ⟨x, y⟩  Var b kC (x, y) ≤ D

(72)

i=1

2p − 1 2p ∥x∥2p (73) 2 ∥y∥2 . D Moreover, the sketches Cx⊗p and Cy⊗p can be computed in O(p (nnz(x) + D log D)) and O(p (nnz(y) + D log D)) time, respectively. ≤

Proof We first outline the structure of the proof. Complex TensorSketch is viewed as a Complex CountSketch applied to the p-fold tensor products x⊗p and y⊗p via suitably defined composite hash and complex random functions. We prove Unbiasedness by expanding the sketched inner product and using properties of the expected value of the random function s(.) to eliminate all cross terms. To analyze the variance, we expand the second moment of the estimator and using the independence between the functions (H, S), the second moment reduces to a scaled second-moment expression involving only the random function s(.). This expression is bounded using a complex AMS moment bound, proved later in Lemma 20. Finally, the variance bound is simplified using the Cauchy-Schwarz inequality, resulting in an O(1/D) bound. We now present the detailed proof. We begin by noting that the TensorSketches Cx⊗p , Cy⊗p are the CountSketches of the tensor product X := x⊗p , Y := y⊗p using the two aggregated functions H : [d]p 7→ [D] and S : [d]p → {1, ω, ω 2 , ω 3 } such that:   p X H(i1 , . . . , ip ) =  hj (ij ) mod D, (74) j=1

S(i1 , . . . , ip ) =

p Y

sj (ij ).

(75)

j=1

Also note that H(.) is 2-wise independent [31]. For further proof, we use u, v ∈ [d]p as the indices of vectors X, Y of dimension dp . Then we expand k̂C (x, y) as, X k̂C (x, y) = ⟨CX, CY ⟩ = Xu Yv S(u) S(v) 1[H(u)=H(v)] , (76) u,v∈[d]p

= ⟨X, Y ⟩ +

Xu Yv S(u) S(v) 1[H(u)=H(v)] .

(77)

u̸=v

As we know, E[S(u) S(v)] = 0, ∀ u ̸= v. Then we have h i E k̂C (x, y) = ⟨X, Y ⟩ = ⟨x, y⟩p . (78) h i For the variance, we first compute E |k̂C (x, y)|2 . Let’s first expand the second moment term, |k̂C (x, y)|2 = ⟨Cx⊗p , Cy⊗p ⟩⟨Cx⊗p , Cy⊗p ⟩

(79)



  X X = ⟨X,Y ⟩+ Xu Yv S(u)S(v)1[H(u)=H(v)] ⟨X,Y ⟩+ Xu Yv S(u)S(v)1[H(u)=H(v)], u̸=v

(80)

u̸=v

 =⟨X, Y ⟩2 + ⟨X, Y ⟩ 

Xu Yv S(u)S(v)1[H(u)=H(v)]

u̸=v

 +

Xu Yv S(u)S(v)1[H(u)=H(v)] 

u̸=v

2

 + 

Xu Yv S(u)S(v)1[H(u)=H(v)] 

(81)

u̸=v

h i h i Now, take the expectation of |k̂C (x, y)|2 and we know that E S(u)S(v) = E S(u)S(v) = 0, ∀u ̸= v. Then,   2 h i X   E |⟨Cx⊗p , Cy⊗p ⟩|2 = ⟨X, Y ⟩2 + E  (82) Xu Yv S(u)S(v)1[H(u)=H(v)]   . u̸=v

Using the fact that functions S and H are independent and Lemma 20 (proved below), we can bound the expectation of the second non-diagonal term in the above equation. " ! 2# " X X Xu Yv S(u)S(v)1[H(u)=H(v)] E Xu1 Yv1 Xu2 Yv2 × · · · =E u̸=v

u1 ̸=v1 u2 ̸=v2

# · · · × S(u1 )S(v1 )S(u2 )S(v2 )1[H(u1 )=H(v1 )] 1[H(u2 )=H(v2 )] ,

(83)

h i E Xu1 Yv1 Xu2 Yv2 S(u1 )S(v1 )S(u2 )S(v2 ) · E[1[H(u1 )=H(v1 )] 1[H(u2 )=H(v2 )] ], (84)

u1 ̸=v1 u2 ̸=v2

≤

i 1 X h E Xu1 Yv1 Xu2 Yv2 S(u1 )S(v1 )S(u2 )S(v2 ) , D

(85)

i 1 X h E |Xu1 | |Yv1 | |Xu2 | |Yv2 | S(u1 )S(v1 )S(u2 )S(v2 ) , D

(86)

u1 ̸=v1 u2 ̸=v2

≤

u1 ̸=v1 u2 ̸=v2

 1  = E D

 2  |Xu | |Yv | S(u)S(v)  .

(87)

u̸=v∈[d]p

We bound the above equation using the second-moment bound of Lemma 20. Therefore, we begin by restating the second-moment bound in the proof of Lemma 20,   2 !p d X X   2 2 2 2 2  E |Xu | |Yv | S(u)S(v) xi yi . (88)  = ⟨x, y⟩ + ∥x∥2 ∥y∥2 − i=1

u,v∈[d]p

Now, we expand the term follows,

u,v∈[d]p |Xu | |Yv | S(u)S(v)

from the above equation as

2

 X

|Xu | |Yv | S(u)S(v)



|Xu | |Yv | +

|Xu | |Yv | S(u)S(v)

(89)

u,v∈[d]p u̸=v

u∈[d]p

u,v∈[d]p





 X = |Xu | |Yu | +  u∈[d]p

X u,v∈[d]p u̸=v

 |Xu | |Yv | S(u)S(v)  × ··· 

  X ··· ×  |Xu | |Yu | +  u∈[d]p



X u,v∈[d]p u̸=v

 |Xu | |Yv | S(u)S(v) ,

(90)



 X = |Xu | |Yu | +  u∈[d]p

X u,v∈[d]p u̸=v

 |Xu | |Yv | S(u)S(v)  × ··· 

  X ··· ×  |Xu | |Yu | +  u∈[d]p

X u,v∈[d]p u̸=v

 |Xu | |Yv | S(u)S(v) .

(91)

By further expanding the RHS of the above equation, we get  2 X X  |Xu | |Yv | S(u)S(v) = |Xu1 | |Yu1 | |Xu2 | |Yu2 | + · · · u,v∈[d]p

··· +

u1 ,u2 ∈[d]p

|Xu1 | |Yu1 | |Xu2 | |Yv2 | S(u2 )S(v2 ) + · · ·

u1 ∈[d]p u2 ,v2 ∈[d]p u2 ̸=v2

··· +

|Xu1 | |Yv1 | |Xu2 | |Yu2 | S(u1 )S(v1 ) + · · ·

u1 ,v1 ∈[d]p u2 ∈[d]p u1 ̸=v1

··· +

|Xu1 | |Yv1 | |Xu2 | |Yv2 | S(u1 )S(v1 ) S(u2 )S(v2 ).

(92)

u1 ,v1 ∈[d]p u2 ,v2 ∈[d]p u1 ̸=v1 u2 ̸=v2 We know that u2 ̸= v2 , ∀u2 , v2 ∈ [d]p ,

|Xu1 | |Yu1 | |Xu2 | |Yv2 | E[S(u2 )S(v2 )] = 0,

(93)

u1 ∈[d]p u2 ,v2 ∈[d]p u2 ̸=v2

h i as E S(u2 )S(v2 ) = 0, ∀ u2 ̸= v2 ∈ [d]p . Similarly, for u1 ̸= v1 , ∀u1 , v1 ∈ [d]p , X X |Xu1 | |Yv1 | |Xu2 | |Yu2 | E[S(u1 )S(v1 )] = 0, u1 ,v1 ∈[d]p u2 ∈[d]p u1 ̸=v1

(94)

h i as E S(u1 )S(v1 ) = 0, ∀ u1 ̸= v1 ∈ [d]p . Substituting this into Equation (92) upon computing expectation, we get !2 X X E |Xu | |Yv | S(u)S(v) = |Xu1 | |Yu1 | |Xu2 | |Yu2 | + · · · u,v∈[d]p

u1 ,u2 ∈[d]p

" X

··· + E

|Xu1 | |Yv1 | |Xu2 | |Yv2 | × · · ·

u1 ,v1 ∈[d]p u2 ,v2 ∈[d]p u1 ̸=v1 u2 ̸=v2

# · · · × S(u1 )S(v1 ) S(u2 )S(v2 ) , 2

 = ⟨X, Y ⟩2 + E 

|Xu | |Yv | S(u)S(v)

u̸=v

Now, we conclude that,  2  X X E  |Xu | |Yv | S(u)S(v) = E  u̸=v

2 |Xu | |Yv | S(u)S(v)

− ⟨x, y⟩2p .

(95)

u,v∈[d]p

Now substitute the value of E

u,v∈[d]p |Xu | |Yv | S(u)S(v)

from Equation (88), we

get 2

 E 

|Xu | |Yv | S(u)S(v)

⟨x, y⟩

+ ∥x∥22 ∥y∥22 −

d X

!p x2i yi2

− ⟨x, y⟩2p .

(96)

i=1

u̸=v

Further we substitute this value in Equation (87), we get   2 X   E  Xu Yv S(u)S(v)1[H(u)=H(v)]   u̸=v

  !p d X 1  2 2 2 2 2 2p  ⟨x, y⟩ + ∥x∥2 ∥y∥2 − xi yi − ⟨x, y⟩ . ≤ D

(97)

i=1

Now we can compute the second moment using Equation (82) as follows,   2 h i X   E |k̂C (x, y)|2 = ⟨X, Y ⟩2 + E  Xu Yv S(u)S(v)1[H(u)=H(v)]   ,

(98)

u̸=v

  !p d X 1  ⟨x, y⟩2 + ∥x∥22 ∥y∥22 − ≤ ⟨x, y⟩2p + x2i yi2 − ⟨x, y⟩2p  . D

(99)

i=1

Now, compute variance as follows h h i i2 Var(k̂C (x, y)) = E |k̂C (x, y)|2 − E k̂C (x, y) ,

(100)

≤ ⟨x, y⟩

  !p d X 1  2 2 2 2 2 2p  + ⟨x, y⟩ + ∥x∥2 ∥y∥2 − xi yi − ⟨x, y⟩ − ⟨x, y⟩2p , D i=1

(101)  =

X 2 2 1  ⟨x, y⟩2 + ∥x∥22 ∥y∥22 − xi yi D

 − ⟨x, y⟩2p  .

(102)

i=1

We can upper bound the above equation by using inequality ⟨x, y⟩2 ≤ ∥x∥22 ∥y∥22 , then p 1 2∥x∥22 ∥y∥22 − ∥x∥2p ∥y∥2p , (103) Var(k̂C (x, y)) ≤ 2 2 D p (2 − 1) 2p ∥x∥2p = (104) 2 ∥y∥2 . D □ Remark 2 (Sketching time for Complex TensorSketch) Let x ∈ Rd and let p ≥ 1 be an integer. The Complex TensorSketch of x⊗p with sketch dimension D can be computed in O(p(nnz(x) + D log D)) time. This follows because TensorSketch avoids explicitly forming the tensor x⊗p . Instead, it applies p independent Complex CountSketch to x, each takes O(nnz(x)) time, and combines the resulting p sketches using circular convolution, which is implemented via FFT in O(pD log D) time. Lemma 20 Let x, y ∈ Rd , let p > 1 be an integer, and let s1 , . . . , sp : [d] → {1, ω, ω 2 , ω 3 } be independent functions, each taking values uniformly from the four fourth roots of unity. Define Z =

p Y

Zsj (x) Zsj (y),

(105)

j=1

Where Zsj (x) =

d X

Zsj (y) =

xi sj (i),

d X

yi sj (i).

(106)

i=1

Then, E[Z] = ⟨x, y⟩p , Var[Z] =

(107)

⟨x, y⟩2 + ∥x∥22 ∥y∥22 −

d X

!p x2i yi2

− ⟨x, y⟩2p ,

(108)

i=1 2p ≤ 2p ∥x∥2p 2 ∥y∥2 .

(109)

Proof Following the approach of [32], adapted from [25][Lemma 8], we compute the expectation and variance of Z. First, we consider the expectation. For each j, we note that " d ! !# d h i X X E Zsj (x) Zsj (y) = E xi sj (i) yk sj (k) , (110) i=1

k=1

d X d X

xi yk E[sj (i)sj (k)],

(111)

i=1 k=1

d X

xi yi E[|sj (i)|2 ] +

i=1

xi yk E[sj (i)sj (k)],

(112)

i̸=k

= ⟨x, y⟩,

(113) 2

Where, E[sj (i)sj (k)] = 0, ∀i ̸= k and E[|sj (i)| ] = 1, ∀ i ∈ [d]. Since the functions sj are independent across different j, we have E[Z] =

p Y

E[Zsj (x)Zsj (y)] = ⟨x, y⟩p .

(114)

j=1

Next, to bound the variance, Var(Z) = E[|Z|2 ] − |(E[Z])|2 .

(115)

Because functions is independent across different j, we may write p i h Y E[|Z|2 ] = E | Zsj (x) Zsj (y) |2 .

(116)

j=1

For each j, expanding the square gives h i E | Zsj (x) Zsj (y) |2 " d ! d ! ! d X X X =E xi sj (i) xi sj (i) yk sj (k) i=1

i=1

k=1

d X d X d X d X

d X

!# yk sj (k)

(117)

k=1

h i xi xi′ yk yk′ E sj (i)sj (k)sj (i′ )sj (k′ ) .

(118)

i=1 i′ =1 k=1 k′ =1

Observing that E[sj (i)sj (k)sj (i′ )sj (k′ )] is nonzero only when the indices form pairs (including the possibility that all four are identical), we have  ′ ′   1, if i = k = i = k ,    1, if i = k ̸= i′ = k′ , E[sj (i)sj (k)sj (i′ )sj (k′ )] = (119) ′ ′   1, if i = i ̸= k = k ,    0, otherwise. The contribution from terms with i = k = i′ = k′ is d X x2i yi2 .

(120)

i=1

Terms with i = k ̸= i′ = k′ contribute X

xi yi xi′ yi′ =

i̸=i′

d X i=1

!2 xi yi

−

d X

x2i yi2 = ⟨x, y⟩2 −

i=1

d X

x2i yi2 .

(121)

i=1

Finally, for i = i′ ̸= k = k′ we obtain X

x2i yk2 = ∥x∥22 ∥y∥22 −

d X i=1

i̸=k

x2i yi2 .

(122)

Thus, summing these contributions, we have ! d d h X 2i X 2 2 2 2 2 2 2 xi yi + ⟨x, y⟩ + ∥x∥2 ∥y∥2 − 2 xi yi , E | Zsj (x) Zsj (y) | = i=1

(123)

i=1

= ⟨x, y⟩2 + ∥x∥22 ∥y∥22 −

d X

x2i yi2 .

(124)

i=1

Substituting this bound into Equation (116) yields 2

E[|Z| ] =

⟨x, y⟩

+ ∥x∥22 ∥y∥22 −

d X

!p x2i yi2

(125)

i=1

Which completes the proof since Var(Z) = E[|Z|2 ] − ⟨x, y⟩2p , =

⟨x, y⟩2 + ∥x∥22 ∥y∥22 −

(126) d X

!p x2i yi2

− ⟨x, y⟩2p .

(127)

i=1

Using the Cauchy–Schwarz inequality, ⟨x, y⟩2 ≤ ∥x∥22 ∥y∥22 , and noting that it follows that 2p Var(Z) ≤ 2p ∥x∥2p 2 ∥y∥2 .

2 2 i=1 xi yi ≥ 0,

(128)

□

References [1] Girard, D.: Un algorithme simple et rapide pour la validation croisée généralisée sur des problèmes de grande taille: applications à la restauration d’image. Rapport de recherche 669, IMAG, Grenoble, France (1987) [2] Hutchinson, M.F.: A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines. Communication in Statistics- Simulation and Computation 18, 1059–1076 (1989) https://doi.org/10.1080/03610919008812866 [3] Meyer, R.A., Avron, H.: Hutchinson’s estimator is bad at kroneckertrace-estimation. SIAM Journal on Matrix Analysis and Applications 47(1), 353–387 (2026) https://doi.org/10.1137/24M1720895 https://doi.org/10.1137/24M1720895 [4] Sun, X., Woodruff, D.P., Yang, G., Zhang, J.: Querying a matrix through matrixvector products. ACM Transactions on Algorithms (TALG) 17(4), 1–19 (2021) [5] Chen, T., Hallman, E.: Krylov-aware stochastic trace estimation. SIAM Journal on Matrix Analysis and Applications 44(3), 1218–1244 (2023) [6] Halikias, D., Townsend, A.: Structured matrix recovery from matrix-vector products. Numer. Linear Algebra Appl. 31(1) (2024) https://doi.org/10.1002/NLA. 2531 32

[7] Bakshi, A., Clarkson, K.L., Woodruff, D.P.: Low-rank approximation with 1/ϵ1/3 matrix-vector products. In: Leonardi, S., Gupta, A. (eds.) STOC ’22: 54th Annual ACM SIGACT Symposium on Theory of Computing, pp. 1130– 1143. ACM, Rome, Italy (2022). https://doi.org/10.1145/3519935.3519988 . https://doi.org/10.1145/3519935.3519988 [8] Persson, D., Cortinovis, A., Kressner, D.: Improved variants of the hutch++ algorithm for trace estimation. SIAM J. Matrix Anal. Appl. 43(3), 1162–1185 (2022) https://doi.org/10.1137/21M1447623 [9] Meyer, R.A., Musco, C., Musco, C., Woodruff, D.P.: Hutch++: Optimal stochastic trace estimation. In: Le, H.V., King, V. (eds.) 4th Symposium on Simplicity in Algorithms, SOSA 2021, January 11-12, 2021, pp. 142–155. SIAM, Virtual Conference (2021). https://doi.org/10.1137/1.9781611976496.16 . https://doi.org/10.1137/1.9781611976496.16 [10] Roosta-Khorasani, F., Ascher, U.M.: Improved bounds on sample size for implicit matrix trace estimators. Foundations of Computational Mathematics 15(5), 1187–1212 (2015) https://doi.org/10.1007/s10208-014-9220-1 [11] Wimmer, K., Wu, Y., Zhang, P.: Optimal query complexity for estimating the trace of a matrix. In: Esparza, J., Fraigniaud, P., Husfeldt, T., Koutsoupias, E. (eds.) Automata, Languages, and Programming, pp. 1051–1062. Springer, Berlin, Heidelberg (2014) [12] Jiang, S., Pham, H., Woodruff, D.P., Zhang, Q.R.: Optimal sketching for trace estimation. In: Proceedings of the 35th International Conference on Neural Information Processing Systems. NIPS ’21. Curran Associates Inc., Red Hook, NY, USA (2021) [13] Bujanovic, Z., Kressner, D.: Norm and trace estimation with random rank-one vectors. SIAM Journal on Matrix Analysis and Applications 42(1), 202–223 (2021) https://doi.org/10.1137/20M1331718 https://doi.org/10.1137/20M1331718 [14] Ahle, T.D., Kapralov, M., Knudsen, J.B.T., Pagh, R., Velingker, A., Woodruff, D.P., Zandieh, A.: Oblivious sketching of high-degree polynomial kernels. In: Chawla, S. (ed.) Proceedings of the 2020 ACM-SIAM Symposium on Discrete Algorithms, SODA 2020, Salt Lake City, UT, USA, January 5-8, 2020, pp. 141–160. SIAM, Philadelphia, PA, USA (2020). https://doi.org/10.1137/1. 9781611975994.9 . https://doi.org/10.1137/1.9781611975994.9 [15] Carter, L., Wegman, M.N.: Universal classes of hash functions (extended abstract). In: Hopcroft, J.E., Friedman, E.P., Harrison, M.A. (eds.) Proceedings of the 9th Annual ACM Symposium on Theory of Computing, May 4-6, 1977, Boulder, Colorado, USA, pp. 106–112. ACM, New York, NY, USA (1977).

https://doi.org/10.1145/800105.803400 . https://doi.org/10.1145/800105.803400 [16] Patrascu, M., Thorup, M.: The power of simple tabulation hashing. In: Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing. STOC ’11, pp. 1–10. Association for Computing Machinery, New York, NY, USA (2011). https://doi.org/10.1145/1993636.1993638 . https://doi.org/10.1145/1993636.1993638 [17] Thorup, M.: Simple tabulation, fast expanders, double tabulation, and high independence. CoRR abs/1311.3121 (2013) 1311.3121 [18] Braverman, V., Krishnan, A., Musco, C.: Sublinear time spectral density estimation. In: Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing. STOC 2022, pp. 1144–1157. Association for Computing Machinery, New York, NY, USA (2022). https://doi.org/10.1145/3519935.3520009 . https://doi.org/10.1145/3519935.3520009 [19] Swartworth, W., Woodruff, D.P.: Optimal eigenvalue approximation via sketching. In: Proceedings of the 55th Annual ACM Symposium on Theory of Computing. STOC 2023, pp. 145–155. Association for Computing Machinery, New York, NY, USA (2023). https://doi.org/10.1145/3564246.3585102 . https://doi.org/10.1145/3564246.3585102 [20] Tsourakakis, C.E.: Fast counting of triangles in large real networks without counting: Algorithms and laws. In: Proceedings of the 2008 Eighth IEEE International Conference on Data Mining. ICDM ’08, pp. 608–617. IEEE Computer Society, USA (2008). https://doi.org/10.1109/ICDM.2008.72 . https://doi.org/10.1109/ICDM.2008.72 [21] Avron, H.: Counting triangles in large graphs using randomized matrix trace estimation. In: Workshop on Large-scale Data Mining: Theory and Applications, vol. 10, p. 9 (2010) [22] Han, I., Malioutov, D., Avron, H., Shin, J.: Approximating spectral sums of largescale matrices using stochastic chebyshev approximations. SIAM J. Sci. Comput. 39(4) (2017) https://doi.org/10.1137/16M1078148 [23] Avron, H., Toledo, S.: Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. J. ACM 58(2) (2011) https: //doi.org/10.1145/1944345.1944349 [24] Charikar, M., Chen, K.C., Farach-Colton, M.: Finding frequent items in data streams. Theor. Comput. Sci. 312(1), 3–15 (2004) https://doi.org/10.1016/ S0304-3975(03)00400-6 [25] Pham, N., Pagh, R.: Fast and scalable polynomial kernels via explicit feature maps. In: Dhillon, I.S., Koren, Y., Ghani, R., Senator, T.E., Bradley, P., Parekh,

R., He, J., Grossman, R.L., Uthurusamy, R. (eds.) The 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD 2013, pp. 239–247. ACM, New York, NY, USA (2013). https://doi.org/10.1145/2487575. 2487591 . https://doi.org/10.1145/2487575.2487591 [26] Pham, N., Pagh, R.: Tensor sketch: Fast and scalable polynomial kernel approximation. CoRR abs/2505.08146 (2025) https://doi.org/10.48550/ARXIV.2505. 08146 2505.08146 [27] Ubaru, S., Chen, J., Saad, Y.: Fast estimation of tr(f(a)) via stochastic lanczos quadrature. SIAM J. Matrix Anal. Appl. 38(4), 1075–1099 (2017) https://doi. org/10.1137/16M1104974 [28] Mor-Yosef, L., Ubaru, S., Horesh, L., Avron, H.: Multivariate trace estimation using quantum state space linear algebra. SIAM J. Matrix Anal. Appl. 46(1), 172–209 (2025) https://doi.org/10.1137/24M1654749 [29] Chen, T., Chen, R., Li, K., Nzeuton, S., Pan, Y., Wang, Y.: Faster randomized partial trace estimation. SIAM Journal on Scientific Computing 46(6), 3427–3447 (2024) https://doi.org/10.1137/23M1620399 https://doi.org/10.1137/23M1620399 [30] Verma, B.D., Pratap, R., Kang, K.: Stochastic trace and diagonal estimator for tensors. CoRR abs/2510.22157 (2025) https://doi.org/10.48550/ARXIV.2510. 22157 2510.22157 [31] Pǎtraşcu, M., Thorup, M.: The power of simple tabulation hashing. Journal of the ACM 59(3), 1–50 (2012) https://doi.org/10.1145/2213556.2213560 [32] Braverman, V., Chung, K.-M., Liu, Z., Mitzenmacher, M., Ostrovsky, R.: AMS Without 4-Wise Independence on Product Domains (2010). https://arxiv.org/ abs/0806.4790

Record · ID 965442 · SHA-256 ecb0c5918a9a4fe8
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.