{"id":"e199ec72-7282-4571-9213-d4a4049b65bc","arxiv_id":"1908.08845","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":4,"one_line_summary":"A stabilized Chebyshev Langevin sampler, proximal SK-ROCK, accelerates Bayesian imaging by raising effective sample sizes 20-40x over MYULA at equal gradient cost.","lead":"This paper swaps the small Euler step in a Bayesian Langevin sampler for a multi-stage Chebyshev integrator, letting the sampler take much larger stable steps in ill-conditioned imaging problems. In deconvolution, hyperspectral unmixing, and tomography tests, the new sampler reaches effective sample sizes 20 to 40 times larger at the same gradient budget.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The printed O(√κ) complexity proof is internally inconsistent: the C in §3.1.1 conflicts with (3.3), (3.7), (3.8), and Appendix B's squared-triangle inequality is invalid; a direct Gaussian gradient-evaluation count is required.","rationale":"I agree with the reader that the complexity statement in §3.1.1 is the fragile point and that the displayed C formula is not derived; that is exactly where the argument is least secure. However, my reading differs on the conclusion. If one uses the C that actually follows from (3.3), (3.7) and (3.8), the per-outer-iteration contraction is roughly constant in κ, so the total number of gradient evaluations is O(s log(1/ε))=O(√κ log(1/ε)) rather than O(κ); the inconsistency in the printed formula may therefore be corrigible. The additional flaw I found, the invalid squared-triangle inequality in Appendix B, is a genuine proof error but a corrected bound with Young's inequality would preserve the κ-scaling because only constant factors change. The empirical evidence (equal-gradient-budget comparisons, ESS speed-ups of 21.8 and 37.9, P-MALA bias checks, long runs) is independent support and should be credited; the authors also candidly state the limitations: Gaussian-only analysis, no general complexity theory, and no acceleration when asymptotic bias dominates. These self-reported limitations were weighed. For these reasons the concern is substantial but not fatal: the right action is to keep the manuscript conditional, requiring a corrected derivation of C and a direct Gaussian gradient-evaluation count, plus ideally code release and repeated-run error bars. I therefore leave the reader's verdict unchanged rather than moving to reject.","tokens_in":25309,"tokens_out":22856,"duration_ms":243189,"concrete_test":"Use the exact expression (3.4) with R1,R2 from (3.3). For κ=10²,10⁴,10⁶,10⁸, set s=ceil(√(η/2)(κ−1)) and δ=(ω0−1)/(ℓω1) as in (3.7)–(3.8), compute the minimal number of gradient evaluations N (outer iterations times s) to reach W2(π,Q_N)<ε with ε=0.1·W2(π,Q_0) and a fixed deterministic X0≠0, and plot log N against log κ. The slope distinguishes O(√κ) from O(κ). Also tabulate C=max_i R1(z_i)² and compare with ((√κ−1)/(√κ+1))² and with 1/T_s(ω0)² for the same κ values.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The load-bearing element of the paper's theoretical claim is the gradient-evaluation count in §3.1.1, and it is not established by the text. First, the displayed contraction constant C≈((√κ−1)/(√κ+1))² cannot be obtained from the algorithm's own ingredients: with δ=(ω0−1)/(ℓω1) and z_i=−δ/σ_i², the slowest coordinate has ω0+ω1 z_i=1, so R1(z_i)=T_s(1)/T_s(ω0)=1/T_s(ω0). Thus C=max_i R1(z_i)²=1/T_s(ω0)²≈cosh⁻²(√(2η))≈0.907 for large κ, independent of κ. The paper cites [20] for the other formula but does not derive it in Appendices A–B. Second, the proof of (3.5) in Appendix B starts from W2(P,Q_{n+1})² ≤ W2(P,π̃)² + W2(π̃,Q_{n+1})², which is not a valid use of the triangle inequality for squared 2-Wasserstein distances: the positive cross term is dropped. A correct bound would have extra constants (e.g., via Young's inequality). So no correct non-asymptotic complexity count follows from the written proof. This does not automatically falsify the acceleration claim: with the corrected C above, the per-gradient contraction is exp(log C/s)=1−O(κ^{-1/2}), so O(√κ) total gradient evaluations could still hold; a direct numerical evaluation of the exact Gaussian formula is needed to decide.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a proximal Markov chain Monte Carlo method, SK-ROCK, that replaces the Euler-Maruyama discretization in the Moreau-Yosida unadjusted Langevin algorithm (MYULA) with an orthogonal Runge-Kutta-Chebyshev stochastic integrator. For Gaussian targets the authors derive an explicit Wasserstein-distance bound and claim that, under a contraction constant C approximately equal to ((sqrt(kappa)-1)/(sqrt(kappa)+1))^2, the number of gradient evaluations needed to reach a prescribed accuracy is O(sqrt(kappa)) rather than O(kappa). The method is then tested on one-dimensional Laplace and uniform targets, image deconvolution, hyperspectral unmixing, and tomographic reconstruction, with equal numbers of gradient/proximal evaluations for MYULA and SK-ROCK. The empirical sections report substantially faster autocorrelation decay, larger effective sample sizes, and lower MMSE estimation error for SK-ROCK at equal computational budget.","tokens_in":25648,"tokens_out":3116,"duration_ms":31018,"significance":"If the acceleration claim holds, the paper is a useful and potentially influential contribution to Bayesian imaging, since it offers a proximal MCMC algorithm that behaves like an accelerated optimization method on ill-conditioned problems. The experimental methodology is careful: comparisons are made at equal gradient-evaluation budgets, long runs are used for autocorrelation and ESS estimates, and a Metropolized P-MALA baseline is included to track asymptotic bias. The reported speedups in ESS for the slowest mixing component (21.77x and 37.93x in the cameraman and hyperspectral experiments) are striking and worth taking seriously. The Gaussian Wasserstein computation in Appendix A is self-contained and useful. However, the central theoretical complexity statement is not established by the written proof, and the paper's own formulas point to an internal inconsistency in the O(sqrt(kappa)) claim.","major_comments":[{"comment":"The claim that SK-ROCK requires O(sqrt(kappa)) gradient evaluations is not supported by the displayed formulas. With delta = (omega_0 - 1)/(ell s omega_1), the slowest coordinate satisfies omega_0 + omega_1 z = 1, so R_1(z) = T_s(1)/T_s(omega_0) = 1/T_s(omega_0), giving C = max_i R_1(z_i)^2 approximately equal to cosh^{-2}(sqrt(2 eta)) = 0.907 for large kappa, which is independent of kappa and contradicts the quoted C approximately ((sqrt(kappa)-1)/(sqrt(kappa)+1))^2. Moreover, even under the quoted formula, each outer iteration costs s approximately sqrt(kappa) gradient evaluations, so O(log(epsilon^{-1}) sqrt(kappa)) outer iterations would translate to O(kappa) gradient evaluations, not O(sqrt(kappa)). Please provide a direct gradient-evaluation count from the exact Gaussian formula in Proposition 3.1, or state the complexity claim as an empirical observation.","section":"Section 3.1.1, Eqs. (3.7)-(3.8)"},{"comment":"The proof of the key bound (3.5) uses W_2(P; Q_{n+1})^2 <= W_2(P; Qtilde)^2 + W_2(Qtilde; Q_{n+1})^2, which is not a valid use of the triangle inequality for squared 2-Wasserstein distances because the positive cross term is omitted. A correct derivation via Young's inequality would introduce an additional factor (such as 1+rho or 1+1/rho), so the displayed constant C as the sole multiplicative factor is not justified. This invalidates the non-asymptotic convergence proof as written, and therefore the complexity conclusion that rests on it.","section":"Appendix B, inequality (B.1)"},{"comment":"The paper correctly notes that no general complexity results for Runge-Kutta-Chebyshev methods exist, and the Gaussian analysis is restricted to diagonal covariance with deterministic initial condition. The conclusion that the O(sqrt(kappa)) behavior carries over to the non-Gaussian, nonsmooth imaging models in Section 4 is therefore heuristic. Please state explicitly which parts of the theoretical claim are proved and which are conjectural, and consider adding a direct numerical check of the exact Gaussian Wasserstein formula for a range of kappa to determine whether the acceleration is genuine before the O(sqrt(kappa)) statement is presented as a theorem-level result.","section":"Section 3.1.1, first paragraph and Proposition 3.1"}],"minor_comments":[{"comment":"The reference to 'Appedix A' in the text is a typo for 'Appendix A'.","section":"Section 3.1.1 and Appendix A"},{"comment":"The notation is inconsistent: the invariant measure is written as tilde pi in the main text and as tilde Q in Appendix B; please unify the notation.","section":"Eq. (3.5) and Appendix B"},{"comment":"The caption of Figure 3 mentions W_2(P; Q_0), but the precise value of Q_0 and the exact numerical procedure used to produce the curves are not given; please add enough detail to reproduce the figure.","section":"Figure 3 and Section 3.1.1"},{"comment":"For the tomography experiment the speed-up for the fastest component is 0.81, below 1, which is not discussed; a sentence explaining why SK-ROCK can be slower on the fastest component would improve the presentation.","section":"Table 5"}],"recommendation":"major_revision","confidential_remarks":"The main concern is the gap between the paper's central theoretical claim and the proof provided; this is a load-bearing issue because the O(sqrt(kappa)) complexity statement is presented as the paper's main theoretical contribution. I do not see evidence of misconduct; the empirical work appears carefully conducted. One small observation for the editor: the contraction formula is attributed to reference [20], which shares an author with the present paper; this is not a problem by itself, but the derivation should be made self-contained or clearly credited at each step. The manuscript fits the journal's scope and could become publishable after a major revision that either proves or carefully qualifies the complexity claim, ideally with a direct numerical evaluation of the exact Gaussian formula."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a solid engineering contribution with a shaky theory section. The SK-ROCK + MYULA combination is genuinely new and the experiments are convincing; the complexity claim is not proven as written.\n\nWhat's new: replacing Euler-Maruyama in MYULA with SK-ROCK, Chebyshev stages, plus the Gaussian Wasserstein analysis and three imaging studies. The empirical comparisons are fair—same gradient/proximal budget, thinning for autocorrelation, P-MALA benchmarks. The speed-ups in ESS (21.7× cameraman, 37.9× hyperspectral, 20.2× tomography) are large and plausible. They also state limitations honestly: no general complexity theory, Gaussian-only analysis, and acceleration disappears when the bias term dominates.\n\nThe weak spot is the theoretical section. The paper claims O(√κ) gradient evaluations for SK-ROCK. The displayed contraction constant C ≈ ((√κ−1)/(√κ+1))² is attributed to [20] and not derived in Appendices A–B. If you instead use the algorithm's own R₁(zᵢ), with δ = (ω₀−1)/(ℓω₁), the slowest coordinate gives C = 1/T_s(ω₀)², which is roughly cosh⁻²(√(2η)) ≈ 0.907, independent of κ. That does not give O(√κ) from the stated formulas; with s ≈ √κ stages, the total gradient count would be different. Also, Appendix B applies the triangle inequality to squared Wasserstein distances, dropping the cross term. That is not a valid bound. These issues don't automatically kill the acceleration mechanism—per-gradient contraction exp(log C/s) could still yield O(√κ)—but the paper needs a direct numeric evaluation of the exact Gaussian formula or a correct analytic argument.\n\nMinor notes: the complexity constant is attributed to a paper co-authored by Zygalakis; that is worth flagging but is not itself a flaw. The experimental section would be stronger with released code and repeated-run error bars, but I do not see a critical flaw in the methodology.\n\nVerdict: this deserves a serious referee. I would send it out with a request to rewrite Section 3.1.1 and Appendix B. The method and experiments are worth publishing; the theory section should not be accepted as is.","headline":"Genuinely useful proximal MCMC method with convincing experiments, but the O(√κ) theory in Section 3.1.1 is not proven as written and needs a correction.","tokens_in":26179,"tokens_out":1694,"would_cite":true,"duration_ms":18569,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","65C40","65L06"],"pacs":[],"model":"deepseek-v4-flash","headline":"This paper claims that replacing the Euler-Maruyama discretisation in proximal MCMC with an explicit stabilised Runge-Kutta-Chebyshev (SK-ROCK) approximation speeds up sampling from ill-conditioned log-concave imaging posteriors: for…","keywords":["proximal Markov chain Monte Carlo","SK-ROCK","explicit stabilised methods","Langevin diffusion","Moreau-Yosida regularization","imaging inverse problems","log-concave posterior sampling","effective sample size"],"falsifier":"Run a controlled Gaussian experiment with known condition number $\\kappa$ over several values of $\\kappa$, count gradient evaluations until the empirical 2-Wasserstein distance to the exact target falls below $\\varepsilon$, and compare the growth rate with $\\sqrt{\\kappa}$. If the count grows like $\\kappa$ rather than $\\sqrt{\\kappa}$, or if direct evaluation of $\\max_i R_1(z_i)^2$ for the recommended $s$ from equation (3.7) disagrees with $((\\sqrt{\\kappa}-1)/(\\sqrt{\\kappa}+1))^2$, the paper's central complexity claim is wrong.","tokens_in":25021,"feed_emoji":"🧮","tokens_out":8713,"duration_ms":85440,"temperature":0.7,"pith_summary":"This paper sets out to remove the main practical bottleneck of proximal Markov chain Monte Carlo for Bayesian imaging: the tiny stepsize that Euler-Maruyama discretisation forces on ill-conditioned posteriors. The authors replace that discretisation with the stochastic orthogonal Runge-Kutta-Chebyshev (SK-ROCK) scheme, which performs several gradient evaluations per step along Chebyshev extrapolants and can therefore take far larger stable steps. For Gaussian targets they argue that reaching a prescribed Wasserstein accuracy needs $O(\\sqrt{\\kappa})$ gradient evaluations rather than $O(\\kappa)$, matching the behaviour of accelerated optimisation methods. In experiments on image deconvolution, hyperspectral unmixing, and tomography, the resulting chains mix faster, produce larger effective sample sizes, and estimate posterior means and variances more accurately at equal computational cost. If correct, this gives imaging scientists an MCMC option whose cost scaling with conditioning is much closer to optimisation than before.","feed_headline":"Bayesian image sampling gets a √κ speed-up","feed_subtitle":"Chebyshev-staged sampling takes bigger stable steps and lifts effective sample sizes for deconvolution, unmixing, and tomography.","key_machinery":"The central object is the stochastic orthogonal Runge-Kutta-Chebyshev (SK-ROCK) integrator: an explicit $s$-stage Runge-Kutta method whose coefficients are generated from Chebyshev polynomials of the first and second kind, applied to the Moreau-Yosida smoothed Langevin diffusion. Each stage evaluates the gradient $\\nabla \\log \\pi_\\lambda(x) = -\\nabla f(x) - \\frac{1}{\\lambda}(x - \\mathrm{prox}^\\lambda_g(x))$ at extrapolated points, and the scalar Gaussian analysis reduces the integrator to the recurrence $X_{n+1} = R_1(z)X_n + \\sqrt{2\\delta}R_2(z)\\xi_{n+1}$ with $R_1(z) = T_s(\\omega_0+\\omega_1 z)/T_s(\\omega_0)$ and $R_2(z) = U_{s-1}(\\omega_0+\\omega_1 z)/U_{s-1}(\\omega_0)(1 + \\omega_1 z/2)$. The stability region of the method grows like $s^2$, so the maximum stable stepsize $\\delta^{\\max}_s = l_s/(L_f + 1/\\lambda)$ is much larger than MYULA's; choosing $s \\approx \\sqrt{\\eta/2}(\\kappa-1)$ stages makes the contraction constant $C = \\max_i R_1(z_i)^2$ approximately $((\\sqrt{\\kappa}-1)/(\\sqrt{\\kappa}+1))^2$. This explicit stabilisation is what converts the $O(\\kappa)$ gradient-evaluation complexity of Euler-Maruyama into $O(\\sqrt{\\kappa})$.","core_discovery":"The discovery is that the stability limitation of the Moreau-Yosida unadjusted Langevin algorithm (MYULA) is not intrinsic to proximal MCMC: it is an artefact of the Euler-Maruyama discretisation. By using the SK-ROCK integrator on the same smoothed Langevin diffusion, the admissible stepsize grows roughly quadratically with the number of stages $s$, so the chain can move meaningfully along slow directions without becoming unstable along fast ones. For a Gaussian target with condition number $\\kappa$, the paper derives an explicit bound $W_2(\\pi, Q_n)^2 \\leq W_2(\\pi, \\tilde{\\pi})^2 + C^n W_2(\\tilde{\\pi}, Q_0)^2$ for the numerical chain, with $C \\approx ((\\sqrt{\\kappa}-1)/(\\sqrt{\\kappa}+1))^2$, and concludes that the gradient-evaluation count to reach accuracy $\\varepsilon$ is $O(\\sqrt{\\kappa})$ when the bias term is dominated by the contraction term. Empirically, at equal numbers of gradient and proximal evaluations, SK-ROCK with $s=10$ to $15$ outperforms MYULA by factors of $21.77$ and $37.93$ in the effective sample size of the slowest mixing component in the cameraman deconvolution and hyperspectral unmixing experiments, and reaches the typical set of the posterior far more quickly.","pith_inferences":["The paper proves the $\\sqrt{\\kappa}$ bound only for Gaussian targets; if the mechanism is the quadratic growth of the stability domain rather than Gaussian-specific algebra, the same qualitative acceleration should appear in other log-concave targets, but that remains a conjecture to test.","The Wasserstein bound separates a contracted term and a bias term that is not accelerated, so the $\\sqrt{\\kappa}$ advantage is for moderate accuracy. A natural follow-up is a Metropolis-Hastings variant whose proposal uses the SK-ROCK step, which would remove the bias while keeping the large step.","The method pays $s$ gradient evaluations per iteration, so the practical win depends on the cost of those evaluations relative to the proximal map; for models with expensive proximal operators, the observed effective-sample-size gains may shrink. This trade-off is not quantified in the paper."],"forward_implications":["For strongly log-concave imaging models, the method permits stepsizes far above MYULA's stability limit at the same number of gradient evaluations, so chains mix much faster and posterior variance estimates become usable at practical budgets.","Under the paper's Gaussian analysis, the gradient-evaluation count for moderate accuracy scales as $O(\\sqrt{\\kappa})$ instead of $O(\\kappa)$, matching the scaling of accelerated first-order optimisation.","In the reported experiments, equal-budget effective sample sizes for the slowest mixing component improve by factors of roughly 14–38 relative to MYULA, which translates directly into lower Monte Carlo variance for estimators.","Because only the gradient and proximal oracle of MYULA is needed, the same integrator can be inserted into other proximal MCMC constructions that use MYULA as a kernel, as the authors note for future work.","The paper recommends $s=15$ stages for strongly log-concave problems and $s=10$ otherwise, giving practitioners a simple tuning rule for imaging models."],"supporting_citations":[{"why":"Defines the SK-ROCK integrator and its stability properties, the scheme the paper adopts.","marker":"[2]"},{"why":"Provides the $O(\\kappa)$ gradient-evaluation complexity for Euler-Maruyama Langevin that serves as the comparison baseline.","marker":"[17]"},{"why":"Introduces MYULA and the Moreau-Yosida smoothed target that SK-ROCK discretises.","marker":"[19]"},{"why":"Supplies the contraction constant $C \\approx ((\\sqrt{\\kappa}-1)/(\\sqrt{\\kappa}+1))^2$ and the Chebyshev-based stabilised-gradient analysis the complexity claim relies on.","marker":"[20]"},{"why":"Defines the Metropolis-corrected proximal benchmark used to assess bias and provides the proximal MCMC formulation.","marker":"[38]"}],"fun_headline_variants":["Proximal MCMC achieves √κ gradient count","Chebyshev-stabilised steps accelerate proximal MCMC","SK-ROCK method speeds up proximal MCMC for imaging","Bigger stable steps for faster Bayesian image sampling","Proximal MCMC hits √κ complexity via stabilised steps"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The claimed $\\sqrt{\\kappa}$ speed-up rests on a particular formula for how fast the Chebyshev chain contracts on Gaussian targets, and on the assumption that the bias from the numerical invariant distribution is negligible compared with the target accuracy; if that formula or that negligibility fails, the $\\sqrt{\\kappa}$ bound fails with it.","fun_headline_variants_meta":{"raw":{"variants":["Proximal MCMC achieves √κ gradient count","Chebyshev-stabilised steps accelerate proximal MCMC","SK-ROCK method speeds up proximal MCMC for imaging","Bigger stable steps for faster Bayesian image sampling","Proximal MCMC hits √κ complexity via stabilised steps"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.001112,"raw_usage":{"total_tokens":4658,"prompt_tokens":999,"completion_tokens":3659,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":615,"completion_tokens_details":{"reasoning_tokens":3578}},"tokens_in":615,"tokens_out":3659,"duration_ms":25900,"temperature":1.0,"reasoning_tokens":3578,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T11:29:48.605300+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run a controlled Gaussian experiment with known condition number $\\kappa$ over several values of $\\kappa$, count gradient evaluations until the empirical 2-Wasserstein distance to the exact target falls below $\\varepsilon$, and compare the growth rate with $\\sqrt{\\kappa}$. If the count grows like $\\kappa$ rather than $\\sqrt{\\kappa}$, or if direct evaluation of $\\max_i R_1(z_i)^2$ for the recommended $s$ from equation (3.7) disagrees with $((\\sqrt{\\kappa}-1)/(\\sqrt{\\kappa}+1))^2$, the paper's central complexity claim is wrong.","supporting_citations":[{"cited_title":"A BDULLE , I","cited_arxiv_id":null,"evidence_quote":"Defines the SK-ROCK integrator and its stability properties, the scheme the paper adopts."},{"cited_title":"D URMUS , S","cited_arxiv_id":null,"evidence_quote":"Provides the $O(\\kappa)$ gradient-evaluation complexity for Euler-Maruyama Langevin that serves as the comparison baseline."},{"cited_title":"D URMUS , E","cited_arxiv_id":null,"evidence_quote":"Introduces MYULA and the Moreau-Yosida smoothed target that SK-ROCK discretises."},{"cited_title":"E FTEKHARI , B","cited_arxiv_id":null,"evidence_quote":"Supplies the contraction constant $C \\approx ((\\sqrt{\\kappa}-1)/(\\sqrt{\\kappa}+1))^2$ and the Chebyshev-based stabilised-gradient analysis the complexity claim relies on."},{"cited_title":"P EREYRA , Proximal Markov chain Monte Carlo algorithms , Statistics and Computing, 26 (2016), pp","cited_arxiv_id":null,"evidence_quote":"Defines the Metropolis-corrected proximal benchmark used to assess bias and provides the proximal MCMC formulation."}],"review_version":1}