Pith. sign in

REVIEW 3 major objections 6 minor 25 references

Stochastic Diagonal Estimation Based on Matrix Quadratic Form Oracles

T0 review · 3 major / 6 minor · reviewed 2026-08-15 · deepseek-v4-flash

Pith's one-line read Diagonal entries of any matrix can be recovered by scalar quadratic-form queries.

desk verdict First sample complexities for diagonal estimation from quadratic-form oracles; element-wise results are correct, but the norm-wise theorem has a real bug (wrong coefficient and a Markov step that uses the squared norm), so Theorem 2 is false as stated. read the letter →

arxiv 2506.15360 v1 pith:PJU34EJM submitted 2025-06-18 math.NA cs.NA

classification math.NAcs.NA MSC 65C0565F99
keywords diagonalestimationmatrixquadraticformoraclestochasticestimatorunbiasedMonteCarlosamplecomplexityHessianzeroth-orderoptimizationnonsymmetricmatrices
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper establishes that the diagonal of a matrix—symmetric or not—can be estimated using only scalar quadratic-form queries $u^\top A u$, without ever forming $A$ or computing a matrix–vector product. The estimator averages $\tfrac12 (u^\top A u)(u_p^2-1)$ over independent standard Gaussian vectors $u$; each term is an unbiased observation of the diagonal entry $A_{p,p}$. The authors compute the exact variance of this estimator and from it obtain sample complexities: for entrywise accuracy $\varepsilon$ with probability $1-\delta$, roughly $O\big((\operatorname{tr}(A)+A_{p,p})^2+\lVert A+A^\top\rVert_F^2\big)/(\delta\varepsilon^2)$ queries suffice, and for the whole diagonal the required sample size is linear in the dimension $d$. This matters because in zeroth-order optimization only function values are available, and a Hessian quadratic form can be approximated by three function evaluations, so the result opens a route to diagonal preconditioning without gradients or Hessian–vector products.

What carries the argument

The load-bearing object is the random variable $\tfrac12 u^\top A u\,(u_p^2-1)$ with $u\sim\mathcal{N}(0,I_d)$; its expectation is exactly $A_{p,p}$, so a single quadratic-form query becomes one unbiased sample of a diagonal entry. The mechanism that makes the sample complexities explicit is the exact variance computation in Lemma 4: using Gaussian moments up to order eight, the paper reduces the variance to $2(\operatorname{tr}(A)+4A_{p,p})^2+\lVert A+A^\top\rVert_F^2+8\lVert A_{p,:}+A_{:,p}\rVert^2-12A_{p,p}^2$. Chebyshev's inequality converts that variance into the entrywise bound of Theorem 1; summing the per-entry variances and applying Markov's inequality yields the dimension-linear norm-wise bound of Theorem 2; and a median-of-means argument with Hoeffding's inequality turns the $1/\delta$ dependence into $\log(1/\delta)$ in Theorem 3.

What would settle it

Take a fixed nonsymmetric matrix with $A_{p,p}=0$ and large off-diagonal entries, so the variance formula predicts $V=\lVert A+A^\top\rVert_F^2+8\lVert A_{p,:}+A_{:,p}\rVert^2$; draw many Gaussian vectors, compute the sample variance of $\tfrac12 u^\top A u\,(u_p^2-1)$, and check whether it matches $V$ to within Monte Carlo error. If it does not, the Gaussian-moment computation behind Lemma 4 is wrong; if it does, run Algorithm 1 with $N=V/(4\delta\varepsilon^2)$ and check that the empirical failure rate is no larger than $\delta$.

Watch

Extended reading notes

Core claim

The central claim is that for any real matrix $A$ and every index $p$, the identity $\mathbb{E}[\tfrac12 u^\top A u\,(u_p^2-1)] = A_{p,p}$ holds for $u\sim\mathcal{N}(0,I_d)$; averaging $N$ independent copies yields an unbiased estimator of the entire diagonal. The paper goes beyond unbiasedness by computing the estimator's variance in closed form, expressed in terms of $(\operatorname{tr}(A)+4A_{p,p})^2$, $\lVert A+A^\top\rVert_F^2$, and $\lVert A_{p,:}+A_{:,p}\rVert^2$. From that variance it derives an element-wise sample complexity (Theorem 1), a norm-wise sample complexity that grows linearly in $d$ (Theorem 2), and a median-based version whose failure probability scales as $\log(1/\delta)$ rather than $1/\delta$ (Theorem 3). The bounds hold for nonsymmetric matrices and depend only on scalar summaries of $A$; experiments on dense random matrices and a large sparse positive-definite matrix are reported in support of the claimed tightness.

