{"id":"978c9ed6-a96e-4b36-b4ca-dfe30671ded1","arxiv_id":"1909.01238","paper_version":1,"verdict":"REJECT","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"high","formal_verification":"none","parameter_count":6,"one_line_summary":"A stochastic quasi-Newton method that learns the Hessian via Gaussian process regression from noisy gradient differences, combined with a stochastic Armijo line search, is demonstrated on nonlinear system identification benchmarks.","lead":"This paper proposes a stochastic quasi-Newton optimizer that models the Hessian with a Gaussian process and adds a stochastic backtracking line search. It targets maximum likelihood identification of nonlinear state space models, where particle filters provide noisy likelihood and gradient estimates.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Equation (36b) has the wrong sign for the GP posterior mean; as written, Algorithm 2's Hessian estimate is inverted, so the central algorithm is not the one claimed.","rationale":"The reader's verdict is REJECT with high correctness risk, and I agree. But the single most load-bearing defect is not, in my reading, the unspecified R: it is the sign error in the GP posterior mean at Eq. (36b). The paper's central claim is that Algorithm 2 learns the Hessian from noisy gradients via a GP and uses it in a quasi-Newton direction. That claim requires Eq. (36b) to be the conditional mean of the joint Gaussian in Eq. (27). It is not: the standard identity has a plus sign, while the paper prints a minus. Because H_k inherits this sign through Eq. (48b), Algorithm 2 as written is not the algorithm described in words. This is an internal consistency failure, not a disagreement with any external consensus, and it is checkable from the equations alone. The reader's chosen weakest assumption (known i.i.d. R) is real but secondary: even with perfectly known R, the GP posterior is wrong as printed. I therefore keep the REJECT/UNCHANGED verdict: the manuscript needs a corrected derivation and a reproducible specification of R and c before the central claim can be accepted.","tokens_in":17609,"tokens_out":7074,"duration_ms":69725,"concrete_test":"Re-derive Eq. (36b) from Eq. (27) by the conditional Gaussian identity; then implement the one-dimensional example of Section 4 (Eq. 38, Fig. 1) twice, once with Eq. (36b) exactly as printed and once with the plus sign. The published minus-sign version should reproduce the dashed blue curve in Fig. 1; if it does not, the paper's Algorithm 2 as written is not the algorithm whose results are reported.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The load-bearing issue is internal, not a matter of consensus. For the joint Gaussian in Eq. (27), the conditional mean of h(x) given y_lk is standard: E[h(x)|y] = mu(x) + K_{x,lk} K_{lk,lk}^{-1}(y_lk - m_lk). Eq. (36b) instead gives phi = mu(x) - K_{x,lk} K^{-1}(y_lk - m_lk). The minus sign is not a notational choice: K_{x,j} is defined in Eq. (35) as a nonnegative kernel integral times Dbar_j^T, so the correction term has the wrong sign. Section 4.1 explicitly reuses Eq. (36) for the simplified update, and Eq. (48b) forms H_k = D phi(x_k). Hence the Hessian estimate used by Algorithm 2 is the negative of the GP regression estimate whenever mu is not already the dominant term. If mu = 0, H_k is exactly the negative of the intended regression estimate. The search direction and the stochastic line-search regularizer are then not the quasi-Newton quantities the paper claims. The experiments in Section 7 cannot be reproduced from the text as written; either Eq. (36b) is a typo and the paper needs a correction, or the algorithm actually run in Section 7 is different from Algorithm 2.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a stochastic quasi-Newton algorithm, called QN-GP, that models the Hessian as a Gaussian process and learns it from noisy gradient differences. A stochastic Armijo line search is introduced and the combined algorithm is applied to maximum likelihood identification of nonlinear state-space models using particle-filter estimates of the likelihood and gradient. The manuscript reports numerical experiments on a linear Gaussian model, a nonlinear benchmark, an interferometry problem, and a MIMO Hammerstein-Wiener system, comparing favorably to an EM-based method in the last example.","tokens_in":17971,"tokens_out":7511,"duration_ms":71425,"significance":"If the technical issues are resolved, the paper would contribute a genuinely novel way to bring second-order information into stochastic optimization via GP regression, with a plausible application to SMC-based system identification where it avoids particle smoothers. The paper also formulates a stochastic Armijo condition and a backtracking line-search that transitions to a classical step-length schedule. However, the current manuscript contains several load-bearing errors and gaps that undermine the claims as written.","major_comments":[{"comment":"The conditional mean of h(x) given y_ell_k is written as phi(x) = mu(x) - K_{x,ell} K_{ell,ell}^{-1} (y_ell - m_ell). The standard Gaussian conditioning formula has a plus sign, not a minus sign. Since K_{x,j} in Eq. (35) is a nonnegative kernel integral times Dbar_j^T, the sign error flips the learned Hessian estimate whenever the prior mean is not the dominant term; in the prior-mean-zero case H_k in Eq. (48b) is the negative of the intended regression estimate. This is not a cosmetic typo: it changes the search direction in Algorithm 2 and makes the algorithm described in the text different from the one that must have been run to produce the experimental results.","section":"Section 4, Eq. (36b)"},{"comment":"The derivation of Theorem 1 contains an inequality sign error. The correct algebra gives alpha (gamma - c(gamma+beta)) >= -E[alpha^2/2 p^T H p], not <=. The paper's displayed inequality is reversed. The stated condition 0 < c < gamma/(gamma+beta) is consistent with the corrected inequality for small alpha, but as printed the proof and the theorem do not match. This makes the stochastic Armijo condition in Section 5 insufficiently supported.","section":"Appendix A, Eq. (A.4)"},{"comment":"No convergence result is proved for the full Algorithm 2. The authors give a descent-direction lemma and a local Armijo condition, and then claim that the line search eventually reduces to alpha_k = xi/k, which 'affords standard convergence results.' But the search direction p_k = -B_k g_k with learned B_k is not the standard stochastic gradient direction, and the paper does not show that the required conditions for stochastic approximation theorems (e.g., on the direction and step sizes) hold for this B_k. A convergence theorem for the proposed algorithm is missing.","section":"Sections 5--6"},{"comment":"The algorithm requires the gradient noise covariance R as an input (Eq. (14)) and uses it both in the GP posterior (Eq. (31)) and in the stochastic Armijo bound (Theorem 1). In the SMC system identification setting the particle-filter gradient estimates have unknown, state-dependent covariance, and the manuscript never states how R was set in the experiments. If R is misspecified, both the Hessian estimate and the line-search acceptance condition are miscalibrated, so the reported results cannot be reproduced or interpreted from the text.","section":"Section 7, experiments"},{"comment":"The numerical comparisons do not isolate the effect of the GP Hessian from that of the line search; there is no baseline that uses the same stochastic line search with a plain gradient direction. In Table 1, the EM results are reported after removing 8/100 runs that 'were trapped in minima far from the global minimum,' while no such removal is done for QNGP; this asymmetric protocol weakens the comparison.","section":"Section 7, Tables 1 and 2"}],"minor_comments":[{"comment":"The relation vech(nabla^2 f) = D h(r_k(tau)) is incorrect; the standard identity is vec(nabla^2 f) = D vech(nabla^2 f). The subsequent equations appear to use the correct identity, but the displayed equation is misleading.","section":"Section 3.2, Eq. (22)"},{"comment":"In the proof of Lemma 1, the term '- E[-v_k]^T B nabla f' contains an erroneous sign; it should be '+ E[v_k]^T B nabla f' or simply omitted since E[v_k]=0.","section":"Proof of Lemma 1"},{"comment":"The comparison methodology in the first experiment adds artificial noise only to the QNGP evaluations; this is a profiling choice, but it should be explicitly framed as such and not as a head-to-head stochastic optimization comparison.","section":"Section 7.1"},{"comment":"The computation of the integrals in Eqs. (31) and (35) is deferred to the authors' technical report [19]. The manuscript would benefit from a brief self-contained description or at least a statement of the main formulas used.","section":"Section 4, after Eq. (37)"},{"comment":"There are multiple typos and grammatical slips, e.g., 'indecies', 'conevient', 'psuedo-code', and an incomplete reference entry [54]. These should be corrected in a revision.","section":"Throughout"}],"recommendation":"major_revision","confidential_remarks":"The paper is not yet ready for publication because the central algorithm as written is affected by a sign error in the GP posterior mean and the appendix contains an algebraic sign error in the main theorem's proof. However, if these are simple typos, a corrected version with a proper convergence analysis and more careful experimental baselines could be a viable contribution. The referee report is based on the manuscript as submitted; I have not assumed that the authors' implementation silently corrected the signs."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things to know. The paper combines a GP model of the Hessian with a stochastic backtracking line search for stochastic optimization, aimed at SMC-based maximum likelihood system identification. That combination is new. Second, the printed mathematics has a load-bearing sign error: Eq. (36b) gives the GP conditional mean with a minus instead of a plus. Since H_k is formed as D times that mean, and the cross-covariance terms are nonnegative kernel integrals, the Hessian estimate is the negative of the intended regression estimate whenever the prior mean is not dominant. As written, Algorithm 2 is not the quasi-Newton algorithm the paper claims, and the Section 7 experiments cannot be reproduced from the text unless this is a typo.\n\nWhat is genuinely good: the integral formulation of the quasi-Newton condition as a line integral of the Hessian, the use of a GP prior over the unique Hessian elements, and the stochastic Armijo line search that eventually decays like 1/k. The application to SMC identification is well motivated, and the comparison with PSEM on the toy and Hammerstein-Wiener problems is a reasonable demonstration. The paper is clearly organized and honestly cites the prior work, including the authors' own CDC 2017 paper.\n\nThe soft spots are real but uneven. The sign error in (36b) is fatal as printed. Eq. (22) also confuses vech and vec — it should be vec(A) = D vech(A), not vech(A) = D h. The Appendix proof of Theorem 1 flips a sign in the reduction; the final bound c < gamma/(gamma+beta) may be correct, but the proof as written does not establish it. The algorithm assumes the gradient-noise covariance R is known, but in the SMC examples it is neither known nor specified, which affects both the Hessian estimate and the line-search acceptance condition. Convergence of the full algorithm is not proven — only the Armijo condition in expectation. Experimental details such as kernel hyperparameters, R, and line-search parameters are missing, so the results are not reproducible from the text.\n\nBottom line: the paper deserves peer review because the core idea is substantive and the main flaw is a specific, checkable sign error. I would not accept this version, but I would engage with a revision that fixes the signs and supplies reproducible experimental settings. The right move is to send it to referees rather than desk-reject.","headline":"A genuinely new GP-Hessian stochastic quasi-Newton idea, but a sign error in Eq. (36b) inverts the Hessian estimate as printed, so the paper needs serious correction before it can be trusted.","tokens_in":18399,"tokens_out":5056,"would_cite":false,"duration_ms":50176,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"Gaussian-process Hessian learning makes quasi-Newton work on noisy gradients.","keywords":["stochastic optimization","quasi-Newton methods","Gaussian process","stochastic line search","system identification","sequential Monte Carlo","particle filter","maximum likelihood estimation"],"falsifier":"Run Algorithm 2 on a nonlinear state-space benchmark while replacing the assumed gradient-noise covariance $R$ in the GP posterior and line search with a deliberately wrong value, say $10R$; if the final parameter estimates and accepted step lengths barely change, the known-$R$ assumption is not load-bearing, whereas a sharp degradation would confirm it is.","tokens_in":17442,"feed_emoji":"📐","tokens_out":11367,"duration_ms":106192,"temperature":0.7,"pith_summary":"The paper proposes a quasi-Newton method for stochastic optimization when only noisy evaluations of the cost and its gradient are available. Instead of assuming the Hessian is constant between iterates, as classical secant equations do, it places a Gaussian-process prior on the Hessian and conditions it on noisy gradient differences, producing a closed-form posterior whose mean supplies curvature for the search direction. It also introduces a stochastic backtracking line search that enforces an Armijo condition in expectation and eventually settles into a $\\alpha_k = \\xi/k$ step-length schedule, so standard stochastic-gradient convergence conditions apply. The combination is demonstrated on maximum-likelihood identification of nonlinear state-space models using particle-filter estimates of likelihood and gradient, where the new algorithm matches a particle-smoother EM baseline while using only forward filtering. If the claims hold, this gives nonlinear system identification a parameter-estimation route that avoids the often difficult particle smoother.","feed_headline":"Noisy gradients can teach a quasi-Newton method the Hessian","feed_subtitle":"A Gaussian-process Hessian plus stochastic line search finds parameters without a particle smoother.","key_machinery":"The load-bearing object is the stochastic quasi-Newton integral, Eq. (17), viewed as a line-integral observation model $y_k = \\bar D_k \\int_0^1 h(r_k(\\tau))\\,d\\tau + w_k$ for the half-vectorized Hessian $h$. Because Gaussian processes are closed under linear operators, conditioning this GP on stacked differences $y_{\\ell_k}$ and $s_{\\ell_k}$ gives the closed-form posterior (36); its mean $\\varphi_{\\ell_k}(x)$ supplies the curvature matrix $H_k$ through the duplication matrix $D$. A second mechanism is the regularized scaling $B_k = (H_k + \\lambda_k I)^{-1}$ with $\\lambda_k = \\epsilon - \\min\\{0, \\eta_k\\}$, where $\\eta_k$ is the minimum eigenvalue of $H_k$, ensuring $B_k \\succ 0$ and hence an expected descent direction. The third mechanism is the stochastic Armijo line search from Algorithm 1, which backtracks while the observed cost fails a sufficient-decrease test and, because backtracking is limited, eventually follows the deterministic schedule $\\alpha_k = \\xi/k$ needed for stochastic-gradient convergence.","core_discovery":"The paper claims that a working quasi-Newton method can be built directly from noisy gradients by treating the Hessian as a Gaussian process and conditioning it on line integrals of gradient differences. Starting from the identity $\\nabla f(x_{k+1}) - \\nabla f(x_k) = [\\int_0^1 \\nabla^2 f(r_k(\\tau))\\,d\\tau] s_k$, the authors derive a stochastic observation model $$y_k = \\bar D_k \\$int_0^{1}$ h(r_k(\\tau))\\,d\\tau + w_k,$$ where $h$ is the half-vectorized Hessian and $\\bar D_k = (s_k^T \\otimes I)D$. With a GP prior on $h$, the closed-form posterior mean $\\varphi_{\\ell_k}(x)$ is obtained from the most recent $p+1$ gradient differences. The search direction is $p_k = -B_k g_k$ with $B_k = (H_k + \\lambda_k I)^{-1}$ and $H_k = D\\varphi_{\\ell_{k-2}}(x_k)$, the regularization chosen so that $B_k \\succ 0$ and therefore the direction is a descent direction in expectation. Step lengths come from a stochastic backtracking line search, Algorithm 1, which enforces an Armijo condition in expectation; Theorem 1 bounds the Armijo constant by $\\bar c = \\gamma/(\\gamma + \\beta)$ with $\\gamma = \\nabla f(x_k)^T B_k \\nabla f(x_k)$ and $\\beta = \\operatorname{Tr}(B_k R)$. On maximum-likelihood identification of nonlinear state-space models with particle-filter likelihood and gradient estimates, the resulting Algorithm 2 matches a particle-smoother EM baseline while using only a forward filter.","pith_inferences":["The paper keeps the gradient-noise covariance $R$ fixed and known, but particle-filter gradient noise is typically state-dependent; estimating $R$ online or modeling it with a second GP would be a direct test of how much the method's calibration depends on that assumption.","The GP construction over the half-vectorized Hessian involves a duplication matrix and squared-exponential kernel, so its cost grows with parameter dimension; extending it to large-scale problems would require sparse or inducing-point approximations.","The stochastic Armijo line search is developed for any search direction, not only the GP-scaled one; composing it with minibatch or adaptive diagonal direction methods would test whether the expectation-form Armijo condition alone improves step-length selection."],"forward_implications":["Algorithm 2 can be run with only a forward particle filter for the likelihood and its gradient, so maximum-likelihood parameter estimation becomes available in settings where particle smoothing is difficult or prohibitive.","Because backtracking is limited to $\\max\\{0, \\tau-k\\}$ steps, the line search eventually produces steps $\\alpha_k = \\xi/k$, placing the iterates under the standard stochastic-approximation convergence conditions.","The GP Hessian posterior uses only the most recent $p+1$ gradient and iterate differences, so per-iteration memory does not grow with the number of optimization steps.","On the reported benchmarks, the new method reaches parameter accuracy comparable to the PSEM baseline; in the nonlinear toy problem, 8 of 100 PSEM runs were trapped in poor local minima while none of the QNGP runs were.","The Hessian is learned continuously during optimization, which in principle lets the method adapt to non-quadratic cost surfaces rather than assuming constant curvature between iterates."],"supporting_citations":[{"why":"Supplies the Gaussian-process regression framework used to model and update the unknown Hessian.","marker":"[41]"},{"why":"Provides the probabilistic interpretation of quasi-Newton updates that motivates representing the Hessian as an object inferred from gradient observations.","marker":"[20]"},{"why":"Sets out the stochastic-gradient convergence conditions and step-length requirements that the line search must eventually satisfy.","marker":"[7]"},{"why":"Supplies the Armijo sufficient-decrease condition that the paper adapts to expectation form.","marker":"[2]"},{"why":"Gives the standard quasi-Newton and line-search theory against which the new algorithm is positioned.","marker":"[38]"},{"why":"Guarantees unbiasedness of particle-filter likelihood estimates, the noisy objective used in the application.","marker":"[11]"},{"why":"Provides particle approximations of the score, producing the noisy gradients that Algorithm 2 consumes.","marker":"[40]"},{"why":"Defines the PSEM particle-smoother EM baseline and the Hammerstein-Wiener benchmark used for comparison.","marker":"[51]"},{"why":"Contributes the numerical method for the squared-exponential covariance integrals appearing in the GP posterior.","marker":"[19]"},{"why":"Supplies the duplication and elimination matrix machinery that handles Hessian symmetry in the observation model.","marker":"[33]"}],"fun_headline_variants":["GP-inferred Hessian makes noisy-gradient quasi-Newton work","Quasi-Newton with GP Hessian learned from noisy gradients","Stochastic line search regularizes GP-Hessian quasi-Newton","No particle smoother needed: GP quasi-Newton matches EM"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the noise in each gradient estimate is independent, Gaussian, and has a known covariance $R$; if the real particle-filter gradient noise is state-dependent or misspecified, both the learned Hessian and the stochastic Armijo test are miscalibrated.","fun_headline_variants_meta":{"raw":{"variants":["GP-inferred Hessian makes noisy-gradient quasi-Newton work","Quasi-Newton with GP Hessian learned from noisy gradients","Stochastic line search regularizes GP-Hessian quasi-Newton","No particle smoother needed: GP quasi-Newton matches EM"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001112,"raw_usage":{"total_tokens":4658,"prompt_tokens":995,"completion_tokens":3663,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":611,"completion_tokens_details":{"reasoning_tokens":3592}},"tokens_in":611,"tokens_out":3663,"duration_ms":26013,"temperature":1.0,"reasoning_tokens":3592,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T05:24:38.346466+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run Algorithm 2 on a nonlinear state-space benchmark while replacing the assumed gradient-noise covariance $R$ in the GP posterior and line search with a deliberately wrong value, say $10R$; if the final parameter estimates and accepted step lengths barely change, the known-$R$ assumption is not load-bearing, whereas a sharp degradation would confirm it is.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the Gaussian-process regression framework used to model and update the unknown Hessian."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the probabilistic interpretation of quasi-Newton updates that motivates representing the Hessian as an object inferred from gradient observations."},{"cited_title":"Bottou, F","cited_arxiv_id":null,"evidence_quote":"Sets out the stochastic-gradient convergence conditions and step-length requirements that the line search must eventually satisfy."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the Armijo sufficient-decrease condition that the paper adapts to expectation form."},{"cited_title":"Nocedal and S","cited_arxiv_id":null,"evidence_quote":"Gives the standard quasi-Newton and line-search theory against which the new algorithm is positioned."},{"cited_title":"Del Moral","cited_arxiv_id":null,"evidence_quote":"Guarantees unbiasedness of particle-filter likelihood estimates, the noisy objective used in the application."},{"cited_title":"Poyiadjis, A","cited_arxiv_id":null,"evidence_quote":"Provides particle approximations of the score, producing the noisy gradients that Algorithm 2 consumes."},{"cited_title":"Wills, T","cited_arxiv_id":null,"evidence_quote":"Defines the PSEM particle-smoother EM baseline and the Hammerstein-Wiener benchmark used for comparison."},{"cited_title":"Evaluating the squared-exponential covariance function in Gaussian processes with integral observations","cited_arxiv_id":"1812.07319","evidence_quote":"Contributes the numerical method for the squared-exponential covariance integrals appearing in the GP posterior."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the duplication and elimination matrix machinery that handles Hessian symmetry in the observation model."}],"review_version":1}