Frictionless Hamiltonian Descent with Discretization and Parallel Optimization
Authors: Heng Zhu
Jun-Kun Wang
Venue: Accepted by TMLR ()

Abstract: Frictionless Hamiltonian Descent is a recently proposed optimization method that leverages a fundamental principle from classical mechanics. The algorithm is based on energy conservation of the Hamiltonian Flow, with resetting the kinetic energy at each iteration, and is shown to be a descent method. However, the idealized frictionless Hamiltonian Descent requires access to the oracle of the Hamiltonian Flow, while exactly implementing the Hamiltonian Flow becomes elusive when the underlying function is not quadratic. Motivated from considerable popularity of Hamiltonian dynamics in sampling, where a geometric numerical integrator is used to simulate the idealized Hamiltonian Monte Carlo, we consider Hamiltonian Descent with two kinds of integrator, which results in some new optimization dynamics. Moreover, we extend the original framework by introducing various forms of kinetic energy. This expansion yields a broad class of optimization algorithms and provides a fresh perspective of algorithm design. We further propose a novel parallelization technique for parallelizing the inherently sequential updates of the proposed optimization algorithms, where gradients at different points are computed simultaneously. The parallelization technique improves the actual running time by 2-3x in practice for multinomial logistic regression across a range of datasets when 4 GPUs is used, compared to approximating the Hamiltonian Flow in the standard sequential fashion by a single GPU.

