{"id":"29f020f0-a877-404e-a856-89ecb014f9cf","arxiv_id":"2608.01003","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":7,"one_line_summary":"A Cauchy-integral framework that decouples preconditioner construction from Krylov iteration for f(A)b, with rational and polynomial branches, but the reported gains omit deflation eigensolver costs.","lead":"This paper proposes a framework for speeding up computations of matrix functions applied to a vector, f(A)b, by embedding preconditioners directly into the Cauchy integral that defines f(A). It gives closed-form and quadrature-based variants and tests them on the 2D Laplacian and a lattice QCD operator; the main open question is whether the claimed speedups hold once all overhead costs are counted.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"LQCD convergence claim rests on an unvalidated matrix-free stagnation proxy; no theorem links epsilon_m to ||f_m - f_exact||, and Figure 13 plots proxies rather than true errors.","rationale":"The reader's weakest assumption is correct: the matrix-free stagnation proxy in Section 6.7 is the least secure load-bearing element of the central efficiency claim. The mathematical derivations in Sections 4-6 are internally consistent: the shift-invariance operator equivalence (5.2-5.3), the closed-form extraction (5.7), and the Newton-basis quadrature decoupling (6.8-6.13) all check out. However, the experimental section validates the stopping criterion on only one small Hermitian example and then applies it without further validation to the hard non-Hermitian LQCD case, where the true error is never computed. Because the paper's headline is that these methods compute f(A)b with far fewer iterations than standard Arnoldi, an unverified stopping rule could invalidate the quantitative convergence claims. The concern is not that the proxy is known to fail, but that the paper provides no theorem, no residual bound, and no independent full-space check for the regime where it matters. A concrete full-space comparison at the stopping iteration would settle this. I do not see a more central threat: the derivations are sound, the experiments are honest about their setup, and the abstract's mismatch between f=sign/non-Hermitian and the tested A^{-1/2}/HPD cases is a scope caveat rather than a correctness flaw. The present conditional verdict is appropriate; if the proposed test passes, the efficiency claim would be substantially stronger, and if it fails, the LQCD section would need revision or re-interpretation.","tokens_in":22231,"tokens_out":3699,"duration_ms":46108,"concrete_test":"Re-run the N=32768 LQCD experiment with p=400, d=40, stopping when epsilon_m < 1e-10, and then assemble the full-space vector f_m using the Newton-basis Horner recursion (cost: d-1 extra SpMVs). Compare it against an independent reference computed by uncontracted full Arnoldi with m=2000 (or, if memory is limiting, against a reference solve on a smaller 4096-dimensional sub-block of the same operator where a dense eigendecomposition of A=DD^dagger is feasible). Accept the stopping criterion only if ||f_m - f_reference|| / ||f_reference|| is also below 1e-10 at the same iteration; otherwise the reported LQCD convergence is a proxy artifact and the efficiency comparison needs to be redone using true error.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The decisive efficiency claim for the LQCD benchmark is that POLY-DEF reaches a 1e-10 tolerance in roughly 15 Krylov iterations with p=400 deflated modes (Figure 13b). The convergence indicator used to stop the algorithm is not the true error but the matrix-free stagnation measure epsilon_m (eq. 6.19), which tracks only the relative change of the base polynomial component in the small Schur space. The derivation in Section 6.7 establishes an isometry between the small-space coefficient vector and the full-space component, but only for that component's change, not for the distance from f_m to f_exact. The sole validation is Figure 4, a single 64x64 Laplacian example, where the proxy 'shadows the descent rate' but is not shown to agree quantitatively with the true error near the stopping tolerance. In the LQCD run, no full-space reference error is computed at all; Figure 13's y-axis is 'relative stagnation' for both methods, so the 1e-10 claim is about epsilon_m, not ||f_m - f_exact||. If epsilon_m stagnates before the full-space solution has converged, the reported speedup is illusory; if it lags, the method is stopped late and the SpMV count is understated. The paper's own conclusion also concedes that deflating too many modes can lose to standard Arnoldi, and the SpMV comparison omits the cost of computing the 400 deflated eigenmodes. These issues do not invalidate the framework's derivation, but they do mean the headline LQCD result is currently conditional on a stopping criterion that has not been certified for that problem.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes an integral-based framework for preconditioning f(A)b. Starting from the Cauchy integral representation, the author uses shift-invariance to construct a single Krylov subspace for a preconditioned operator and evaluates f(A)b through a mapped integral. Two branches are developed: (i) rational shift-and-invert preconditioning, yielding the closed-form extraction in Eq. (5.7) with Double MGS reorthogonalization to suppress phantom poles; and (ii) polynomial preconditioning, where Schur decomposition and dual-regime scaling stabilize contour quadrature, and LR-deflation removes critical low modes. Numerical experiments on the 2D Laplacian with f=exp and on a coarse Wilson-Dirac operator from lattice QCD with f=A^{-1/2} are presented. The algebraic derivations in Sections 4-6 are internally consistent, and the shift-invariance operator equivalence is a useful unification of existing techniques.","tokens_in":22609,"tokens_out":6280,"duration_ms":64171,"significance":"If the experimental claims are fully established, the paper would be a valuable contribution: it extends shift-invariance ideas from multi-shift linear systems to the continuous contour-integral setting, provides a closed-form extraction that avoids explicit contour quadrature in the rational case, and proposes a polynomial branch with a matrix-free stopping criterion. The derivations in Sections 4-6 are a strength, as is the availability of Octave code on GitHub. However, the headline LQCD efficiency result is conditional on an unvalidated stagnation proxy and an incomplete cost accounting, so the significance of the numerical claims is currently limited.","major_comments":[{"comment":"The matrix-free stagnation measure epsilon_m is not a bound or reliable estimate of the true error ||f_m - f_exact||. It tracks only the relative change of the base polynomial component in the small Schur space; no theorem connects it to the full-space error. The only validation, Figure 4, is a single small Laplacian case, and the caption contradicts the text on the grid size (32x32 vs 64x64). In Figure 13, the y-axis is 'relative stagnation' for both methods, so the claimed 1e-10 convergence for the LQCD run is not a true-error statement. This is load-bearing for the LQCD speedup claim. The paper should either prove a bound relating epsilon_m to ||f_m - f_exact||, or present reference-error comparisons for the LQCD cases, and clearly label what is plotted.","section":"Section 6.7, Eq. (6.19)"},{"comment":"The SpMV comparison for POLY-DEF with p=400 omits the cost of computing the 400 deflated eigenmodes of the preconditioned operator. The text states these are computed to 1e-10 relative tolerance, but no cost is included. Since m is only about 15 in this run, the eigensolve likely dominates the total cost. The comparison is fair only if the eigensolve cost is amortized over many right-hand sides, as suggested in the introduction. The paper should state this assumption explicitly and provide either the eigensolve cost or an amortized break-even analysis.","section":"Section 7.3.1, Figure 13"},{"comment":"The SpMV-based complexity metric ignores the cost of the contour quadrature itself. For the non-deflated LQCD case (p=0), the paper reports 1e6 quadrature nodes; each node requires solving an m x m triangular system (O(m^2) work). With m=47, this is about 2e9 flops, which is not negligible relative to the reported ~1900 SpMVs. The paper should report wall-clock time or total flops, or at least quantify the quadrature cost, before claiming that the p=0 POLY-DEF run is advantageous in overall efficiency.","section":"Sections 6.4 and 7.3.1"},{"comment":"The theoretical bound in Section 5.2 concerns only eigenvalues: it requires the largest eigenvalue of W to be bounded away from 1.0. The matrix (I_m - \\tilde H_m) whose inverse appears in the closed-form extraction (5.7) is non-normal (a Hessenberg matrix), so eigenvalue bounds do not control ||(I_m - \\tilde H_m)^{-1}||. The claim that the extraction is shielded from singularity is therefore not fully supported by the asymptotic analysis. The paper should either prove a norm bound, or explicitly state that the empirical conditioning history in Figure 7 is the primary evidence for stability.","section":"Section 5.2, Eqs. (5.8)-(5.9), Figure 7"}],"minor_comments":[{"comment":"The text in Section 7.2.1 states that validation of the stopping criterion uses a 64x64 grid (N=4096), but the Figure 4 caption says '32x32 grid'. Correct the inconsistency.","section":"Figure 4"},{"comment":"The legend describes 'STD-ARN (True Error Proxy)' while the y-axis is 'Relative Stagnation'. Clarify what proxy is used for the baseline and whether it is the true error or a residual estimate.","section":"Figure 13"},{"comment":"After Eq. (6.19), the text says the criterion demands 'absolutely zero SpMVs'. This is true per iteration, but the Arnoldi process itself still generates the basis; clarify that the saving is relative to reconstructing f_m at each iteration, not relative to the overall algorithm.","section":"Section 6.7"},{"comment":"The operator D is called non-Hermitian, but the matrix actually used, A=DD^dagger, is Hermitian positive definite. Distinguish D from A to avoid confusion.","section":"Section 7.3"},{"comment":"The legend references 'Eq. 31', which does not match the paper's equation numbering; it should refer to Eq. (6.19).","section":"Figure 4 legend"}],"recommendation":"major_revision","confidential_remarks":"The algebraic core of the paper is solid, and the author acknowledges some limitations (e.g., deflation can lose to standard Arnoldi). The main concern is that the experimental section overstates the reliability of the convergence claims, particularly the LQCD results. The referees found that the stagnation proxy is validated on a single small example and that the SpMV counts omit the eigensolve cost; these are exactly the load-bearing points that need rework. The paper is otherwise a good fit for math.NA."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The paper delivers a real framework, not just tweaks. Two branches: a rational shift-and-invert method with a closed-form extraction (eq. 5.7) and a polynomial contour method with Schur decomposition, dual-regime quadrature, and LR-deflation. The derivations in Sections 4–6 check out; the closed form is genuinely new and neat, and the Newton-basis formulation for high-degree Ritz polynomials is a useful contribution. Code is provided, which helps reproducibility.\n\nThe Laplacian experiments are careful, and the conditioning history in Figure 7 supports the phantom-pole story. The paper is also honest about some limitations: the conclusion explicitly concedes that deflating too many modes can lose to standard Arnoldi.\n\nNow the soft spots, in proportion. First, the abstract advertises f=sign and non-Hermitian Wilson-Dirac, but the experiments actually use f=A^{-1/2} with A=DD†, which is Hermitian positive definite. That mismatch matters for how readers weigh the claims. Second, and more important, the LQCD speedup (about 15 Krylov iterations in Figure 13b) is measured by the matrix-free stagnation proxy epsilon_m, which is validated only on a 64×64 Laplacian. There is no theorem linking epsilon_m to ||f_m - f_exact||, and Figure 13 plots the proxy, not true error. If the proxy stagnates before the full-space solution has converged, the reported speedup is overstated; if it lags, the SpMV count is understated. Third, the SpMV comparison omits the cost of computing the 400 deflated eigenmodes, which is not negligible. Together these make the headline LQCD claim conditional, not established.\n\nNone of this kills the derivation. The algebra stands. But the experiments need revision: report true full-space errors on a moderate LQCD-sized example, include the deflation cost, and align the abstract with what was actually run.\n\nWho is this for? Someone working on matrix functions, Krylov methods, or lattice QCD will get value from the framework and the derivations. It deserves a serious referee, with a request for honest experimental accounting. I would not cite the claimed speedups until the stopping criterion is validated, but I would cite the closed-form extraction and the Newton-basis decoupling.","headline":"A coherent preconditioning framework for f(A)b with real algebraic novelty, but the headline LQCD speedup rests on an unvalidated stopping proxy and an omitted eigensolve cost — worth serious refereeing, not desk rejection.","tokens_in":23145,"tokens_out":1667,"would_cite":true,"duration_ms":18745,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65F60","65F10","65N22"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper shows that f(A)b can be preconditioned through a Cauchy contour integral by building one Krylov basis from a shifted operator, with a rational closed form and a matrix-free polynomial branch.","keywords":["matrix functions","polynomial preconditioning","rational preconditioning","Krylov subspace methods","Cauchy integral representation","shift-invariance","LR-deflation","matrix inverse square root"],"falsifier":"Compute the polynomial branch on a non-Hermitian matrix with explicitly known spectrum and compare its reported $\\varepsilon_m$ with the true relative error $\\|f_m - f_{\\mathrm{exact}}\\|/\\|f_{\\mathrm{exact}}\\|$; if $\\varepsilon_m$ reaches the stopping tolerance while the true error remains orders of magnitude larger, the matrix-free stopping criterion is falsified and the reported LQCD convergence is not certified.","tokens_in":22047,"feed_emoji":"🧮","tokens_out":24484,"duration_ms":202536,"temperature":0.7,"pith_summary":"The paper tries to establish that the expensive action $f(A)b$ can be preconditioned without leaving the Cauchy integral representation: a single shift-invariance condition lets one generate one Krylov basis for a well-conditioned shifted operator and then recover $f(A)b$ through a mapped contour integral. In the rational branch this produces an exact closed-form extraction, $f_m = \\mu(\\mu I_N + A)^{-1} V_m (I_m - \\tilde H_m)^{-1} f(\\hat H_m)(\\beta e_1)$, so no numerical contour quadrature is needed, and double reorthogonalization suppresses the phantom poles that small shifts otherwise create. In the polynomial branch the same identity yields a matrix-free method in which a Schur decomposition reduces each contour integrand to triangular back-substitution, a Newton basis keeps high-degree coefficients stable, and LR-deflation protects the quadrature from near-singular low modes. The payoff, if correct, is that ill-conditioned evaluations—stiff diffusion and the Wilson-Dirac operator from lattice QCD—can be computed with far fewer iterations and a much smaller Krylov basis than unpreconditioned Arnoldi.","feed_headline":"Preconditioned f(A)b with one Krylov basis via shifted contour","feed_subtitle":"A contour-integral shift lets one Krylov basis cut iterations on stiff diffusion and Wilson-Dirac problems.","key_machinery":"The load-bearing mechanism is shift-invariance of Krylov subspaces fused with the Cauchy integral representation of $f(A)b$. For every contour node $z$, the paper asks for a scalar $\\eta_z$ and a shifted preconditioner $M_z^{-1}$ satisfying $\\eta_z I_N - A M_0^{-1} = (z I_N - A) M_z^{-1}$; then $\\mathcal{K}_m(A M_0^{-1}, b) = \\mathcal{K}_m((z I_N - A)M_z^{-1}, b)$, so one orthonormal basis $V_m$ and its projected Hessenberg matrix $\\tilde H_m$ represent every shifted solve. The rational specialization of this identity factors all $z$-dependence into a small matrix $\\hat H_m = \\mu(I_m - \\tilde H_m)^{-1}\\tilde H_m$ and delivers the closed-form extraction $f_m = \\mu(\\mu I_N + A)^{-1} V_m (I_m -","core_discovery":"The central claim is that preconditioning $f(A)b$ reduces to solving, for each contour node $z$, the operator equivalence $\\eta_z I_N - A M_0^{-1} = (z I_N - A) M_z^{-1}$; when it holds, the Krylov subspace generated by $A M_0^{-1}$ is identical to the one generated by every shifted operator $(z I_N - A) M_z^{-1}$, so one basis serves the whole integral. With the rational choice $M_0^{-1} = (\\mu I_N + A)^{-1}$, the shifted systems collapse to the closed form $f_m = \\mu(\\mu I_N + A)^{-1} V_m (I_m - \\tilde H_m)^{-1} f(\\hat H_m)(\\beta e_1)$, where $\\hat H_m = \\mu(I_m - \\tilde H_m)^{-1}\\tilde H_m$; the paper proves a lower bound on $\\mu$ keeps $(I_m - \\tilde H_m)$ nonsingular and uses double mod","pith_inferences":["An extension the paper leaves implicit: with many right-hand sides $b_i$, the Krylov basis, Schur factorization, and deflated subspace are constructed once and reused, so the per-vector cost of the polynomial branch should approach the cost of the polynomial assembly loops alone; a timed benchmark varying the number of right-hand sides would quantify this.","A testable extension of the same framework is to apply the polynomial branch to other analytic functions with branch cuts, such as $f(A)=\\log(A)$ or a different contour for $A^{-1/2}$, on the same twisted-mass operator; the machinery does not depend on the sign function's specific structure.","Since the paper's equivalence condition is sufficient rather than necessary, an approximate version of $M_z^{-1}$ may still give an acceptable single-basis approximation; measuring how the contour quadrature degrades as the equivalence is relaxed would map the boundary of the framework's applicability.","A natural strengthening would be a formal bound linking the cheap stagnation measure $\\varepsilon_m$ to the full-space error; such a bound would certify the LQCD results without relying on the proxy's observed behavior."],"forward_implications":["On the 2D Laplacian on a 512×512 grid with $f(A)=\\exp(-cA)$, the polynomial branch at degree $d=5$ reaches the precision floor in about 88 Krylov iterations (about 527 SpMVs), while standard Arnoldi needs 600 iterations and a memory-heavy basis of dimension 600; increasing stiffness is absorbed by tuning the rational shift $\\mu$.","For the 32768-dimensional twisted-mass Wilson-Dirac operator with $f(A)b = A^{-1/2}b$, a degree-40 harmonic-Ritz polynomial with 400 deflated modes reaches a $10^{-10}$ tolerance in about 15 Krylov iterations (fewer than 600 SpMVs), whereas the unpreconditioned baseline requires a Krylov dimension around 1300.","The rational branch's closed form $f_m = \\mu(\\mu I_N+A)^{-1}V_m(I_m-\\tilde H_m)^{-1}f(\\hat H_m)(\\beta e_1)$ removes contour quadrature entirely, and the lower bound $\\mu > \\lambda_{\\max}\\tau$ is the tunable safety margin protecting the extraction from singularity.","The matrix-free stopping criterion $\\varepsilon_m$ allows the iteration to be halted from small-space coefficient vectors only, with zero additional SpMVs; the paper demonstrates on the Laplacian that it shadows the true error's descent.","Because the polynomial branch never algebraically isolates the function from the squared polynomial, it avoids the sign ambiguity that earlier polynomial preconditioners for $A^{-1/2}$ had to manage."],"supporting_citations":[{"why":"Supplies the operator-equivalence condition and shift-invariance for families of shifted linear systems that the paper lifts into the contour integral.","marker":"[1]"},{"why":"Gives the LR-deflation construction and the exact invariant-subspace component that the polynomial branch adapts to the preconditioned operator.","marker":"[4]"},{"why":"Supplies the Schur decomposition used to reduce the small dense contour integrals to triangular back-substitution.","marker":"[7]"},{"why":"Dynamic Ritz-based polynomial preconditioning combined with exact LR-deflation on ill-conditioned operators; the paper's deflated polynomial method is inspired by it.","marker":"[9]"},{"why":"Establishes that shifting a matrix leaves the Krylov subspace unchanged, the property that lets one basis serve every contour node.","marker":"[10]"},{"why":"Previous polynomial preconditioning for the matrix square root and inverse square root whose algebraic sign ambiguity the contour formulation removes.","marker":"[12]"},{"why":"Documents the phantom-pole phenomenon in finite-precision rational approximations that motivates the double reorthogonalization.","marker":"[15]"},{"why":"Supplies recent evidence that high-degree polynomial preconditioning accelerates convergence for difficult matrices, used to justify the Ritz-based polynomial degrees.","marker":"[19]"},{"why":"Provides the foundational Arnoldi approximation $f(A)b \\approx V_m f(H_m) V_m^H b$ that the integral framework generalizes.","marker":"[22]"},{"why":"Provides the contour-integral evaluation and the parabolic contour choice used for the continuous quadrature.","marker":"[24]"}],"fun_headline_variants":["One Krylov basis for shifted contour: precondition f(A)b efficiently","Shift-invariant f(A)b: one Krylov basis, any contour point","Cut f(A)b iterations on stiff problems with one shifted basis","Precondition f(A)b for many vectors: reuse one Krylov basis","Shifted contour precondition: one Krylov subspace for all nodes"],"cache_read_input_tokens":2816,"weakest_assumption_plain":"The stopping guarantee rests on the assumption that the measured lack of change in the base polynomial component inside the small projected space, $\\varepsilon_m$, faithfully tracks the true error of the full-space solution; this was validated on one 64×64 Laplace example and then used to stop the LQCD runs without further support.","fun_headline_variants_meta":{"raw":{"variants":["One Krylov basis for shifted contour: precondition f(A)b efficiently","Shift-invariant f(A)b: one Krylov basis, any contour point","Cut f(A)b iterations on stiff problems with one shifted basis","Precondition f(A)b for many vectors: reuse one Krylov basis","Shifted contour precondition: one Krylov subspace for all nodes"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000908,"raw_usage":{"total_tokens":3823,"prompt_tokens":906,"completion_tokens":2917,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":650,"completion_tokens_details":{"reasoning_tokens":2825}},"tokens_in":650,"tokens_out":2917,"duration_ms":21006,"temperature":1.0,"reasoning_tokens":2825,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-06T00:37:00.644420+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Compute the polynomial branch on a non-Hermitian matrix with explicitly known spectrum and compare its reported $\\varepsilon_m$ with the true relative error $\\|f_m - f_{\\mathrm{exact}}\\|/\\|f_{\\mathrm{exact}}\\|$; if $\\varepsilon_m$ reaches the stopping tolerance while the true error remains orders of magnitude larger, the matrix-free stopping criterion is falsified and the reported LQCD convergence is not certified.","supporting_citations":[{"cited_title":"Analysis of some","cited_arxiv_id":null,"evidence_quote":"Supplies the operator-equivalence condition and shift-invariance for families of shifted linear systems that the paper lifts into the contour integral."},{"cited_title":"Toward efficient polynomial preconditioning for","cited_arxiv_id":null,"evidence_quote":"Supplies the Schur decomposition used to reduce the small dense contour integrals to triangular back-substitution."},{"cited_title":"An iterative method to compute the sign function of a non-","cited_arxiv_id":null,"evidence_quote":"Dynamic Ritz-based polynomial preconditioning combined with exact LR-deflation on ill-conditioned operators; the paper's deflated polynomial method is inspired by it."},{"cited_title":"PoS , pages=","cited_arxiv_id":null,"evidence_quote":"Establishes that shifting a matrix leaves the Krylov subspace unchanged, the property that lets one basis serve every contour node."},{"cited_title":"Parabolic and hyperbolic contours for computing the","cited_arxiv_id":null,"evidence_quote":"Documents the phantom-pole phenomenon in finite-precision rational approximations that motivates the double reorthogonalization."},{"cited_title":"2019 , school=","cited_arxiv_id":null,"evidence_quote":"Supplies recent evidence that high-degree polynomial preconditioning accelerates convergence for difficult matrices, used to justify the Ritz-based polynomial degrees."},{"cited_title":"Coarsest-level improvements in multigrid for lattice","cited_arxiv_id":null,"evidence_quote":"Provides the foundational Arnoldi approximation $f(A)b \\approx V_m f(H_m) V_m^H b$ that the integral framework generalizes."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Provides the contour-integral evaluation and the parabolic contour choice used for the continuous quadrature."}],"review_version":1}