Load-bearing premise

The sample-complexity theorems assume that each query returns the exact value $u^\top A u$; in the motivating Hessian application the query is instead a finite-difference approximation whose error is bounded separately in Proposition 1 and never enters the stated sample sizes.

Editorial extensions

If this is right

  • For any matrix $A$, including nonsymmetric and indefinite ones, each diagonal entry $A_{p,p}$ can be recovered to additive accuracy $\varepsilon$ with probability $1-\delta$ using $V/(4\delta\varepsilon^2)$ quadratic-form queries, where $V$ is the closed-form variance parameter in Eq. (2.4).
  • Estimating the whole diagonal to norm-wise accuracy requires a sample size proportional to $d$, with a proportionality constant controlled by $\operatorname{tr}(A)$ and the off-diagonal part of $A+A^\top$.
  • In zeroth-order optimization, the quadratic-form oracle can be implemented by three function evaluations, so Hessian-diagonal preconditioning becomes available in settings where gradients and Hessian–vector products are too memory-expensive to compute.
  • Because the element-wise bound contains $A_{p,p}^2$ in the denominator when written as a relative error, diagonal entries of small magnitude are intrinsically harder to estimate; the paper's plots confirm this.
  • The median-of-means variant (Algorithm 2) achieves the same entrywise guarantee with total sample complexity proportional to $\log(1/\delta)$ instead of $1/\delta$.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Beyond the paper: the unbiasedness identity does not require symmetry, so the estimator could be applied to implicit nonsymmetric linear operators (for example Jacobians in derivative-free settings) whenever a scalar quadratic query is available.
  • Beyond the paper: the finite-difference perturbation bound in Proposition 1 suggests a combined sample complexity that folds the approximation error $\alpha(d+\log(1/\delta))^3$ into the noiseless bound; making that dependence explicit would close the gap between the theorem and the zeroth-order application.
  • Beyond the paper: replacing the Gaussian law by lighter-tailed zero-mean unit-variance distributions is testable; the proof uses Gaussian moments up to eighth order, and the variance constant $V$ may change with the law.
Share X Bluesky LinkedIn Reddit HN

Signed reviews

No signed human review yet.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 6 minor

Summary. The paper proposes a stochastic estimator for the diagonal of an implicitly given d×d matrix A using only quadratic-form queries u^T A u. The estimator is g = (1/(2N)) Σ_j u_j^T A u_j (u_j^2 - 1_d) for standard Gaussian vectors u_j, and the authors prove an unbiasedness identity and provide element-wise (Theorem 1), norm-wise (Theorem 2), and high-probability median-based (Theorem 3) sample complexities. The paper also gives Gaussian moment computations in Section 3 and numerical experiments on Gaussian, uniform, and a real sparse matrix. The main contribution claimed is that diagonal estimation is possible with a weaker oracle than the matrix-vector products used in prior work, with a norm-wise sample complexity linear in d.

Significance. The element-wise estimation result is a genuine and useful contribution: the estimator is unbiased for arbitrary (not necessarily symmetric) matrices, the variance computation in Theorem 1 checks out, and the comparison with matrix-vector-product estimators is informative. The paper also provides full moment derivations rather than black-box concentration arguments, which is a strength. However, the norm-wise theorem as stated is not proven and is in fact false, because the proof bounds the squared norm while the theorem concerns the norm itself. This is a load-bearing issue since the d-linear norm-wise sample complexity is presented as a central contribution. The paper is potentially publishable after a substantive revision that repairs Theorem 2 and the related experimental validation.

