Heavy-Ball Momentum Method in Continuous Time and Discretization Error Analysis

We thank the reviewer very much for the constructive suggestions and appreciation of our work. Below we answer your questions.

Weaknesses Part

1. Comment: "It would be...HB method."

   Response: We thank the reviewer for this insightful question. To make the HBF a valid approximation, there are two necessary conditions:

   • 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$.

     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.

   • 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.

2. Comment: "I am concerned..."

   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.

3. Comment: "Additionally, I...reproduction."

   Response: Thanks for the suggestion. We will upload the code, including the newly supplemented numerical experiments, in the revision.

Minor Comments

• We thank the reviewer for this nice suggestion. We will formally define $\gamma_k$ as the counter term in the revision.
• We have changed the colors so that these curves are now more discernible.
• We thank the reviewer for pointing the typos and grammatical error, and we will correct them accordingly.

Questions Part

1. Comment: "When does the ... deteriorate?"

   Response: We thank the reviewer for the series of interesting questions. We answer these questions point by point below.

   • 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.

   • 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).

     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.

   • 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.

   • 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.

     • 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$.
     • 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}.$

2. Comment: "(Eq. (17) (18))... case?"

   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.

3. 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.

   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:

   • find the equations for the counter term $\gamma_k$ (similarly, $G_k$ should degenerate to sign gradient, the lowest-order continuous approximation of Adam);
   • solve the equation obtained in the last step in the series form.

   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.

4. 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.$

5. Yes. This is because HB utilizes the update history while GD does not.

6. Please refer to our response to the second point of the Weaknesses part.

7. We run our experiments on a Linux platform with an Intel(R) Xeon 626 CPU E5-2683 and an NVIDIA Tesla A100.

   • For Fig.1, the running time is about 1 hour and 54 minutes.
   • For Fig.2, the running time depends on $\eta/\eta_{\text{euler}}$. For $\eta/\eta_{\text{euler}}=10$, the experiments take about 2 minutes.

Reference

Leclerc and Madry. The two regimes of network training, 2020.

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:

1. Comment: "The numerical experiments ... in Miyagawa (2023)."

   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 $`Linear`(784\times128)\to `SiLU` \to `Linear`(128\times128) \to `BatchNormalization`\to `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.

2. Comment: "The mathematical ... are limited."

   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.

   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.

   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.

3. Comment: "I have some difficulties ... well-organized."

   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.

   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, \dot{\beta} = - G_k - \eta \gamma_k, \tag{a.1} such that the discretization error can be controlled arbitrarily low, where the design of the ODE Eq.(a.1) follows two basic principles:

   • 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)$.
   • 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.

   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:

   • 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.
   • 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$."

4. Thanks for pointing the typos out. We will correct this in the revision.

Questions Part

Comment: "How can the ... optimizer design?"

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.

Reference

Kovachki and Stuart. Continuous time analysis of momentum methods.

Ghosh et al. Implicit regularization in heavy-ball momentum accelerated stochastic gradient descent.

Miyagawa. Toward equation of motion for deep neural networks: Continuous-time gradient descent and discretization error analysis.

We thank the reviewer very much for the valuable comments. Below we address your concerns.

Our Improvements

We would like to first highlight that we have made a series of important improvements in this version:

1. 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.
2. We have improved the statement of Theorem 2.1, where now we:
   • first state what the formulation of the proposed ODE with discretization error to arbitrary order of $\eta$ will be;
   • then state what the explicit formulations of the components of the ODE will be by defining necessary notations.
3. 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.

Weaknesses Part

Below we address your concerns mentioned in the Weaknesses part.

1. Comment: "I found the main theorem ... mathematically sophisticated."

   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.

   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.

   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.

   Intuition: Overall, for the discrete HB method, our goal is to find a counter term $\gamma_k$ for the modified ODE proposed by us $\dot{\beta} = - G_k - \eta \gamma_k,$ 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:

   • 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.
   • 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$.

2. Comment: "Some notations ... Similar situation for $\gamma_k$"

   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.

Questions Part

1. Comment: "In Figure ..."

   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)).

2. Comment: "This analysis can be generalized ... future direction."

   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.

We thank the reviewer very much for the insightful comments. We appreciate this opportunity to address your concerns.

Weaknesses Part

Weaknesses 1

We address concerns in Weaknesses 1 point by point.

• Comment: "While I...Appendix."

  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$.

• Comment: “First of...term of $\gamma$”

  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.

  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:

  \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\\\\ \boxed{\text{Step 5}}\qquad\qquad\quad&\text{Solution of }\gamma_k \end{aligned}$$ Specifically, - 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) - 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}$. - 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). - 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 $$ I_k^++\mu I_k^--\gamma_k+\mu\gamma_{k-1}=O(\eta^{\alpha-1})\tag{a.2} $$ in Eq.(34) (line 454-455, page 13). - 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)}, $$ then Lemma A.1 provides us the corresponding $I_k^{\pm}$, which is also in the series form: $$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., $$\gamma_k = \sum_{\sigma}^{\alpha - 2} \eta^{\sigma} \gamma_k^{(\sigma)},$$ guarantees Eq.(a.2).

• Comment: “Besides,...calculation.”

  Response: We first confirm that Eq.(47) still introduces $n$-th order differential calculation. The advantage of Eq.(47) is as follows.

  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.

• Comment: "In addition,...equation (44)?"

  Response: We first confirm that the domain of $G_k$ is the entire space. Below we answer the questions separately.

  • 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).

  • 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)$.

  • 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$.

Weaknesses 2

• Comment: “The use...proof.”

  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.

  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.

  • 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, 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$.
  • 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}.$

  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.

Finally, 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.

Questions Part

We thank the reviewer for suggesting these nice related works. We will discuss these related works in the revision.