REVIEW 3 major objections 5 minor 2 cited by
Operator Splitting with Hamilton-Jacobi-based Proximals
T0 review · 3 major / 5 minor · reviewed 2026-08-15 · deepseek-v4-flash
Pith's one-line read Replacing exact proximal operators with Hamilton-Jacobi-based Monte Carlo approximations preserves almost-sure convergence in PPM, PGD, DRS, DYS, and PDHG.
desk verdict The central convergence theorems outrun the assumptions: Theorem 3.3's Lipschitz requirement is missing from every splitting theorem, but the flaws are identifiable and the framework is worth refereeing. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The HJ-Prox operator approximates $\mathrm{prox}_{tf}(x)$ by the ratio $\mathbb{E}_{y\sim N(x,\delta t I)}[y e^{-f(y)/\delta}] / \mathbb{E}_{y\sim N(x,\delta t I)}[e^{-f(y)/\delta}]$, a smoothed version of the proximal map evaluated from function values alone. The load-bearing bound is Theorem 3.3: for convex, LSC, $L$-Lipschitz $f$, the Monte Carlo error falls below $\sqrt{8J^*M^*/(\alpha N)} + \sqrt{n t \delta}$ with probability at least $1-\alpha$, where $J^*=\exp(2L^2 t/\delta)$. This bound converts into almost-sure summability, and Theorem 3.1 (perturbed Krasnosel'skiĭ–Mann) then supplies convergence once each algorithm map is known to be averaged.
What would settle it
Run HJ-Prox PPM with the Assumption 3.4 schedule on the one-dimensional problem $\min_x x^2$ (so $g=0$, $f(x)=x^2$): the functions are proper, LSC, and convex but not globally Lipschitz, making $J^*$ and $M^*$ infinite. Divergence, stagnation, or an error sequence that fails to be summable would show the stated assumptions are insufficient; convergence to $0$ would suggest the Lipschitz requirement can be relaxed.
Extended reading notes
Core claim
The central claim is that convergence theory for operator splitting does not require exact proximals; it only requires errors that vanish fast enough. The paper proves that HJ-Prox-based PPM (Theorem 3.5), PGD (Theorem 3.6), DRS (Theorem 3.8), DYS (Theorem 3.9), and PDHG (Theorem 3.10) all generate iterates that converge almost surely to a global minimizer, provided $\delta_k$, $\alpha_k$, $N_k$, and $t_k$ follow Assumptions 3.4 or 3.7. The proof strategy is to write each HJ-Prox iterate as an exact split iterate plus a random error, then apply a perturbed Krasnosel'skiĭ–Mann theorem pathwise after Borel–Cantelli establishes almost-sure summability.
Load-bearing premise
The Monte Carlo error bound that drives the whole argument assumes the approximated function is globally $L$-Lipschitz, which yields the finite constant $J^*=\exp(2L^2t/\delta)$; the convergence theorems state only proper, lower-semicontinuous, convex functions, so without adding global Lipschitzness the key summability bound does not follow.
Editorial extensions
If this is right
- Following Assumption 3.4, HJ-Prox PPM and PGD converge almost surely for convex objectives with no closed-form proximal at all; only calls to the objective are needed.
- Under Assumption 3.7, DRS, DYS, and PDHG keep a fixed step size $t$ while $\delta_k$ decays, so existing splitting implementations can swap in HJ-Prox for whichever term lacks a proximal formula.
- Because the Monte Carlo constant $J^*$ depends exponentially on $L^2 t/\delta$, splitting beats applying HJ-Prox to the whole sum: the Lipschitz constants of separate terms enter additively rather than through the square of their sum.
- Hybrid schemes that apply exact proximals where available and HJ-Prox only to the non-proximable component have smaller fixed-point residuals than fully approximate schemes, as the non-negative LASSO experiment demonstrates.
Reading between the lines
- A testable extension beyond the paper: the practical fixed-sample-size runs ($N=1000$) contradict the sufficient condition $N_k\to\infty$; a plausible explanation is that the $J^*$ bound is far from tight, and an adaptive $N_k$ rule could be benchmarked against the paper's fixed-$N$ performance.
- If the global Lipschitz assumption required by Theorem 3.3 is genuinely necessary, then a simple quadratic $f(x)=x^2$ would lie outside the theorem's scope despite satisfying the stated convexity assumptions; checking convergence on such a function would reveal whether the stated assumptions are sufficient.
- The same perturbed-KM route could be applied to other inexact proximal schemes, for instance random block-coordinate or randomized proximal updates, by replacing the Monte Carlo error there with the corresponding sampling error.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes replacing exact proximal operators with a Monte Carlo approximation called HJ-Prox in five standard operator splitting algorithms: proximal point (PPM), proximal gradient descent (PGD), Douglas–Rachford splitting (DRS), Davis–Yin splitting (DYS), and primal–dual hybrid gradient (PDHG). The main theoretical contribution is a set of almost-sure convergence theorems (Theorems 3.5–3.10) that model the approximation error as a summable perturbation of a perturbed Krasnoselskii–Mann iteration, using a new probabilistic bound on the HJ-Prox error (Theorem 3.3) together with a deterministic smoothing error bound (Theorem 3.2). The paper also reports numerical experiments on LASSO, sparse group LASSO, trend filtering, multitask learning, total variation denoising, non-negative LASSO, and overlapping group LASSO, claiming that HJ-Prox matches analytical solvers.
Significance. The idea of embedding a zeroth-order, derivative-free proximal approximation into splitting algorithms is timely and potentially impactful: it could remove the need for closed-form proximal operators in a broad class of convex optimization problems. The Monte Carlo error bound in Theorem 3.3 is a genuinely new technical result, and the paper correctly identifies splitting as a way to reduce the exponential-in-Lipschitz constant sample complexity compared with applying HJ-Prox to the composite objective. If the convergence theorems were valid under the stated assumptions, this would be a strong contribution. However, as detailed below, the central theoretical claim is currently unsupported because the error bound used throughout requires a global Lipschitz condition that is never assumed in the convergence theorems, and in the PDHG case the required assumption is not even satisfiable for standard problem data. The numerical experiments are encouraging but are run with fixed sample size N, which explicitly violates the theoretical sufficient conditions, so they do not compensate for the theoretical gap.
major comments (3)
- [Section 3.2, Assumptions 3.4/3.7 and Theorems 3.5–3.10] The convergence proofs for all five algorithms invoke Theorem 3.3, whose hypotheses require the function whose proximal is being approximated to be convex, lower-semicontinuous, and globally L-Lipschitz. Assumptions 3.4 and 3.7 reference the constants Jk and Mk of Theorem 3.3, but they never add a Lipschitz assumption on f, g, or g*. The theorem statements likewise assume only proper, LSC, and convex (plus smoothness for f or h in PGD and DYS). For example, f(x)=||x||^2 and g=0 satisfy the hypotheses of Theorem 3.5, yet for this data L=infinity, so J*=exp(2L^2 t/delta) and M* are infinite, the required sample size Nk in Assumption 3.4 does not exist, and the error bound in (7) is vacuous. This is not a minor technicality: every one of Theorems 3.5, 3.6, 3.8, 3.9, and 3.10 relies on bounding the HJ-Prox error by expressions that are undefined for non-Lipschitz functions. Adding a global Lipschitz assumption to the theorems would substantially narrow their scope (excluding, e.g., quadratic data terms) and still would not resolve the PDHG issue described below.
- [Theorem 3.10 and Appendix G (PDHG)] The PDHG iteration (12) applies HJ-Prox to prox_{sigma g*}, the proximal operator of the Fenchel conjugate of g. For typical proper LSC convex g such as g(u)=||u||_1 or an indicator function, the conjugate g* is extended-valued, taking the value +infinity outside a bounded set. Theorem 3.3 requires the approximated function to be real-valued and globally L-Lipschitz, so it cannot be applied to g* at all. The proof in Appendix G does not introduce any additional condition on g* (or a substitute error bound), and indeed the only stated assumptions on g in Theorem 3.10 are proper, convex, and LSC. Consequently, the PDHG convergence claim is unsupported by the paper's own error analysis. This flaw cannot be repaired by a simple added Lipschitz assumption, because g* is not Lipschitz for the standard functions that motivate PDHG.
- [Section 4 and Assumptions 3.4/3.7 (sample-size schedule)] The convergence theory requires the Monte Carlo sample size Nk to grow according to Nk >= 8Jk/alpha_k in Assumptions 3.4 and 3.7, yet all experiments in Section 4 use a fixed N=1000, with the text stating that 'this is not the setting of our theorem' implicitly. This is not itself an error, but it means the experiments do not provide evidence for the claimed almost-sure convergence regime; they demonstrate only empirical stability under much weaker conditions. That gap is acknowledged in Section 5, but the paper should state more explicitly that the theorems do not explain the reported fixed-N performance.
minor comments (5)
- [Theorem 3.3 and Appendix B] The notation in Theorem 3.3 is inconsistent: the statement reads 'Let [prox^delta_t f(x) denote...' with an opening bracket instead of a hat; the hat notation \widehat{prox} is used in Appendix B. Please fix the typographical error.
- [Section 3.2, first paragraph] The phrase 'Theorem 3.2assumes exact integral evaluation' is missing a space before 'assumes'; please correct.
- [Lemma C.1 in Appendix C] Lemma C.1 states the step sizes as t_k = 1/(k+1), but the proof then refers to a general divergent step-size assumption (sum t_k = infinity). Either state the lemma for general t_k with sum t_k = infinity and lsc convergence, or align the proof with the specific choice t_k=1/(k+1).
- [Section 4.4 and Appendix H.7] The number of overlapping groups is reported as 298 in Section 4.4 but as 289 in equation (175) in Appendix H.7. Please make these numbers consistent.
- [References] The reference list contains both 'Mckenzie, D.' and 'McKenzie, D.' entries; the capitalization should be uniform.
Circularity Check
No circular reductions in the derivation chain; the convergence theorems are conditional on explicit summability and sampling assumptions and are anchored in external results. The score reflects one minor non-load-bearing self-citation, not a circular step.
full rationale
The paper's central claim is a conditional convergence theorem: if the parameter sequences satisfy the summability and sampling conditions in Assumptions 3.4 and 3.7, then the HJ-Prox perturbations are summable almost surely and Combettes' perturbed Krasnosel'skiĭ–Mann theorem applies pathwise. The deterministic HJ-Prox bias bound is attributed to Crandall and Lions (1983) and is re-proved in Appendix A from strong convexity and integration by parts; the Monte Carlo bound in Theorem 3.3 is derived in Appendix B from self-normalized importance sampling, the Tsirelson–Ibragimov–Sudakov inequality, and a Poincaré inequality. These arguments do not assume the convergence conclusions. The only self-citation of note is the parenthetical '(Di et al., 2025)' in Section 3.1, which points to the authors' own prior workshop paper for the exact-integral case; the Monte Carlo theorems in Appendices C–G are proved independently from Theorem 3.3 and Theorem 3.1, so this citation is not load-bearing. A genuine weakness is that Theorem 3.3 assumes the approximated function is globally L-Lipschitz, while Theorems 3.5–3.10 and Assumptions 3.4 and 3.7 do not state a Lipschitz hypothesis, so J_k and M_k can be undefined for non-Lipschitz convex objectives such as a squared norm. That is a hypothesis mismatch and a correctness risk, not a circular reduction: the proofs do not define their conclusion into their assumptions, and no fitted parameter is relabeled as a prediction. Since no circular step can be exhibited, the score is kept low and is raised only slightly for the non-load-bearing self-citation.
Assumptions & free parameters
free parameters (2)
- Monte Carlo sample size N =
N=1000 in all experiments
- Smoothing schedule exponent p =
p>0, with delta_k = O(1/k^(2+p))
assumptions (5)
- domain assumption The objective functions f and g are proper, lower semicontinuous, and convex.
- domain assumption Every function whose proximal is approximated is globally L-Lipschitz.
- standard math Combettes' perturbed Krasnoselskii-Mann theorem applies to the constructed algorithm maps.
- standard math Crandall and Lions' deterministic error bound sup_x ||prox^delta_t f(x) - prox_t f(x)|| <= sqrt(n t delta).
- standard math Tsirelson-Ibragimov-Sudakov concentration and Poincare inequalities hold for the importance-sampling analysis.
Cite this review
Pith. "Pith review of Operator Splitting with Hamilton-Jacobi-based Proximals." pith.science (2026). https://pith.science/paper/SICNB4GS
@misc{pith2026260122370,
author = {Pith},
title = {Pith review of: Operator Splitting with Hamilton-Jacobi-based Proximals},
year = {2026},
howpublished = {\url{https://pith.science/paper/SICNB4GS}},
note = {Machine review of arXiv:2601.22370}
}
read the original abstract
Operator splitting algorithms are a cornerstone of modern first-order optimization, decomposing complex problems into simpler subproblems solved via proximal operators. However, most functions lack closed-form proximal operators, which has long restricted these methods to a narrow set of problems. Hamilton-Jacobi-based proximal operator (HJ-Prox) is a recent derivative-free Monte Carlo technique based on Hamilton-Jacobi PDE theory, that approximates proximal operators numerically. In this work, we introduce a unified framework for operator splitting via HJ-Prox, which allows for deployment of operator splitting even when functions are not proximable. We prove that replacing exact proximal steps with HJ-Prox in algorithms such as proximal point, proximal gradient descent, Douglas-Rachford splitting, Davis-Yin splitting, and primal-dual hybrid gradient preserves convergence guarantees under mild assumptions. Numerical experiments demonstrate HJ-Prox is competitive and effective on a wide variety of statistical learning tasks.
Figures
Figures from the paper (2 more)
Forward citations
Cited by 2 Pith papers
-
Implicit Neural Optimal Transport via Fixed-Point Optimization
A single-network implicit neural optimal transport method that solves the c-transform via proximal fixed-point iteration for stable, non-adversarial training.
-
Probabilistic Gaussian Homotopy: A Probability-Space Continuation Framework for Nonconvex Optimization
PGH deforms Boltzmann densities of nonconvex objectives into a soft-min homotopy whose Monte Carlo gradients bias descent toward low-energy regions and recovers a posterior-mean Moreau envelope.
Reference graph
Works this paper leans on
-
[2012]
The Graphical Lasso: New Insights and Alternatives
doi: 10.48550/arXiv.1111.5479. URL https: //arxiv.org/abs/1111.5479. McKenzie, D., Heaton, H., and Fung, S. W. Differentiating through integer linear programs with quadratic regular- ization and davis-yin splitting.Transactions on Machine Learning Research, 2024. Mckenzie, D., Heaton, H., Li, Q., Wu Fung, S., Osher, S., and Yin, W. Three-operator splittin...
-
[2014]
The PGD iterates are computed by applying the mappingTk(x) = proxtkg(x−t k∇f(x))
section 4.2. The PGD iterates are computed by applying the mappingTk(x) = proxtkg(x−t k∇f(x)). By Lemma D.2 and D.1, all Tk is an averaged operator andxk→x∗∈T k≥0 FixTk =X∗. The HJ-PGD iterates can be written as ˆxk+1 =dproxδk tg (ˆxk−tk∇f(ˆxk)) =T k(ˆxk) +εk,(119) where εk =dproxδk tkg(ˆxk−tk∇f(ˆxk))−prox tkg(ˆxk−tk∇f(ˆxk)).(120) From Theorem 3.3, we hav...
work page 1979
-
[2018]
with the HJ-Prox called on the overlapping group lasso. 27
-
[2020]
URL https://optimization-online. org/2020/02/7644/. Boucheron, S., Lugosi, G., and Massart, P.Concentration Inequalities: A Nonasymptotic Theory of Independence. Oxford University Press, 2013. Brice˜no-Arias, L. M., Chierchia, G., Chouzenoux, E., and Pesquet, J.-C. A random block-coordinate Douglas–Rachford splitting method with low com- putational comple...
arXiv 2020
-
[2024]
URL https: //arxiv.org/abs/2412.11485
doi: 10.48550/arXiv.2412.11485. URL https: //arxiv.org/abs/2412.11485. Zhang, Z., Fung, S. W., Kyrillidis, A., Osher, S., and Vardi, M. Y . Thinking out of the box: Hybrid sat solving by unconstrained continuous optimization.arXiv preprint arXiv:2506.00674, 2025. 10 HJ-Prox-based Operator Splitting A. Proof of HJ-Prox Error Bound For completeness and ease...
-
[4094]
PMLR, 2018. URL https://proceedings. mlr.press/v80/pedregosa18a.html. Pe˜na, V . H., Lai, T. L., and Shao, Q.-M.Self-Normalized Processes: Limit Theory and Statistical Applications. Probability and Its Applications. Springer, Berlin, Hei- delberg, 1 edition, 2009. ISBN 978-3-540-85636-8. doi: 10.1007/978-3-540-85636-8. URL https:// doi.org/10.1007/978-3-5...
arXiv 2018
Reviewed August 15, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.