Heavy-Ball Momentum Method in Continuous Time and Discretization Error Analysis
Authors: Bochen Lyu
Xiaojing Zhang
Fangyi Zheng
He Wang
Zheng Wang
Zhanxing Zhu
Venue: NeurIPS 2025 poster ()
Keywords: continuous approximation
momentum
HB
Primary area: theory
TL;DR: We establish a continuous time approximation for Heavy-Ball momentum method with explicit discretization error.

Abstract: This paper establishes a continuous time approximation, a piece-wise continuous differential equation, for the discrete Heavy-Ball (HB) momentum method with explicit discretization error. Investigating continuous differential equations has been a promising approach for studying the discrete optimization methods.  Despite the crucial role of momentum in gradient-based optimization methods, the gap between the original dynamics and the continuous time approximations due to the discretization error has not been comprehensively bridged yet. In this work, we study the HB momentum method in continuous time while putting more focus on the discretization error to provide additional theoretical tools to this area. In particular, we design a first-order piece-wise continuous differential equation, where we add a number of counter terms to account for the discretization error explicitly. As a result, we provide a continuous time model for the HB momentum method that allows the control of discretization error to arbitrary order of the learning rate.  As an application, we leverage it to find a new implicit regularization of the directional smoothness and investigate the implicit bias of HB for diagonal linear networks, indicating how our results can be used in deep learning. Our theoretical findings are further supported by numerical experiments.

Heavy-Ball Momentum Method in Continuous Time and Discretization Error Analysis
Bochen Lyu
Xiaojing Zhang
Fangyi Zheng
He Wang"
Zheng Wang"
Zhanxing Zhu"*
University of Southampton
DataCanvas 'Pony.ii University College London
Uniiversity of Leeds
{bochen. lyu, z zhutoooooon.ccc.
Abstract
This paper establishes continuous time approximation, piece- wise continuous differential equation, for the discrete Heavy-Baie (HB) momentum method with explicit discretization error. Investigating continuous differential equations has been promising approach for studying the discrete optimization methods. Despite the crucial role of momentum in gradient-based optimization methods, the gap between the original discrete dynamics and the continuous time approximations due to the discretization error has not been comprehensively bridged yet. In this work, we study the HB momentum method in continuous time while putting more focus the discretization error to provide additional theoretical tools to this area. In particular, we design a first-order piece-wise continuous differential equation, where we add a number of counter terms account for the discretization error explicitly. As result, we provide continuous time model for the HB momentum method that allows the control of discretization error to arbitrary order of the step size. As an application, we leverage it to find a new implicit regularization of the directional smoothness and investigate the implicit bias of HB for diagonal linear networks, indicating how our results can be used in deep learning. Our theoretical findings are further supported by numerical experiments.
Introduction
Gradient descent (GD) and its variants momentum methods, such as Polyak's Heavy-Ball momentum
(HB) (Polyak] 1964), Nesterov's method accelcrated gradients (NAG) (Nesterov) 1983). and Adam (Kingma and Ba 2017), are at the of the success of training deep neural networks. This leads to the importance of understanding the gradient-based optimization methods, whereas the direct analysis for their discrete learning dynamics is challenging. Hence optional approach is generally applied in this area: leveraging tools from dynamical systems to study the learning dynamics in continuous time shed light on the dynamical behaviors of the discrete updates, e.g., BBarrtt and Dherin (002)); y Lyu and Li (2020); Ji and Telgarsky (2019); Li et al. ):22); Lyu and Zhu (2023); Miyagawa (203)) Rosca al. (203)) However, fundamental gap- the discrepancy between the discrete optimizations and their continuous time models named as the discretization error- becomes
manifest when employing this approach.
While recent works have made significant efforts to fill the aforementioned gap for GD (Brrrttt and Dherin 2022 Miyagawa 2023 Rosca et al. 2023), the analysis for the discretization error of the important variants of GD momentum methods in continuous time still far from comprehensive. Along this line of research, Kovachki and Stuart (2020) demonstrated that the continuous time approximation of HB and NAG, which solves ming L(3), can be expressed as rescaled gradient
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

