{"id":"bdc4c39a-535c-46ca-8acc-3e93ce8dd8ee","arxiv_id":"2506.15360","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"A randomized estimator using standard Gaussian vectors recovers any diagonal entry of a matrix from scalar quadratic-form queries, with explicit sample complexity bounds for element-wise and norm-wise accuracy.","lead":"This paper gives a method to estimate the diagonal entries of a matrix using only scalar quadratic form queries, useful when a matrix is too big to store or to multiply by vectors. The method is simple and comes with explicit formulas for how many random queries are needed, and it is tested on random matrices and a 10,480 by 10,480 engineering matrix.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Theorem 2's claimed norm-wise bound is not what its proof establishes: the proof bounds ||g-diag(A)||_2^2 via Markov, not ||g-diag(A)||_2, so the stated sample complexity in Eq. (2.8) can fail for the actual norm event.","rationale":"I read the paper as claiming an unbiased quadratic-form-only diagonal estimator with tight element-wise and d-linear norm-wise sample complexities. The element-wise analysis (Lemma 1, 3, 4 first formula, Theorem 1) checks out; the median trick in Theorem 3 is also sound (the '1/4' in the text should read '3/4', but the Hoeffding step uses 3/4). The load-bearing failure is in the norm-wise guarantee: Theorem 2's proof and experiments compare squared-norm variance to a non-squared norm, so the headline 'norm-wise accuracy linear in d' is not established for the stated metric. This is not an extra-consensus disagreement but an internal inconsistency, and it is directly checkable. The reader's weakest_assumption (unmodeled finite-difference noise in the Hessian application) is a legitimate limitation, but I do not treat it as the most load-bearing concern because the core theorems are stated for an exact quadratic-form oracle and the perturbation can in principle be controlled by choosing α small. My verdict remains CONDITIONAL: the element-wise claims are solid, and the norm-wise claim is repairable by either restating the guarantee for squared norms or changing the sample complexity to include 1/ε^2 and 1/δ^2; the experiments should then compare matched quantities.","tokens_in":24098,"tokens_out":10384,"duration_ms":97731,"concrete_test":"Take A=I_1 (d=1), ε=0.1, δ=0.1. Compute N=1850 from the corrected Eq. (2.8): ((2d+16)tr^2+(d+8)||A+A^T||^2+20∑A_ii^2)/(4δεS)=74/(0.04)=1850. For this N, simulate 10^5 independent runs of Algorithm 1 with standard Gaussian vectors and record whether |g-A_11|>0.1; the observed failure probability should be well above 0.1 (for the stated norm event, the required N is 74/(4δε^2)=18500). This direct Monte Carlo check settles whether Theorem 2's stated guarantee holds as written.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"In the proof of Theorem 2 (Section 3.2), the target event is ||g-diag(A)||_2 > ε ∑ A_ii^2, but the Markov step 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], which is E[||g-diag(A)||_2^2]. The derived sample size Eq. (2.8) therefore guarantees P(||g-diag(A)||_2^2 > ε S) ≤ δ, not P(||g-diag(A)||_2 > ε S) ≤ δ with S=∑ A_ii^2. For A=I_1, ε=δ=0.1, Eq. (2.8) (even with corrected coefficient) gives N=1850, while the stated norm event needs N≥18500 via Var=74/(4N ε^2). The experiments inherit this mismatch: Eq. (4.3) defines relative error as ||g-diag||_2/S, but the theoretical curve Eq. (4.4) is derived from the squared-norm event. Separately, the norm-wise variance sum in Lemma 4 has a coefficient error: summing V_p gives (2d+16)tr(A)^2, not (4d+16); the d+8 coefficient for ||A+A^T||_F^2 is correct once ∑_p||A_p,:+A_:,p||^2 is included. This makes the published Theorem 2 false as stated, though a squared-norm variant with the same d-linear scaling is recoverable.","agreement_with_reader":"disagree"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":24454,"tokens_out":15348,"duration_ms":139927,"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":[{"comment":"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":"Section 3.2"},{"comment":"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":"Section 3.1, Lemma 4"},{"comment":"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.","section":"Section 4.2"}],"minor_comments":[{"comment":"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.","section":"Section 2.1"},{"comment":"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":"Algorithm 2"},{"comment":"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.","section":"Section 4.2"},{"comment":"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.","section":"Proposition 1"},{"comment":"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.","section":"Appendix C"},{"comment":"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.","section":"General"}],"recommendation":"major_revision","confidential_remarks":"The element-wise Theorem 1 and the supporting Gaussian moment computations appear sound and worth publishing. The norm-wise Theorem 2 is not correct as stated; the error is concrete and reproducible, so this is not a matter of taste. I recommend major revision: the authors should repair Theorem 2 (either by changing the event to a squared-norm guarantee or by increasing the sample complexity to 1/(ε^2δ)), correct the coefficient in Lemma 4, and regenerate the corresponding experimental curves. The paper should not be accepted in its current form."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Quick take: this is the first paper I know that gives sample complexities for diagonal estimation from u^T A u queries. The element-wise results are right, and the estimator is nice. But the norm-wise part has a real bug: the variance sum in Lemma 4 is miscomputed, and the proof of Theorem 2 applies Markov to the wrong quantity. Theorem 2 is false as stated, and the norm-wise experiments inherit the error.\n\nWhat's new: the estimator g = 1/(2N) \\sum (u^T A u)(u^2 - 1) is a natural Gaussian-identity extension of Hutchinson, but nobody had worked out non-asymptotic element-wise bounds for the quadratic-form oracle. The unbiasedness for nonsymmetric A is a plus. Lemma 3's moment computations check out; Theorem 1 and the median version (Theorem 3) are correct. The comparison with Bekas/Hallman/Baston-Nakatsukasa is honest; the d factor in sample count is a genuine information difference, not an artifact.\n\nSoft spots: First, Lemma 4's sum over p: \\sum 2(tr(A)+4A_pp)^2 gives (2d+16)tr(A)^2, not (4d+16). The d+8 coefficient on ||A+A^T||_F^2 is fine once you include the row/column terms. Second, Theorem 2's proof bounds P(||g-diag||_2 > \\epsilon S) via E[||g-diag||_2], but then plugs in E[||g-diag||_2^2] (the variance sum) without a square root. That's a category error. The consequence is that the stated sample complexity depends on 1/\\epsilon instead of 1/\\epsilon^2; for A = I_1, \\epsilon=\\delta=0.1, the paper's formula gives N\\approx 1900, but a correct Chebyshev bound needs about 18500. You can rescue a squared-norm guarantee with the same d-linear structure, but it's a different theorem. Third, the experimental validation of tightness uses \\delta=1 (not a probabilistic bound), reports means of 10 runs only, and the norm-wise theoretical curves come from the defective Eq (4.4). So the \"tightness\" claims are overstated. Minor: Proposition 1's finite-difference error is never folded into the oracle model, so the zeroth-order application is motivational rather than covered by Theorems 1-3.\n\nOverall: the element-wise analysis is a solid contribution to the trace/diagonal estimation literature, and the paper is clearly written. The norm-wise theorem needs to be restated and reproved. I'd send it to a referee, but with instructions that Theorem 2 and the norm-wise experiments need substantial revision before acceptance.","headline":"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.","tokens_in":25028,"tokens_out":5002,"would_cite":false,"duration_ms":46639,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65C05","65F99"],"pacs":[],"model":"deepseek-v4-flash","headline":"Diagonal entries of any matrix can be recovered by scalar quadratic-form queries.","keywords":["diagonal estimation","matrix quadratic form oracle","stochastic estimator","unbiased Monte Carlo","sample complexity","Hessian diagonal","zeroth-order optimization","nonsymmetric matrices"],"falsifier":"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$.","tokens_in":23861,"feed_emoji":"🎲","tokens_out":13613,"duration_ms":119430,"temperature":0.7,"pith_summary":"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.","feed_headline":"Quadratic-form probes recover any matrix's diagonal","feed_subtitle":"The estimator needs only scalar queries, so Hessian diagonals become reachable in gradient-free optimization.","key_machinery":"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.","core_discovery":"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.","pith_inferences":["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."],"forward_implications":["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$."],"supporting_citations":[{"why":"introduces Monte Carlo diagonal estimation via matrix–vector products, the baseline setting this paper contrasts with the quadratic-form oracle.","marker":"Bekas et al., 2007"},{"why":"supplies the sharper matrix–vector-product sample complexities (Eq. (2.5) and Corollary 2) that the paper compares against to quantify the information gap.","marker":"Baston & Nakatsukasa, 2022"},{"why":"provides the Monte Carlo diagonal estimation methods whose matrix–vector-product query model the new estimator bypasses.","marker":"Hallman et al., 2023"},{"why":"its cubic-Lipschitz bound is used in Proposition 1 to control the finite-difference error when a quadratic form is approximated by function values.","marker":"Nesterov & Polyak (2006)"},{"why":"supplies the Gaussian norm concentration inequality used in Proposition 1's failure-probability bound.","marker":"Ye et al. (2025)"},{"why":"the probability inequality that underlies the median-of-means argument in Theorem 3.","marker":"Hoeffding (1994)"}],"fun_headline_variants":["Scalar queries unlock every diagonal entry","Gaussian quadratic probes estimate diagonals","Diagonal from quadratic form: unbiased and tight","Probing diagonals with a single scalar oracle","Quadratic oracle recovers full matrix diagonal"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Scalar queries unlock every diagonal entry","Gaussian quadratic probes estimate diagonals","Diagonal from quadratic form: unbiased and tight","Probing diagonals with a single scalar oracle","Quadratic oracle recovers full matrix diagonal"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000258,"raw_usage":{"total_tokens":1539,"prompt_tokens":860,"completion_tokens":679,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":476,"completion_tokens_details":{"reasoning_tokens":614}},"tokens_in":476,"tokens_out":679,"duration_ms":7716,"temperature":1.0,"reasoning_tokens":614,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-15T19:37:16.324992+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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$.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"introduces Monte Carlo diagonal estimation via matrix–vector products, the baseline setting this paper contrasts with the quadratic-form oracle."},{"cited_title":"C., & Saibaba, A","cited_arxiv_id":null,"evidence_quote":"provides the Monte Carlo diagonal estimation methods whose matrix–vector-product query model the new estimator bypasses."},{"cited_title":"& Polyak, B","cited_arxiv_id":null,"evidence_quote":"its cubic-Lipschitz bound is used in Proposition 1 to control the finite-difference error when a quadratic form is approximated by function values."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"the probability inequality that underlies the median-of-means argument in Theorem 3."}],"review_version":2}