Published in Transactions on Machine Learning Research (02/2026)
Frictionless Hamiltonian Descent with Discretization and Parallel Optimization
Heng Zhu
hez0070ucsd. edu
HDSI, University of Colifornia San Diego
Jun-Kun Wang ECE & HDSI, University of California San Diego
jkw005 @ucsd. edu
Reviewed on OpenReview: https: // openreview. net/ ptt/formm? id= =d=1110Q33Je
Abstract
Frictionless Hamiltonian Descent is a recently proposed optimization method that lever- ages fundamental principle from classical mechanics. The algorithm is based on energy conservation of the Hamiltonian Flow, with resetting the kinetic energy each iteration, and is shown to be descent method. However, the idealized frictionless Hamiltonian De- scent requires access to the oracle of the Hamiltonian Flow, while exactly implementing the Hamiltonian Flow becomes clusive when the underlying function is not quadratic. Moti- vated from considerable popularity of Hamiltonian dynamics in sampling, where a geometric numerical integrator is used to simulate the idealized Hamiltonian Monte Carlo, we consider Hamiltonian Descent with two kinds integrator, which results some new optimization dynamics. Moreover, we extend the original framework by introducing various forms of ki- netic energy. This expansion yields broad class of optimization algorithms and provides fresh perspective of algorithm design. We further propose novel parallelization technique for parallelizing the inherently sequential updates of the proposed optimization algorithms, where gradients at different points are computed simultaneously. The parallelization tech- nique improves the actual running time by 2-3x in practice for multinomial logistic regression across a range of datasets when 4 GPUs is used, compared to approximating the Hamiltonian Flow in the standard sequential fashion by single GPU,
1 Introduction
Hamiltonian mechanics a fundamental subject classical mechanics that provides a mechanism to predict the behaviors of many physical systems by unfolding the dynamics over time. Specifically, the Hamiltonian dynamic offers powerful mathematical framework for describing the dynamic of a particle in phase space that follows the physical laws in Nature. The nice and elegant Hamiltonian perspective exposes important features of the systems such as energy conservation, and it also allows for deeper study because of its rich mathematical structure. In machine learning and computational statistics, the applications of the Hamiltonian dynamics principle lead to the classical Hamiltonian Monte Carlo algorithm (HMC) (Duane
et al.
1987:
Neal 2011 Durmus et. al.] 2017 Wang & Wibisono, 2023 and its numerous variants (Kook
ct al.
2022
Lee al. 2021 Noble et al. 2024; Mangoubi
Vishnoi 2018 Ver Stoog Galstyan 2021
Hoffman & Gelman| 2014
Chaari al. 2016), which are for the task of sampling and have been found useful
in many other fields such systems biology (Kook et al. 2022).
Recently the principle of Hamiltonian mechanics has inspired designing optimization algorithms for solving an optimization problem, i.e., solving minxers [(x) (see e.g., Wang (2024); Diakonikolas & Jordan (2021): Teel et al. (0019)). Specifically, Wang (2024) proposes new optimization algorithm called frictionless Hamiltonian Descent. They define Hamiltonian H(x, v) as H(x,v) == f(x) 1 where is the position of the update, and is its velocity. The function f(-) in their definition of the Hamiltonian is the underlying function being optimized. From the perspective of Hamiltonian mechanics, the function value

---

Published in Transactions on Machine Learning Research (02/2026)
Algorithm Idealized Prictiontess HAMILTONIAN DESCENT (IDEALIZED-HD) (Wang 2024)
1: Input: an initial point X1 E R", number of iterations K, and a scheme of integration time (m7))=1-
for k = 1 to K do
3: 1/ Execute Hamiltomian Flow (1) for integration time The from the initial position X E IRd with the initial
velocity 0 Rd to obtain Xa+1, L.., the position of the particle time ?/-
4: (X++1, (i+1))Ve41) HF,, (Xx,0) 0).
end for
6: Output: XK+1-
(x) is the polential energy of one particle with its position x E R while the term 2 in the Hamiltonian is the kinetic energy. In other words, the total energy H(x, v) of a particle consists of its potential energy and its kinetic energy. With a specified Hamiltonian, the corresponding Hamiltonian Flow is described by a
system of differential equations:
dx 0H dv Ov dt
and 0H Ox
(Hamiltonian Flow)
(1)
dt
Wang (2024) observes that by energy conservation of the Hamiltonian flow, i.e., H Hx,,, vo) (xn, vo)
any time t, which is well-known principle in classical mechanics (Arnold 1989 Greiner 2003), it follows
f(x:) - 1 (xx,, vo) = vo) = f(zo),
when the initial velocity of the flow is zero, i.c., Vo 0 E Rd. Hence, the mechanism can be used to design a descenl method, which is Idealized Frills Hamiltonian Descent (Idealized-HD) (Algorithm []. Idealized HD is direct counterpart of Idealized- HMC in optimization, the initial velocity for executing Hamiltonian Flow is set to be a random vector from the normal distri- bution in HMCs. Wang (2024) analyzes Idealized-HD for the class of strongly convex quadratic problems, where exactly simulating the Hamiltonian Flow is feasible, as the Hamiltonian system ① has a. closed-form solution for this specific class of functions. also proposes Idealized Frictionless Coordinate Hamiltonian Descent (CHD), which updates a coordinate in each iteration, and shows that classical algorithms such as Gauss-Seidel methods, Successive Over-relaxation, Jacobi method, and weighted Jacobi method for solving lincar systems of equations are all instances of CHD. Precisely, they are CHD with different schemes of the integration tine for executing the Hamiltonian Flow. The optimization mechanism of rumning friction- less Hamiltomian Flow with resets by Wang (024) is novel in the sense that it challenges the conventional wisdom in optimization, which believes that incorporating a friction term into designing the optimization dynamics can facilitate the optimization progress. Specifically, a common approach in the optimization literature is to include a friction term yv into the update, e.g., di 3x - yv, where 7 > 0 serves as a damping factor to dissipate energy during the flow. This dissipation underpins various momentum- based methods, c.g., Nesterov (1983); Polyak (1964): (0644) O'Donoghue & Maddison (2019); Wang & Aberncthy (2018):
the only difference between them is that
Su et al. (2016): Maddison et al. (0018) However, Wang (2024) demonstrates that even without the friction term, the Tdealized Hamiltonian Descent (Idealized-HD) achieves an accelerated convergence rate as those of friction-based accelerated methods (Nesterov) 1983: Polyak, 1964 Wang & Abernethy 2018 Su et al.
2016) for strongly convex quadratic functions. We further note that recently, Fu & Wibisono (2025) studied Trictionless Hamiltonian Descent with randomized integration time and demonstrated that. their proposed scheme enables Hamiltonian Descent to achieve the provable accelcrated rate general strongly convex optimization.
neat, the Idealized-HD requires the exact Hamiltonian Flow, and hence whether it is realizable for minimizing general non-quadratic functions is not clear, since the differential equation system (1) does not have a closed-form solution in general. Motivated from the area
While the results presented in Wang 2024
of sampling (Mangoubi & Smith 2017 Chen & Gatmiry 2023 Chen et al. 2020 Durmus et al. 2017 Mangoubi Vishnoi [2088) Monmarché 2022 Camrud et al. 2023 Bou-Rabee & Sanz-Sernal (018)) we consider two geometric integrators numerically approximate the Hamiltonian Flow; one is the Störmer-
Verlet integrator (Hairer et al. 2003 Ernst. Hairer 2006 (a.k.a. the Leapfrog integrator) and the other is
2

---

Published in Transactions on Machine Learning Research (02/2026)
recently proposed stratified Monte Carlo integrator [Bou-Rabee & Marsden 2022]. When applying these integrators to approximate the Hamiltonian Flow in Hamiltonian Descent, the resulting updates exhibit novel momentum-liks dynamics, where the velocity v of the update x is an accumulation of gradients along the discretized Hamiltomian Flow, with a reset after cach flow. Interestingly, the use the Leapfrog leads to an optimization dynamic that compules the gradients at the points along the flow, while the straiffied Monte Carlo integrator generates dynamic where the gradient is computed at a point located a certain distance from the update, determined by random step in the direction of the current velocity.
We also extend the framework of optimization via Hamiltonian mechanics in Wang 2024 to more general casc. More preciscly, we consider a more general form of the Hamiltonian in our work:
H(x, H(x,v)== == f(x+ Φ(v),
where the kinetic energy d(-) (to be specified) satisfies <(v) 2 0, Vv E R' and Φ(v) = 0 if and only the velocity v = We found that under different forms of the kinctic energy Φ(v), the resulting dynamics of Hamiltonian Descent with the Leapfrog integrator lead different known and new updates, including Gradient Descent (GD); Normalized Gradient Descent, (Normalized GD), Sign Gradient Descent (Sign GD), Coordinate Descent (CD) well as their counterparts that incorporate certain novel momentun-liss dynam- ics. Specifically, while the use of a single integrator step to simulate the Hamiltonian Flow at each iteration of Hamiltonian Descent produces these known algorithms, executing multiple steps of the integrator induces some new dynamics that update the current position of the iterate via the gradients collected since the last reset of the velocity.
We further propose a method speed up simulating the Hamiltonian Flow via simultaneously computing the gradients at some points along the discretized flow. Specifically, given access to a computational environment with multiple GPUs, we aim to accelerate the process of numerically simulating Hamiltonian Flow by having each GPU compute a gradient at different steps of the flow in parallel and in an iterative fashion. Although the Hamiltonian Flow is sequential by nature, i.e., a gradient at the current position is used to update the iterate before computing the gradient at the next step, we found that via the underlying mechanism of Picard Tteration (Teschl, 2012 Khalil 2002 Andrade et al. 2023), computing gradients at different steps in parallel becomes feasible, which in turn leads to a faster optimization process in practice. We conduct experiments the proposed Hamiltonian Descent via. Parallel Picard Iteration for solving multinomial logistic regression across multiple datasets. We found that the improvement of the speedup in running time is strikingly consistent across all the datasets. For example, compared to the case of using single GPU, the parallel optimization with 4 GPUs reduces the actual running time of simulating a Hamiltonian Flow by 2.13x on a8a dataset and by 3.23x on MNIST dataset.
2 Frictionless Hamiltonian Descent with Geometric Numerical Integrators
We begin by observing that the idea of optimization via Hamiltonian mechanics in Wang 2024 be more applicable when considering the Hamiltonian in its general form 3 IIndedd, the energy conservation along the Hamiltonian Flow 1 still holds.
Lemma 1. (See e-g-, Arnold 1989). Greiner is 0 along the Hamiltonian Flow i.e., dt
2003 dF O.
The time derivative of lime-indepenndnnt Hamiltonian
Proof.
dH de
00 0x
2H Dv
dv
HH
> +
0H
2H 8x
0.
8x
DE
dE
0x dv
Ov
Using Lemma and that the kinetic energy is zero when the velocity is zero, i.c., 6(0) 0, it becomes evident that Idealized-HD (Algorithm ! with the general notion of the Hamiltonian satisfies:
f(xx+1) + ((((=11 =H(x va+1) =H((x, 0) f(xx), Vk,
(4)
3

---

Published in Transactions on Machine Learning Research (02/2026)
Algorithm 2 FRICTIONLESS HAMILTONIAN DESCENT WITH GEOMETRIC NUMERICAL INTEGRATOR.
Input: an initial point X E RR, number of iterations K, scheme of integration time (nx ) and the step size used in geometric numerical integrator 0. Specify the Hamiltonian H (x, v). 3: for k= == 1 to K do
(K)
Obtain (xX41, V++1) = INTEGRATORR, 0), where S₂ == // Numerically simalate the Hamillonion Flow, where INTEGRATOR can be, for example, LeapFrog (Algo- rithm or Stratified (Algorithm 4
5: end for
6: Output: Xk+1.
where the kinetic energy nonnegative, i.e, Φ(v++1) 2 Equation (4) reduces to (2) when x(-) = 1 11. On the other hand, it also encapsulates the cases when the kinctic energy <(-) is in different form. For example, in the following, we will also consider <bv) || [[vll, and ||v||
As described in the introduction, the differential equation system ( does not have a closed-form solution in general, and one may need to use a geometric numerical integrator to simulate Hamiltonian flow. Wo hence propose Algorithm 2 which utilizes an integrator to simulate the contimuous- dynamic in discrete time. We note that approximating the flow via a geometric numerical integrator standard in sampling (Mangoubi & Smith 2017 Chen & Gatmiry 2023 Chen et al. 2020
Durmus et al. 2017 Mangoubi & Vishnoi 2018 Monmarché 2022 Camrud et al. 2023 Bou-Rabee & Sanz-Serna 2018). and perhaps one of
the most popular integrators is the Leapfrog integrator (Hairer et al. 2003 Ernst Hairer 2006). We thus
consider Hamiltonian Descent (HD) (Algorithm 2 with the Leapfrog integrator first, i.c., let INTEGRATOR
on Line 4 in Algorithmm be LeapFrog (Algorithmn 3). We note that since Leapfrog integrator is symplectic integrator, thus not exactly conserving the Hamiltonian or total energy. However, it is still widely used in sampling area due to its bounded energy crror. The energy oscillates near the true Hamiltonian at the order O(6"), where is the time step size (Neal 2011 Ernst Hairer 2006). Unlike other non-symplectic integrators which brings systematic energy drift, the energy error of Leapfrog integrator grows slowly with time, making it good choice to approximate the continuous Hamiltonian dynamic.
In this paper we use x and v to denote the variables inside one iteration. It turns out that when a single Leapfrog step (i.e., S = 1) is used to simulate the Hamiltonian Flow at each iteration of HD, the resulting dynamics recover Gradient Descent, Normalized Gradient Descent, Sign Gradient Descent, and Coordinate Descent, under different forms of the kinetic energy. We summarize the results in Lemma 2
Algorithm 3 THE LEAPFROG INTEGRATOR LeapFrogs, (Xo, Vu).
Algorithm 4 THE STRATIFIED INTEGRATOR
Stratified (X., Vo). (same input as Alg. 3
1: Input: an initial point (xo, Vo), # of Leapfrog steps S, and the step size 0 > 0.
1: for 0 to S - 1 do *+4 = x, + C. 0,H V.), where So
2: X
for s to
do V V. H Vs).
Uniform(0, 0)
3: V V, $ V [x++3] V.).
4: X.+1 = X. + + 0.0 Xx, Vo+
4: 8,+1 = x, + 0 H0 X,, V
V++1=
Vo+.)
5: 7+41=V.41 2
0
H x+ v.).
6: end for
6: end for
7: Output: LeapFrogs (Xo, vo) = (Xs, Vs)-
7: Output: Stratified / (Xo, vo) = (Xs, Vg)-
Lemma 2. Consider Frictionless Hamiltonion Descent with the Leapfrog integrator (HD-LF), i.e., Algo-
rithm 2 with INTEGRATOR LeapFrog. Set a single Leapfrog step Sk at each iteration k. Then, we have
4

---

Published in Transactions on Machine Learning Research (02/2026)
1. Setting Φ(v) = 2/ recovers GRADIENT DESCENT (GD):
0
XK+1 Xx V)((x))
2. Selling Φ(v) = ||vi|2 recovers NORMALIZED GD:
Vx(x))
Xk+1 = X
Vf(x))]
3. Setting Φ(v) = ||vi[a recovers SIGN GD:
Xk+1 = X - signn (V)()))
4. Setting ø(v) = |ov recovers COORDINATE DESCENT (CD):
X++1 = X - sign (V (Vf(x)) [max] eswax!
where mmax arg ma [xx))l]
Proof. We provide the proof of the first item here and defer the proof of the rest of the items to Appendix A
When ((v) = 1 we have Hx, v) = vx(x) and VvH(x, v) = From Algorithm 3 with the initial
condition Xo f- Xx and vo + 0, one can see that the dynamic with a single Leapfrog step 18:
f4 vo - (((o) = 0
0
2 fx(x))
8²
X1 = x + v = Xx 2 vf(x)).
(5)
With only a single step, Leapfrog outputs x1, which is the position of the update X at the next ileration k ++1 1. We hence obtain the update of GD, i.e,, X6+1 Xx. T (x(x).
On the other hand, when more than a single iteration executed in Leapfrog, the dynamics differ from the known ones. For example, under the form of the kinetic energy d(v) ] [|vv the resulting dynamic with S Leapfrog steps
S02 Vf(x))) 0² S fx,)
Xg =
(6)
In addition, when the form of the kinctic energy Φ(v) ||v|22 the resulting dynamic with Leapfrog steps is
S-1 1 ( + Σ j=1 Vf(x,) [x(xx)>||
(x(x))
Xs Xo
0
V (x,) + E Vf(x))
Furthermore, when the form of the kinetic energy Φ(v) || the resulting dynamic with S Leapfrog steps
is
X3 X
Osign V f(xo))
8-1
Σ V
0
sign
(8)
2
5