major comments (3)
  1. [Section 3.2] The Markov step in the proof of Theorem 2 is incorrect. The proof writes P{||g−diag(A)||_2 > ε Σ A_ii^2} ≤ E[||g−diag(A)||_2]/(ε Σ A_ii^2) and then replaces E[||g−diag(A)||_2] with (1/(4N)) E[Σ_p (u^T A u u_p^2 − u^T A u − 2A_pp)^2]. The latter expression is E[||g−diag(A)||_2^2], not E[||g−diag(A)||_2]. Consequently, Markov applied to the squared norm gives P(||g−diag(A)||_2^2 > ε Σ A_ii^2) ≤ δ, while the theorem claims P(||g−diag(A)||_2 > ε Σ A_ii^2) ≤ δ. This makes Theorem 2 and Eq. (2.8) false as stated. A concrete check: for A = I_1 and ε = δ = 0.1, the corrected Eq. (2.8) yields N = 1850, but the stated norm event requires N ≥ 18500 (since the variance is 74/(4N) and the threshold is ε^2 S^2 = 0.01 with S = 1). The theorem can be repaired either by changing the guarantee to ||g−diag(A)||_2^2 ≤ ε (Σ A_ii^2)^2 with the same 1/(εδ) rate, or by keeping the stated norm and using a sample complexity proportional to 1/(δ ε^2). This must be fixed because the claimed d-linear norm-wise sample complexity is a central contribution.
  2. [Section 3.1, Lemma 4] The summation of the element-wise variance in Lemma 4 contains an algebra error. Summing Eq. (3.5) over p gives 2Σ_p (tr(A) + 4A_pp)^2 = 2d(tr A)^2 + 16(tr A)Σ_p A_pp + 32Σ_p A_pp^2 = (2d + 16)(tr A)^2 + 32Σ_p A_pp^2, so the total variance is (2d+16)(tr A)^2 + (d+8)||A+A^T||_F^2 + 20Σ_i A_ii^2. The coefficient (4d+16) reported in Eq. (3.6) is therefore wrong and propagates to Eq. (2.8) in Theorem 2 and to Eq. (4.4) in the experiments. Because the error overestimates the variance, it does not by itself invalidate a repaired squared-norm guarantee, but the coefficient must be corrected.
  3. [Section 4.2] The numerical validation of the norm-wise bound is affected by the same issue. The norm-wise relative error defined in Eq. (4.3) is ε = ||g−diag(A)||_2 / Σ A_ii^2, but the theoretical curve in Eq. (4.4) is derived from the squared-norm event P(||g−diag(A)||_2^2 > ε S) ≤ δ, not from the norm event. Therefore the figures do not validate the stated norm-wise sample complexity; they validate a different (weaker) event. After Theorem 2 is repaired, the theoretical curve should be re-derived from the corrected event (likely requiring a square root of the variance bound) and the figures should be regenerated. As written, the experimental 'tightness' claim for the norm-wise bound is not supported.
minor comments (6)
  1. [Section 2.1] The notation ||A^T_{p,:} + A_{:,p}|| in Theorem 1 is inconsistent with the later expression A_{p,:} + A_{:,p} in Lemma 4; please unify the notation, for instance by writing A_{p,:} + (A_{:,p})^T or by explicitly defining the vector expression.
  2. [Algorithm 2] The number of repetitions T = 8 log(1/δ) should be taken as an integer, e.g. T = ceil(8 log(1/δ)); otherwise the Hoeffding argument and the median step are not well-defined for non-integer T.
  3. [Section 4.2] Setting δ = 1 in Eqs. (4.2) and (4.4) makes the probability guarantee vacuous (the confidence level becomes 0). If the theoretical curves are intended as variance-based heuristics rather than probability bounds, this should be stated explicitly; otherwise a concrete value of δ, such as 0.1, should be used.
  4. [Proposition 1] The paper should explicitly state that Theorems 1–3 assume an exact quadratic-form oracle, and that the zeroth-order Hessian application requires choosing α small enough that the O(α (d + log(1/δ))^3) perturbation in Proposition 1 is below the target tolerance. Without such a statement, the guarantees do not directly apply to the motivating Hessian-diagonal setting.
  5. [Appendix C] In the first displayed equation of the proof of Lemma 13, the factor u^2_p appears where u^n_p is clearly intended (the expression is for general n). This typo should be corrected for readability.
  6. [General] The loose bound ||u|| ≤ 2d + 3 log(1/δ) in the proof of Proposition 1 is extremely pessimistic and makes the α-dependence appear much worse than a standard Gaussian norm concentration result would give; the authors may want to replace it with the usual O(sqrt(d + log(1/δ))) bound to give a tighter perturbation estimate.

Circularity Check

0 steps flagged · score 0.0 of 10

No circularity: the estimator is unbiased by direct Gaussian moment computation, and the sample-complexity bounds are derived from computed variances, not from fitted inputs.

full rationale

The derivation chain is not circular. The estimator (2.2) is constructed from a direct identity, Lemma 1, which computes E[u^T A u (u_p^2 - 1)] = 2 A_p,p from standard Gaussian moments; there is no parameter fitted to the target diagonal and no definition of A_p,p embedded in the construction. The element-wise sample complexity (Theorem 1) follows from the exact variance formula (Lemma 4) via Chebyshev, and the high-probability version (Theorem 3) is a standard median argument. The norm-wise claim (Theorem 2) is derived from the displayed variance sum; whatever the correctness of the Markov step (it appears to bound the squared norm rather than the norm), that is a proof error, not circularity, because the bound is not assumed as its own input. Self-citations (Ye 2023; Ye et al. 2025; Zhao et al. 2025) are used for motivation and for the standard Gaussian norm concentration lemma in Proposition 1; they do not supply the main sample-complexity results, and the Gaussian tail bound is an externally checkable fact. No fitted parameters are renamed as predictions, and no uniqueness theorem or ansatz is imported from the authors' prior work.

Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

