{"id":"74ecd0d0-9c49-49c1-8e93-d5a361e981f6","arxiv_id":"2412.12598","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Differentiating a conformally mapped contour integral gives fast, accurate gradients of matrix functions, applied here to free-energy saddle searches for defects in silicon.","lead":"This paper computes derivatives of matrix functions, such as logarithms and square roots, by differentiating a contour-integral formula with automatic differentiation. The method is demonstrated on silicon defect simulations, where it makes free-energy gradients with vibrational entropy practical.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The gradient of the quadrature approximation is only validated against finite differences of itself, leaving unexamined the dependence of the contour nodes and weights on the spectral bounds m, M and the convergence of d(f_N)/du to df/du.","rationale":"I read the central claim as the accuracy and efficiency of AD through the Hale--Higham--Trefethen contour quadrature for matrix functions, with the silicon free-energy gradients as the application. The most load-bearing condition is that the derivative of the numerical approximation converges to the derivative of the true matrix function. The paper provides convergence results from Hale et al. for f_N(X) itself, but the derivative of an approximation need not converge at the same rate or even to the same limit if the approximation family depends on parameters such as the spectral bounds m and M. The experiments compare AD gradients to finite differences of f_N, so they validate only the internal consistency of the implementation, not its agreement with the exact derivative. Given this gap, the reasonable outcome is a conditional acceptance: the method is plausible, reproducible code is provided, and the application is meaningful, but the derivative accuracy claim needs an additional convergence test against a spectral reference. This matches the reader's verdict, so no adjustment is recommended.","tokens_in":20503,"tokens_out":2242,"duration_ms":26058,"concrete_test":"Use a small symmetric positive-definite example, e.g. n = 10 and u a scalar, with X(u) = B(u)^T B(u) + eps I, and compute d/d u Tr log(X(u)) by direct eigendecomposition as a reference. Then compute the same derivative with the contour quadrature and AD for N = 10, 20, 40, 80 quadrature points, in two configurations: (a) m and M fixed at their values at the base point, and (b) m and M recomputed from X(u) and differentiated through. Plot the relative error against the eigendecomposition reference as a function of N. If the error does not decay to the quadrature convergence level, the gradient is differentiating the wrong contour and the central claim fails. Also record the finite-difference error of f_N itself to confirm that comparing against FD of f_N is insufficient.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim requires that differentiating the conformally mapped contour quadrature gives accurate gradients of the true matrix function. The exact derivative in Lemma 2.1 (Section 2.3) is derived for a fixed contour C enclosing the spectrum. However, the numerical method in Eq. (11) treats m and M as inputs to the conformal map and quadrature nodes z(t_j). If an AD implementation treats m and M as constants when differentiating, the resulting gradient omits terms involving dm/du and dM/du. If instead m and M are recomputed from X(u) and differentiated through, the AD trace must differentiate through eigenvalue-extremum estimates or proxies for them, and this is not analyzed anywhere in the paper. No theorem establishes uniform convergence of d(f_N)/du to df/du, nor is there a numerical test that isolates this quantity. The experiments in Section 4.2 and Section 5.2 compare AD gradients to central finite differences of the same quadrature approximation f_N; finite differences of f_N approximate d(f_N)/du, not df/du, so they cannot detect an error in the quadrature derivative. Since the silicon application in Section 5 needs gradients of Trace log^+(F H(u) F), the accuracy of the actual quantity of interest under quadrature truncation is untested. This is the most load-bearing gap in the paper's correctness argument.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a method for computing Jacobians and gradients of analytic matrix functions f(X(u)) by applying automatic differentiation through the conformally mapped contour-integral quadrature of Hale, Higham, and Trefethen. The authors present an exact contour-integral formula for derivatives (Lemma 2.1), derive complexity estimates for dense, banded, and 2D/3D sparse matrices, and implement forward and reverse mode AD in Julia. They validate the approach on toy problems and use it to compute free-energy barriers for vacancy and interstitial migration in silicon, where the quantity of interest is the gradient of Trace log^+(F H(u) F). The central claim is that differentiating the truncated quadrature yields accurate derivatives at a cost that exploits sparsity, e.g., O(ℓ b p n^2) in reverse mode for banded systems.","tokens_in":20904,"tokens_out":6778,"duration_ms":60580,"significance":"If fully established, the method is practically useful: it provides a way to differentiate matrix functions such as the logarithm and p-th roots without dense eigendecompositions, with complexity that exploits structural sparsity. The derivative formula is exact and standard, the complexity tables are plausible, no fitted parameters enter the migration-barrier calculations, and the companion code is provided. However, the numerical evidence currently validates only the derivative of the truncated quadrature, not the convergence of that derivative to the true derivative, and the treatment of spectral bounds under differentiation is unspecified. These gaps are load-bearing for the paper's central accuracy claim.","major_comments":[{"comment":"Figures 4b and 6b compare the AD gradient of the truncated quadrature f_N(X(u)) with centered finite differences of the same f_N. This only demonstrates that differentiation through the quadrature is consistent with the quadrature itself; it does not test whether d f_N/du approximates d f/du. Please add either a theorem bounding ||d/du f(X(u)) - d/du f_N(X(u))|| in terms of N and the spectral condition (for example, by applying the Hale–Higham–Trefethen convergence analysis to the derivative integrand (zI-X)^(-1) X_u (zI-X)^(-1), whose singularities are also confined to the spectrum of X), or a numerical test against an independent high-accuracy reference such as an eigendecomposition on small matrices. This is load-bearing because the silicon application in Section 5 needs the gradient of Trace log^+(F H(u) F), not merely the gradient of its quadrature approximation.","section":"§2.3, Figs. 4b/6b"},{"comment":"The quadrature nodes z(t_j) and weights are defined through the conformal map, which depends on the spectral bounds m and M. The manuscript nowhere states whether m and M are treated as constants when differentiating f_N(X(u)). If they are constant, the derivative omits dm/du and dM/du terms, and the resulting gradient is valid only if the same fixed contour encloses σ(X(u)) for every u in the parameter set considered; no such spectral-inclusion condition is given. If instead m and M are recomputed for each u, then the AD trace must differentiate through the eigenvalue-extremum estimation (or its proxy), which is neither described nor analyzed. Please specify the convention used in the code and experiments, and add the required condition or analysis.","section":"§2.2–2.3, Eqs. (9),(11)"},{"comment":"The VTST rate expression contains an inconsistent sign. The text states E_ξξ(ξ_saddle) < 0 and concludes that HTST overestimates the rate, which requires the variational correction to reduce k relative to k_HTST; however, the printed equation shows exp(β T^2 S_ξ^2/(2E_ξξ)) (equivalently exp(β T^2 T_e)) with the opposite sign, and the following display then writes exp(-β T^2/(2 T_e)). Please correct the signs and define T_e unambiguously, since the direction of the entropy correction is one of the reported application results.","section":"§5.3.2, Eq. (29)"}],"minor_comments":[{"comment":"The derivative du/dt is stated as sn(t); the derivative of sn(t) is cn(t)dn(t). The subsequent formulas use cn(t)dn(t), so this is a typo but should be corrected.","section":"Eq. (6)"},{"comment":"The sentence 'centered finite difference or automatic differentiation [58]' cites [58], which is the ComplexElliptic.jl package; a differentiation reference such as [51] or [36] is intended.","section":"Appendix 8.2"},{"comment":"The toy potential in Eq. (16) contains the term (1/2)|u_i-u_j|^3 while the appendix version in Eq. (54) writes δ|u_i-u_j|^3; the coefficient and the parameter δ should be harmonized.","section":"Eqs. (16) and (54)"},{"comment":"The sentence 'Suppose that f is analytic on a closed set D, where D is an open set containing the spectrum' is self-contradictory; rephrase as analytic on an open set containing a closed neighborhood of the spectrum.","section":"§2.1"},{"comment":"The notation T_e is introduced immediately after its use in the exponential; define it before the formula and ensure it is nonnegative as an 'effective temperature'.","section":"Eq. (29)"}],"recommendation":"major_revision","confidential_remarks":"The paper fits the journal's scope and the method is likely correct, but the missing derivative-convergence analysis and the unspecified treatment of spectral bounds are central to the claimed accuracy. I would not require fundamentally new algorithmic work; a derivative-error theorem or a decisive numerical test, together with clarification of the m/M convention, would be sufficient."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Tim, quick read on arXiv:2412.12598.\n\nThe thing to know: this is not a new mathematical identity. Lemma 2.1 is the standard derivative of the Cauchy integral, and the quadrature is Hale–Higham–Trefethen 2008. What is genuinely new is the engineering: a clean AD-through-the-conformal-map implementation, a sensible complexity accounting (Table 3, with the reverse-mode banded claim of O(ell b p n^2)), and a realistic demonstrator in which gradients of Trace log^+ are used to run NEB on the free-energy surface for vacancy and interstitial migration in silicon. The 0K vacancy barrier comes out 0.515 eV against 0.52 eV from the cited DFT work, which is the right ballpark for a Stillinger-Weber model. Code and data are shipped, so the numerical claims are reproducible. That deserves credit.\n\nNow the soft spots.\n\nThe main one is the one the stress-test flags. All derivative validation in Sections 4.2, 4.3 and 5.2 compares AD against central finite differences of the same truncated quadrature f_N. That checks the implementation differentiates the discrete sum correctly; it does not check that d(f_N)/du converges to df/du. The conformal map depends on m and M, the spectral bounds of X(u). If those are frozen, the gradient silently drops dm/du and dM/du terms; if they are recomputed, AD must differentiate through eigenvalue estimates, which is not analyzed. For the silicon application this is probably a mild error in practice — the spectrum of the preconditioned Hessian should not move much along a migration path — but the paper does not even say which convention is used, and there is no numerical isolation of the effect. This is a genuine gap, not a manufactured one, and it is load-bearing for the correctness argument.\n\nTwo smaller things. The complexity tables ignore the cost of recomputing LU factorizations at each quadrature node, which is fine since they are counted per solve; the selected-inversion aside is just an aside. And the VTST section has a sign/clarity stumble around Eq. (29) — the effective temperature Te ends up negative for a saddle, making the rate correction term look odd until you track the signs through. It is a presentation issue, not a fundamental one.\n\nWho should read it: anyone working on free-energy barriers for defects who needs gradients of matrix functions at scale, and anyone teaching AD through numerical quadrature. It deserves a serious referee. I would ask the authors to (i) state whether m and M are constants or differentiated, (ii) add a derivative-convergence test against eigendecomposition or a very fine reference quadrature on a small example, and (iii) clean up the VTST notation. Without the first two, the accuracy claim is under-supported.\n\nRecommendation: send it out, conditional on those revisions.","headline":"A practical AD-through-contour-integral package with plausible complexity gains and a working silicon application, but derivative accuracy is tested only against finite differences of the same quadrature, leaving the m/M dependence and convergence of d(f_N)/du untested.","tokens_in":21253,"tokens_out":2450,"would_cite":true,"duration_ms":22302,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65F60"],"pacs":[],"model":"deepseek-v4-flash","headline":"Differentiating a conformally mapped contour-integral quadrature yields accurate gradients and Jacobians of matrix functions such as the logarithm and $p$-th roots, at costs that scale with matrix structure, and makes free-energy barrier…","keywords":["matrix functions","automatic differentiation","contour integration","conformal maps","matrix logarithm","vibrational entropy","free energy surface","silicon defects"],"falsifier":"Take a small symmetric matrix $X(u)$ with a known eigendecomposition and spectral endpoints that move strongly with $u$, and compare $\\partial f_N(X(u))/\\partial u$ from forward-mode differentiation of Eq. (11) with the exact derivative obtained from the closed-form derivative of $X^{1/p}$ or $\\log X$; if the difference grows with $dM/du$ in a way not bounded by the function error $\\|f(X)-f_N(X)\\|$, the central derivative-accuracy claim would be refuted.","tokens_in":20308,"feed_emoji":"⚛️","tokens_out":14181,"duration_ms":116038,"temperature":0.7,"pith_summary":"This paper aims to show that derivatives of analytic matrix functions $f(X(u))$—including the matrix logarithm, $p$-th roots, and trace-logarithms—can be computed accurately and cheaply by differentiating a conformally mapped contour-integral quadrature instead of using finite differences or eigendecompositions. The motivating case is the vibrational entropy $S(u)=-\\frac{1}{2}\\operatorname{Trace}\\log^+(F H(u) F)$ of a crystal, whose gradient is required to search for free-energy barriers in defect migration. By applying forward or reverse automatic differentiation to the quadrature, the paper obtains gradients and Jacobians whose accuracy is controlled by the number of quadrature points $\\ell$ and whose cost scales with matrix structure, for example $O(\\ell b p n^2)$ in reverse mode for banded matrices with bandwidth $b$ and $p$ outputs. If the method holds, finite-temperature free-energy-surface calculations become practical at system sizes where direct differentiation of the matrix function is prohibitive, and the paper demonstrates the point on vacancy and interstitial migration in silicon.","feed_headline":"Differentiating contour quadrature gives fast, accurate gradients","feed_subtitle":"Trace-log gradients for silicon free-energy barriers become cheap enough for nudged-elastic-band searches.","key_machinery":"The load-bearing object is the Cauchy integral representation $f(X)=\\frac{1}{2\\pi i}\\oint_C f(z)(zI-X)^{-1}\\,dz$ together with the conformal map that sends the doubly connected slit domain to an annulus through logarithmic, Jacobi-elliptic, and Möbius transformations. After the change of variable $w=\\sqrt{z}$ for functions with a branch cut on $(-\\infty,0]$, the trapezoidal rule yields the quadratures (9) and (11), whose convergence is bounded by $O(e^{-\\pi^2 N/\\log(M/m+3)})$ where $m$ and $M$ are the spectral endpoints of $X$. The paper differentiates these quadrature sums through automatic differentiation, with Lemma 2.1 providing the resolvent-integral form of the parameter derivative as justification; the computational mechanism is simply that the conformal map, the resolvent solves, and the quadrature weights are all composed from differentiable elementary operations, so the same code path that evaluates $f_N(X(u))$ also evaluates its gradient or Jacobian in forward or reverse mode.","core_discovery":"The paper claims that differentiating the quadrature approximation $f_N(X(u))$ reproduces the derivative of $f(X(u))$ to the accuracy needed in applications, and supports that claim with comparisons against second-order centered finite differences. The theoretical starting point is the Cauchy integral formula and its parameter derivative $\\frac{\\partial f(X)}{\\partial u_k} = \\frac{1}{2\\pi i}\\oint_C f(z)(zI-X)^{-1}\\frac{\\partial X}{\\partial u_k}(zI-X)^{-1}\\,dz$; in practice the paper differentiates the explicit trapezoidal conformal-map quadrature (Eqs. (9) and (11)) through automatic differentiation, so forward mode propagates directional derivatives and reverse mode propagates adjoints through the same resolvent solves. This avoids forming the full $n\\times n\\times m$ Jacobian tensor $\\partial X/\\partial u$ when only a scalar output is needed. In the silicon application the method produces a 0K vacancy migration energy of $0.515$ eV, matching the literature value of $0.52$ eV, and temperature-dependent free-energy barriers that decrease as entropic contributions are included.","pith_inferences":["A natural extension would be to prove that $\\partial f_N/\\partial u$ inherits the quadrature's exponential convergence; the paper's numerical validation compares against finite differences, and a closed-form derivative bound would make the accuracy guarantee independent of finite-difference step-size tuning.","Because the quadrature nodes depend on the spectral endpoints $m$ and $M$, and those endpoints are only estimated in practice, an untested source of error is whether $dm/du$ and $dM/du$ contribute to the gradient; differentiating through the eigenvalue estimates and comparing with treating them as constants would settle the size of that term.","The same differentiate-the-quadrature strategy should extend to second derivatives, giving Hessians of the vibrational entropy and hence variational transition-state-theory prefactors without finite differences, even though the paper only demonstrates first derivatives."],"forward_implications":["For a scalar output and a banded matrix, reverse mode costs $O(\\ell b n^2)$, so gradient-based searches on a free-energy surface can scale to much larger supercells than the 64-atom cells shown.","The convergence estimate $O(e^{-\\pi^2 N/\\log(M/m+3)})$ means the number of quadrature points needed for a target accuracy grows only logarithmically with the condition number $M/m$, so ill-conditioned spectra remain tractable.","Because the 0K vacancy migration energy comes out at $0.515$ eV against a literature value of $0.52$ eV, and the barriers shift with temperature, free-energy-surface nudged-elastic-band paths can be computed with the same workflow as ordinary potential-energy paths.","For full Jacobians of $f$, reverse mode costs $O(\\ell n^4)$ in the banded one-dimensional case, so the paper's tables imply that forward mode remains preferable for many-output problems and full Jacobians are only competitive when the number of outputs $p$ is small."],"supporting_citations":[{"why":"Supplies the conformal-map contour-integral quadrature for matrix functions whose evaluation and differentiation the paper extends.","marker":"[31]"},{"why":"Provides the forward-mode dual-number automatic differentiation engine used for pushforward derivative computations.","marker":"[51]"},{"why":"Provides the reverse-mode automatic differentiation engine used to propagate adjoints through the quadrature code.","marker":"[36]"},{"why":"Supplies the chain-rule adjoint interface that makes reverse-mode differentiation of the quadrature possible.","marker":"[68]"},{"why":"Supplies the nudged-elastic-band method used to locate minimum energy paths on the free energy surface in the silicon examples.","marker":"[33]"},{"why":"Provides the literature vacancy migration energy value used to validate the computed 0K barrier.","marker":"[53]"}],"fun_headline_variants":["Fast differentiation of matrix functions via contour AD","Contour integrals give fast gradients for matrix log","Skip the Jacobian tensor for matrix function gradients","Efficient matrix-function gradients via contour AD","Fast trace-log gradients for silicon defects"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that differentiating the finite $N$-point quadrature sum gives the true derivative of the matrix function, including when the spectral endpoints $m$ and $M$ that set the quadrature nodes themselves depend on $u$; the paper relies on this without proving a separate convergence bound for the derivative.","fun_headline_variants_meta":{"raw":{"variants":["Fast differentiation of matrix functions via contour AD","Contour integrals give fast gradients for matrix log","Skip the Jacobian tensor for matrix function gradients","Efficient matrix-function gradients via contour AD","Fast trace-log gradients for silicon defects"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000857,"raw_usage":{"total_tokens":3683,"prompt_tokens":865,"completion_tokens":2818,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":481,"completion_tokens_details":{"reasoning_tokens":2751}},"tokens_in":481,"tokens_out":2818,"duration_ms":21480,"temperature":1.0,"reasoning_tokens":2751,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-11T13:55:37.855391+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Take a small symmetric matrix $X(u)$ with a known eigendecomposition and spectral endpoints that move strongly with $u$, and compare $\\partial f_N(X(u))/\\partial u$ from forward-mode differentiation of Eq. (11) with the exact derivative obtained from the closed-form derivative of $X^{1/p}$ or $\\log X$; if the difference grows with $dM/du$ in a way not bounded by the function error $\\|f(X)-f_N(X)\\|$, the central derivative-accuracy claim would be refuted.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the chain-rule adjoint interface that makes reverse-mode differentiation of the quadrature possible."},{"cited_title":"Formation and migration energies of the vacancy in Si calculated using the HSE06 range-separated hybrid functional","cited_arxiv_id":null,"evidence_quote":"Provides the literature vacancy migration energy value used to validate the computed 0K barrier."}],"review_version":1}