---

Published in Transactions on Machine Learning Research (02/2026)
Under the form of the kinctic energy Φ(v) // the resulting dynamic with S Leapfrog steps is
X5 X Osign (Vf(xo)) max))
C.p
8-1 ( +
0
Σ sign
V((x)) [
(9)
where iaax is the index corresponding to the largest absolute value of f((o). and tix is the index corre- sponding to the largest absolute value of vector (fV ビ VS(X;))
In optimization literature, Heavy Ball Polyak unfolding its dynamic over steps, one finds
1964 is one of the popular ummentuu methods. After
ys yo
fly,),
(10)
where n is the step size and /3 the momentum um parameter. By comparing the dynamics, (6) and 10p, we see that simulating single Hamiltonian Flow in HD via Leapfrog yields a certain momentum-lin dynamic with the momentum parameter 3 = 1 (i.e., without roducing the weights of the accumulated gradients). On the other hand, at the end of the flow in HD, the velocity v and hence the kinetic energy are reset to 0. The reset of the kinetic energy in HD is necessary to make the dynamical system to be dissipative in order to converge. Without the reset, by energy conservation of the Hamiltonian Flow, the update might continue moving past a desired point and hence fail to converge. This mechanism of HD differs from existing friction- based momentum methods [Maddison et al.] 2018 O'Donoghue & Maddison 2019: De Luca & Silverstein, 2022 França et al., 2020 Mai & Johansson 2020 Wibisono et al.] 2016 Willon et al., [202)), which typically involve movement with friction to dissipate kinetic energy, e.g., via specifying the value of the momentum parameter B to be $ < 1.
We further consider Hamiltonian Descent (HD) (Algorithm 2 with the recently proposed Stratified Monte Carlo integrator (Bou-Rabee & Marsden 2022) (HD-ST), i.e., instantiate INTEGRATOR on Line in Algo- rithm _ as Stratified (Algorithm [ The Stratified Monte Carlo integrator was shown to outperform the
Leapfrog integrator for HMC in sampling (Bou-Rabee & Marsden
2022). When the kinctic energy is in the
form Φ(v) = v the dynamie of HD-ST becomes:
X =x, VoV.
Vs+ =v. - fxxx++)))
8,+1 X, + +V.+
Vot1 -V++
2 v ( 8+)))
(11)
Therefore, different from HD-LF, the gradient at each intogrator step is computed at a different point than x, in HD-ST. More precisely, the point where the gradient is computed is determined by a random step Ç, in the direction of the current velocity v, from x,. Interestingly, when only a single integrator step is executed in each iteration k, the resulting dynamics recover the known ones Lemma albeit with random step size Ce.
3 Speeding Up Hamiltonian Descent via Parallel Optimization
One challenge in simulating continuous-time Hamiltonian Flow via discretization is that while reducing the step size 0 brings the approximated flow closer to the exact Hamiltonian flow, it comes at the expense of requiring more steps in each itcration, as the number of steps is S₁ |n/00] at & given a fixed intogration time Tik thereby might increase the lime needed lo complete an iteration of HD. To tackle this issue, we aim to answer the following question: can we run the steps in an approximated Hamiltonian Flow in parallel, given access to multiple GPUs that allow computing gradients at diffcrent points simultancously? The answer is affirmative, and our method to achieve the goal based on parallel implementation of Picard Iteration.

---