The main theorems rest only on standard Gaussian moment identities and classical concentration inequalities. The finite-difference connection in Proposition 1 assumes a gamma-Lipschitz Hessian and a Gaussian norm concentration lemma from prior work, but this is motivational rather than used in the main sample complexity proofs. No free parameters are fit to data and no new entities are postulated.

assumptions (4)
  • standard math Standard Gaussian moment identities: E[u^2]=1, E[u^4]=3, E[u^6]=15, E[u^8]=105, and odd moments vanish.
    Used throughout Lemma 3 and Appendix C for variance computations.
  • standard math Classical concentration inequalities: Chebyshev, Markov, and Hoeffding.
    Used to convert variance bounds into sample complexity with probability guarantees in Theorems 1, 2, and 3.
  • domain assumption Exact quadratic-form oracle QA(u) = u^T A u.
    The main theorems are stated for this oracle; Proposition 1 shows only a finite-difference approximation for Hessians, whose error is not incorporated into the sample complexities.
  • domain assumption For Proposition 1, the Hessian is gamma-Lipschitz and the Gaussian norm bound ||u|| <= 2d + 3 log(1/delta) holds with probability 1-delta (Lemma 4 of Ye et al. 2025).
    Needed to bound the finite-difference error in the motivational zeroth-order application; not used in the main matrix-oracle theorems.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Stochastic Diagonal Estimation Based on Matrix Quadratic Form Oracles." pith.science (2026). https://pith.science/paper/PJU34EJM

@misc{pith2026250615360,
  author       = {Pith},
  title        = {Pith review of: Stochastic Diagonal Estimation Based on Matrix Quadratic Form Oracles},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/PJU34EJM}},
  note         = {Machine review of arXiv:2506.15360}
}
abstract

We study the problem of estimating the diagonal of an implicitly given matrix $\Ab$. For such a matrix we have access to an oracle that allows us to evaluate the matrix quadratic form $ \ub^\top \Ab \ub$. Based on this query oracle, we propose a stochastic diagonal estimation method with random variable $\ub$ drawn from the standard Gaussian distribution. We provide the element-wise and norm-wise sample complexities of the proposed method. Our numerical experiments on different types and dimensions matrices demonstrate the effectiveness of our method and validate the tightness of theoretical results.

Figures

Figures reproduced from arXiv: 2506.15360 by the authors.