flow (RGF):
B VI(8) 1
where L is the objective function and μ is the momentum factor. Despite its convenience for studying discrete HB and NAG, this approximation has one drawback: the insufficient consideration for discretization error when approaching the continuous time limit. As a result, it is insufficient to differentiate between HB and GD in continuous time. For example, the solutions of Eq. ①I and that of the gradient flow (GF) B = -VL(B), the simplest continuous time model of GD, are not very different: if 3* (t) solves RGF and (t) solves GF, then 38t (1) = (a - p))t), which suggests that GD and HB will converge similar points. However, this is inconsistent with the actual behavior, e.g., in Fig. 2(a) GD and HB converge to different points, which reveals that the RGF cannot fully capture the difference between HB and GD and motivates us to study continuous approximations for HB with smaller discretization error.
Recently, backward error analysis (Hairer et al. 2006) has been successfully applied to construct continuous time approximations for GD (Barrett and Dherin 2022 Miyagawal [0233 Rosca et al., 2023) with lower discretization error than GF. Inspired by these works, Ghosh et al. (2023) provided a continuous differential equation for HB, which achieves smaller discretization error compared to
RGF (Eq. (T) While promising, it still does not fully capture the actual discrete learning dynamics since the discretization error is only to the second order of the step size.
To this end, we also study the HB momentum method in continuous time while putting more focus on closing the gap due the discretization error. In particular, we propose HB Flow (HBF), pice-wise continuous differential equation as a novel continuous time approximation for the HB momentum method. We add number of counter terms to the improved version of Eq. (I to cancel the discretization error, obtaining a continuous time differential equation that can be arbitrarily close to HB, i.e., the discretization error can be controlled to arbitrary order of the step size. Our HBF provides a more reliable foundation for studying momentum methods in continuous time when the direct study of discrete learning dynamics cumbersome. As a case study, we examine the implicit bias of HB, the preference for certain kind of solution without explicit regularization, by employing the HBF for the popular diagonal linear networks where there are already abundant results for GD (Azulay et al. 202 Even et al. 2023 Pesme et al. 2021 Pillaud- Vivien et al. 2020
Woodworth et al. 2020 Yun et al. 2021).
Contributions. We study HB one of the most important gradient- based optimization algorithms with momentum-nn continuous time. In particular,
1. We propose piece-wise continuous differential equation (HBF) by adding counter term (Theorem 2.1 y that can be expanded in a series formulation to different orders of the step size the RGF (Eq. (T), such that the HBF can precisely capture the learning dynamics of the discrete HB.
2. We explicitly show the leading order of the discretization error for the HBF and indicate the way how one can obtain HBF with discretization error to any order of the step size (Section 2.2h. revealing the explicit precision of HBF for approximating the discrete HB. These results are firstly revealed in this work to the best of our knowledge, bridging the gap between the discrete HB and the its precise characterization in continuous time.
3. We leverage our HBF to examine learning dynamics of HB. For example, we reveal that, as HBF implicitly has a regularization term for the directional smoothness, the learning dynamics of HB exhibits smaller directional smoothness compared to GD (. In addition, for the application of HBF in deep learning, we investigate the implicit bias of HB for the diagonal linear network (Thcorem 3.1) through the lens of HBF as a case study, revealing its difference compared to that of GF which cannot be obtained by the RGF.
1.1 Related Works
Continuous approximation Prior works (Shi et al. 2022, Su et al. 2016 Wilson et al. 2016) have constructed second-order ODEs to study the convergence properties of HB/NAG in continuous time, where the momentum factor typically depends on the step size and the iteration count. this paper, we study HB with constant momentum factor that is independent of step size and iteration
2

---

count, a setting that is more consistent with the practical case, e.g., PyTorch Paszke et al. (211)). and we focus more on the order of the discretization error of the continuous time model. Towards this direction, besides Kovachki and Stuart (2020); Ghosh et al. (2023). Cattaneo et al. (2023) focused on the continuous limit of more general adaptive gradient- based methods with momentum, Adam, achieving discretization error to the second order of the step size. As comparison, we provide continuous differential equation that can be arbitrarily close to the discrete HB.
Implicit bias of optimizers The implicit bias of GD for various deep neural networks has been widely studied in, e.g., Ji and Telgarsky (2019); Soudry et al. (2018); Yun et al. (202T) for linear networks and Chizat and Bach (2020): Lyu and [.| (2020) for homogeneous networks. For diagonal linear networks studied this paper, Azulay et al. (202T): Pesme et al. (202T): PPillaud Vivien et al.) (2020); Woodworth et al. (2020) revealed the interesting transition from kernel to rich regime by altering the scale of initialization. Papazov et al. (2024) studied HB using a second-order ordinary differential equation (ODE) with discretization error to the second order of step size. For linearly separable data and linear model, Waan et al. 22022) showed that momentum method converges to the na-max-marghh solution, which is the same as GD. As comparison, Zhang et al. (2024) revealed that Adam with negligible stability constant exhibits the preference of & -max-mmaaggin solution. Furthermore, Wang al. (202T) studied the implicit bias of adaptive optimization algorithms on homogeneous deep neural networks.
1.2 Preliminaries
Notations For a vector 3 E Rd that depends on time t, we use 3 and / to denote its first and second derivative with respect to time i, respectively. We use 3, to denote its j-th component and ||3||p for its lp norm. We use a to denote the inner product and @ to denote elementwise product, e.g.. Rd gdu,ds B Rd
a. A.AAB B denotes a A3 for E 1di A
We use [N] for integers between [0, N].
Heavy-Ball momentum method HB Polyykk 1964)
employs a two-step updating scheme (Sutskever et al., 2013), rather than the single-step manner of GD. Particularly, HB first accumulates the history of past iterations before updating the model parameter 3 E IR i.e.,
m jum;- - 7/ (3xx) Bk+1 8x + m2+1 where μ E (0, 1) is the momentum factor, v) is the step size, & is the iteration number, and m E Rd is the momentum, which can be further written in single equation
Discrete HB: 3x-1 - Bx L(ab (Ba Bk-
2 HB Momentum Methods in Continuous Time and Discretization Error
In this section, we will propose an ODE
B(t = -((3
(3)
that be arbitrarily close the discrete learning dynamics of HB Eq. (3) gap between the discrete learning dynamics and the ODE Eq. given N discretization error of Eq. 3
To characterize the 0, we define the
Wk E [N] & = B(kr) - Bk,
where (kk)) solves Eq. 3 and we will denote ta = kn for convenience in the rest of this paper. An investigation of Euler forward method will identify as the step size. When 6 O (17%), we say that the discretization error of the continuous time model Eq. (3) is to the a-th order of the step size. We aim to find the formulation of ((3) such that Eq. (3) can be arbitrarily close to the discrete HB, i..., E O (17²) ) for any given a.
Intuition of our approach Given the discrete HB Eq. 33). Kovachki and Stuart (2000) showed that & O (7) if ((3) = LL/(1- μ1) as in Eq. (I. To coincide with such observation, our overall goal is to find the formulation of a modified ODE can be controlled arbitrarily low, where its design should follow two basic principles: (i). it should coincide with RGF, the simplest continuous approximation of HB RGF, when v) is very small, thus, G' should degenerate to LL p);; (ii). by adding the counter term 7/, it should allow us to further
- 77 such that the discretization error
3

---

decrease the discretization error of RGF to get better continuous approximations for HB. To achieve this, we need derive the exact forms of G and y under the condition E (17)), ). And there will be two key steps:
1. find the equations that G and y must satisfy if we require the discretization error EY = (17)) for any a > O;
2. solve these equations to give the formulations of and 7.
Overview of our approach We now apply the above intuition to establish such continuous time approximation of HB. Adding the counter term y directly to the RGF is, however, problematic due to the fact that each iteration of momentum methods exploits the history of previous iterations. This renders the local error analysis unreliable since it ignores previous updates. To perform a global analysis, instead of directly utilizing the backward error analysis in BBrrett and Dherin (2022), propose a piece-wise continuous differential equation by decomposing ((3) into two parts, named as HB Flow (HBF),
E [tx, tx+1): - == -
(5)
This was previously discussed in Ghosh et al. (2023). Eq. (5). k denotes the iteration count for the discrete updates, f4 krj, and G5 depends on k and should degenerate to VL/(1 - (4) as in Eq. |I for small step size. Additionally, the counter term "A is designed to cancel the further discretization error brought by GB)) such that Et: can be controlled to be arbitrarily small, hence the name counter term. As discussed in the above intuition, deriving the formulation y: needs set of equations for it to satisfy. We establish such equations as follows. When approximating ty from & >> tx and < ty. respectively, Taylor expansion for 3 in Eq. (5) provides us
B(1x+1) - 3(x)) = m)((4))+ 1 G/ - 11 7k + 17221 8(tx) - ((xx-1) = - -
(6)
G 72---1 11 1
where t and tx mean that we approximate ta from t > to and t < tx. respectively, " plats ((T) - dT, and we apply Eq. (5) for a(t)). A subtraction of the discrete HB update Eq. (2 the first equality of Eq. (6) now allows us to construct the relation between EK and ykk
Ek+1 - E = n(5x 1 - 1 |Gx pG (L(33))
+ v/ [T + +1l 7h + 1 (7)
As - can be expressed by Yk. requiring Ek+1 E 017 in Eq. (7) 77mimmddiatee builds an equation for yk which can be solved to give the form of 76. Hence, the discretization error of the
continuous time model Eq. 5 for the discrete update Eq. 2 can be controlled to the a-th order of
the step size yj. The formulation of Eq. (5) is presented in Theorem 2.1 and the detailed technical
proofs are deferred to Appendix A
2.1 HBF with Discretization Error to Arbitrary Order of the Step Size
Following the intuition and overview of our approach, we now solve Eq. 7 to derive Yo and Gg, and reveal that the obtained HBF by doing so is indeed an (11²)) (piece-wise) continuous time model of
HB. As discussed earlier in Eq. (7). we need to find ?/k to ensure that Ek+1 - & = O (1)+)), hence the following integral functional equation
(1 + pl5 - Tk + w7k-1) = O (17+1)
(8)
must be satisfied. We solve this equation following three steps shown below: (1), write the solution of Eq. (8), Yk. as a series form
7
(9)
(2). derive F explicitly by employing Eq. 9 in the learning dynamics Eq. 5 in the series form
+ Σ (o
(10)
We will show in Appendix that this condition is sufficient for proving ex (17)) by induction.

---

(3), match the terms of Eq. ( with that of Eq. 1O to each order of the step size 1), i.c., Vo E N
7k (0) (o))),,, As result, a simple truncation of " will automatically lead Eq. 7 to give us
Ek+1 Ex = (1att) as desired. Below formalize the aforementioned discussion.
Theorem 2.1 (HBF with the discretization error O (n Let k [N] be the iteration count,
be the step size, and ix = kn, the piece-wise contimuous time differential equation HB Flow (HBF) for the discrete update Eq. (33)
HBF:
- (3) t [tx, tx+1)
(11)
has discretization error that satisfies
Eb+1 Ek - O (1)**)
(12)
and, as result, EK = O (17)) for a 2 1, if G. and The has the following formulations:
G₁ =
(13)
a-2 2
Σ
(0)
(14)
0=0
In particular, for ease of notation, let L (k)o)
=7k
(a-1),
V be a differential operator and
7 -Gki
(-1) (-1)=
Sm, ((01111000000 ..00:....000
Σ
q.
m + 2, E EN
==1
x+2
(0)
(-1)
(f.01)
()0m (0m
ull (j- (j-1,0s) ... (j 10-10m-1 m)) 1) (0m-1)
Xj
L 8--1
mn!
Sm..
then each term of in Eq.
14 can be simply written as
E E 7
Σ E 1=0
(15)
Remark Eq. 12 is equality obtained from solving the integral functional equation Eq. 8 rather than an inequality bound. In addition, Eq. (12) is established to arbitrary order of the step size for any given a > 1. Our approach is different from the previous continuous time model for HB method (Ghosh et al. 2023) where the discretization error is specifically constructed to the second order of the step size.
the formulation of G intuitively resembles the update of momentum the discrete learning dynamics of HB m& 1mm--11 VL. Interestingly, G& can be further simplified as
Remark 2
Eq. 13
1 - k+1 VL
largo VL
G
(16)
P
p
which is exactly the R.H.S of the RGF (Eq. (T). Moreover, Eq. (15) indicates that each iteration of
HB depends on the history of previous iterations as (a) incorporates information of all previous
(o) with S k; and such dependence decays very fast due to the coefficient y.² < 1 for small j. By letting μ = 0, all the dependence on l will disappear and our results can recover those of GD. Interestingly, it is worth mentioning that the difference between HBF and the continuous approximations of GD is closely related to the powers of 711 + p1/(1 - 2)2 as we will show in Section 2.2
Xj
Remark The formulations of 7/. are obtained in a recursive manner, (x()) is obtained from
(o'); with o' <)), hence Theorem 2.1 allows us to always build continuous time models for HB method with smaller discretization errors given those with larger discretization errors, instead of performing the whole set of analysis. For example, to build the HBF with EK = (7²) given HBF
The
with E = 0 (172, we only need to derive ( (since "Yk is truncated to the order 2 1), which
can be obtained by applying 7/
O (17²). In addition, the recursive manner in Theorem 2.T provides a possibility to calculate the involved terms automatically by using
(0) provided by the HBF with
software for symbolic mathematics such SymPy, which could be interesting future work.
5

---

Ft (i())
GD
HB
VL
Kovachki Sturrt 2020
3 -V Pl
Barrett and Dberin
2022
3
and Ghosh et 2023
-
8
VL-
1)
3
V.
Ten [T-
?
Miyagown
2023
Rosea
2023
d' B-pa)
120
of this work
Arbirary
Miyagawa
2023
Rosca ot al.
2025
Theoremn 11 of this work
Discrete
Aott 0u (VM al A6-1)
Table 1: Continuous approximations for GD and HB up to differrnnt orders of discretization error.
2.2 HBF with Discretization Error E₁ = O (17²) and E O (11*
In this section, we derive HBF with discretization error the second and third order of the step size to indicate how our approach works. There are basically three steps for finding HBF with E O (11*) for 1:
1. truncate "h the order a Eq. (14 i.e, 7h = -a=0 (o)
7
(0) o a 2. find all Xj with k by identifying the
2. from the smallest
corresponding Sn with m ( . 2} for each
3. derive the expression of (c) for all s a 2 in recursive manner using the relation 7
7 = Σグ= o (a)
Below we discuss the cases for 2, 3. We also summarize these results in Table | Note that the case for states that HBF is RGF, i.e., 3 VL/(1 n),, which might not fully characterize the difference between momentum methods and vanilla GD.
2.2.1 HBF with E7 0 (m²*
According to Theorem 2.1 there only one term in the series of 7/k i.e., 7k Recall that L L
Gş and there only one element the set Sm=2,0=0- i.c., Sm=2,0=0 obtain for j > 1:
{ = 0,0g = 0)},
00)
L 1..0) (-1) (1) and 7i
μL (1- (1-1.0) (-1)
(17)
Thus, using the definition of L (1.0) derive that
according to Eq.
15 in Theorem 2.1
we can immediately
3
(0) VL L 2(1
7 = 7
K-J [1
+11 + u(1 nº)3].
(18)
1 1=
Interestingly, all the dependence on the iteration count & exists the form of n² then for large iteration count k, the form of 7/ can be largely simplified as 7h se 1t VL.VLL This gives us [1-a)3
HBF with E O (17²):
VL
1+ 1 VL V 1(1 3
(19)
which is consistent with the O (1,² continuous approximation of HB in Ghosh et al. 2023 while our derivation of HBF is in different approach that does not depend inequality bounds. worth to mention that when f 0, HBF recovers the (??) continuous approximation of GD as expected.
2.2.2 HBF with ε = O (n)²)
As shown in Theorem
2.1 in this case, the series of th should be truncated to the order a - = 1,
(0) (I) Since have already derived (0) for HBF with Ek O (n²) in
hence 7
% 77k

---

Section 2.2.1 we only need to find Xx which will give us % We first find the collection of sets
(1)
Sm.a. where m = {2.3] given o = 1 as m can taken as integers between [2, a + 2] according to
Theorem 2.1 Specifically, we have
2 = {(0) = (, 0,00 00 1)};
(20)
S11 {(01 0, 03 0)).
We defer the rest of the detailed calculation to Appendix A and directly present the results for large below:
(1) (1+ 12)²
1+ 10p 2²
7k 4(1 14)5
W1 +
W2
(21)
3(1 + 10)2
where we let
41 (VL. T.). F21,
w2 VL. (VL. v2r.
As result, the HBF with E O (17²) is
VL 1 1
VL V of2 (1 + 11² (1 11)5
1 10 n²2
8 8= 1-M 2 - 3
W1 +
3(1 + 14)2 wg
(22)
Implicit regularization of HBF
According Eq. (1)). HBF with E O (17²) in Section 2.2.1
indicates that momentum induces a stronger implicit gradient regularization (IGR, Barrett and Dhenn TGD where THE is the implicit regularization HB while
(2022)) i.e., 7HB (1 + p1/(1 TGD is that of GD. For HBF with & = O (17²) we can conclude that the difference between HB
MLP on CIFARIO 5k Subset (full-batch)
and GD more complicated since HBF now re- lies more on w2 that primarily depends on L L. Interestingly, the formulation of HBF with a 3 suggests that HB will implicitly impose regular-
Directional Smoothness
15.0
ization effect of directional smoothness, which
12.5
is not the case for GD. In particular, for ze 1,
10.0
the third term of Eq.
22 is close to
7.5
VL (VLLV))),
(23)
5.0
GD
HB
which an approximation for the directional smoothness (Ann et al. 2022
0.0
250
500
750
1000
1250
1500
1750
2000
Epoch
VL(3) (VLB)) VL(B-
Figure 1: Comparison of directional smoothness for HB and GD for MLP CIFAR-10 with full- batch GD and HB (μx = 0.9) with v = 0.1.
VL(3))2
by expanding VT V(()) around The directional smoothness measures the extent of oscillating behavior of optimization algorithms, i.e., the discrepancy between two adjacent iterates. For GD, Ahn et al. 2022) revealed that it exhibits oscillatory behavior such that its directional smoothness, gon would saturate around 2/ /n. As comparison, due the implicit regularization for the directional smoothness of our HBF, HB prefers learning dynamics with smaller directional smoothness HB GGD implying an oscillatory behavior to less extent compared GD. This verified in Fig. where @GD & 2/11 20 while HB < gGD and keeps decreasing. The numerical experimental details can be found in Appendix C
3 Implicit Bias of Momentum Methods through HBF
The HBF proposed Theorem 2.1 provides a reliable mathematical tool for analyzing wide variety of properties of HB. One crucial aspect is its implicit bias in deep learning. To demonstrate the significance of HBF, we characterize the implicit bias of HBF specifically for the two layer diagonal study.
linear network (Wodwooth et al. 2020)

---

The formulation of 2-layer diagonal linear networks A 2-layer diagonal linear network (DLN) with parameter w (W+, W-) where W F. is equivalent a linear predictor
f(; w) x' T (w O w_ w_ OW--,,
(24)
where we the parameterization model a popular proxy model for deep neural networks as it shares many interesting phenomena with more complex architectures, e.g., the transition from kernel to rich regime. In this section, focus on this model with (0) = (0).
+
w (Woodworth al. 2020). This
Leaning task For our task, given a dataset ((z;, y:) with n samples where I; E Rd and Ui E IR, we assume that <, d and consider the regression problem with quadratic loss T( w-) D (o) w w1)2/((n)) We use X E iea to denote the data matrix and y (y1,..., yr)" E I2".
Implicit bias of GF for diagonal linear networks For GF, Azulay et al. (2021) Woodworth et al. (2020) showed that the limit point of w is equivalent to the solution of the constrained optimization problem w(o)) = arg minw A AGF (w; x(0)), s.l. Xw y where (0) W+. (0)w-. (0), the potential function A GF
(w; x()) = j== A (w; (0)), and
GP (w;))()) 1 W arcsinh
Wj 2k,(0))
V /4(k((0))² + w3 2x,
2x, (0)
(25)
Note that x00 controls the transition from rich regime to kernel regime, i.e., A GF [w; k()) - ||
for small k(0) while A AGP (w; + ||w||2 large k(0) (Woodworth et al. 2020
3.1 Implicit Bias of HBF for Diagonal Linear Networks
According to Theorem 2.1 the learning dynamics of the diagonal linear networks f(x; w) can be written
Vw. L
V m
W+
"
(26)
1
1
where we use > E Rd and > E IRd for HBF of w. and HBF of w., respectively, and we use Yij -w to denote its j-th component. Compared to RGF (Eq. (.). Eq. (26) has one extra term that
accounts for the high-order discretization error. The implicit bias of w under the RGF is similar to
that of GF, which, however, is not the case for Eq. (66).
Theorem 3.1 (Implicit bias of HBF for diagonal linear networks). If the dynamics of diagonal linear network f(z; w) = x= w where w = W+ C W+ - w_ 3 w 2w--follde follows HBF defined in Theorem 2.1 and if (co) converges to an interpolation solution, let (t) W+)(0) (0) me; (f) where e,(1) (t) = fg ds (7%
w. (s)/w+ss (+:((s))))))))))))))))>)) + 7
then w(x) satisfies that
w(20) = arg aig A(W;R) s.t. Xw
(27)
w
where A(w;k)= Σj (w 0o; x(;) with
As(w,t;((t)) = AGF((; wwys((1) p,(() = / ds
Ys (s)
To (s)
(28)
W+y
w-y(s)
Comparison with the implicit bias of GF Compared the implicit bias of GF in Eq. 25 there are two differences brought by the high-order correction terms of HBF: (1). the potential function GF (w;x(0) x()) for GF becomes A (w; k(oo)) )) for HBF where x(00) is different from x(0). meaning that HBF induces an effect equivalent to a rescaling of the initialization; (2). A,(w, oo; oo; additionally depends on the product (oc). A similar term will appear in the potential function of GF A GF the initialization no longer satisfies w (0) = w w-(0). In this sense, HBF also brings an effect that is equivalent to breaking the symmetry of the initialization. Theorem 3. also applies to the case for higher-order continuous approximation of GD by setting μ O, suggesting an effect equivalent to the rescaling of the initialization that has been verified in GD (Even et al. 2023). This further reveals the reliability and usefulness of high-order continuous approximations.
8

---

The comprehensive characterization for the HBF requires a detailed investigation for the formulations of y- specifically for the diagonal linear networks, which will be an open problem, while below we focus on O (11²) as an example.
Corollary 3.2 (Implicit bias of HBF for diagonal linear networks with EN = O (172). Under conditions of Theorem 3.1 if we use HBF with Ey O (vp²,, then
(xtXa),
7 + 11) (1 ni)
(1( 1- +
kg(1) (0) exp
(29)
where (t) 4 So dww L)2 >,>>, 0, q E Rd with = w/ 4m - (0.
For the exponent of Kyy (oo), when $ (x0) 0 dominates, e.g., (0 (co), we will conclude that (o)) < (O), hence the rescaling effect brought by HBF equivalently reduces the initialization. Therefore, compared to A AGF A(w; x) will closer to the nnoom and the solution w(oo) will enjoy better sparsity. This finding is consistent with parts of results in Papazov et al. 2024), which analyzed the implicit bias of HB also with a continuous time differential equation that is second-ord -odder ODE,
while HBF is first-order ODE and can also cover the case for GD simply by letting μ 0.
4 Numerical Experiments
Discretization error of HBF
15
HEF
2.5
|[cxil
10.2
2.0
RGF
HBr.0=
10
HBF.0
0.0
0.5
100 1.5
2.0
23
30
10
12
19'
Iteration
(a)
(b)
Figure 2: (a). Trajectories for learning dynamics of GD, HB, RGF, and HBF with discretization error O (172) and O (n²) in 2-d model. All dynamics start from the same point (a1 2.8, a2 3.5). The convergence point of HB is denoted as black star. The black dotted line denotes the set of all global minima. (b). Discretization errors for different continuous approximations HB during training in (a).
In this section, we show numerical experiments on simple 2-d model to verify our theoretical claims, and we present numerical experiments details and more experiments for diagonal linear networks in Appendix C
Our simple 2-d model has the formulation f(z; a1, a2) aj where 01 02 E are the model parameters and x, E IR is the training data. The loss function is L (f(x; a1, 02) ) y)² All are global minima. To show that higher-order HBFs with discretization eтror O (17*) are better approximations for HB, we visualize trajectories for different learning dynamics, i.e., GD, HB, RGF, HBF with a 2, and HBF with a 3, in Fig. 2(a) The trajectory of HBF with a 3 is closer to that of HB than both RGF and HBF with = 2. Furthermore, Fig. ((a) also reveals that RGF is more similar to GD and it cannot capture the discrete learning dynamics of HB well. We also plot the norm of discretization errors || for these continuous approximations during training in Fig. 2bb where HBF with a = has the lowest discretization error after several steps. These results validate the reliability of HBF proxy of HB.
parameters a1 a2 satisfying d d2T

---

5 Conclusion
In this paper, we have established new continuous time model for the discrete HB method (Eq. 2 namely HBF, with explicit discretization error that can be controlled to arbitrary order of the step size 1). In particular, our approach constructs relation, which is a functional integral equation, between discretization errors of adjacent iterates for any step and can be solved to arbitrary order of the step size. This is a different approach compared to prior works (Kovachki and Stuart, 2020 Ghosh et al. ) Our results provide a reliable foundation for analyzing the momentum methods in the continuous time limit. We leverage our HBF to shed lights on a newly observed implicit regularization effect of the HB method: the preference for small directional smoothness compared to GD. In addition, as another interesting application of our HBF in deep learning, we study the implicit bias of HBF for the popular proxy model diagonal linear networks, and reveal the difference between the implicit bias of HB and that of GD which cannot be captured by RGF.
Limitation and future directions The framework in this paper does not consider optimization methods with adaptive learning rate, e.g., Adam (Kingma and Ba 2017). A generalization of our framework to such case would be an interesting future direction. In addition, our analysis can be generalized to the stochastic case by replacing V L with VL. the approximate stochastic gradient, by following Li et al. (2018); 8); |atz (202T). Finally, we only consider the simple diagonal linear networks in this paper, and future works can explore more complex deep learning models with our HBF to study the implicit bias of HB.
Acknowledgments and Disclosure of Funding
B.L. is funded by a studentship provided by the School of Electronics and Computer Science, Univer- sity of Southampton. The authors acknowledge DataCanvas AlayaNeW for providing computational resources. The authors thank the insightful and constructive feedback from the anonymous reviewers.
References
Kwangjun Ahn, Jingzhao Zhang, and Suvrit Sra. Understanding the unstable convergence of gradient descent, 2022.
org/ org/
2204. 01050
Shahar Azulay, Edward Moroshko, Mor Shpigel Nacson, Blake E Woodworth, Nathan Srebro, Amir Globerson, and Daniel Soudry. On the implicit bias of initialization shape: Beyond infinitesimal mirror descent. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, puges
PMLR, 18-24 Jul 2021. URL
mlr www.mm/..........mm press/vv3//
azulay21a.htmr
David G. T. Barrett and Benoit Dherin. Implicit gradient regularization. In International Conference Machine Learning, 2022.
Matias D. Cattaneo, Jason M. Klusowski, and Boris Shigida. On the implicit bias of adam, 2023.
Lenaic Chizat and Francis Bach. Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss. In Conference on Learning Theory, 2020.
Mathieu Even. Scott Pesme, Suriya Gunasekar, and Nicolas Flammarion. (s)gd over diagonal linear networks: Implicit regularisation, large stepsizes and edge of stability, 2023.
Avrajit Ghosh, He Lyu, Xitong Zhang, and Rongrong Wang. Implicit regularization in heavy-ball
momentum accelerated stochastic gradient descent, 2023.
Emstt Hairer, Gerhard Wanner, and Christian Lubich. Geometric Numerical Integration: Structure- Preserving Algorithms for Ordinary Differential Equations. Springer Berlin, Heidelberg, 2006.
Ziwei and Matus Telgarsky. Gradient descent aligns the layers of deep linear networks. In International Conference Learning Representations, 2019.
2019.
10

---

Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In aar.ii:: 2017.
Nikola B. Kovachki and Andrew M. Staart, Continuous time analysis momentum methods. In Journal of Machine Learning Research, 2020.
Alex Krizhevsky, Geoffrey Hinton, al. Learning multiple layers of features from tiny images, 2009.
Jonas Latz. Analysis of stochastic gradient descent in continuous time. Statistics and Computing, 31(4). May 2021. ISSN 1573-1375, doi: 10.1007/s1 1ttp:/s1122222221111 0rg/10.1007/ 10111000007/111222-112 022-122222-00---------
URL http:///t..x... doi
Guillaume Leclerc and Aleksander Madry. The two regimes of deep network training, 2020. URL
10376
Qianxiao Li, Cheng Tai, and Weinan E. Stochastic modified equations and dynamics of stochastic gradient algorithms i: Mathematical foundations, 2018. URL httpps://axiii.orrtt wtwp:/vww.r.g/////////aaa 01558
Zhiyuan Li, Tianhao Wang, and Sanjeev Arora. What happens after SGD reaches zero loss? -a
mathematical framework. In International Conference on Learning Representations, 2022. URL
httpps://pperrevier httpeeeviiiieeeeefeeeeeeeeeeee
Bochen Lyu and Zhanxing Zhu. Implicit bias of (stochastic) gradient descent for rank-1 linear neural network. In Thirty-sevnnth Conference on Neural Information Processing Systems, 2023.
Kaifeng Lyu and Jian Li. Gradient descent maximizes the margin of homogeneous neural networks. In International Conference on Learning Representations, 2020.
Taiki Miyagawa. Toward equation of motion for deep neural networks: Continuous-time gradient descent and discretization error analysis. In 36th Conference Neural Information Processing Systems (NeurIPS 2022), 2023.
Yurii Nesterov. A method of solving convex programming problem with convergence rate o(1/k?). In Soviet Mathematics Doklady, 1983.
Hristo Papazov, Scott Pesme, and Nicolas Flammarion. Leveraging continuous time to understand momentum when training diagonal linear networks. In International Conference on Artificial Intelligence and Statistics, 2024.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Ed- ward Z. Yang. Zach Detito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, pirf-performmnce deep 01703
learning library. CoRR, Abs//912.00173. 2019. URL http://aaaa
Scott Pesme, Loucas Pillaud-Vivve Vivien, and Nicolas Flammarion. Implicit bias of SGD for diagonal linear networks: a provable benefit of stochasticity. A. Beygelzimer, Y. Dauphin, Liang, and J. Wortman Vaughan, editors. Advances in Neural Information Processing Systems, 2021. URL https://openrevier net www H
Loucas Pillaud- Julien Reygner, and Nicolas Flammarion. Label noise (stochastic) gradient descent implicitly solves the lasso for quadratic parametrisation. In Conference Learning Theory, 2020.
Boris Polyak.
Some methods of speeding up the convergence of iteration methods. Ussr Computa- tional Mathematics and Mathematical Physics, 1964.
Mihaela Rosca, Yan Wu, Chongli Qin, and Benoit Dherin. On continuous time model of gradient descent dynamics and instability in deep learning. Transactions on Machine Learning Research, 2023. ISSN 2835-8856. URL httpt:tt://///t /ttp::////peerevieeee wwt.:/nnett.oooommeeeeeetttt
Bin Shi, Simon S. Du, Michael I. Jordan, and Weijie J. Su. Understanding the acceleration phe- nomenon via high-resolution differential equations. Mathematical Programming, 195, 2022.

---

Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro. The implicit bias of gradient descent on spearable data. Journal of Machine Learning Research, 2018.
Weijie Su, Stephen Boyd, and Emmanuel J. Candès. differential equation for modeling nesterov's accelerated gradient method: Theory and insights. Journal of Machine Learning Research, 17 (153):1-43 43, 2016. URL http::/ jmlr. wwtppappes/vv///22222 htm.
Ilya Sutskever, James Martens, George Dahl, and Geoffrey Hinton. and momentum in deep learning. In Sanjoy Dasgupta and David McAllester, editors, Proceedings of the 30th International Conference on Machine Learning, volume 28 of Proceedings of Machine Learning Research, pages 1139-1147, Atlanta, Georgia, USA, 17-19 Jun 2013. PMLR. URL https::/procceedngss mlr. wtt.rpereas/////ussseeeeeereee htmm
On the importance initialization
Bohan Wang, Qi Meng, Wei Chen, and Tin-Yan Liu. The implicit bias for adaptive optimization algorithms on homogeneous neural networks, 2021. URL /ttps://arxxv... wrg.g/r.iiio.........ooo..o abs//bs/11 06244
Bohan Wang, Qi Meng, Huishuai Zhang, Ruoyu Sun, Wei Chen, Zhi-Ming Ma. and Tie- Yan Liu. Does momentum change the implicit regularization on separable data? In Alice Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL jttpss/ts///opeerrre
net/forumiid-i-Be -uquurrji
Ashia C. Wilson, Benjamin Recht, and Michael Jordan. A Iyapunov analysis of momentum methods in optimization. In arXiv: 2016.
Blake Woodworth, Suriya Gunasekar, Jason D Lee, Edward Moroshko, Pedro Savarese, Itay Golan, Daniel Soudry, and Nathan Srebro. Kernel and rich regimes in overparametrized models. In Conference on Learning Theory, 2020.
Chulhee Yun, Shankar Krishnan, and Hossein Mobahi. unifying view on implicit bias in training linear neural networks. In International Conference on Learning Representations, 2021.
Chenyang Zhang, Difan Zou, and Yuan Cao. The implicit bias of adam separable data, 2024. URL butps://arxts org/ abs/ 2406. 10650
12

---

Appendix
In Appendix A
we provide proof' for Section 2
In Appendix B
we present proof' of Section
In Appendix we show details of numerical experiments in Section 2.2.2 and and present more related numerical experiments support our theoretical claims.
A Proofs for Section
We prove Theorem 2.1 in Appendix A.2 and give the details for the first several orders HBF in A.3 We first discuss the conditions that guarantee HBF effective approximation of HB.
A.1 Conditions of the Effectiveness of HBF
To make the HBF a valid continuous approximation of HB, there are two necessary conditions:
1. It is crucial control the ratio between y and i avoid r) > 1 n., which might lead the Taylor series to diverge. More interestingly, we conjecture that it is the magnitude of special composite quantity
71
(30)
(1- 1)²
that matters for the effectiveness the HBF. This quantity spontaneously appears both HBF with = 2 and 3 but not in the RGF, i.e., for HBF with the counter term is proportional to v while for HBF with a = a new counter term proportional to e will appear. If v is too large, then results would longer hold.
Hence, we need to fix the value of and treat only 1) as the variable to denote the higher- order terms (17²) while hide y in the expansion. And it would be interesting for future works to study the case when both μ and v) are treated as variables such that higher-order terms are denoted (1)aa) for 1. In addition, the regime of large μ and large 7, the model might not be trained properly either: the update direction coming from the gradient and that from the momentum will jointly affect the training direction significantly, while these two directions can be very different due to the large value of μu and hence cannot give a consistent updating direction.
In addition, the dependence of HBF on the special composite quantity &b is consistent with the empirical observation in Leclere and Madry (2020), where the optimization curves for different momentum values can be recovered by corresponding change in the learning rate. The dependence of HB on n and μ at the same time further indicates the advantage of HBF with > and that RGF, which only depends on μ., is not sufficient to reflect the optimization properties of HB.
2. Given the continuous approximations include derivatives of L up to the a-th order, hence should least be a-times continuously differentiable and || [| should be upper bounded.
A.2 Proof of Theorem 2.1
Given the HBF for l E [N]
& (441)
/=
m)x(3)
(31)
with unknown G3 and 7k. we expect that the counter term 7/ could cancel higher-order discretization errors and G2 should degenerate to rescaled gradient, i.c., VL/(1 μ). Hence, 7 and G& should be designed in such way that 3(4x) is close to 3g in the sense that the discretization enror
Eg B(tx) Bk
(32)
13

---

is small. Inspired by Miyagawa (2023), we first present the outline of our three main steps for deriving their formulations below:
Step I
Unknown Gk. 7k
determine
Taylor expansion residual 1
Taylor Expansion of ((ta11) 35
(36)
constructs
Step
Expression of Ex+1 - EK
+ required to be 11 for EK = O17)²) (Lemmal A.2
Equalities for Gk.
solved by matching to each order of
Step III
Solution of Gg and "
We now discuss the detailed proof following this outline.
Proof. We start with Step which deals with how F is determined by Ge and 7k-
Step I Recall that the discrete learning dynamics of HB is
B 3x VB(b) + p(Bx Bx-1),
(33)
where y is the momentum factor and lt is the iteration number. Based our discussion in Section the continuous differential equation for HB is
2
B Gn(B) (3)
(34)
fort E [ta, 1x++) where tu kn and the solution is B(t). For arbitrary unknown "ik in Eq. (34) + is determined, which is also unknown depends on ?/k- Specifically, the first-order Taylor expansion with the remainder term in the integral form gives us
((t++1) - 3(tx) = n3(t + 12 = G 12m0 + +,223
(35)
where i means we approximate 4 from 1. L..
r - $ 3 (n(k + s)) (1 s)ds,
and Eq. (34) in the second equality. Similarly, when approximating t, from t tg. we obtain that
B(4x) - B(xx-))) = - - 1
12
(36)
Now we construct the dependence 1 on 7 in the series form, as shown in Lemma. A.1 (proof can be found in Appendix A.2.I).
Lemma A.1. Given the series form
0C m
Σ
ry
(a)
(37)
o=0
and the continuous time differential equation Eq. 34 series forms:
in Eq.
35
and Eq.
(36) have the following
+ Σ
((
Σ
(-1)
I
(k.o.)
L (R.o, 1)1(00
(a.-1)
(38)
q!
0=
p=0 Ed
op=
I
Σ
Tp (()))
ΣΣ
Σ
q!
1
(k- (k-1,00))
I
(k-1og
1)_
(0,00)1)))) (0g-1)
7/ (39)
p=0g=2 E /ー=ーー ++
14

---

Step II Given the dependence of F on 7k we are now able to write the Taylor expansion of ((t explicitly. This allows us to construct relation between Ek+1 and Eg by subtracting Eq. (33) from both sides of Eq. (35):
ck+1
-1 (G+(B(ta)) VL(3)))() - y + I
1(5x Bk-1))
(40)
Note that
3K - = =((x)) - B(tx - (Ex - =k-1)
-Tf 7/ 7/-1 1/ (ex Ek-1);
we obtain the expression of Ex+1 E as in the Step Il of our outline:
Ek+1 Es = pl(x 1 17 G PG - VI((((+) ex)
+ [ t V + pYk-1]
(41)
Eq. (41 builds the connection between the discretization error ε. and the counter term 7h in the continuous time differential equation Eq. (5). We can now construct the equalities for Gg and 7k under the constraint of low discretization error by following the lemma below (proof can be found in Appendix(A.22.
Lemma A.2. For the continuous differential equation Eq. Eq. 33 if
34 and the discrete sequence given by
Gx(B(tx)) (3Gt-((((tt) + VL((tt))
(42)
with G-1 0 and
It + pll Y + 1Mk-1 0 (17a
(43)
as in Eq.
49
then we have
Ek Ek- 0 (1)+1)
and, result,
8k = 0
As shown in Lemma. A.2 to ensure that the leading order of the L.H.S of Eq. 41 is to the order of 1. we only need to require ER+11 E
O (17*) which can be guaranteed by Eq. (42) and the
functional integral equation Eq. (43). As Eq. (42 can be easily solved by induction, we only need to
solve Eq. (43). To achieve this, we build a stronger functional equation below:
It + p
7 17
(44)
which is the final equation that we aim to solve to derive -
Step III The Step III of our outline is solving the functional integral equation Eq. (44). Our
idea is simple: as the series form of I has already derived in Lemma. A.1 we make Eq. 44
satisfied by matching both sides of it for each order of 17. In particular, given the series forms of T in Lemma A.l and that of
7% (Eq. (37 we require
E T( (p) 12 == = (19))00 (1)00++($55) (1++)(+ (())))00))
(45)
which, let
(0) (g)) + 1(
(46)
gives the recursive relation of %/ for EN
(7)
Th
(a)
(0) (a) = wo + X
(47)
(o') (0) only depends on % with o' < a according to Lemma A.1 because X
Now given E z we
can truncate 7/ to the order a - such that
a-2 2
(a)
7
"
Tk
(48)
0=0
15

---

then the functional integral equation is solved to the (a 2)-th order of the step size
a-2 (a)
1 + ulg 7h
1
1
Xk
Σn 10 (11
w
(p)
0=0
a=0
(49)
Σ
= 0 (1p-1))
0=a
which is exactly the condition Lemma.
A.2 Therefore, by constructing * following Eq. 47
and truncating "/k preserver its first terms, we can prove that the discretization error of the continuous time differential equation Eq. (34 is to the order O (17*)
A.2.1 Proof for Lemma A.1
Proof. We first rewrite 1 as follows:
1 1 bm
((T) = d)
bn
(kr))
(tn 7')dd'
J0
n!
du
n
(bb))" n d"
3(t²)
(50)
(n + 2) den
0=
where we use (n - T')dT'= need the expression of d' 3/dt" and we start with t :
12
12
12
in the last equality, To continue, we
ntI n+2
(x+1)(m+2)
d"
dt ( ( B((x))
d
(d"-1
3((
dr
dtn
(da-1
b(tt) V
B(tt))
dln
dd=-1
= -- + mik) V
B((+)
dim
= (-1)" (1))m
n-1 (G& + 177k)
(51)
where we denote the differential operator L # Now suppose that Y' can be written as a series
= (Gx + 7 and use Eq. 334 in the third equality.
7
Σ
n
7k 7
(0) (-1)
Gk:
0=0
then Eq.
5I becomes
den 8(1) B(1)= = (-1)"
Σ
(01-1)
01=0
In- p
(an (0n-1-1) v
Σ
y g (1n
m==1=0
on-1=0
= (-1)"
Σラ a Lg
(k(01)
L
(ki
(R.on=)) (0n-1)
(52)
0 00000=00
16

---

Combined with Eq.
50
we obtain the form of 1 for E (tg,tb+1) as
1f
d
1
=
2)
dt (tt)
==0
20
(-1)
n+ En Lg L
Lg ( (0n42-1)
(n + 2)! 11
n=0 01mmp04200
(-1)at2
*
(k:01) 1(,0))
L (k.onti)
(0n+2-1)
(n + 2)!
2)!
R= m=0 La=1 By m
#+2
Σ
(-1)²
q! L
(k.01)
L 11 (0.-1)
(53)
P= 21=1 -
where we let - n +
12
j= Jj. f n + 2, in the last equality. Similarly, when L. we have
dm (())) (-1) (LS 1 n- (Gk-1 77k-1)
which implics that
p+2
- --ΣΣ
Σ
9!" 71 1k 1,-1,0)) ... ... L kk-1
L
(k-1,00-1) (00 (00-1)
(54)
jj
3
A.2.2
Proof for Lemma A.2
Proof. We first present several useful relations. As Eq. 49 is established by solving the functional equation for any iteration count k we can write the relation between Ek+1 and Ey. Eq. (41) as
Ek+1 - Eb μ(Ek Ek-1) (VL()(())) (L((((4) ew)) + O (1a+)).
(55)
If E O(7)), then Eq. (55) implies
||e411 Exll S E Ek-1 (F((tt)))) VL (3(te) -x) |l c
(56)
S μl|ER - E=- + q>|l=>/ + 01
for some constant (j where we use G - pG = VL and let X = max (2I(3)]
Denoting
2
2c1 k
Vk: : (k)
cs(k)
f
now prove by induction.
For the first step (A = 0), by definition we have
Eo = S c2(0))²
Note that G_1 0 and 7-1 0 by definition, then we have
2c1
||sı sall S c17jat1
a0))
1- /
since < 1.
2. Suppose that for the k-th step the following relations hold:
| |l s e2 (
(57)
|| - Exll < ca (k))))
17

---

Then for the (A: + 1)-th step, we have
||+++||| = ||et11 Ek + Ek
$ ||Ex11 Ekll + |
(k)
S c2(h)
1 (2(k)) 17a
S 02 (b)ee 17
= c2 + 1)m²
where the last inequality is because e"' > + x for x > 0. Similarly,
|| Ek+1ll S MllEx41 - = + m>| k + c1n
< (k0g(k) + Ae2 (k + 1) + 1))
24n Me
--
1 -
+
1 1-
m (k+1)
cg(k 1)n
2
2
1
1 1-
1
+
1- H
ca(k 1)matt
2
2
ca(k + 1)mat1.
A.3 O (17)) HBF for specific
In this section, we derive the form of O (ai)-close HBF for given a specific a. There are basically three steps to find HBF that is O (17*t) to HB:
1. truncate Yk to the desired order i.e, Tk = Eo= (o)
7
2. from the smallest a, find all (a) with j S k by finding the corresponding Sm,a with (2..... 2} for each o;
(0) 3. derive the expression of 7k for all a 2 in recursive manner using the relation (a)
=
In the following, we give the cases for a = 2 and 3 as examples. With this approach, one can in fact find HBF with arbitrary order of closeness to HB.
A.3.1 = 2.
According Theorem 2.1 the series of 7 is truncated to the first term, i.e., 7k 77/7// where
(0)
T = Ee=o - X2 Thus the first step is to find x)² which can be given by first identifying the
(0)
set S:
Sm=2,0=0 = ((01 = 0, = 0)}
(58)
(0)
therefore there is only one term in X,
Xj (0)
1 (-10)) (-1)
+ n 70, 7 j-1
1.0 (-1)
Recall that
2- Gj - 1i+1 VL, r
(59)
which, according our definition in Theorem 2.1 leads to
1.0 L 3
7
(-1)
v G5 V.
18

---

we obtain that
1 2
x3
G. VG, + #G VGj-1]
1
(1 2 + Tpl 1²) VL L.L.
(60)
2(1 11/2
Thus
1 k
(0)
Σ
1 -
[G, VG, + 1, -G--I VG
7k
2
1=0
VL V
11 x [11 222 pll - p2))
2(1 /1) j=0
VL V
k
[(1 + p) + + (1 + 1)) - 4]]
(61)
2(1 12)
j=0
When & is larege, the above expression can be simplified as
(1+ Σ 12
1+ 1 2(1 - pt)
VL L2L
VL V2L.
211-w)² 22²
A.3.2 a = 3.
Similarly, in this case we first truncate the series of Yt to the desired order, i.e., 7 = 7k + 77k where we have already obtained 7/ in the last section, thus we only need to find 7 and X which can be done by first finding the set 2m=2,0=1 and Sm=3,0=1:
(1)
1 {(0) 02 0), (01 0, 1)},
Sa. = ((0) = 0, = 0, = 0))-
Therefore there are three terms of X;
(1).
x (1)
1
L
1,1 (-1)
+ 2
1
LB 75
(0)
+
1-10 (0)
2
(-1)
11.100 (-1) 3/
(62)
6
Recall that 7/3
(111) Gj, 1 = GG V. and L #
= ) 50)
V, the first line of Eq. 622) is
2
1 [ 00)
?
VG, 17, 00)
VG,-1 + Gj V (0) + pGJ-1 V
(63)
while the second line
G, V VG,) - pCy-1- VG,-1)
(64)
To simplify these terms, we can either replace all 73 (0) with the expression in Eq. 6 and write
Gj explicitly, or notice the recursive relation between G, and Gý-1 1 in Theorem 2. i.e. ,G
jGG-1 1 VL., then Eq. (33 becomes
2 ?
(0) 2 +
VL
5 (0)
+
1 (1 + +,, - V (2 (0) + 21)]
and Eq. 664 is now
VL VL. VG,)
| Gj-1 VG VG, Gy-1 VG,-1).
(65)
Summing over these terms gives us x;
(1)
x = 1 + p
(66)
19

---

where
// 1 (0) V + VL v),-0) 1 VL V(G VG;)
2 (/
1
(
(0) + 00
VGj-1 + G,-1 (2. + +1/4,))
2
G--1 V(G VG, GJ-1 VG
We now find (1) 7
through its definition:
k
h (1)
& + 11
Σ
Σ
(67)
j=0
1=0
j=0
In the following, derive the form of
(1)
When & is large. According to Eq.
61
have
VL- L
k A- (0)
(1 + + (a (1+11)++))) 11)-44)))
2(1 14) i=0
VL V
(1+ 11) 1(11)
1 (1 + (1 (1-t)) 1- 11
4(j + 1)
2(1 11-11 11)2
1-
VL L
(1 + (H
k- ~- - -
1 (1 + 1n((1 - -mettt)
44j + 1)
2(1 -1)²
1
1 1
VL L
(1+ p))
a8t5+1(1+ 11)
4(j 11)
2(1 - 11)3
1 1
1-
(1 +
VI. 1
(68)
2(1 11)3
and, according to Eq. (59)
1 (1 2 + ((++1)))
GG VG,
VL VLL= (1- VL V
(69)
(1 n1)
Combining Eq.
68 and 699
gives the form of ( when k is large:
+ 1))
[ V2L) V'L+VL-V v-L)]
4(1 1
ドー
VL V (VL ALL)
6(1 n2)
which immediately leads to
v (1)
1=
(1 u))
77 V (VL- V²L)
[VLL - V2L) - VL V (VL. A-D)]
4(1 -1)1
6(1 4)3
(1+1)VVZ V'L) VLL+ + (1
5p)
VL v (VL. voL)
(70)
4(1 nr
3
The left part is now deriving the form of , ()) (1) which can be done by first finding
1 (0) VG
(1 + μ1)
( V2L) VG,-1
2(1 p)3
2
(1 11+1 μ)
(VL V²L) VLL 10
1 VG,-1
(0)
(71)
2(1 p
20

---

and
2 ドJ
1 M G VG, VG, - G--1-
1L V (VL.
(VL.V2)),
(72)
(1 nn
thus
(1 (1+11) 14)
Σ 1
[VV V2 vLL V
VL V VL]]-
(73)
2(1 1.)s
1=
Combing this equation with Eq.
(70.. we can now conclude the form of 7/ when k: large:
(1)
(1)
k
1 (1)
)
7k
1=0
(1 (1+H)(VL (V PD VLL+
(1+ 5p)
VL (VL VoL)
4(1 44
3
+ m11 + 11) +11+
(VL L) VLL VL V2
2(1 11)5
(1 + 22
VV L) V 1 10p + 2
VL V (VL FLL)
(74)
4(1 - 4²)
3(1 + 112
Note that when μ = 0 we recover the result of GD, i.e.
(1) 7
(VLV'L)--V' VLA(VVAV))
12
B Proofs for Section 3
Given data (ma, y:), the architecture of 2-llayer diagonal linear network is
d
f(x; = x (w+ o ww+ w O w-
Tij (w w2
(75)
j=1
and the empirical loss function is
I(w)
1 2n
Σ (s(s::) y1)
We let (r1, (ri,.... TT R" be the residual where Vi = f(x;; w) y4- According to Theorem 2.1 the HBF learning dynamics of model parameters W+ and w_ will be
V L
V L
Th
m
(76)
W d and % where we use % respectively, and the gradients are
E IRd to represent the error terms for HBF of w4 and w
X ww, = 2w+ O L. V L= 2w 2 V.
V w -x
(77)
71
(78)
Using the expressions above, can be easily verified that
w Vwr I. + w. @ V=_
(79)
and we will frequently use this relation later. Recall the definition of K. W W-gj we now present useful lemmas before proving Theorem 3.1
Lemma B.I. Let Kz = W+:s(1) (l) -1), Vkej denote the j-th component of 7 wt and
(()) f ds
W Tk:s
(s)
Yku (8)
(80)
(s)
w
w
then we have
kt(t) (0)())))))))
70(())
(81)
21

---

Proof. The proof applies the dynamics of w and that of w-:
dx
W
① W- + W-
3 W+
dt
V +
V
17k
W+
W
177
w-
1 P
p
W+
W+
%
(82)
where use Eq. (26) in the second equality and Eq. j-th component of K, we have
(79)
in the third equality. As result, for the
W
Tkyj
Kg = -1pK
W+i
W
k((() (0))0))))))))
(83)
is also interesting to investigate the dynamics of w as shown below.
Lemma B.2. If W+ is run with HBF, then the dynamics of w satisfies that
V
w 4v
(84)
where we les
W+ W+ +W
IE 2 (TH w W+
7 C w w-).
(85)
Proof. Using the dynamics of W Eq. 26)..
we can show that
w = 2w
② 2w
① w
Vw_L
= 2 1 1
Mi
W+
O W+
-2 2
w w
1 1- 4
1
-4 (w O W + W
W
2n
(ik
W
7 o w-).
(86)
L
To show the implicit bias of HBF, we need to first explore the dynamics of w, which is present in the following lemma.
Lemma B.3 (Dynamics of w for diagonal linear networks under HBF). Under conditions of Theo-
rem 3.1
if the diagonal linear network f (a; w) is trained with HBF (Theorem 2.1 let
(W))) (L)
w, (1(
wt(()
w3(()
aresinh
+
(1)(()
(1y(1)
4x3(1) (1)
f ds
Tk: (s)
Tki (s)
()(1)
4
(+)(()) W
W
w_)
1,(w G (w;x(1) + (1),,
(87)
then the learning dynamics of the parameter w satisfies that
0w, Lb
- Ow
+
0.
(88)
di
1-
The proof of this lemma can be found in Appendix of Theorem 3.1
B.2 In the following we first focus on the proof
22

---

B.1 Proof of Theorem 3.1
Now we can prove Theorem 3.1 with above helper lemmas.
Proof. Recall the definition of Ay in Lemma B.3 and we further define
A(w, tw,
(89)
then Lemma B.3 gives
d v (w,(w,;;))=
d dt
A1 (i),) Owa Aaw,,(W
W
dt
X²1
n11 - 11)
= ow 0o; - Ww-((((0),0;s) = 2e. So r((T)dr
Σ a4 (90)
n11
where we let G.
_ ra(r)d
Let (w and recall the definition of A(w; x)
a11 ())
in Theorem 3.T then Eq. equivalent
72
Σ
324
which is exactly the KKT condition of arg minw: A(w; proposed in Theorem 3.1 Therefore, we finish the proof.
B.2 Proof of Lemma B.3
In this section we present the proof of Lemma B.3
Proof. For simplicity, in the following we write the subscripts explicitly. According to Lemma B.2 the dynamics of wy can be written as
w =
v
vjOws
T. - nTT kij
(91)
#
Note that
v3 w/ = 4w
4.j
w -13
v3 V w/ + 4n²
(92)
then Eq.
911 can be written
W
Ow: L
nw
(93)
4 w3 + 4x²
w/3
11
4 /w3 + 4m²
We now define a function
(ww,t:k) = Ay(w, + (1)
(94)
for some Ag (w,, f; x) and (t) such that
w, n k
dt Ow (sW,
(95)
4v //w3 w3 + 4x3
the we can prove this lemma. Now we continue to find the A, (w, t; k) and (t).). By definition,
d dt Bw, As (w, = Ol, A,W, + 0,8w, A + bb.
(96)
23

---

which, when compared with Eq. (95). implies that
Se, A
(97)
av ww3 4n3
Solving this equation gives us
In W/3 + 4mg + w))
Ow Aw A3 /
dw
+
(98)
4k
where is a constant and can be determined by requiring Ow, Agle= (2;(0)
c
10()k)00))/4. Thus Eq. (98) becomes
w3 + 4xz(1) + wj
(1)
Ow,A, Ag In
2xy(t)
4
where we have used the definition of (t(t) Lemma B.1 the form of A
Solving the above equation will give us
A, = ʃ dw, arcsinh
wj
(t,(t)
(1,((1)
4
= 4 w aresinh
W
V w/ + 4ng(t) + 2R2((t)
(t)(t),,
2k, (t)
4
G (w: x()) ())s(()w
(99)
4
where we use the definition of A GF in Eq. requires that
25
Comparing the rest parts of Eq. (96) with Eq.
(95
RJ
Brdw, A, + +3 =
Av w + 4t)
n (t)
= (t((t) =
+
kis
(100)
w, + w/ (t) / ww + 4x)(t)
1v w/ + 4x3(t)
When combined with the form of Ag. we can find the form of Ay:
ds
k))ss)
Ay(w, (:) A,TH (w; k()) nww
w, + v w/ + 4x2(s) éj
/w3 (x)
w3 + 4ng(s)é
Tw
kid
ds
Wjéj + kij
- (w;)(tt) + nww
w/ 4x3(()) (s)
4
w Tkuj
AGF (w; k()) + 7w,
ds
(101)
W+jj
W
where we use the definition of é (Lemma B. and Eq.
92 in the last equality.
B.3 Implicit Bias of HBF for Diagonal Linear Networks when a = 2
In this case, the correction term Wt will be
1 1 V w. L. V2 Vw L.
2(1 11)
24

---

We need to first find the Hessian V wy derive the Hessian by writing the subscripts explicitly. We start with w+.
L. Due to the element- product, it will be convenient to
2
Sw I. Switi
( (W++( r);)
n
2
()((XT)) +
Σ
W S (Fagle) w yc))
11
2
6xy(X²' r) + Σ
w W+15 2tj Fea W W
(102)
n
c==
where we use the delta symbol diy if j j otherwise dig = O. Therefore, we can conclude that
diag( (Tr) + Σ
(w+ 2 (w+) 3 xx)
(103)
W
cml
Following a similar approach, we obtain that for w-
2 dw
Ow Ow
(w (W-s(XTT))) r)))
71
n
5ig (1 r)j + Σ
w Neg w
(104)
n
2
diag(X (X r) + Σ (w_ C [w)(w- 2 2x,)²
(105)
72
It is now left for us to find the form of V L. V L. Again, it is convenient to write the subscripts explicitly:
d
(Vw+L
V.L.VV...
L.
wx 1)
Ow+n Own.s Ldw L
i=1
(ij(X* r) +2 Σ w Fe (+
n
i=1
4
W W+;i((XT r), +2
tyjeey (Te ) (106)
n
Similarly,
(v L. LW L.
(x r),) -2 - Ww (107)
n
c=1
Using Eq. (106) and (177). we can derive that
W
2(1 + H))
((XTr)))² + Σ
Boy (
(108)
(1 nn
c=1
which further gives us the integral ej:
E +9
4(1 + (4)
Tey (e (wi,. - -w
(1 (1- (1-1)3n2 n2
c=1 tml
4(1 + +1)) [1 (1- nn
((xT)) r)
Tenso w (w ( (XTr))
4(1 Jr)
(Vw ((X(wWVwL)))
(109)
(1 11)3
25

---

On the other hand, according to Lemma
B.2
Ow. L can be written as
W;
11)
111 p2)
(110)
4v; Av
which further gives us that
7 1 ds
wi(f)
w:(s)
=-1(1 11)
dw:
+0 (17))
w.(())
4v:(8)
will)
ws(s)
=-1(1 - p))
Tey
dw:
4K7(s) +07))
w.(()
4
(1 14)
[w + 4x2 v w? (0) + 4x²(0)
4
where we use Lemma B.2 in the first equality and Eq. and Lemma
(2 in the second equality. Since w(0) =
B.1 obtain
1 do
1711 - 11)
(vw2(1) + 4x²(0) - 2 (x((())
4
c.l s=1
1111 11)
(X (x(())())
(111)
4
where we let Rd with
(t) = V (1) + 4x² 2nz, 2 O.
Now combining Eq. (109) and Eq. (111) we can derive
(t) = 4n(1 + (1) ds ds(Ow, L)²
m(1 1) (1
(x" Xq) + +0 (1²)
(112)
(1 11)3
2Pn
To obtain the full potential function, we still need to find the form of yoy- According to the definition of v and 3j and Eq. 108). we can derive
W+
W Vku
27kg W W+s
YTk W
Wjēz= =
W
4(1 + 1) (1
Ve, L,
(113)
3n
i=
which, when combined with the definition of 2j Lemma B.3
further gives us
(1 p.)
n
c=1 1=1 4 = 1/ (1 - n1
Tej V L
(11 + p.) 4(1 13 c=1 t=1
Tey w/ + (
(114)
where we use Eq. (110) in the second equality. As result,
7/(1 + 11)
71( + μ) 1
())))
(xT (xxXw);
(115)
4(1 11) n
One interesting thing aspect of ;0 if w converges to an interpolation solution where Xw((o) = = is
17(1 + (4) 4(1 - p1)² Ow, L(0).
(,(00)=
(116)
In summary, the potential function O (172)-close is
k)(( (0) exp
47/(1 + 11)
Jo
ds(Ow; 1)2 + 71(1 + 1))
(xf xx(xo))))
(1 1133
(1 - u22n
As(w,oo; = A (G x(0)) + 1 + 11)
w, w,Ww, L(0).
(117)
4(1 - 11
26

---

C Details for Numerical experiments
The experiments conducted CentOS Linux 7.9.2 platform equipped with an Intel(R) Xeon CPU E5-2683 at 3.00 GHz, 256GB of RAM, and an NVIDIA Tesla A100 graphics card.
C.1 Details for Section 2.2.2
For the experiment of Figure I we conduct observation on the comparison of directional smoothness for HB and GD on the CIFAR-10 dataset Krizhevsky et al. (2009). A multilayer perceptron with two hidden layers (each of width 200) is trained for 2000 epochs using full-batch GD and HB (μ 0.9). and the step size set to 0.1.
C.2 Details for Section 4
For the discrete learning dynamics of HB and GD, we set the step size as 1 and the momentum factor For the continuous approximations, we Tiller 1010 as the Euler step sizes approximate the dynamics. These hyper-parameters are listed in Table
E. y Starting point
0.6
a1
2.8, 02 = 3.5
x 10-3
/ Mualer
0.7
x 10
Table 2: Hyper- parameters for 2-d model.
We let the model parameter be (21, E R.2. For RGF, the ODE
V Bx+1 Ba 1/Euler /
VaL
Formulations of HBFs with a = 2,3 are denoted in Table We denote 1, = 1===1,... (1,...,1) 1)" E Rd. For the dataset [x = we set n 40, 100. The data point follows a Gaussian distribution (O,, Ia). To make the ground truth solution w sparse, let 5 components of it be nonzero. Recall that the initialization is x(0)= = s² s²1, where s controls the initialization scale. In Fig. 6(b) we make the initialization as W w_ s1, with 0.01. We set the step size for HB as I0 For RGF and HBF, we let the Euler step size 1/Ealer 10-4 to simulate the continuous dynamics. In Fig. 6(a) and 6(b) set 7) 10 For the initialization, make the task slightly harder, we let w 1sl. and w sla/0 with 9 = 0.9 such that we still have x0) = s2 1, while the initialization symmetry slightly broken.
Addition experiments for different Baler We additionally run the experiments Fig. for each of T/ukerr {1/10,, 10, 1/ 1/000, n/ 2/1000J, and confirm that the observations and conclusions hold in all different 7/suter Fig.
C.3 Additional Numerical Experiments for Non-Linear Networks
We conduct experiments of Fig. 2(b) in the MNIST dataset, where we now train a three-layer fully-connected neural networks (FCNN). The FCNN has structure of Linear (784x 128) - SiLU->Linear (128 x for the loss function. The batch size is 60,000 and the momentum factor (0.7, 0.8). The learning rate is = 0.01. The results (reported in Fig. well align with the results for the toy model in Fig. 2b) HBF with a has lower discretization error compared HBF with a and both of them are better than the RGF.
(128 10). Cross-entropy is used
C.4 Additional Numerical Experiments for Diagonal Linear Networks
We now investigate the implicit bias of HB for 2-layer diagonal linear networks
f(x; w) w a (w+ w 3 w_
(118)
27

---

25.
HB= Oower 6/10
HBF.a 2. 0.7, 0010
1115.a.3,jp00...0mm (aa 1/10
RGF, H- (por 6/100
10F.e=2.yv007,,0
HBF.a=3,u-0.33 0.7, m/100
RGF, p= Otver g/000
HBF,a 2. Power
HBF.a 3, 0.7, g/00
0.0
0.5
2.0
2.5
a1
Figure 3: Trajectories of continuous approximations with different Muuler
Discretization error of HBF in MNIST
102, grw D.T
1.73
HeN, are 2.p
10%, prm Du
1.58
1.26
113
1.00
0.25
Iteration K
Figure 4: Discretization error HBF in MNIST when training a three non-linear networks.
for a dataset (e))) where x IRd y E R. The empirical loss is
I(w+. w-) ()(:::))- - y;)
(119)
We let n < and denote the ground truth solution as w such that w = y. We let w be sparse. For given scale we
Ki (0) w + (0) c w- (0)
1)T Rdd
(120)
C.4.1 Discretization Error for Difffrrnt Approximations
Our first experiment explores the discretization error, where let & denote the iteration count and first obtain wHB by training [z; w) with HB. In addition, we also train [a; w) with RGF (Eq. (1) and HBF (Corollary 3.2. respectively. We calculate the discretization error as
|
w(t+))l
(121)
for w(x)) obtained from HBF or RGF and present the results in Fig. 5, where HBF enjoys smaller discretization error than RGF for different μu, supporting theoretical claims.
C.4.2 Implication for difference of implicit bias between HB and GD
We first discuss the implication obtained from our HBF, Setting f = Corollary gives us the implicit bias of O (cl)) continuous approximation of GD, i.c., IGR (Barrett and Dherin 2022
28

---

10
1| ((G
10
10
HBF, # #
RGF, #=066
HBF, 0=
RGF, 0
HBF, 0.8
10
RGF, 0.8
2000
4000
6000
8000
10000
Iteration
Figure 5: Discretization errors l| HB
(1
2 for HBF (dotted lines) and RGF (solid lines),
respectively, when training the 2-layer diagonal linear networks.
Flow (IGRF). Both IGRF and HBF have the initialization rescaling effect. The difference between them is closely connected with the composed parameter
7 + 1)
(122)
(1 14)2
and the value of (/(1-) 14). And the discrepancy between the implicit bias of HB and that of GD will be more obvious for large value of μ. These observations stand in contrast to the case for O (1)see RGF, which cannot distinguish the implicit bias of HB from that of GD.
In addition, we note the following observations:
(i). the value of
1741+ +
( ds(O)))) (13
exp(
(123)
is related to the speed of convergence since it depends on the integral of gradient along the training trajectory, implying that faster speed of convergence would possibly lead a smaller fo daL))) which then leads to a larger (ii). smaller x implies a better sparsity and generalization performance for sparse regression, because the objective function Corollary 3.2 will be closer to the In-ormm
Now let k HBF and x IGRF be obtained from HBF and IGRF, respectively. As HB converges faster than GD in practice, then ÖL becomes neglectable very quickly for HBF if it converges too fast, and, according to our observation (i) above, we conclude that K HBF > x IGRF, Then according to our observation (ii) above, IGRF will generalize better than HBF. On the other hand, if the speed of convergence of HBF and that of IGRF are similar (e.g., blue line in Fig. 6(b)). then HBF and IGRF will have similar values of 1+
fo dsL)) while HBF additionally depends on coefficient TI-A]ø > 1, thus it is possible that in this cae K HBF < IGRF which implies that HBF will generalize better in this case. In summary, there might exist a tension between the speed of HB, i.e., if HB converges too fast, K. would be larger for HB hence solutions of GD would enjoy better generalization properties.
convergence and the generalization
Below we conduct experiments for the above claims. In particular, we compare the implicit bias of HB with that of GD. Given we train f(z; w) with GD and HB, respectively. We calculate the distance between the returned solution w(20) and the ground truth solution w' i.c., || ||w((o) as measure of generalization performance and report the results in Fig. 6(a) It can be seen that, when the initialization scale s small, solutions GD generalize better than those of HB. This can be explained by Corollary 3.2 compared to GD, when is small, L(w) decreases much faster for HB (green lines Fig. 6(b) which leads to a smaller ddsL(w) and weaker initialization mitigation effect, thus the solutions of HB generalize worse than GD solutions. Recall that in Corollary 3.2, as ky (0) = 22 increases, & determines the generalization performances for HB and GD since controls the extent of the initialization mitigation effect. Furthermore, L(w) does not decrease much faster
l|2-
29

---

25
H8.
p- 0.8
12
12
-10
1.0
2.0
W
14
1.5
A(W
0.4
0.5
10742
0.5 Initialization scale
0.4
0.7
Iteration K
(a)
(b)
Figure 6: (a). Generalization performances |[() - ww for different initialization scales s when fæ; w) is trained by GD and HB with different values of μ.. (b). L(w) during training processes of HB (μ = 0.9) and GD for different s.
for HB than for GD (blue lines in Fig. 6(b) thus GD and HB have similar value of t which is further enhanced by a factor of (1 + / 1 for HB according to Corollary 3.2 As a result, HB solutions will generalize better than GD and the discrepancy between them is more significant for large (large (1 + p1/(1- - =) as shown in Fig. 6(a)
30

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: This paper provides new continuous time model for HB method. HBF, and investigates several properties of the discrete HB method through lens of the HBF, These contributions exactly match those in the abstract and introduction.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Please see the separate paragraph "Limitation" Section
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-specfiiaatinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g. if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
31

---

Answer:
[Yes]
Justification:
Proofs for all theorems are presented in Appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: We have presented numerical experiments in Section 2.2.2 and Section 4 with the corresponding experimental details in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
32

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [NA]
Justification: We use synthetic data and open-source data, which are discussed in Ap-
pendix C
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.. cc/
Polooyooollyyooooo for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
//nips.cc/pbble Sodi:///uuuddee/deeuuuubbiiii
(CodeSubiisss wwtplubmmmsssiioo ssiiiiolcccy for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Please Appendix C
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer:
[NA]
Justification:
We do not have error bars to report.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi-
dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for
example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
33

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative rates).
show in tables or
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Please see Appendix C
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't 't make it into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics http://eeerrpp..cccoubbbttttbbba
Answer:
[Yes]
Justification:
The research conform with the NeurIPS Code of Ethis in every respect.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- cration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: This theoretical paper, which dose not have direct potential societal impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
34

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This is theoretical paper, which does not provide data models.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We use the open-source data CIFAR-10, which are properly credited in this paper.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
35

---

If assets are released, the license, copyright information, and terms of use in the thcode has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
latasets.paperssit
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: This paper does not provide new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
These topics are not covered by this paper.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [Yes]
Justification: This paper does not involve any of these subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
36

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM for what should or should not be described.
37


Reviews:
[{"review_id": "tNT5GJAg7b", "reviewer": "Reviewer_nQN7", "summary": "In this paper, the authors propose a continuous ode approximation of the popular Heavy-ball momentum optimization algorithm. Compared to the previous studies, the framework proposed in this paper can achieve an error as small as arbitrary order of the small interpolation, i.e. the learning rate $\\eta$. Based on this conclusion, authors propose a new implicit bias of diagonal linear networks. While simple and idealized, this model usually serves as the first step for multiple theoretical studies.", "questions": "While not a major issue, as a study including the implicit bias of momentum, I suggest the author could consider discussing [1, 2, 3] in the related works. \n\n[1]. Wang et al. The implicit bias for adaptive optimization algorithms on homogeneous neural networks. ICML.\n\n[2]. Wang et al. Does Momentum Change the Implicit Regularization on Separable Data? Neurips.\n\n[3]. Zhang et al. The implicit bias of Adam on separable data. Neurips.", "limitations": "The potential incorrectness in the proof by using the notation of $O(\\cdot)$ in induction.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strength:**\n\n1. Compared to the previous studies, the framework proposed in this work can better capture the optimization trajectories of the Heavy-ball momentum, as it can achieve an arbitrarily small error, which, also claimed as the major technical contribution of this work.\n\n2. The study of the implicit bias of momentum implies the potential benefits of studying the the higher-oder approximations.\n\n**Weaknesses:**\n\n1. While I found the approximation w.r.t. an arbitrary order of the interpolation size $\\eta$ is relatively interesting, the intuition behind such an approach is hard to follow, especially after reading the proof details in the Appendix. First of all, from my understanding, the presentation of equation (6) is quite different from the standard form of second-order Taylor expansions. Indeed, after checking the definition of $\\eta^2 I_k^+$ in line 140, it is exactly $\\beta(t_{k+1}) - \\beta(t_{k}) - \\eta\\dot{\\beta}(t_{k}^+)$. Therefore, from my perspective, term $I_k^+$ is entirely determined by $\\beta(\\cdot)$ and is fixed through the entire stage. However, from the proof and the discussion of the proof, it seems that $I_k^+$ is **specifically designed** to match the term of $\\gamma$. Besides, what is the advantage to utilize the calculation of equation (47)? It seems that it still introduces the $n$-th order differential calculations here. I wonder more intuitive explanation of such a calculation. In addition, I hope that authors can also clarify some of the notations. For example, is the domain of $G_k(\\beta)$ only from $\\beta_k$ to $\\beta_{k+1}$, or the entire space? If it is the former, then $G_k(\\beta)$ would always be a multiple of $\\nabla L(\\beta)$, in which case it is unclear what is gained by defining distinct $G_k(\\beta)$ functions. If it is the latter, then how should we interpret the expression $G_{k-1}(\\beta_{t_k})$ in equation (44)? Moreover, when using the shorthand $G_k$, what is the intended independent variable?\n\n2. The use of the notation $O(\\cdot)$ can lead to significant incorrectness in induction and should therefore be avoided in such contexts. To illustrate what I claim, let's consider a classical paradox when using the $O(\\cdot)$ in induction. For any integer $K$, it is obvious that we can write $K=\\sum_{k=1}^K a_k$, with $a_k=1$. Clearly, $a_1 = O(1)$ holds. Next, assume that at the $k'$-th iteration $\\sum_{k=1}^{k'} a_k = O(1)$, which leads to the conclusion that $\\sum_{k=1}^{k'+1} a_k = O(1) + 1 = O(1)$. Consequently, for any integer, no matter how large it is, we can always claim that it is $O(1)$, which, is obviously incorrect. To make induction mathematically rigorous, authors should explicitly specify the universal absolute constant hided in $O(\\cdot)$, and keep such a constant fixed during the proof. \n\n\nI would like to increase my score if the authors could provide further illustration and modify their proof (especially the $O(\\cdot)$ notation used in their induction proof of Lemma A.2.)", "was_revised": true, "final_justification": "The rebuttal addresses my concerns. I have updated my score to 4.", "rebuttal": "We thank the reviewer very much for the insightful comments. We appreciate this opportunity to address your concerns.\n\n---\n## Weaknesses Part\n\n### Weaknesses 1\nWe address concerns in Weaknesses 1 point by point.\n- **Comment**: *\"While I...Appendix.\"*\n\n   **Response**: To better illustrate the intuition of our approach, we first summarize it below. Overall, our goal is to find the formulation of the counter term $\\gamma_k$ for the ODE $$\\dot{\\beta}=-G_k-\\eta\\gamma_k\\tag{a.1}$$ such that the discretization error can be controlled as $\\varepsilon_k = O(\\eta^{\\alpha})$. To achieve this, our approach aims to find a functional equation for $\\gamma_k$ that must be satisfied under the condition $\\varepsilon_k = O(\\eta^{\\alpha})$, and this functional equation can be solved to give the desired $\\gamma_k$.\n- **Comment**: *“First of...term of $\\gamma$”*\n\n   **Response**: Eq.(6) is the first-order Taylor expansion with the remainder term $I_k$ in the integral form. $I_k$ is determined by $\\gamma_k$ as noted by the reviewer, and it is not specifically designed to match $\\gamma_k$. Instead, it is $\\gamma_k$ that is specifically designed to solve the functional equation mentioned above by matching terms.  \n \n   To illustrate this in detail, below we indicate how we establish and solve such a functional equation. Given the unknown $\\gamma_k$ in Eq.(a.1), our method follows a general flow:\n    $$\\begin{aligned}\\boxed{\\text{Step 1}}\\quad\\text{unknown }&\\gamma_k\\overset{\\text{determines}}{\\longrightarrow}I_k^{\\pm}\\\\\\\\&\\ \\searrow\\qquad\\ \\swarrow\\\\\\\\\\boxed{\\text{Step 2}}\\quad\\ \\ \\ \\text{Taylor }&\\text{Expansion of }\\beta(t_{k\\pm1})\\\\\\\\&\\qquad\\ \\downarrow\\\\\\\\ \\boxed{\\text{Step 3}}\\quad\\quad\\text{Expres}&\\text{sion of }\\varepsilon_{k+1}-\\varepsilon_k\\\\\\\\&\\qquad\\ \\downarrow\\text{\\small required to be }O(\\eta^{\\alpha+1})\\text{ \\small to guarantee }\\varepsilon_k=O(\\eta^{\\alpha})\\\\\\\\\n   \\boxed{\\text{Step 4}}\\quad\\ \\text{Function}&\\text{al equation for }\\gamma_k\\\\\\\\&\\qquad\\ \\downarrow\\text{\\small solved by matching to each order of }\\eta\\\\\\\\\n    \\boxed{\\text{Step 5}}\\qquad\\qquad\\quad&\\text{Solution of }\\gamma_k\n    \\end{aligned}$$\n   Specifically,\n    - In Step 1, as noted by the reviewer, given arbitrary unknown $\\gamma_k$ in Eq.(a.1), $I_k^{\\pm}$ is then determined, which is also unknown but depends on $\\gamma_k$ (e.g., their relation is constructed in Lemma A.1 (line 463-466, page 13) if $\\gamma_k$ is a series)\n    - In Step 2, given the unknown $\\gamma_k$ and the corresponding $I_k^{\\pm}$, we can derive the Taylor expansion for $\\beta(t_{k\\pm1})$ (see Eq.(6) in line 138-139), which depends on both $\\gamma_k$ and $I_k^{\\pm}$.\n    - In Step 3, given the Taylor expansions in Step 2, we can construct a set of expressions for the discretization error $\\varepsilon_{k+1}-\\varepsilon_k$ (Eq.(7) in line 141-142 page 4). \n    - In Step 4, to guarantee $\\varepsilon_k = O\\eta^{\\alpha})$, the condition $\\varepsilon_{k + 1} - \\varepsilon_k = O(\\eta^{\\alpha + 1})$ is sufficient (this statement will be formally proved in Lemma A.2 (line 477-481 page 14)), which can be satisfied by letting $$\n    I_k^++\\mu I_k^--\\gamma_k+\\mu\\gamma_{k-1}=O(\\eta^{\\alpha-1})\\tag{a.2} $$ \n    in Eq.(34) (line 454-455, page 13). \n    - In Step 5, as $I_k^{\\pm}$ depends on $\\gamma_k$, Eq.(a.2) is in fact a functional equation for $\\gamma_k$. To solve this equation, we write $\\gamma_k$ in the series form $$\\gamma_k=\\sum_{\\sigma}^{\\infty}\\eta^{\\sigma}\\gamma_k^{(\\sigma)}, $$\n    then Lemma A.1 provides us the corresponding $I_k^{\\pm}$, which is also in the series form: \n    $$I_k^{\\pm}=\\sum_{\\sigma}^{\\infty}\\eta^{\\sigma}(\\mathscr{I}_k^{\\pm})^{(\\sigma)},$$ where each $(\\mathscr{I}_k^{\\pm})^{(\\sigma)}$ is determined by the lower-order terms of $\\gamma_k$. If we match $\\gamma_k-\\mu\\gamma _ {k-1}$ with $I_k^{+}+\\mu I_k^{-}$ to each order of $\\eta$, then we obtain a set of conditions that must be satisfied by $\\gamma_k^{(\\sigma)}$ for any $\\sigma$, which can be solved to give the forms of $\\gamma_k^{(\\sigma)}$.Then a simple truncation of $\\gamma_k$ to the $(\\alpha-2)$-th order, i.e., \n  \n      $$\\gamma_k = \\sum_{\\sigma}^{\\alpha - 2} \\eta^{\\sigma} \\gamma_k^{(\\sigma)},$$ guarantees Eq.(a.2).\n- **Comment**: *“Besides,...calculation.”*\n    \n   **Response**: We first confirm that Eq.(47) still introduces $n$-th order differential calculation. The advantage of Eq.(47) is as follows. \n\n   Intuitively, we aim to derive the formulation of $I_k$ in the series form such that the functional equation Eq.(a.2) can be solved by matching terms to different orders of $\\eta$. Hence, we need to group together terms of the same order in $\\eta$ in $d^n\\beta / dt^n$ (as $I_k$ is determined by $d^n\\beta / dt^n$). The advantage of Eq.(47) is that the operator, $(\\sum_{\\sigma = 0}^{\\infty} \\eta^{\\sigma}\\gamma_k^{(\\sigma)})\\cdot\\nabla$, depends on $\\eta$ explicitly, allowing us to gather terms in the same order easily.\n- **Comment**: *\"In addition,...equation (44)?\"*\n    \n   **Response**: We first confirm that the domain of $G_k$ is the entire space. Below we answer the questions separately. \n   - The reason why we use the notation $G_k$ is to indicate that we expect it to be dependent on the iteration count $k$ explicitly in the modified equation Eq.(a.1).\n    - Interpret Eq.(44): We would like to first mention that $G_{k-1}(\\beta_{t_k})$ in Eq.(44) is a typo, and it should be $G_{k-1}(\\beta(t_k))$. We apologize for potential misunderstanding due to this typo. As mentioned earlier, the condition $\\varepsilon_{k + 1} - \\varepsilon_k=O(\\eta^{\\alpha+1})$ is sufficient to guarantee $\\varepsilon_k = O(\\eta^{\\alpha})$. To make this condition satisfied at $\\beta(t_k)$, besides the functional equation Eq.(a.2), we require a general recursive relation for $G_k$ in the entire space of $\\beta$: $$\\forall k\\in [N]:G_k(\\beta)=\\mu G_{k-1}(\\beta)+\\nabla L(\\beta),\\tag{a.3}$$which implies $$\\forall k\\in [N]:G_k(\\beta)=\\frac{1-\\mu^{k+1}}{1-\\mu}\\nabla L(\\beta)\\tag{a.4} $$in Eq.(16) (line179-180, page 5). As the domain of $G_k(\\beta)$ is the entire space, Eq.(a.3) is also satisfied at $\\beta(t_k)$, which is exactly Eq.(44). Thus, along with the functional equation Eq.(a.2), $\\varepsilon_{k + 1} - \\varepsilon_k = O(\\eta^{\\alpha + 1})$ can be guaranteed. For large $k$, $G_k$ becomes the rescaled gradient $\\nabla L/(1-\\mu)$.\n\n   - The use of $G_k$: When using the shorthand $G_k$, the intended independent variable is $\\beta$, i.e., it is a shorthand of $G_k(\\beta)$ and closely related to gradient via Eq.(a.4). The subscript $k$ only indicates that $G_k$ will change for different iteration count $k$, and this dependence on $k$ will be negligible for large $k$.\n### Weaknesses 2\n- **Comment**: *“The use...proof.”*\n    \n   **Response**: We thank the reviewer for this insightful suggestion. We provide our proof below, where we specify the universal absolute constant in the $O$ notation explicitly as suggested by the reviewer. \n    \n    We first present several useful relations. As Eq.(43) is established by solving the functional equation for any iteration count $k$, we can write the relation between $\\varepsilon_{k + 1}$ and $\\varepsilon_{k}$ (Eq.(34), line454-455 in page 13) as $$\\varepsilon_{k+1}-\\varepsilon_{k}=\\mu(\\varepsilon_{k}-\\varepsilon_{k-1})-\\eta[\\nabla L(\\beta(t_k))-\\nabla L(\\beta(t_k)-\\varepsilon_{k})]+O(\\eta^{\\alpha + 1}).\\tag{a.5}$$If $\\varepsilon_{k}=O(\\eta^{\\alpha})$, then Eq.(a.5) implies $$||\\varepsilon_{k+1}-\\varepsilon_{k}||\\leq\\mu||\\varepsilon_{k}-\\varepsilon_{k-1}||+\\eta\\lambda||\\varepsilon_{k}||+c_1\\eta^{\\alpha+1}$$for some constant $c_1$ where we use $G_k-\\mu G_{k-1}=\\nabla L$ and let $\\lambda=\\max_{\\beta}||\\nabla^2L(\\beta)||$. Denoting$$\\forall k:c_2(k)=\\frac{c_1}{\\lambda}e^{\\frac{2\\lambda \\eta}{1-\\mu}k},\\ c_3(k)=\\frac{2c_1}{1-\\mu}e^{\\frac{2\\lambda\\eta}{1-\\mu}k},$$we now prove by induction.  \n    - For the first step ($k = 0$), by definition we have $$\\varepsilon_0=0\\leq c_2(0)\\eta^{\\alpha}.$$\n    Note that $G_{-1}=0$ and $\\gamma_{-1}=0$ by definition, we have $$||\\varepsilon_1-\\varepsilon_0||\\leq c_1\\eta^{\\alpha+1}\\leq\\frac{2c_1}{1-\\mu}=c_3(0)\\eta^{\\alpha+1}$$since $\\mu<1$.\n    - Suppose that for the $k$-th step $$||\\varepsilon_k||\\leq c_2(k)\\eta^{\\alpha},\\ ||\\varepsilon_{k+1}-\\varepsilon_k||\\leq c_3(k)\\eta^{\\alpha+1}.$$ Then for the $(k + 1)$-th step, we have $$\\begin{aligned}||\\varepsilon_{k+1}||&\\leq||\\varepsilon_{k+1}-\\varepsilon_k||+||\\varepsilon_k||\\\\\\\\&\\leq c_2(k)\\left(1+\\eta\\frac{c_3(k)}{c_2(k)}\\right)\\eta^{\\alpha}\\\\\\\\&\\leq c_2(k)e^{\\frac{2\\lambda \\eta}{1-\\mu}}\\eta^{\\alpha}=c_2(k+1)\\eta^{\\alpha}\\end{aligned}$$where the last inequality is because $e^x>1+x $ for $x>0$. Similarly, according to Eq.(a.4), $$\\begin{aligned}||\\varepsilon_{k+2}-\\varepsilon_{k+1}||&\\leq\\mu||\\varepsilon_{k+1}-\\varepsilon_{k}||+\\eta\\lambda||\\varepsilon_{k+1}||+c_1\\eta^{\\alpha+1}\\\\\\\\&\\leq\\left[\\mu c_3(k)+\\lambda c_2(k+1)+c_1\\right]\\eta^{\\alpha+1}\\\\\\\\&=\\left[\\mu e^{-\\frac{2\\lambda\\eta}{1-\\mu}}+\\frac{1-\\mu}{2}+\\frac{1-\\mu}{2}e^{-\\frac{2\\lambda\\eta}{1-\\mu}(k+1)}\\right]c_3(k+1)\\eta^{\\alpha+1}\\\\\\\\&\\leq\\left[\\mu+\\frac{1-\\mu}{2}+\\frac{1-\\mu}{2}\\right]c_3(k+1)\\eta^{\\alpha+1}=c_3(k+1)\\eta^{\\alpha+1}.\\end{aligned}$$This now proves our claim: for a fixed time $T = N\\eta$ such that $N = \\lfloor\\frac{T}{\\eta}\\rfloor$,  we have $$\\forall k\\in [N]: ||\\varepsilon_k||\\leq\\frac{c_1}{\\lambda}e^{\\frac{2\\lambda}{1-\\mu}k\\eta}\\eta^{\\alpha}.$$\n\n     Note that we need the total time $T = N\\eta = O(1)$ to control $\\varepsilon_k$ globally. This is a typical requirement for studying continuous approximations with bounded global discretization error. We will clearly state this condition in the revision.  \n\nFinally, we thank the reviewer for these insightful questions and suggestions, and we will update the proof and provide a more clear sketch of the corresponding intuition in the revision to improve the readability.\n\n----\n\n## Questions Part\nWe thank the reviewer for suggesting these nice related works. We will discuss these related works in the revision."}, {"review_id": "bf4IlEVPSa", "reviewer": "Reviewer_YrV8", "summary": "This work proposes a continuous time framework--HB Flow (HBF)--for the heavy-ball momentum method (HB). The authors show HBF approximate HB up to arbitrary order. They also study the implicit bias of HBF on diagonal linear networks which brings insights into the implicit bias problem in momentum methods. Experiments on simple neural networks are done to verify their results.", "questions": "1. In Figure 2(a), I noticed that the $\\alpha=2$ case gives a better approximation away from the convergence point while $\\alpha=3$ is the opposite. I wonder if the authors have some intuition for this behavior.\n\n2. This analysis can be generalized to the stochastic case by for example replacing $\\nabla L$\n by $\\tilde{\\nabla} L$\n where $\\tilde{\\nabla}$\n denotes the approximate gradient. The approximation can be done by either introducing a diffusion term [1] or an index switching process [2]. I think this might be an interesting future direction.\n\n[1] Stochastic modified equations and dynamics of stochastic gradient algorithms I: Mathematical foundations, Li et al.\n\n[2] Analysis of Stochastic Gradient Descent in Continuous Time, Latz.", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "*I have reviewed this paper last year. I don't see many significant changes. So I am borrowing some of my previous review here, which the authors had not yet addressed.*\n\nStrengths.\n\n1. The authors propose HB flow--a continuous-time approximation of the heavy ball momentum method. By adding back the discretization error in each time interval, the HB flow achieves an arbitrarily close distance to the discrete-time dynamic. This is justified rigorously in theorem 2.1\n\n2. More precisely, the distance between two dynamics is measured using $L^\\infty[0, T]$ on the trajectory space. The $L^\\infty$ is qualified by an upper bound in the learning rate $\\eta^\\alpha$. The authors give the specific form for the case $\\alpha = 1, 2, 3$ as examples, which help the readers to understand the construction of the continuous dynamic.\n\n3. The authors computed the implicit bias of HBF for diagonal linear networks in theorem 3.1 which makes this framework useful in studying the generalization properties of HB. This type of behavior cannot be captured by lower order models like RGF.\n\n4. In numerical experiments, Figure 2 illustrate the trajectories in a clear way and verifies main results from the paper.\n\n5. I see the authors added one experiments about the implicit regularization. \n\nWeaknesses.\n\n1. I found the main theorem hard to follow. Specifically, Theorem 2.1 is stated with many notations without enough intuition/explanations. I am trying to understand this at a high level. So basically find the difference between continuous and discrete dynamics and use Taylor expansion for it. The order $\\alpha$ in $\\eta^\\alpha$ just means how many times we do Taylor expansion in $t$ because the time interval is of length $\\eta$. Although the computation needed to find the exact expressions is not trivial, one could argue this is not mathematically sophisticated.\n\n2. Some notations are not consistently defined. For example, on page 4, the authors stated $G_k \\in R^d$, i.e. a d dimensional vector, but $G_k$ is a function. Similar situation for $\\gamma_k$.", "was_revised": true, "final_justification": "all concerns are addressed.", "rebuttal": "We thank the reviewer very much for the valuable comments. Below we address your concerns.\n\n---\n\n## Our Improvements \nWe would like to first highlight that we have made **a series of important improvements** in this version:\n\n1. We have added some paragraphs to help the readers to understand the merit of our approach, e.g., the overview of our approach in line 128-147, page 3.\n2. We have improved the statement of Theorem 2.1, where now we: \n   - first state what the formulation of the proposed ODE with discretization error to arbitrary order of $\\eta$ will be; \n   - then state what the explicit formulations of the components of the ODE will be by defining necessary notations.\n3. As noted by the reviewer, we have also added new experiments about the implicit regularization for the directional smoothness of HB, which is a novel result and reveals the difference compared GD. This clearly shows the significance of studying continuous time model with lower discretization order.\n\n---\n\n## Weaknesses Part\nBelow we address your concerns mentioned in the Weaknesses part.\n\n1. **Comment**: *\"I found the main theorem ... mathematically sophisticated.\"*\n\n   **Response**: We would like to first clarify that our results  are not obtained by performing Taylor expansion to the specific order of $\\eta$. Instead, our approach performs Taylor expansion with remainder term to the first order. Meanwhile, our approach involves a delicate construction of the modified ODE such that we are always able to build continuous time models with discretization errors to arbitrarily high order of $\\eta$ by directly employing existing lower order ones, rather than performing the whole set of analysis. This cannot be achieved by performing Taylor expansion to a specific order in advance. \n    \n    In addition, our result builds a continuous time model for HB with discretization error to arbitrary order of $\\eta$ for the first time, making it a significant one. Moreover, our approach might also be generalized to establish higher-order continuous time models for more optimization methods, e.g., Adam. This is because the intuition behind our approach is natural  and effective, and we do not think this is a downside of our approach to make it less mathematically sophisticated.\n\n    We now address your concerns. On one hand, the notations in Theorem 2.1 are unavoidable due to the fact that higher-order continuous approximations typically require a more delicate construction. On the other hand, we will further add more explanations and intuitive illustration for our approach at a high level (please find the explanation of our intuition below) as suggested by the reviewer in the revision.\n\n    **Intuition:** Overall, for the discrete HB method, our goal is to find a counter term $\\gamma_k$ for the modified ODE proposed by us $$\n    \\dot{\\beta} = - G_k - \\eta \\gamma_k, $$\n    such that the discretization error can be controlled as $\\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha})$. To achieve this, mathematically, we need to derive two sets of equations for $\\gamma_k$ and $G_k$, respectively, that must be satisfied under the condition $\\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha})$. Then these equations can be solved to give us the desired formulations for $\\gamma_k$ and $G_k$, respectively. Therefore, there are two crucial parts of our approach:\n    \n    - Find the equations that $G_k$ and $\\gamma_k$ must satisfy if we require $\\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha})$: this is achieved by (i). deriving the formulation of $\\varepsilon_{k + 1} - \\varepsilon_k$ by Taylor expanding of $\\beta(t_{k\\pm 1})$ at $t = t_k$ to the first order of $\\eta$ with the remainder term in the integral form; (ii). requiring $\\varepsilon_{k + 1} - \\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha + 1})$, which will be sufficient for us to prove $\\varepsilon_k = \\mathcal{O}(\\eta)$ hence can give us the corresponding equations for $G_k$ and $\\gamma_k$, respectively.\n    - Solve these equations to give the formulations of $G_k$ and $\\gamma_k$: deriving $G_k$ is straightforward, and the construction of $\\gamma_k$ is done by carefully deriving $I_k$ as an explicit functional for $\\gamma_k$ such that both of them will be in the series form, then the corresponding equation can be solved by matching terms in each order of $\\eta$.\n2. **Comment**: *\"Some notations ... Similar situation for $\\gamma_k$\"*\n\n   **Response**: We thank the reviewer for pointing this. For general $\\beta\\in \\mathbb{R}^{d}$, $G_k$ and $\\gamma_k$ are functions (vector fields). We will correct this in Page 4. At a specific point of $\\beta$, $G_k$ and $\\gamma_k$ will be $d$-dimensional vectors.\n\n---\n\n## Questions Part\n1. **Comment**: *\"In Figure ...\"*\n\n   **Response**: Intuitively, continuous approximations with discretization error $\\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha})$ with $\\alpha > 1$ are obtained by adding counter terms to that with $\\alpha = 1$, hence these counter terms will lead their trajectories to intersect the trajectory of discrete HB during the early stage. In particular, the intersection for HBF with $\\alpha = 2$ occurs in the early stage, making its performance appear to be better. In addition, although the trajectory is long in the early stage, the points are sparsely distributed as the gradient is large at this stage. Therefore, when evaluating the overall discretization error across all points, HBF with $\\alpha = 3$ tracks the discrete HB better, e.g., it demonstrates the lowest discretization error very quickly (after only 10 steps according to Fig.2(b)).\n2. **Comment**: *\"This analysis can be generalized ... future direction.\"*\n\n   **Response**: We thank the reviewer for the insightful suggestion. We agree with the reviewer that incorporating stochasticity by replacing $\\nabla L$ with $\\tilde{\\nabla} L$ is an interesting direction, and we believe it might be more interesting if our approach can be generalized to include the stochasticity for the more advanced optimizer Adam. We will add a discussion of the extension to stochastic HB by applying techniques from these suggested works in the revision."}, {"review_id": "m23vTTz2Fo", "reviewer": "Reviewer_PYzD", "summary": "This paper proposes a new continuous ODE model that can reduce the discretization error from the discrete HB method.\nBy ultilizing the Taylor expansion and discrete matching on the learning rate, the discretization error can be arbittrary order of the learning rate. This result provides a reliable foundation for analyzing the momentum methods in the continuous time limit.", "questions": "1.see the weakness.\n\n2.How can the theoretical insights from Theorem 2.1 be used to improve practical optimizer design?", "limitations": "yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The theoretical results are very attractive, and can be considered as a good extension of previous works. \nWhile (Kovachki and Stuart, 2020) and (Ghosh et al., 2023) explicitly anaylzed the discretization error to $\\mathcal{O}(\\eta)$ and $\\mathcal{O}(\\eta^2)$, this paper extended the result to $\\mathcal{O}(\\eta^{\\alpha})$ with arbitrary positive integer $\\alpha>0$.\nThis continuous model is applied to analyze the regularization term for the directional smoothness and the implicit bias of diagonal linear networks.\nThis bridge the gap between the discrete HB and the its precise characterization in continuous time.\nAlthough I did not check all the proofs in detail, I believe the overview of the proof is correct.\nThe following parts could be improved a bit more:\n\n1. The numerical experiments are far from enough to reflect the correctness of the main Theorem 2.1. This paper only conduct results on 2d-simple model and 2-layer diagonal linear networks for simple regression problems. More experiments on deep learning taskes should be explored to examing the discretization errors. What I expect to see is the discretization errors vs iterations on different $\\alpha$, like the Figure 2(b) in your paper (but in different tasks), or Figure3-4 in Miyagawa (2023).\n\n2. The mathematical results Theorem 2.1 are got by ultilizing the Taylor expansion and discrete matching on the learning rate used in (Miyagawa ,2023), (Rosca et al., 2023) for GD, and extended the techniques to HB, which covers the theoretical results of previous works  (Kovachki and Stuart, 2020) and (Ghosh et al., 2023). From this point of view, the theoretical contributions are limited.\n\n3. I have some difficulties when I first read the paper and before searching the cited works. I think the reason is that the \"overview of our approch\" in page 3 full fill of mathematical symbols without explicit explaning, and the following proofs are based  on the overview. The motivation and main technniques in this overview should be well-organized.\n\n4. Some typos: lacking $\\eta$ for the discrete GD in Table 1, lacking $\\frac{\\nabla L \\nabla^2 L}{2}$ in equation (20).", "was_revised": true, "final_justification": "Most of my questions have been resolved. I would like to see more experiments to support the theoretical validity of the work. So I  maintain the positive score.", "rebuttal": "We thank the reviewer very much for the constructive comments and suggestions. We will make the corresponding improvements as suggested by the reviewer in the revision. In particular:\n1. **Comment**: *\"The numerical experiments ... in Miyagawa (2023).\"*\n\n   **Response**: To address this concern, during this rebuttal period, we have conducted experiments similar to those in Fig.2(b) in the MNIST dataset, where we now train a three-layer fully-connected neural networks (FCNN) similar to that in Miyagawa (2023). The FCNN has a structure of $\\texttt{Linear}(784\\times128)\\to \\texttt{SiLU} \\to \\texttt{Linear}(128\\times128) \\to \\texttt{BatchNormalization}\\to \\texttt{Linear}(128\\times10)$. Cross-entropy is used for the loss function. The batch size is 60,000, i.e., we perform full-batch training. We let the momentum factor $\\mu\\in\\\\{0.7, 0.8\\\\}$ and take the learning rate as $\\eta=0.01$. The results well align with that for the toy model in Fig.2(b): HBF with $\\alpha = 3$ has lower discretization error compared to HBF with $\\alpha = 2$, and both of them are better than the RGF. We will add these experiments in the revision, as in this stage we are not allowed to upload such results. \n2. **Comment**: *\"The mathematical ... are limited.\"*\n\n   **Response**: We would like to highlight that our approach, while inspired partly by  prior works, bridge an important gap of establishing arbitrary order approximations of momentum methods and can be further developed as a general framework to include more optimization algorithms. In particular, compared to results for GD, our theoretical results successfully considered one important aspect of modern optimization methods--momentum--in the continuous time model. This extension is nontrivial. For example, methods with momentum explicitly incorporate effects of historical iterations while GD does not, hence the extension to momentum methods needs to carefully consider how to characterize such effects. \n\n    In addition, approaches in Kovachki and Stuart, (2020) and Ghosh et al., (2023) require specifying the desired order of the discretization error before performing the analysis, hence are not suitable to be generalized to the case with $\\mathcal{O}(\\eta^{\\alpha})$ for arbitrary $\\alpha > 0$. As a comparison, our approach bridge this gap by considering an arbitrary $\\alpha$ at the first place. \n   \n   Therefore, our results provide a viable way to develop a general framework for studying high-order continuous time models for a variety of advanced optimization methods that also utilize momentum such as Adam. Hence our contributions could be significant.\n\n3. **Comment**: *\"I have some difficulties ... well-organized.\"*\n\n   **Response**: We thank the reviewer for the suggestion. In the revision before presenting the discussion of \"overview of our approach\" in page 3, we will add a new paragraph to provide a more clear and concise verbal explanation of the motivation for our approach. This organization allows readers to first understand the motivation behind our approach before exploring the technical aspects.\n\n   **Paragraph that will be added to the revision** \"Given the discrete HB method, our overall goal is to find the formulation of a modified ODE, $$\n    \\dot{\\beta} = - G_k - \\eta \\gamma_k, \\tag{a.1} $$\n    such that the discretization error can be controlled arbitrarily low, where the design of the ODE Eq.(a.1) follows two basic principles:\n    - When $\\eta$ is very small, Eq.(a.1) should coincide with the simplest continuous approximation of HB, the rescaled gradient flow. Hence, $G_k$ should be able to degenerate to the rescaled gradient $\\nabla L / (1 - \\mu)$.\n    - By adding the counter term $\\gamma_k$, Eq.(a.1) should allow us to further decrease the discretization error of the rescaled gradient flow to $\\mathcal{O}(\\eta^{\\alpha})$ for any $\\alpha > 0$ to get better continuous approximations for HB.\n\n    With Eq.(a.1), it is left for us to derive the exact forms of $G_k$ and $\\gamma_k$. For this purpose, there will be two key steps:\n    - Find the equations that $G_k$ and $\\gamma_k$ must satisfy if we require the discretization error $\\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha})$ for any $\\alpha > 0$: this is achieved by (i). deriving the formulation of $\\varepsilon_{k + 1} - \\varepsilon_k$; (ii). requiring $\\varepsilon_{k + 1} - \\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha + 1})$, as this condition will be sufficient for us to prove $\\varepsilon_k = \\mathcal{O}(\\eta^{\\alpha})$ hence can give us the corresponding equations for $G_k$ and $\\gamma_k$, respectively.\n    - Solve these equations to give the formulations of $G_k$ and $\\gamma_k$: deriving $G_k$ will be straightforward, and the construction of $\\gamma_k$ will be done by carefully expressing  the corresponding equation obtained in the first step as an explicit functional for $\\gamma_k$ using the series form, which can then be solved by matching terms in each order of $\\eta$.\"\n4. Thanks for pointing the typos out. We will correct this in the revision.\n\n---\n\n### Questions Part\n\n**Comment**: *\"How can the ... optimizer design?\"*\n\n   **Response**: One interesting observation according to HBF with $\\alpha = 3$ is its implicit preference of low directional smoothness, and this preference can be explicitly added to optimizers as a regularization to improve the performance, since low directional smoothness typically is related to low sharpness and better generalization performance.\n\n---\n\nReference\n\nKovachki and Stuart. Continuous time analysis of momentum methods.\n\nGhosh et al. Implicit regularization in heavy-ball momentum accelerated stochastic gradient descent.\n\nMiyagawa. Toward equation of motion for deep neural networks: Continuous-time gradient descent and discretization error analysis."}, {"review_id": "Fx6Ni6BszE", "reviewer": "Reviewer_iEY4", "summary": "The paper highlights the differences between the discrete and continuous time heavy-ball (HB) momentum methods.\nThe authors explore the continuous-time HB method and report notable findings that distinguish their work from previous studies.\n\n- Research Question\n    \n    - The discretization error between gradient flow (GF) and gradient descent (GD) has been studied in [Barrett and Dherin, 2022; Miyagawa, 2023; Rosca et al., 2023]. Rather than GD, the authors focus on the HB method, which has yet to be explored.\n        \n    - More specifically, the continuous time HB as a surrogate model of the discrete HB has been formalized as a rescaled gradient flow (RGF: Eq. (1)) [Kovachki and Stuart (2020)], but the authors point out that it cannot differentiate gradient descent (GD) and HB (Sections 3 and 4). Another continuous time HB has been proposed by [Ghosh et al. (2023)], which reduces the discretization error compared to RGF, but it is limited to the second order of the step size. [Cattaneo et al. (2023)] has analyzed Momentum and Adam but is also limited to the second order.\n        \n- Solution\n    \n    - To address this issue, the authors propose HB Flow (HBF), a first-order piece-wise continuous differential equation built upon RGF, that approximates the discrete HB momentum method and can control discretization error to *arbitrary order* of the step size (learning rate).\n        \n    - The idea is to add a number of counter terms to RGF order-by-order (Theorem 2.1), akin to backward error analysis [Hairer et a., 2006; Miyagawa, 2023; Rosca et al., 2023].\n        \n- Applications\n    \n    - Examining the directional smoothness and the diagonal linear network, the authors demonstrate that HBF captures the differences between HB and GD, which has been impossible by RGF.\n        \n        - The authors show that the learning dynamics of HB exhibits smaller directional smoothness compared to GD (Fig. 1), through the lens of HBF.\n            \n        - In addition, they find the implicit bias of HB for the diagonal linear network (Theorem 3.1) that cannot be obtained by RGF.\n            \n    \n- Overall, the paper argues that, despite its popularity, continuous-time analysis fails to fully capture the learning dynamics of discrete-time GD methods when neglecting discretization error. Although this limitation has been noted in prior work, such as [Miyagawa, 2023], it is worth emphasizing that the issue is general and extends to various optimization methods such as the HB momentum method.", "questions": "- When does the approximation accuracy of HBF for the discrete HB method deteriorate?\n    \n    - I suspect it occurs in regions where the derivatives of $L$ are large and where $1 - \\mu \\ll \\eta$, causing the Taylor series to diverge.\n        \n    - Does this happen in practical settings?\n        \n    - Could you provide rough bounds on the derivatives required to keep the discretization error small?\n        \n    - It would also be helpful to provide a rough estimate of the $\\alpha$\\-th order discretization error of a $d$\\-dimensional DNN at the $k$\\-th step from the results in Section 2.2, for example.\n        \n- (Eq. (17) & (18)) This question relates to the one above. When $\\mu \\sim 1$, it appears that the dynamics are dominated by higher-order derivatives or even a constant term, rather than the first-order gradient. I am concerned that this does not reflect the correct learning dynamics. Why is that the case?\n    \n- How can we extend the authors’ approach to adaptive learning rate, e.g., Adam? What is the challenge specifically? The authors mentioned the extension to Adam in the Limitation section, and I agree this is an interesting and important step toward practically relevant scenarios in view of the prevalence of Adam.\n    \n    - To my knowledge, existing continuous-time formulations of Adam are quite complex. How might one begin to approach this problem?\n        \n- (Theorem 2.1) Can the result be extended from $\\alpha \\geq 1$ to $\\alpha > 0$?\n    \n- To my understanding, a difference between the paper and [Miyagawa, 2023] is, aside from HB and GF, that the former proposes a piecewise continuous differential equation, whereas the latter provides a globally continuous one. Is it correct?\n    \n- (Appendix C.2) Was $\\eta_\\mathrm{Euler} = \\eta/10$ sufficiently small? The authors can check this, for example, by changing $\\eta/10, \\eta/100,... $ and confirming the learning trajectory does not change dramatically or the change is much smaller than the discretization error, the quantity is of interest.\n    \n- How long did it take to run the experiments?", "limitations": "Yes, but additional discussions can be added: under what conditions the proposed HBF fails to serve as a valid and accurate approximation of the discrete-time HB method?", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "### Strengths\n\nPlease see also the Summary field:\n> - (Novelty) The authors explore the continuous-time heavy-ball (HB) method and report notable findings that distinguish their work from previous studies.\n> - (Novelty & Technical Contribution) The authors show that the learning dynamics of HB exhibits smaller directional smoothness compared to GD (Fig. 1), through the lens of HBF.\n> - (Novelty & Technical Contribution) They find the implicit bias of HB for the diagonal linear network (Theorem 3.1) that cannot be obtained by RGF.\n> - (Broader Impact) The paper argues that, despite its popularity, continuous-time analysis fails to fully capture the learning dynamics of discrete-time GD methods when neglecting discretization error. Although this limitation has been noted in prior work, such as [Miyagawa, 2023], it is worth emphasizing that the issue is general and extends to various optimization methods such as the HB momentum method.\n\n- The proposed HBF might be seen as a natural progression of [Kovachki and Stuart (2020); Miyagawa, 2023; Rosca et al., 2023]. In fact, the the authors’ approach and contents in Section 2 in the paper are similar to  Sections 3 & 4 in [Miyagawa, 2023]. However, I do not think that this is a weakness because:\n    \n    - (Clarity, Novelty, & Technical Contribution) The contribution is clear$\\textemdash$as summarized in the Summary field above$\\textemdash$and significant for the field of optimization, as it helps bridge the gap between conventional continuous-time theory and practical discrete-time methods.\n        \n    - (Technical Contribution) In addition, the derivation of HBF is nontrivial$\\textemdash$performing all-order backward error analysis involves careful calculations, including the use of Lie derivatives$\\textemdash$demonstrating the paper's technical soundness.\n        \n    - (Novelty & Technical Contribution) All analyses in Section 3 and 4 are built upon the proposed HBF and thus are novel, indicating HBF can distinguish GD and HB unlike RGF.\n        \n- (Soundness) The newly-observed implicit regularization of HBF discussed at the end of Section 2 matches our intuition: HB prefers smaller directional smoothness to GD because of the momentum (inertia).\n    \n- (Empirical Validation) (If my concern below is addressed) The experiment in Figure 2 is convincing: larger $\\alpha$ leads to smaller discretization error. Although the model is small, this choice is acceptable given that computing higher-order derivatives becomes prohibitively expensive for large models.\n\n\n### Weaknesses\n\n- (See Questions below) It would be helpful to clarify under what conditions the proposed HBF fails to serve as a valid and accurate approximation of the discrete time HB method. \n\n- (See Questions below) I am concerned about $\\eta_\\mathrm{Euler}$.\n\n- Additionally, I recommend that the authors provide their code for reproduction.\n\n\n## Minor Comments\n\n- [Comment (minor)] I recommend that the authors explicitly refer to $\\gamma_k$ as the *counter term* in line 134. Currently, the term “counter term” appears only in the Abstract, Introduction, and Appendix A.1, which may lead to confusion for readers.\n    \n- [Comment (minor)] The colors in Figure 2(a) are hard to distinguish (three blues).\n\n- [Typo] The numerator of $\\frac{ }{q!}$ in Eq. (50) is missing.\n    \n- [Typo] (Lines 188, 193, 219, and 223) Theorem 3.1 should be Theorem 2.1?\n\n- [Typo] (Line 831) Please Appendix C. -> Please see Appendix C.\n    \n- [Grammatical error] (Line 184) … it is worth to mention… -> … it is worth mentioning…\n    \n- [Comment (minor)] (Lines 192--194)\n    \n    > … the recursive manner in Theorem 3.1 provides a possibility to calculate the involved terms automatically by using software for symbolic mathematics such as SymPy, ….\n    \n    Totally agree $\\textemdash$ combining SymPy with AI coding tools makes tedious analytical math way easier.\n\n\n## Review Summary\nThe paper presents strong technical contributions, offers novel insights, and is clearly written. While it raises a few concerns noted above, I believe the strengths outweigh them. I recommend acceptance and look forward to the discussion with the authors.", "was_revised": true, "final_justification": "The paper presents strong technical contributions, offers novel insights, and is clearly written. \n\n> - (Novelty) The authors explore the continuous-time heavy-ball (HB) method and report notable findings that distinguish their work from previous studies.\n> - (Novelty & Technical Contribution) The authors show that the learning dynamics of HB exhibits smaller directional smoothness compared to GD (Fig. 1), through the lens of HBF.\n> - (Novelty & Technical Contribution) They find the implicit bias of HB for the diagonal linear network (Theorem 3.1) that cannot be obtained by RGF.\n> - (Broader Impact) The paper argues that, despite its popularity, continuous-time analysis fails to fully capture the learning dynamics of discrete-time GD methods when neglecting discretization error. Although this limitation has been noted in prior work, such as [Miyagawa, 2023], it is worth emphasizing that the issue is general and extends to various optimization methods such as the HB momentum method.   \n\nSpecifically,\n> - (Clarity, Novelty, & Technical Contribution) The contribution is clearly stated and significant for the field of optimization, as it helps bridge the gap between conventional continuous-time theory and practical discrete-time methods.      \n> - (Technical Contribution) The derivation of HBF is nontrivial$\\textemdash$performing all-order backward error analysis involves careful calculations, including the use of Lie derivatives$\\textemdash$demonstrating the paper's technical soundness.\n> - (Novelty & Technical Contribution) All analyses in Section 3 and 4 are built upon the proposed HBF and thus are novel, indicating HBF can distinguish GD and HB unlike RGF.        \n> - (Soundness) The newly-observed implicit regularization of HBF discussed at the end of Section 2 matches our intuition: HB prefers smaller directional smoothness to GD because of the momentum (inertia).\n> - (Empirical Validation) The experiment in Figure 2 is convincing: larger $\\alpha$ leads to smaller discretization error. Although the model is small, this choice is acceptable given that computing higher-order derivatives becomes prohibitively expensive for large models.\n\nIn the rebuttal, most of the reviewers’ concerns have been addressed, and no critical issues remain.\nThe authors have provided additional experimental results on the dynamics of discretization error with a larger model and additional theoretical analysis on the bound of discretization error. These results have made the paper further convincing, thereby highlighting the contributions.\nMoreover, the proof in the paper has been improved in the rebuttal, particularly regarding the $\\mathcal{O}$ notation.\n\nTherefore, I support the acceptance of this paper.", "rebuttal": "We thank the reviewer very much for the constructive suggestions and appreciation of our work. Below we answer your questions.\n\n---\n\n## Weaknesses Part\n1. **Comment**: *\"It would be...HB method.\"*\n\n   **Response**: We thank the reviewer for this insightful question. To make the HBF a valid approximation, there are two necessary conditions:\n    - It is crucial to control the ratio between $\\eta $ and $1-\\mu$ to avoid $\\eta\\gg1-\\mu$, which might lead the Taylor series to diverge as suggested by the reviewer. More interestingly, we believe it is the magnitude of a special composite quantity $$\\psi:=\\frac{\\eta}{(1-\\mu)^2}\\tag{a.1}$$that matters for the effectiveness of the HBF. This quantity spontaneously appears in both HBF with $\\alpha=2$ and $\\alpha=3$ but not in the RGF, i.e., for HBF with $\\alpha=2$ the counter term is proportional to $\\psi$ while for HBF with $\\alpha=3$ it is proportional to $\\psi^2$. If $\\psi$ is too large, then our results would no longer hold. Hence, we need to fix the value of $\\mu$ and treat only $\\eta$ as the variable to denote the higher-order terms as $\\mathcal{O}(\\eta^{\\alpha})$ while hide $\\mu$ in the expansion. And it would be interesting for future works to study the case when both $\\mu$ and $\\eta$ are treated as variables such that higher-order terms can be denoted as $\\mathcal{O}(\\psi^{\\alpha})$ (Eq.(a.1)) for $\\alpha\\geq1$. \n\n      In addition, the dependence of HBF on the special composite quantity $\\psi$ is consistent with the empirical observation in Leclerc and Madry (2020), where the optimization curves for different momentum values can be recovered by a corresponding change in the learning rate. The dependence of HB on $\\eta$ and $\\mu$ at the same time further indicates the advantage of HBF with $\\alpha > 1$ and that RGF, which only depends on $\\mu$, is not sufficient to reflect the optimization properties of HB.\n    -  Given $\\alpha$, the continuous approximations include derivatives of $L$ up to the $\\alpha$-th order, hence $L$ should at least be $\\alpha$-times continuously differentiable and $||\\nabla^{\\alpha}L||$ should be upper bounded.\n2. **Comment**: *\"I am concerned...\"*\n\n   **Response**: We thank the reviewer for the suggestion. To validate the influence of the value of $\\eta_{\\text{euler}}$, we run the experiments in Fig.2 for each of $\\eta_{\\text{euler}}\\in\\\\{\\eta/10,\\eta/100,\\eta/1000\\\\}$, and confirm that the observations and conclusions hold in all different $\\eta_{\\text{euler}}$. As we are currently not allowed to upload any materials or links, we will add results of these new experiments in the revision.\n3. **Comment**: *\"Additionally, I...reproduction.\"*\n\n   **Response**: Thanks for the suggestion. We will upload the code, including the newly supplemented  numerical experiments, in the revision.\n### Minor Comments\n- We thank the reviewer for this nice suggestion. We will formally define $\\gamma_k$ as the counter term in the revision.\n- We have changed the colors so that these curves are now more discernible. \n- We thank the reviewer for pointing the typos and grammatical error, and we will correct them accordingly.\n\n---\n\n## Questions Part\n1. **Comment**: *\"When does the ... deteriorate?\"*\n\n   **Response**: We thank the reviewer for the series of interesting questions. We answer these questions point by point below.\n   - We agree with the reviewer that the derivatives of $L$ are not large and that $\\eta\\ll1-\\mu$ are two necessary conditions for HBF to be accurate approximations. Please see our response to the first point of the Weaknesses part.\n   - In practice, using large $\\mu$ and small $\\eta$ simultaneously can still reproduce the results of the case with small $\\mu$ and large $\\eta$, and both parameter sets correspond to a small $\\psi$. This phenomenon has also been observed by Leclerc and Madry (2020). \n\n     However, the regime where both $\\eta$ and $\\mu$  are large (hence large $\\psi$) will lead to failures of HBF due to the divergence of the counter terms, as suggested by the reviewer. In addition, in this regime, the model might not be trained properly either: the update direction coming from the gradient and that from the momentum will jointly affect the training direction significantly, while these two directions can be very different due to the large value of $\\eta$ and $\\mu$ hence cannot give a consistent updating direction.\n   - As the arbitrary order approximations involve complex combinations of higher-order derivatives, it would be intractable to derive an explicit bound. However, under the condition that the Taylor expansion is eligible, we can express the discretization error at the $k$-th step (please see below) with an explicit bound of the second-order derivative.\n   - We thank the reviewer for the nice suggestion. Below we give an estimate of the discretization error of the $\\alpha$-th order continuous approximation at the $k$-th step. We first present several useful relations. As Eq.(43) is established by solving the functional equation for any iteration count $k$, we can write the relation between $\\varepsilon_{k + 1}$ and $\\varepsilon_{k}$ (Eq.(34), line 454-455 in page 13) as $$\\varepsilon_{k+1}-\\varepsilon_{k}=\\mu(\\varepsilon_{k}-\\varepsilon_{k-1})-\\eta [\\nabla L(\\beta(t_k)) - \\nabla L(\\beta(t_k) - \\varepsilon_{k})]+\\mathcal{O}(\\eta^{\\alpha+1}). \\tag{a.2}$$If $\\varepsilon_{k}=\\mathcal{O}(\\eta^{\\alpha})$, then Eq.(a.2) implies $$||\\varepsilon_{k+1}-\\varepsilon_{k}||\\leq\\mu||\\varepsilon_{k}-\\varepsilon_{k-1}||+\\eta\\lambda||\\varepsilon_{k}||+c_1\\eta^{\\alpha+1}$$for some constant $c_1$ where we use $G_k-\\mu G_{k-1}=\\nabla L$ and let $\\lambda=\\max_{\\beta}||\\nabla^2L(\\beta)||$ be the upper bound of the second-order derivative. Denoting $$\\forall k: c_2(k)=\\frac{c_1}{\\lambda}e^{\\frac{2\\lambda\\eta}{1-\\mu}k},\\ c_3(k)=\\frac{2c_1}{1-\\mu}e^{\\frac{2\\lambda \\eta}{1-\\mu}k},$$we now prove by induction. \n     - For the first step ($k = 0$), by definition we have $$\\varepsilon_0=0\\leq c_2(0)\\eta^{\\alpha}.$$Note that $G_{-1}=0$ and $\\gamma_{-1}=0$ by definition, then we have $$||\\varepsilon_1-\\varepsilon_0||\\leq c_1\\eta^{\\alpha+1}\\leq\\frac{2c_1}{1-\\mu}=c_3(0)\\eta^{\\alpha+1}$$since $\\mu<1$.\n     - Suppose that for the $k$-th step $$||\\varepsilon_k||\\leq c_2(k)\\eta^{\\alpha},\\ ||\\varepsilon_{k+1}-\\varepsilon_k||\\leq c_3(k)\\eta^{\\alpha+1}.$$Then for the $(k+1)$-th step, we have $$\\begin{aligned}||\\varepsilon_{k+1}||&\\leq||\\varepsilon_{k+1}-\\varepsilon_k ||+||\\varepsilon_k||\\\\\\\\ & \\leq c_2(k)\\left(1+\\eta\\frac{c_3(k)}{c_2(k)}\\right)\\eta^{\\alpha}\\\\\\\\ &\\leq c_2(k)e^{\\frac{2\\lambda \\eta}{1-\\mu}}\\eta^{\\alpha}=c_2(k+1)\\eta^{\\alpha}\\end{aligned}$$where the last inequality is because $e^x>1+x$ for $x>0$. Similarly, according to Eq.(a.4), $$\\begin{aligned}||\\varepsilon_{k+2}-\\varepsilon_{k+1}||& \\leq\\mu||\\varepsilon_{k+1}-\\varepsilon_{k}||+\\eta\\lambda||\\varepsilon_{k+1}||+c_1\\eta^{\\alpha+1}\\\\\\\\ &\\leq\\left[\\mu c_3(k)+\\lambda c_2(k+1)+c_1\\right]\\eta^{\\alpha + 1}\\\\\\\\ &=\\left[\\mu e^{-\\frac{2\\lambda\\eta}{1-\\mu}}+\\frac{1-\\mu}{2}+\\frac{1-\\mu}{2}e^{-\\frac{2\\lambda\\eta}{1-\\mu}(k +1)}\\right]c_3(k+1)\\eta^{\\alpha+1}\\\\\\\\&\\leq\\left[\\mu+\\frac{1-\\mu}{2}+\\frac{1-\\mu}{2}\\right]c_3(k+1)\\eta^{\\alpha+1}=c_3(k+1)\\eta^{\\alpha+1}.\\end{aligned}$$Hence, for a fixed time $T = N\\eta$ such that $N = \\lfloor\\frac{T}{\\eta}\\rfloor$, at the $k$-th step, the discretization error is bounded by $$||\\varepsilon_k||\\leq \\frac{c_1}{\\lambda}e^{\\frac{2\\lambda}{1-\\mu}k\\eta}\\eta^{\\alpha}.$$\n2. **Comment**: *\"(Eq. (17) (18))... case?\"*\n\n   **Response**: As discussed earlier, as long as $\\psi=\\eta/(1-\\mu)^2$ is small, HBF can still track the learning dynamics. However, if $\\mu\\approx1$ while $\\eta$ is also large, though HBF may not be valid, the training dynamics of the discrete HB is not solely governed by the gradient either, as it will be largely affected by the gradient and the momentum in the current step simultaneously and it might even not converge. Hence, we believe the failure of HBF in such case is acceptable.\n3. We agree with the reviewer that existing continuous-time formulations of Adam are quite complex, hence studying the continuous approximations for Adam should be an important future direction. \n    \n    The main challenge, from our perspective, lies in that Adam incorporates momentum and adaptive learning rate at the same time. According to our approach, there will be two key steps for the problem:\n    - find the equations for the counter term $\\gamma_k$ (similarly, $G_k$ should degenerate to sign gradient, the lowest-order continuous approximation of Adam);\n    - solve the equation obtained in the last step in the series form. \n\n    Hence, mathematically, given the first order ODE $$\\dot{\\beta}=-G_k-\\eta\\gamma_k,$$one difficulty lies in establishing the equations that must be satisfied by $\\gamma_k$ from requiring $\\varepsilon_{k+1}-\\varepsilon_k=\\mathcal{O}(\\eta^{\\alpha+1})$ (this condition will be sufficient for us to show $\\varepsilon_k=\\mathcal{O}(\\eta^{\\alpha})$) for Adam, as the formulations of Adam are more complicated compared to HB. In addition, we believe solving the equations obtained in the last step would also be challenging for Adam. Once these two challenges are resolved, the continuous approximations for Adam can also be established.\n4. As the current approach relies on Taylor expansion to the first order of $\\eta$ with remainder term in the integral form, it cannot be directly generalized to the case of arbitrary $\\alpha>0.$ \n5. Yes. This is because HB utilizes the update history while GD does not.\n6. Please refer to our response to the second point of the Weaknesses part.\n7. We run our experiments on a Linux platform with an Intel(R) Xeon 626 CPU E5-2683 and an NVIDIA Tesla A100.\n   - For Fig.1, the running time is about 1 hour and 54 minutes.\n   - For Fig.2, the running time depends on $\\eta/\\eta_{\\text{euler}}$. For $\\eta/\\eta_{\\text{euler}}=10$, the experiments take about 2 minutes.\n\n---\nReference\n\nLeclerc and Madry. The two regimes of network training, 2020."}]

Author rebuttal:

Meta-review: This paper was accepted because it provides a significant technical contribution by developing a more accurate continuous-time framework for the Heavy-Ball momentum method, which allows for a precise analysis of discretization error. All reviewers agreed that this was a major strength and that the resulting insights into implicit bias were valuable. While there were initial concerns about the clarity of the presentation and the sufficiency of the experiments, the authors provided a thorough rebuttal that satisfied most of these concerns, including more formal proof sketches, and new experiments. The positive response from the reviewers, including two who raised their scores, makes this a clear accept.
Decision: Accept (poster)