Published in Transactions on Machine Learning Research (02/2026)
Picard Iteration is classical algorithm that oftentimes serves technique to prove the existence of the dx
unique solution of a differential equation, dr Q(x,tt) where the underlying function Q(:,)) is continuous
and Li-Lischhizz, ie., [Q(x, Qyy,t) S <Lo||x y (see e.g., Khalil 2002 at each iteration m of Picard Iteration, the update is
Teschl 2012 Specifically,
(m+1) + f Q(x(m)) .r)dr. (r)dr
(Picard Iteration)
(12)
However, evaluating the integral might not be possible in general, and one may need to resort to numerical computations. Specifically, a way to numerically approximate the update 12 by computing x; (m++1)
Xo + ((x))") j0), where 0 > 0 is the step size. The update x; (m+1) can be viewed a discrete-time
approximation to x (m+1) of 12 in continuous time, with t == s0.
We now consider Picard Iteration for simulating Hamiltonian flow (1) and present here the case when the kinetic energy is <6-) = 211- From ① and (12). the update can be shown to be:
(m+1) vo J' ( dr
(m+1) xo tvo S ((x))
(13)
Then, by discretizing the update scheme 13 we obtain
(m+1)) vo-
fxx (m)) Vs E [S].
1=0
mm+1) Xg (s0) 0 (s i)0v f(x (m)),
(14)
where 0 is the step size, 17 is the integration time of the Hamiltonian flow, and S = [n/0] is number of steps to simulate the flow. We further note that the initial velocity in HD is zero (i.e., Vn = 0). The update scheme in (14h implies that one can execute the steps in the flow in parallel. More precisely, one would like to simulate the flow for S steps, then all (())s s can be initialized randomly at the beginning (i.e., 0). Furthermore, at each iteration m, gradients at different steps in the flow, i.c., V (x:"") (x. can be computed in parallel,
and obtaining x (1++1) from (") can also be done in a parallel fashion, which involves the gradients at those x (m))s, ss where < We refer the reader to Figure 1 for the illustration of this parallel update scheme. It is known that when m = S, then x5 is guaranteed to be exactly that of sequentially computing the update for S steps, i.e., x = Ks. We formalize this argument. with the following lemma. Fappriiaally we found that the number of iterations does not need to be as large S to achieve competitive performance that of S steps.
Lemma 3. When the number of Picard ileration is S, then the update of Picard ileralion x(5) is evactly the output of sequentially computing the discretized Hamiltonion flow after S steps, i.c., x Xg-
Proof. We prove this argument by induction. The strategy is to show that at step 8, we have x!.) 8, for all j 2
For the base case s = 1, the condition is true since for all j 2 1, we have
x = (( No (x)) = X1.
we have (5) x, for all j 2 Then, at next step + 1, for
Now assume that at particular step all j 2 s we have
x x+1
-x0-
0
88 +1 - ,パ 11 _0-- 0 + (){)) = X..
1=0
1=0

---

Published in Transactions on Machine Learning Research (02/2026)
Picard Iteration
x
(0)
X1
(0)
X 1
(0)
X
(0)
(1)
(1)
(1)
x
1-1
X
(m)
X
(m)
X
(m) 5-1
(m) S
The update ((m) receives gradients from the
Figure 1: Hamiltonian Flow via Parallel Picard Iteration. preceding ones x (m-1) Vi < which can be done in parallel.
since j' s 8. Therefore, the condition holds at step 8 + 1, which completes the induction. That is, for all j 2 S, we have (ff) Xs, which also implies that xg (8) Xs.
Algorithm 5 FRICTIONLESS HAMILTONIAN DESCENT VIA PARALLEL PICARD ITERATION
Input: initial point xı Rdd number of iterations K, scheme of integration time (mx) scheme of step sizes (ma) window size W, and a threshold parameter т.
for l- = to do
3: // /Numerically simulate the Hamiltonian Flow (Xt+1,-) HF., (Xx, 0).
5: Init x (m)
x-+=-1 f- Vi E [W].
Set St = -%%
7:
while < St do
8:
for w = to W do Compute gradients (x(( end for
//ean be executed in porallel
10:
11:
for w = to W do (m++1) (m) 02 (i (x
/ean be creculed in parallel
12:
13:
end for
14:
for w 1 1 to W do
//ean be executed in parallel
15:
Calculate relative error
Cstw = || (m+1) - x ||
16:
end for
17:
Update the stride Ht of the window:
R min(W, mi f). for r= to R do
18:
19:
Init points that window will cover at next inner iteration: x (m++1; end for
** +- (m++1)
*+W"
20:
21:
Update 8 R and +- m 11
22:
Update W 4- min(W, St - s).
(m+1) to initialize the update at the next outer iteration + L. Set X++1 + Xs+
23:
24:
end for
25:
Output: XK+1-
The idea of parallelizing the steps in Picard iteration was recently applied the task of quickly generating images by diffusion models (Shih et al., 2024). where obtaining image from diffusion model involves solving ordinary differential equation. Motivated from their results, consider parallelizing the steps in Hamiltonian Flow for optimization. However, when the number of steps S is large, maintaining multiple (m) at different steps becomes memory consuming. Hence, adapt sliding- approach in Shih

---

Published in Transactions on Machine Learning Research (02/2026)
et al. (2024), where only the steps of the flow in current window are updated, and the window moves forward along the flow at each m. After the updates of the the steps in the current window, the relative error, i.e., [x(=++1) || for all steps in the window is caleulated, which determines the stride of the window for the next iteration. Algorithm 5 details the proposed method.
4 Numerical results
We experiment with the proposed methods of frictionless Hamiltonian Descent via discretization and par- allel optimization and report their performance. Specifically, Subsection 4.1.1 reports the empirical results of HD-LF (Algorithm [ with Leapfrog) and HD-ST (Algorithm / with Stratified Integrator) under four dif- ferent forms of kinetic energy o(-) as well their counterparts such as GD, Normalized GD, Sign GD, and Coordinate Descent, with and without Heavy Ball (HB) momentum for multinomial logistic regression, while Subsection 4.2 provides results of neural network training. Appendix B.1 and B.2 include more experimental results of minimizing non-convex function and matrix sensing.
In subsection 4.3 we evaluate the performance of the proposed HD via Parallel Picard Iteration (Algo- rithm We primarily display results of multinomial logistic regression on different classification datasets. Our results clearly shows consistent speedup across the datasets. Appendix [B3 also displays the speedup across different cases of the kinetic energy (i.e., vv(v 1 |vl22, |v and [vlo))
Table 1: Detailed information of the datasets.
DATASET
ALL SAMPLES
TRAINING SET
TEST SET
# OF FEATURES
A8A
32561
22696
9865
123
A9A
48842
32561
16281
123
MNIST
60000
50000
10000
784
USPS
9298
7291
2007
256
The details of the datasots used for the experiments are described in Table The experiments are conductod on an Amazon EC2 P3 instance with 4 Tesla V100 GPUs. The code for reproducing the experiments avail- able on https://gi thub. com/ hwt.uu..oo///oaaaae howaah/aaaaaaallll and can also be accessed via the link in the supplementary material.
4.1
Hamiltonian Descent with Geometric Numerical Integrator
4.1.1 MNIST dataset
Figure 2 and Figure 3 compare HD-LF, HD-ST, and their counterparts with different forms of the kinetic energy for the task of multinomial logistic regression. When the kinetic energy Φ(v) is the squared l2 norm, the counterparts of HDs GD and HB (Polyak 1964), where the former is without momentum and the latter is momentum method. Similarly, when o(v) [v]]2, the baselines of HDs are Normalized GD and Normalized GD with HB dynamics (Normalized GD-HB); when Φ(v) = ||v|| the comparators are Sign GD and Sign GD with HB momentum (Sign GD-HB); and for the other form of Φ(v) = ||v|loo, Coordinate Descent and Coordinate Descent with HB momenturn (CD-HB) serve as the counterparts of HDs. For HD-LF and HF-ST with four forms of kinetic energy, we report their performance under the specification where the step size 8 is iterations in Algorithm 2 For the counterparts of HDs, the step size is r 0.01, which is the best from the set [0.001, 0.01,,11....0. 0.1], and the value of the momentum parameter is B 0.9. The parameter 8a is set to 0.01, since Lemma 2 shows that HD is equivalent to GD with a step size of order 0(02). Note that in all experiments we use the same initialization for compared algorithms but we plot the results after one round of itcration, not from the initialization.
6² 0.01 and the number of integrator steps is S 10 for all the
All the datasets available http:://www...iatatt.tt..///t/tt

---

Published in Transactions on Machine Learning Research (02/2020)
Test Accuracy
Test Accuracy
80
60
40
Norma Ized GB
HB
Noma feed 6D HB
20
HD-UF
20
HD-UF-12
HD-ST
HD-ST42
101
100
Computations
Computations
(a) (0())
Y 11
(b) Φlv) 1|vll2-
Test Accuracy
Test Accuracy
80
CD
CD-HS HB
HD4F-inf
Test accurec (%)
ID:TTiml
60
70-
50
30
gn
go HB
20
ND-UP-43
HD- HD-S74L1
150
250
Computations
Number Computations
(c) Φ(v) = ||v||
(d) @(v) ||
Figure/ 2: Test accuracy on MNIST dataset. Each line corresponds HD-LF, HD-ST, or their counterparts with and without HB momentum, under four different forms of the kinetic energy o(-).
As HD has inner steps for simulating the Hamiltonian Flow at each iteration k, we report the testing accuracy and training loss vs. the total number of gradient computations in Figure 2 and Figure 3 to have a fair comparison. Specifically, we plot the progress of the algorithms up to 300 number of gradient computations; in the case of HDs, this corresponds to K = 30 iterations with S = 10 steps to simulate the flow at each iteration k. From both figures, HD-LF and HD-ST consistently outperform their counterparts without momentum under different cases of the kinctic energy except when Φ(v) ||v|| (i.c., Sign GD). Furthermore, both HD-LF and HD-ST have competitive performance with their counterparts that incorporate the momentum dynamics. We note that in the literature, Sign GD does not necessarily converge to global optimal point in convex optimization, see c.g., Karimireddy et al. (2019) for an example of provable failure of Sign GD. This is aligned with our observation in Subfigure (c) of both figures, where the algorithms scemingly converge to different points. Interestingly. HD-L.F and HD-ST (and Sign GD with HB momentum) converge faster and achieve higher test accuracy than Sign GD, while having higher training loss compared to Sign GD.
4.2 Training ResNet on CIFAR10 Dataset
We further conducted experiments on training ResNet18 model on CIFAR10 datasct with Hamiltonian Descent. We still use Euler discretization and kinetic function Φ(v) || in HD and compare with GD and HB.
In the experiments we use stochastic gradient rather than the full gradicnt. We actually implement a stochastic version of Hamiltonian Descent and compare it with SGD and SGD with Momentum. To simplify, we still use HD, GD and HB to represent these algorithms. The learning rate of these algorithms is 0.01 and the batch size 64. The momentum of HB 0.9 and the integration time 0.2. Fig. 4aa) and (b) show the testing accuracy and training loss of the HD, GD and HB algorithms. The HD can get similar final performance GD and HB, meaning it can be effectively applied to neural network training. Fig. 4(c) and
10

---

Published in Transactions on Machine Learning Research (02/2020)
Training Loss
Training Loss
6D
Nommaize ized 6D
HB
Norma leed 6D H6
2.0
HD HD-ST
HD LF
-+o-
HD-ST42
1.5
1.0
0.5
Computations
Computations
(a) (0v) 2 llv
(b) Φlv) l|v
Training Loss
Training Loss
20.0
Snn CD
2.50
Sinn CD-HB
HI-IF--))
2.25
HID-STLL]
2.00
12.5
10.0
1.50
1.25
5.0
CD
1.80
CD HB HDLF-IM
2.5
MDSTIN
Computations
Computations
(c) Φ(v) ||v||ı-
(d) Φ(v) ||
Training loss on MNIST dataset. Each line corresponds to HD-LF, HD-ST, or their counterparts with and without HB momentum, under four different forms of the kinetic energy o(-).
Figure 3:
Table 2: Execution time of frictionless HD via Parallel Picard Iteration with 4 GPUs (ø(v) = [
RUNNING TIME
A8A
A9A
MNIST
USPS
1 GPU
93.2s
99.6s
403.1s
83.1s
4 GPUs SPEEDUP
47.7s
50.3s
128.9s
42.7s
1.95x
1.98x
3.13x
1.95x
(d) display the testing accuracy and training loss of HD with different integration time. It seen for HD with stochastic gradient, the impact of integration time is not obvious. We can easily choose an integration time to get an acceptable performance.
4.3 Hamiltonian Descent via Parallel Picard Iteration
Figure and compare the proposed frictionless Hamiltonian Descent with Parallel Picard Iteration (Algo- rithm when using GPU, GPUs, and 4 GPUs, where the kinetic energy Φ(v) = 1 l/v The step size O.I, the number of integrator steps is 20, and the threshold for relative crror т is 0.001. The sliding window size corresponds to the number of GPUs used in our experiments. Choosing the threshold 7 can be delicate process: larger error threshold reduces the accuracy of the Picard iteration approximation, convergence is not fully achieved within each window, whereas a smaller threshold enforces more iterations within a single window, significantly increasing the computational cost. In our experiment, we choose the threshold that is the best. among the set [0.0001, 0.001, 0.01,
Both figures show consistent speedup across all datasets when multiple GPUs are used exccute the proposed algorithm. A significant speedup observed when using GPUs compared to 1 GPU, and similarly, even greater speedup is observed with the of 4 GPUs. In addition, conduct the experiments for
11

---

Published in Transactions on Machine Learning Research (02/2026)
Testing Accuracy
Training Loss
GD
Testing Accuracy
80
HD
Loss
1.0
60
20
30
40
10
20
30
40
50
Number Epochs
Number Epochs
(a) Testing Accuracy.
(b) Training Loss.
Testing Accuracy
Training Loss
1.6
HD-0.1
HD-0.2
Testing Accuracy
80
1.2
ND-0.5
HD-0.3
70
60
HD-0.1
0.6
HD-03
HDD005
Number Ebochs
Number Epochs
(c) Testing Accuracy.
(d) Training Loss.
Figure Testing Accuracy and Training Loss of neural network training with different algorithms and parameters.
Accuraay
1-991 Accuracy
Tesurayy
Accurayy
50.0
725
32
Nuves Time 1em0
L3
(43558
(a) a8a dataset.
(b) a9a dataset.
(c) MNIST dataset.
(d) usps dataset.
Figure 5: Testing accuracy vs. running time (seconds) on the four datasets HD via Parallel Picard Iteration (Algorithm with d(v) 1
Thaning Loss
Training Loss
Training Loss
Training Loss
2.5
2 CTU
0.52
CrU
3.8
2
0.50
3.4
3.45
0.33
6.1
04
6.5
20
15
Suning Time (4200
Auaning C42000
65
(a) a8a dataset.
(b) a9a dataset
(c) MNIST datasct.
(d) usps catasett
Figure 6: Training loss vs. running time (seconds) on the four datasets of HD via Parallel Picard Ttera- (Algorito with @(v) 3
12

---

Published in Transactions on Machine Learning Research (02/2026)
different forms of the kinetic energy and report the results in the supplementary, where also observe that with the use of multiple GPUs, Hamiltonian Descent via Parallel Picard Iteration reduces the actual running time, compared to the case of single GPU. We note that due to the criteria for updating the stride of the sliding window in each inner itcration of Algorithm the actual number of steps when approximating the Hamiltonian flow in the inner loop of Algorithm 5 can vary across different dalasels, resulting in varying speedup via multiple GPUs across different datasets. Though the speedup can vary, we found that the proposed HD via Parallel Picard Iteration is effective in accelerating the process of simulating a flow in the experiments.
5 Conclusions
In this paper, we make the idea of frictionless HD more applicable by presenting different forms of kinetic energy and considering frictionless HD with the Leapfrog integrator and the Stratified Monte Carlo integrator. We found some interesting connections of the resulting dynamics to the existing ones and also some novel update schemes. We further propose speeding up approximately simulating the flow via parallel optimization and show exciting empirical results of the acceleration via multiple GPUs. Nevertheless, couple of directions remain open in this work. First, Frictionless HD with Parallel Picard Iteration (A (Algorithm 5 requires W copies of the model to be loaded into memory in order to perform parallel-in-time updates, where W is the window-size parameter. This results in a significant memory footprint, especially for large models, which in turn may limit the potential applications of the proposed algorithm. How to elegantly reduce this seemingly inherent and substantial memory cost is therefore worthwhile direction for future research. Any concrete and effective approach to reducing the memory footprint could substantially improve the applicability of the algorithm to large-scale model training, such as large language models, which typically rely on large transformer architectures. Second, while we show that discretizing the Hamiltonian flow with multiple integration steps induces different momentum-like dynamics under different choices of kinetic energy, potentially highlighting the modularity of the frictionless Hamiltonian Descent framework. When and how these new optimization dynamics yield tangible performance benefits remain underexplored, especially in the context of non-convex optimization.
Acknowledgements
The authors thank the anonymous reviewers and the action editor for their useful comments. JW also appreciates Andre Wibisono for the discussions. This work is supported by NSF CCF-2403392.
References
Francisco Andrade, Mario AT Figueiredo, and Joao Xavier. Distributed Banach-Picars iteration: Application to distributed parameter estimation and PCA. IEEE Transactions on Signal Processing, 71:17-30, 2023.
V. 1. Arnold. Mathematical methods of classical mechanics. Springer, 1989.
Nawaf Bou-Rabee and Milo Marsden. Unadjusted Hamiltonian MCMC with stratified Monte Carlo time integration. arXiv preprint orXiv:2211.11003 2022.
Nawaf Bou-Rabee and Jesús Maria Senz-Serna. Geometric integrators and the Hamiltonian Monte Carlo method. Acta Numericn, 27:113 206, 2018.
Evan Camrud, Alain Oliviero Durmus, Pierre Monmarché, and Gabriel Stoltz. Second order quantitative bounds for Unadjusted Generalized Hamiltonian Monte Carlo. arXie preprint arXiv:2306.095.1 2023.
Lotfi Chaari, Jean- Yves Tourneret, Caroline Chaux, and Hadj Batatia. A Hamiltonian Monte Carlo method 5594, 2016.
for non-smooth energy sampling. IEEE Transactions on Signal Processing,
Yuansi Chen and Khashayar Gatmiry. When does Metropolized Hamiltonian Monte Carlo provably outper- form Metropolis-adjusted Langevin algorithm? arXiv preprint xxXiw:234404722, 2023.
13

---

Published in Transactions on Machine Learning Research (02/2026)
Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step gradients. The Journal of Machine Learning Research, (1(:::64 3717, 2020.
Damek Davis, Dmitriy Drusvyatskiy, and Liwei Jiang. Gradient descent with adaptive stepsize converges (nearly) linearly under fourth-order growth. Mathematical Programming, pp. 1 66, 2025.
Giuseppe Bruno De Luca and Eva Silverstein. Born-infeld (BI) for AT: energy-conserving descent (ECD) for optimization. In International Conference on Machine Learning, pp. 4918 4936. PMLR, 2022.
Jelena Diakonikolas and Michacl I Jordan. Genoralized momentum-based methods: A Hamiltonian perspoc- tive. SIAM Journal on Oplimization, 31(1):915 944, 2021.
Simon Duane, A. D. Kennedy, Brian J. Pendlcton, and Duncan Roweth. Hybrid Monte Carlo. Physics Lellers B, 1987.
Alain Durmus, Eric Moulines, and Ecro Saksman. On the convergence of Hamiltonian Monte Carlo. arXiv: 705.00166, 2017.
Christian Lubich Ernst Hairer, Gerhard Wanner. Geometric Numerical Integration. Springer, 2006. URL https://link.springge www.woooowo.oooo/ooooooooooo.o 1ttp://00..0003-5--00066-------0
Guilherme França, Jeremias Sulam, Daniel Robinson, and René Vidal. Conformal symplectic and relativistic optimization. Advances Neural Information Processing Systems, 33:16916-16926, 2020.
Qiang Fu and Andre Wibisono. Hamiltonian descent algorithms for optimization: Accelerated rates via randomized integration time. In Annual Conference on Newral Information Processing Systems (NewrIPS), 2025.
Walter Greiner. Classical mechanics: systems of particles and Hamiltonian dynamics. Springer, 2003.
Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric numerical integration illustrated by the Störmer-Verlet method. Acta numerica, 12:399-456 2003.
Matthew D Hoffman and Andrew Gelman. The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15:1351-1381, 2014.
Sai Praneeth Karimireddy, Quentin Rebjock, Sebastian Stich, and Martin Jaggi. Error feedback fixes SignSGD and other gradient compression schemes. In International Conference on Machine Learning, pp. 3252 3261. PMLR, 2019.
Hassan K Khalil. Nonlinear systems. Prentice Hall, 2002.
Yunbum Kook, Yin-Tat Lee, Ruoqi Shen, and Santosh Vempala. Sampling with Riemannian Hamiltonian Monte Carlo in constrained space. Advances in Neural Information Processing Systems, 35:31684-311 31666, 2022.
Yin Tat Lee, Ruoqi Shen, and Kevin Tian. Lower bounds on Metropolized sampling methods for well- conditioned distributions. Advances in Neural Information Processing Systems, 34:18812 18824, 2021.
Chris J Maddison, Daniel Paulin, Yee whye Tch, Brendan O'Donoghue, and Arnaud Doucet. Hamiltonian descent methods. arXio preprint arXiv: 000005042, 2018.
Vien V. Mai and Mikael Johansson. Convergence of stochastic gradient method with momentum for non-smooth non-convex optimization. ICML, 2020.
Oren Mangoubi and Aaron Smith. Rapid mixing of Hamiltonian Monte Carlo strongly log-concave distributions. Xix preprint Xrxi 708. 078.07714, 2017.
Oren Mangoubi and Nishecth Vishnoi. Dimensionally tight bounds for second-order Hamiltonian Monte Carlo. Advances in Neural Information Processing Systems, 31, 2018.
14

---

Published in Transactions on Machine Learning Research (02/2026)
Pierre Monmarché. HMC and Langevin united in the unadjusted and convex case. arXiv preprint arXiv:2202.0000 00977, 2022.
Radford M Neal. MCMC using Hamiltonian dynamics. Handbook of markov chain monte corlo, 2(11):2,
2011.
Yurii Nesterov. A method of solving a convex programming problem with convergence rate O(1/k²). In Doklady AN USSR, volume 269, pp. 543 547, 1983.
Maxence Noble, Valentin De Bortoli, and Alain Durmus. Unbiased constrained sampling with self-concordant barrier Hamiltonian Monte Carlo. Advonces in Neural Information Processing Systems, 36, 2024.
Brendan O'Donoghue and Chris Maddison. Hamiltonian descent for composite objectives. Advances in Neural Information Processing Systems, 32, 2019.
B.T. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1 17, 1964.
Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari. Parallel sampling of diffusion models. Advances in Neural Information Processing Systems, 36, 2024.
Weijie Su, Stephen Boyd, and Fammanul J Candes. A differential equation for modeling Nesterov's ac- celerated gradient method: Theory and insights. Journal of Machine Learning Reseorch, 17(153):1 43, 2016.
Andrew R Teel, Jorge Poveda, and Justin Le. First-order optimization algorithms with resets and Hamil- tonian flows. In 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 5838-5843. IEEE, 2019.
Gerald Teschl. Ordinary differential equations and dynamical systems, volume 140. American Mathematical Soc., 2012.
Greg Ver Steeg and Aram Galstyan. Hamiltonian dynamics with non-Newtonian momentum for rapid sampling. Advances Neural Information Processing Systems, 34:11012 11025, 2021.
Jun-Kun Wang. Hamiltonian descent and coordinate Hamiltonian descent. arXiv preprint Xxxv:2402.13985 2024.
Jun-Kun Wang and Jacob D Abernethy. Acceleration through optimistic no-regret dynamics. Advances in Newral Information Processing Systems, 31, 2018.
Jun-Kun Wang and Andre Wibisono. Accelerating Hamiltonian Monte Carlo via Chebyshev integration time. International Conference on Learning Representations (ICLR), 2023.
Andre Wibisono, Ashia C Wilson, and Michael I Jordan. A variational perspective accelerated methods in optimization. Proceedings of the National Academy of Sciences, 113(47))73351-E7358 2016.
Ashia C. Wilson, Michael Jordan, and Benjamin Recht. A Lyapunov analysis of momentum methods in optimization. JMLR, 2021.
A Proof of Lemma 2
Proof. When d(v) = |/ we have VxHl(x, v) Vx(x) and V.ll(x,v) = M Similarly, From Algo-
rithm 3 with the initial condition Xo +- X and Vo +- 0, one can see that the dynamic with a single Leapfrog step is:
8 V fx(x)
Vo vf(zo) =
V
vf(x+)
x Xo
[
Xx
[V{(x)]]
(15)
15

---

Published in Transactions on Machine Learning Research (02/2026)
With only single step, Leapfrog outputs X1, which is the position of the update X6+1 at the next iteration ()()))
A +++1 We hence obtain the update of Normalizod GD, I.c., X = Xx [[x())]
When Φ(v) = |w]], we have Vx v) (x) and one subgradient is (x, v) = sign(v). Similarly, we can see that the dynamic with single Leapfrog step is:
: Vo 0 ((xn)
0 - - [
x, = Xo + sign(V;;) = X2 Asign( (fign(V((xx))
(16)
With only a single step, Leapfrog outputs X1, which is the position of the update X at the next iteration k ++1 1. We hence obtain the update of Sign GD, i.e, X - X - Osign(V/((xp)
When d(v) |ov we have (x, v) (x() and one subgradient, is (x,, v) sign (V(Fmax)) Omaax: where Aaux arg maxie(d] ivfi]l, since ||v|| is the absolute value of the element in with largest absolute value. And m,a is the one-hot vector with index fmax being non-zero. Similarly, we can see that the dynamic
with a single Leapfrog step is:
: = Vo - = - - (xxx)
0
x = Xy + 0 sign(( (v:) Xx - 0 sign (V (xx)) Cimax
(17)
With only a single step, Leapfrog outputs X1, which is the position of the update x++1 at the next iteration k + 1. We hence obtain the update of Coordinate Descent, i.c., Xx+1 = X 0 sign (Vf(x)) [(max) e where max E arg max [x(xx)]
B Additional Experiments
B.1 Ackley function
We also conduct an experiment for minimizing the Ackley function, which is known to be highly non-convex function with two variables:
F(x) 20 exp 1..2 z5)
eep [0.5 2TT cos 2xT22) 20.
(18)
In this experiment, we consider HD with the kinetic energy being d(v) = |
10a
10f
10
10.27-e-100 1.00
10-1
HD-LF dididi.... DE
10-3
10-0n-dw
10
70
LOSS sunccion
10
20
20-2
10-3
10-3
1D-J---2.365
10
10-2
10-3
10"
200
<00
600
ado
<00
nãc
Number Greddere Computations
Number
Goedent Computaitins
Number of Ceacient Computations
(a) Function value F(-) vs. # gradient (b) HDs with different step sizes 0 (in- (c) HDs with different values of integra- tion time (step size 0 fixed).
computations.
tegration time fixed).
Figure 7: Minimizing the Ackley function.
Figure 7 I(a) compare the proposed HDs with the baselines, GD and HB, starting from the same point. For HD-LF and HD-ST, the integration time is fixed 1j 0.01 and step size is chosen 0 0.001. For GD
16

---

Published in Transactions on Machine Learning Research (02/2026)
and HB, the step size 17 is tuned on the grid [10-1.10-3 ..... 10=0] for the best one. From the figure, we can see that. HB outperforms GD, and our proposed HD-LF and HD-ST can further outperforms HB. We also observe that HD-LF behaves more stably than HB during the optimization process, while the update of HD-ST is less stable due to the randomness of the step size Ç..
.10
We then investigate the impact of parameter choices in HD-LF. Figure 7 (b) shows the performance of HD-LF with different step size 0 with a fixed integration time. The integration time 17 is fixed as 1 = 0.01. Then, the number of integrator steps S varies with A. With a smaller step size, the discretized flow tracks the exact Hamiltonian flow better and can converge to lower value. However, the computational burden incrcases with more integrator steps within one execution of the Hamiltonian flow, suggesting the need to accelerate the process of simulating Hamiltonian flow (and we propose HD via Parallel Optimization in Section 8 to address this issue). Figure [1(c) displays the performance HD-LF with different value of the integration time 7,, while the step size is fixed as 0.001. When rj = 0.01 or 0.05, the performance of HD-LF is better than HD-LF with a smaller value than 0.01 or larger value than 0.05, showing that there is a range of 7) obtain the best performance of HD-LF.
B.2 Matrix Sensing
Low rank matrix sensing problem is typical non-convex problem widely applied to image processing, multi-task regression, and metric embeddings (Davis et al. 2025). It aims to recover symmetric positive semidifinite matrix X Rddd with low rank < d from set of linear measurements y XA,,)), where A; E Rdxd are known measurement matrices. Using low rank faclorization of the matrix X= = BB",, where B E Rd we can formulate the optimization problem as:
min BEROXA
1 4m
Σ (y:- (A,,
(19)
where m is the number of measurements. The rank k of matrix B is usually larger than r since we do not know the exact value of r in practice and we choose an overestimated value k. In the experiments A, is generated from A, = a. ā.ā? where a; and af are standard Gaussian vectors.
In this matrix sensing problem, we test our proposed Hamiltonian Descent (HD) algorithm and compare it with Gradient Descent (GD) and Heavy Ball (HB) methods. To simplify, we Fuler discretization and kinetic function Φ(v) v in HD.
Distance to Optimal Solution
Distance to Optimal Solution
10a
GD
HD-0.5
HD-2.0
HD-5.0
HD-8.0
10
10 -1
10-2
10
1000 2000 3000 4000 5000 6000
1000 2000 3000 4000 5000 6000
Number of Gradient Computations
Number of Gradient Computations
(a) Comparison to GD and HB.
(b) HDs with different values of integration time
Figure 8: Distance to the optimal solution in matrix sensing experiments.
In the experiments, the exact rank r is 2 and the overparameterized rank 4. The dimension of matrix d is 100, and number of measurements m 1000. Fig. (a) shows the performance of HD and the comparison to GD and HB. The learning rate of these algorithms are 0.05. For HD, the integration time is 2, with 40 inner steps per iteration. For HB, the momentum is 0.9. We can see the HD can outperform GD and HB
17

---

Published in Transactions on Machine Learning Research (02/2026)
with lower distance to the optimal solution. Fig. g (b) shows HD performance with different integration time values y as 0.5, 2.0, 5.0, 8.0. When the integration time is larger, the final distance to the optimal solution can be smaller, while the convergence is slower.
B.3 Parallel Optimization with Different Kinetic Energy
We also provide experimental results of Algorithm 5 with different forms of the kinetic energy in Figure 9 where we also observe that with the use of multiple GPUs, Hamiltonian Descent via Parallel Picard Iteration reduces the actual running time, compared to the case of a single GPU.
Test Accuracy
Test Accuracy
Test Accuracy
92
70-
GPU
90
2 GGU
Tost Rocuracy (N4)
GPU
68
50-
84
76
74-
82-
GPU
12-
2 GPU GPU
80
GPU
70
20-
15
Running Time (x1000
Running Time 1x1000
Running Time (x1000 1)
(a) ((v) ||v||2-
(b) Φ(v) = |
(c) d(v) = ||v|
Training Loss
Training
Training Loss
1.0
GFU
GPU
GP
0.9
CPU
2 CPU
GPN
GRU
2.0-
0.7-
1.8-
0.64
0.5
0.4
12-
0.3
Running Time (x
Running (x1000
Running Time (x1000
(d) <(v) = ||vl2
(c) d(v) = liv|12
(f) 6(v) = ||v|
Figure 9: Testing Accuracy vs. running time (scconds) on MNIST of HD via Parallel Picard Iteration (Algorithm under the cases of different forms of kinetic energy Φ(-).
18


Reviews:
[{"review_id": "BGLhI8cZkN", "reviewer": "Reviewer_9Mhf", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper extends the Frictionless Hamiltonian Descent algorithm, originally proposed in Wang (2024), which leverages Hamiltonian flows to minimize unconstrained optimization problems $\\min_{x \\in \\mathbb{R}^d} f(x)$. Unlike traditional momentum-based methods that incorporate friction to dissipate energy, HD relies on energy conservation with periodic resets of kinetic energy, ensuring descent without damping. Key advancements include:\n\n- **Discretization via Geometric Numerical Integrators:** The idealized HD requires exact simulation of the Hamiltonian flow, which is infeasible for non-quadratic functions. The authors address this by employing two geometric integrators to approximate the flow:\nThe Leapfrog integrator, a symplectic method that maintains bounded energy errors.\nThe Stratified Monte Carlo integrator, which introduces randomness and is shown to outperform Leapfrog in related sampling contexts.\nThese lead to novel discrete-time optimization dynamics. For instance, with a single integration step, the resulting algorithms recover established methods such as Gradient Descent, Normalized Gradient Descent, Sign Gradient Descent, and Coordinate Descent, depending on the form of the kinetic energy $\\phi(v)$ in the Hamiltonian $H(x, v) = f(x) + \\phi(v)$. With multiple steps, new momentum-like updates emerge, accumulating gradients along the flow without friction-based decay.\n\n- **Generalization of the Hamiltonian:** The authors broaden the framework by considering various non-negative kinetic energy functions $\\phi(v)$, where $\\phi(v) = 0$ only when $v = 0$. This yields a family of algorithms, providing a unified perspective on optimization design and revealing connections to existing methods.\n\n- **Parallel Optimization Technique:** To accelerate the sequential nature of flow simulation, the paper introduces a parallelization method based on Picard Iteration, a fixed-point technique for solving differential equations. This enables simultaneous gradient computations across multiple GPUs via iterative updates, approximating the discretized flow. Theoretically, after sufficient iterations, it exactly matches the sequential output. Empirically, it achieves 2-3x speedups in runtime for multinomial logistic regression on datasets like a8a and MNIST when using 4 GPUs compared to a single GPU.", "claims_and_evidence_yn": "Yes", "claims_explanation": "Yes, the claims are well-supported by a combination of theoretical derivations and empirical evaluations. The theoretical contributions, such as the generalization of the Hamiltonian framework and the descent properties derived from energy conservation, are substantiated through clear lemmas and mathematical proofs that demonstrate the preservation of key properties like non-negative kinetic energy and equivalence to established algorithms under specific conditions. The discretization using geometric integrators (Leapfrog and Stratified Monte Carlo) is rigorously analyzed, with derivations showing how single-step approximations recover methods like Gradient Descent and multi-step versions yield novel momentum dynamics. Empirically, the parallel optimization technique via Picard Iteration is validated through experiments on multinomial logistic regression across datasets such as a8a and MNIST.", "audience_yn": "Yes", "audience_explanation": "The paper introduces a novel extension of Hamiltonian-based optimization, addressing practical challenges in simulating flows for non-quadratic functions through discretization and parallelization, which holds clear relevance for the TMLR community focused on optimization theory and algorithms.", "requested_changes": "My primary concern is the absence of modern baselines and architectures in the experimental evaluation. The comparisons are primarily against standard gradient descent and heavy ball methods. In the context of deep learning and non-convex optimization, the standard baseline is Adam or AdamW (although they are variants of SignGD). Omitting adaptive gradient methods makes it difficult to assess the algorithm's competitiveness in modern ML workflows. Another point is that optimization papers in 2024/2025 are often expected to demonstrate effectiveness on Transformer architectures (e.g., GPT variants) to demonstrate they aren't limited to toy problems or older architectures like ResNet. Since LLM training is more memory-bounded, maintaining multiple steps for Picard Iteration is memory consuming and proposes a sliding window. An LLM experiment would be the perfect stress test to see if this sliding window approach is actually viable when the model states themselves are gigabytes in size.", "broader_impact_concerns": "N/A", "additional_comments": "This paper provides a valuable contribution to the literature and has theoretical novelty. The numerical results are good. Implementation details are provided. Hence I recommend accept if authors can address my concerns.", "extra_scores": {}, "extra_text": {}}, {"review_id": "zi0N39ZMyp", "reviewer": "Reviewer_nmHB", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper addresses the challenge of applying Hamiltonian Descent to general non-quadratic optimization problems, where the Hamiltonian flow lacks a closed-form solution and existing approaches such as Idealized Frictionless Hamiltonian Descent are therefore inapplicable. To overcome this limitation, the authors propose a discretized formulation of the Hamiltonian flow and employ state-of-the-art numerical integrators to approximate the resulting dynamics.\n\nThe first contribution is a systematic study of alternative kinetic energy formulations, including squared- $\\ell_2, \\ell_2, \\ell_1$, and $\\ell_\\infty $ norms. The authors show that when a single integrator step is used per iteration, the resulting updates resemble those of existing descent-based optimization methods, whereas using multiple integrator steps induces dynamics similar to momentum-based algorithms.\n\nThe second contribution is the use of Picard iteration to simulate the Hamiltonian flow, which, after discretization, enables parallel execution of the flow updates and allows efficient implementation on multiple GPUs.\n\nThe proposed approach is evaluated on certain non-convex functions and machine learning datasets. Experimental results report training loss and test accuracy comparisons against standard optimization methods. The results indicate that the proposed Frictionless Hamiltonian Descent with parallelization achieves performance competitive with momentum-based methods and outperforms non-momentum baselines. In addition, the paper evaluates the quantitative speedup obtained by increasing the number of GPUs during training and testing.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The update schemes corresponding to the different choices of kinetic energy are derived through theoretical analysis that is clear and easy to follow. Experimental results are presented across multiple settings, including non-convex benchmark functions and machine learning datasets. These results show improved performance over non-momentum baseline algorithms and competitive performance relative to momentum-based methods when training neural networks. Additionally, the proposed algorithm outperforms the heavy-ball method in the Ackley function and matrix sensing examples.\n\nThe use of parallel Picard iteration within the proposed framework is also clearly demonstrated, although the performance benefits of parallelization are not consistently evident across all datasets. Overall, the experimental results are largely representative of the claims made in the paper.\n\nThere are a few points where additional clarification would strengthen the presentation. First, providing details of the neural network architecture used for the MNIST experiments would improve reproducibility. Second, stating the final form of the objective function used in the matrix sensing example (Eq. 9) would be helpful.", "audience_yn": "Yes", "audience_explanation": "The methodology and findings presented in the paper are of interest and may be relevant to researchers exploring alternative optimization algorithms for training neural networks.", "requested_changes": "While the methodology is technically sound, the paper would benefit from a clearer motivation for the practical usefulness of the proposed approach relative to well-established momentum-based optimization methods. In particular, the manuscript does not sufficiently articulate the scenarios in which the proposed method is expected to offer advantages. Including a dedicated discussion that clarifies the intended use cases, potential benefits, and scenarios where the method may outperform existing momentum-based algorithms would help place the contributions in a broader optimization context.", "broader_impact_concerns": "This work does not appear to raise any significant broader impact or ethical concerns.", "additional_comments": "Section 4.3- In addition, we conduct the experiments $ \\textbf{for} $ different...", "extra_scores": {}, "extra_text": {}}, {"review_id": "KRj2s9qLma", "reviewer": "Reviewer_oHeX", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes an extension to Frictionless Hamiltonian Descent (HD) by introducing geometric numerical integrators—specifically the Leapfrog integrator and the Stratified Monte Carlo integrator—to simulate Hamiltonian flow, replacing the requirement for an exact flow. The authors generalize the framework by introducing various forms of kinetic energy (e.g., squared $\\ell_2$, $\\ell_1$, and $\\ell_\\infty$ norms).\n\nA key theoretical contribution is presented in Lemma 2, demonstrating that executing a single step of the Leapfrog integrator with these specific kinetic energies recovers standard optimization algorithms such as Gradient Descent, Normalized GD, Sign GD, and Coordinate Descent. Furthermore, to address the computational bottleneck of sequential updates in Hamiltonian flow, the authors propose a parallelization technique using Picard Iteration. This allows for simultaneous gradient computations across multiple GPUs. Experiments are conducted on multinomial logistic regression, neural network training (ResNet18 on CIFAR10), and non-convex problems (Ackley function, Matrix Sensing), demonstrating that the proposed method is competitive and that the parallelization yields significant speedups (approx. 2-3x).", "claims_and_evidence_yn": "Yes", "claims_explanation": "1. Theoretical Evidence: The authors' claim that their framework unifies various optimization algorithms is rigorously supported by Lemma 2.\n2. Empirical Evidence for Optimization Performance: The claim that the proposed Frictionless Hamiltonian Descent (HD) is a competitive optimization method is supported by a diverse set of experiments.\n3. Evidence for Parallelization Speedup: The claim that the proposed Parallel Picard Iteration (Algorithm 5) accelerates training is backed by clear, quantitative data. Table 2 and Figures 5-6 demonstrate consistent wall-clock time reductions across multiple datasets.", "audience_yn": "Yes", "audience_explanation": "There are at least two distinct sub-communities within machine learning that would find this paper of interest: The theoretical contributions of this paper are likely to interest those studying the foundations of optimization algorithms. On the other hand, by proposing a \"Parallel Picard Iteration\" method that leverages multiple GPUs to parallelize these steps, the paper offers a concrete engineering solution that achieves measured speedups (2-3x) on standard hardware. This aspect is relevant to researchers focused on hardware efficiency and distributed training.", "requested_changes": "1. In Section 4.1 (Ackley function), the paper explicitly states that the step size for baselines (GD and HB) was tuned on a grid. However, the tuning protocol for the Multinomial Logistic Regression (Section 4.1.1) and Neural Network (Section 4.2) experiments is less explicit regarding the baselines. To ensure the comparison is rigorous, please add a brief description (either in the main text or Appendix) detailing the tuning range and selection criteria used for the baselines in these sections.\n2. Section 2 provides the update equations for Hamiltonian Descent when using multiple Leapfrog steps ($S > 1$) 2. These equations introduce \"novel momentum-like dynamics\" involving weighted accumulations of gradients. To improve readability and theoretical depth, please expand the discussion in this section to provide more intuition on why these specific dynamics are beneficial. For instance, how does the weighting of gradients in Equation 6 differ qualitatively from the accumulation in standard Heavy Ball momentum?\n3. Algorithm 5 introduces a \"sliding window\" approach with window size $W$ and error threshold $\\tau$. Please add a brief clarification (likely in Section 3 or Appendix) regarding how these parameters were selected for the experiments and how sensitive the algorithm's performance or stability is to choices of $W$ and $\\tau$.\n4. Please proofread the manuscript for minor typos. For example:\n- \"runnning time\": the 2nd line from the bottom of section1\n- \"singe GPU\": the 5th line from the bottom of section 4; the 3rd line of section B.3", "broader_impact_concerns": "N/A", "additional_comments": "N/A", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)