Figure 1
Figure 1. Experiment results on random Gaussian matrix. [PITH_FULL_IMAGE:figures/full_fig_p013_1.png] view at source ↗
Figure 2
Figure 2. Experiment results on random matrix with each entry uniformly drawing from [0 [PITH_FULL_IMAGE:figures/full_fig_p014_2.png] view at source ↗
Figure 3
Figure 3. Experiment results on “Boeing msc10480” matrix. [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

25 extracted references · 19 canonical work pages

  1. [1]

    Baston, R. A. & Nakatsukasa, Y. (2022). Stochastic diagonal estimation: probabilistic bounds and an improved algorithm. arXiv preprint arXiv:2201.10684

  2. [2]

    Bekas, C., Kokiopoulou, E., & Saad, Y. (2007). An estimator for the diagonal of a matrix. Applied numerical mathematics , 57(11-12), 1214--1229

  3. [3]

    Chen, Y., Zhang, Y., Cao, L., Yuan, K., & Wen, Z. (2025). Enhancing zeroth-order fine-tuning for language models with low-rank structures. In Proceedings of the Thirteenth International Conference on Learning Representations (ICLR 2025)

  4. [4]

    Constantine, P. G. & Diaz, P. (2017). Global sensitivity metrics from active subspaces. Reliability Engineering & System Safety , 162, 1--13

  5. [5]

    & Kressner, D

    Cortinovis, A. & Kressner, D. (2022). On randomized trace estimates for indefinite matrices with an application to determinants. Foundations of Computational Mathematics , (pp.\ 1--29)

  6. [6]

    Davis, T. A. & Hu, Y. (2011). The university of florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS) , 38(1), 1--25

  7. [7]

    Dery, L., Kolawole, S., Kagy, J.-F., Smith, V., Neubig, G., & Talwalkar, A. (2024). Everybody prune now: Structured pruning of llms with only forward passes. arXiv preprint arXiv:2402.05406

  8. [8]

    N., Tropp, J

    Epperly, E. N., Tropp, J. A., & Webber, R. J. (2024). Xtrace: Making the most of every sample in stochastic trace estimation. SIAM Journal on Matrix Analysis and Applications , 45(1), 1--23

Show all 25 references
  1. [9]

    Goedecker, S. (1999). Linear scaling electronic structure methods. Reviews of Modern Physics , 71(4), 1085

  2. [10]

    & Teter, M

    Goedecker, S. & Teter, M. (1995). Tight-binding electronic-structure calculations and tight-binding molecular dynamics with localized orbitals. Physical Review B , 51(15), 9455

  3. [11]

    C., & Saibaba, A

    Hallman, E., Ipsen, I. C., & Saibaba, A. K. (2023). Monte carlo methods for estimating the diagonal of a real symmetric matrix. SIAM Journal on Matrix Analysis and Applications , 44(1), 240--269

  4. [12]

    Hoeffding, W. (1994). Probability inequalities for sums of bounded random variables. The collected works of Wassily Hoeffding , (pp.\ 409--426)

  5. [13]

    & Iooss, B

    Kucherenko, S. & Iooss, B. (2015). Derivative-based global sensitivity measures. In Handbook of uncertainty quantification (pp.\ 1--24). Springer

  6. [14]

    & Gy \"o rgy, A

    Lattimore, T. & Gy \"o rgy, A. (2023). A second-order method for stochastic bandit convex optimisation. In The Thirty Sixth Annual Conference on Learning Theory (pp.\ 2067--2094).: PMLR

  7. [15]

    Liu, H., Li, Z., Hall, D., Liang, P., & Ma, T. (2023). Sophia: A scalable stochastic second-order optimizer for language model pre-training. arXiv preprint arXiv:2305.14342

  8. [16]

    & Tsang, I

    Lyu, Y. & Tsang, I. W. (2021). Black-box optimizer with stochastic implicit natural gradient. In Machine Learning and Knowledge Discovery in Databases. Research Track (pp.\ 217--232).: Springer International Publishing

  9. [17]

    D., Chen, D., & Arora, S

    Malladi, S., Gao, T., Nichani, E., Damian, A., Lee, J. D., Chen, D., & Arora, S. (2023). Fine-tuning language models with just forward passes. Advances in Neural Information Processing Systems , 36, 53038--53075

  10. [18]

    & Polyak, B

    Nesterov, Y. & Polyak, B. T. (2006). Cubic regularization of newton method and its global performance. Mathematical programming , 108(1), 177--205

  11. [19]

    Prashanth, L., Bhatnagar, S., et al. (2025). Gradient-based algorithms for zeroth-order optimization. Foundations and Trends in Optimization , 8(1--3), 1--332

  12. [20]

    Wathen, A. J. (2015). Preconditioning. Acta Numerica , 24, 329--376

  13. [21]

    Yao, Z., Gholami, A., Shen, S., Mustafa, M., Keutzer, K., & Mahoney, M. (2021). Adahessian: An adaptive second order optimizer for machine learning. In proceedings of the AAAI conference on artificial intelligence , volume 35 (pp.\ 10665--10673)

  14. [22]

    Ye, H. (2023). Mirror natural evolution strategies. arXiv preprint arXiv:2308.00469

  15. [23]

    J., & Zhang, T

    Ye, H., Huang, Z., Fang, C., Li, C. J., & Zhang, T. (2025). Hessian-aware zeroth-order optimization. IEEE transactions on pattern analysis and machine intelligence

  16. [24]

    D., Yin, W., Hong, M., et al

    Zhang, Y., Li, P., Hong, J., Li, J., Zhang, Y., Zheng, W., Chen, P.-Y., Lee, J. D., Yin, W., Hong, M., et al. (2024). Revisiting zeroth-order optimization for memory-efficient llm fine-tuning: A benchmark. arXiv preprint arXiv:2402.11592

  17. [25]

    Zhao, Y., Dang, S., Ye, H., Dai, G., Qian, Y., & Tsang, I. W. (2025). Second-order fine-tuning without pain for llms: A hessian informed zeroth-order optimizer. In Proceedings of the Thirteenth International Conference on Learning Representations (ICLR 2025)

Pith tools

Reviewed August 15, 2026 · model on record in the stance chip above.