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 →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
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$.
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
- 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.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
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)
- [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.
- [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.
- [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)
- [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.
- [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.
- [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.
- [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.
- [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.
- [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
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
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.
- standard math Classical concentration inequalities: Chebyshev, Markov, and Hoeffding.
- domain assumption Exact quadratic-form oracle QA(u) = u^T A u.
- 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).
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
Reference graph
Works this paper leans on
-
[1]
Baston, R. A. & Nakatsukasa, Y. (2022). Stochastic diagonal estimation: probabilistic bounds and an improved algorithm. arXiv preprint arXiv:2201.10684
arXiv 2022
-
[2]
Bekas, C., Kokiopoulou, E., & Saad, Y. (2007). An estimator for the diagonal of a matrix. Applied numerical mathematics , 57(11-12), 1214--1229
work page 2007
-
[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)
work page 2025
-
[4]
Constantine, P. G. & Diaz, P. (2017). Global sensitivity metrics from active subspaces. Reliability Engineering & System Safety , 162, 1--13
work page 2017
-
[5]
Cortinovis, A. & Kressner, D. (2022). On randomized trace estimates for indefinite matrices with an application to determinants. Foundations of Computational Mathematics , (pp.\ 1--29)
work page 2022
-
[6]
Davis, T. A. & Hu, Y. (2011). The university of florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS) , 38(1), 1--25
work page 2011
-
[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
arXiv 2024
-
[8]
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
work page 2024
Show all 25 references
-
[9]
Goedecker, S. (1999). Linear scaling electronic structure methods. Reviews of Modern Physics , 71(4), 1085
1999
-
[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
1995
-
[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
2023
-
[12]
Hoeffding, W. (1994). Probability inequalities for sums of bounded random variables. The collected works of Wassily Hoeffding , (pp.\ 409--426)
1994
-
[13]
& Iooss, B
Kucherenko, S. & Iooss, B. (2015). Derivative-based global sensitivity measures. In Handbook of uncertainty quantification (pp.\ 1--24). Springer
2015
-
[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
2023
-
[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
2023 arXiv
-
[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
2021
-
[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
2023
-
[18]
& Polyak, B
Nesterov, Y. & Polyak, B. T. (2006). Cubic regularization of newton method and its global performance. Mathematical programming , 108(1), 177--205
2006
-
[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
2025
-
[20]
Wathen, A. J. (2015). Preconditioning. Acta Numerica , 24, 329--376
2015
-
[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)
2021
-
[22]
Ye, H. (2023). Mirror natural evolution strategies. arXiv preprint arXiv:2308.00469
2023 arXiv
-
[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
2025
-
[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
2024 arXiv
-
[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)
2025
Reviewed August 15, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.