A New Accelerated Gradient Method Inspired by Continuous- Time Perspective
Year: 2020

Abstract:
Nesterov’s accelerated method are widely used in problems with machine learning background including deep learning. To give more insight about the acceleration phenomenon, an ordinary differential equation was obtained from Nesterov’s accelerated method by taking step sizes approaching zero, and the relationship between Nesterov’s method and the differential equation is still of research interest. In this work, we give the precise order of the iterations of Nesterov’s accelerated method converging to the solution of derived differential equation as step sizes go to zero. We then present a new accelerated method with higher order. The new method is more stable than ordinary method for large step size and converges faster. We further apply the new method to matrix completion problem and show its better performance through numerical experiments.

heading: 1 Introduction
text: Optimization is a core component of statistic and machine learning problems. Recently, gradientbased algorithms are widely used in such optimization problems due to its simplicity and efficiency for large-scale situations. For solving convex optimization problem
min x∈Rd F (x),
where F (x) is convex and sufficiently smooth, a classical first-order method is gradient descent. We assume that f(x) = ∇F (x) satisfies L-Lipschitz condition, that is, there exists constant L such that ∥f(x)− f(y)∥ ≤ L∥x− y∥, ∀x, y. Under these conditions, gradient descent achieves a convergence rate of O(n−1), i.e., ∥F (xn) − F (x∗)∥ decreases to zero at a rate of O(n−1), where xn denotes the nth iteration and x∗ denotes the minimum point of F (x) in Rd. Nesterov’s accelerated method (Nesterov, 1983) is a more efficient first-order algorithm than gradient descent, of which we will use the following form: starting with x0 = x1,
yn = xn + n− 3 n (xn − xn−1),
xn+1 = yn − sf(yn) (1.1)
for n ≥ 1. It is shown that under abovementioned conditions, Nesterov’s accelerated method converges at a rate of O(n−2). Accelerated gradient method has been successful in training deep and recurrent neural networks (Sutskever et al., 2013) and is widely used in problems with machine learning background to avoid
sophisticated second-order methods (Cotter et al., 2011; Hu et al., 2009; Ji & Ye, 2009). To provide more theorical understanding, an important research topic of Nesterov’s accelerated method is to find an explanation of the acceleration. On this topic, Nesterov’s method was studied via a continuous-time perspective (Su et al., 2014). They considered a curve x(t), introduced the ansatz xn ≈ x(n √ s) and substituted it to (1.1). Letting s → 0, they obtained the following differential equation. ẍ+ 3
t ẋ+ f(x) = 0. (1.2)
The differential equation was used as a tool for analyzing and generalizing Nesterov’s scheme. Furthermore, this idea has been studied from different directions. A class of accelerated methods have been generated in continuous-time (Wibisono et al., 2016). ODE (1.2) can also be discretized directly using Runge-Kutta method to achieve acceleration (Zhang et al., 2018). Although many results have been achieved, the process of obtaining the differential equation (1.2) has not been rigorous, and the method is still time-consuming for large-scale problems. In this work, we give the precise order of the iterations of Nesterov’s accelerated method converging to solution of the differential equation (1.2) with initial conditions
x(0) = x0, ẋ(0) = 0 (1.3) as step size s goes to zero. Inspired from this perspective, we present a new accelerated method to make this convergence faster. As we expected, iterations of the new method are closer to the solution x(t) of differential equation (1.2) than original Nesterov’s method. Moreover, we find the new method is more stable than original Nesterov’s method when step size is large. Based on abovementioned observations, we try to take advantage of the new method in more practical problems. We apply the new method to matrix completion problem. We combine the new method with proximal operator (Parikh & Boyd, 2014) into a new algorithm, which we call modified FISTA. We find that the new method performs better than FISTA (Beck & Teboulle, 2009) and acclerated proximal gradient method (Parikh & Boyd, 2014) because it can work with larger step sizes. This paper is organized as follows. In section 2, we prove that iterations of Nesterov’s accelerated method converge to solution of the differential equation (1.2). In section 3, we present a new method to make the convergence faster and show its better stablity through two simple examples. In section 4, we apply the new method to matrix completion problem.

heading: 2 A strict analysis of the relation between Nesterov’s method and its continuous-time limit
text: We refer to x(t) as the solution of differential equation (1.2) with initial conditions (1.3). Existance and uniqueness of such solutions have been proved (Su et al., 2014). In this section, We give the order of the iterations of Nesterov’s accelerated method converging to x(t) as step sizes go to zero. For convenience, we substitute the first equation in Nesterov’s method (1.1) to the second one to get
xn+1 = xn + n− 3 n
(xn − xn−1)− s · f ( xn +
n− 3 n
(xn − xn−1) ) .
We write s = h2 and rewrite the above recurrence relation as
xn+1 = xn + n− 3 n
(xn − xn−1)− h2 · f ( xn +
n− 3 n
(xn − xn−1) ) . (2.1)
Inspired by the ansatz xn ≈ x(n √ s) (Su et al., 2014), we consider the convergence between xn and x(nh). More precisely, we show that for fixed time t, xn converges to x(t) as h goes to zero, where n = th .
2.1 Truncation error
Firstly, we consider the following ‘truncation error’:
L[x(t);h] =x(t+ h)− 2t− 3h t x(t) + t− 3h t x(t− h)+
h2f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) .
(2.2)
(2.2) is obtained from (2.1) by replacing xn+1, xn, xn−1 with x(t+h), x(t), x(t−h) and substituting the relation n = th . Our first result is the order of truncation error L[x(t);h]. Theorem 1. Assume f satisfies L-Lipschitz condition, and solution x(t) of the derived differential equation (1.2) has a continuous third derivative. For fixed time t, the truncation error (2.2) satisfies
L[x(t);h] = O(h3).
Theorem 1 shows the size of error caused by a single iteration when the starting point is just on x(t). Then we have to add up these errors to prove the convergence proporty we need.
2.2 Convergence theorem
We now come to the convergence theorem. In this theorem, we give the precise order of the iterations of Nesterov’s method converging to solution of the derived differential equation. Theorem 2. Under conditions in Theorem 1, for fixed time t, xt/h converges to x(t) as h goes to zero at a rate of O(h ln 1h ) if x0 = x(0) and x1 = x(h).
Theorem 2 coincides with derivation of ODE (1.2) (Su et al., 2014).

heading: 3 New accelerated method
text: 3.1 Derivation of the new method and analysis of truncation error
Inspired from the continuous-time perspective and our proof of the convergence from iterations of Nesterov’s method to its continuous-time limit, we present a new method to make this convergence faster. Precisely, the new method has a higher truncation order. We need one more step in our scheme than in Nesterov’s method to achieve higher truncation order in the following analysis, so we consider a recurrence relation with form
4∑ i=1 ( αi + βi n + γi n2 ) xn+2−i = −sf ( xn + n− 3 n (xn − xn−1) ) , (3.1)
where {αi}, {βi} and {γi} are to be determined.
Now we expand x(t− h) to first order. Calculation shows that
f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) =− hx(3)(t)− ( 3h
t + 1
) x(2)(t)
+
( 3h
t2 − 3 t
) x(1)(t) +O(h2).
Substitute this expansion to truncation error
L[x(t);h] = 4∑ i=1 ( αi + βih t + γih 2 t2 ) x(t+ (2− i)h)
+ h2f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) ,
and choose parameters appropriately to eliminate low-order terms, we get the following recurrence relation
xn+1 = 10n2 + 9n+ 6
4n2 + 8n xn −
4n2 + 3
2n2 + 4n xn−1 + 2n− 1 4n+ 8 xn−2
− n 2n+ 4 sf
( 2n− 3
n xn − n− 3 n xn−1
) .
(3.2)
Here we rewrite this scheme as Algorithm 1.
Algorithm 1 The new method (3.2) Input: step size s Initial value: X2 = X1 = X0. (k− 1)th iteration (k ≥ 2). Compute
Yk = 10k2 + 9k + 6
4k2 + 8k Xk −
4k2 + 3
2k2 + 4k Xk−1 + 2k − 1 4k + 8 Xk−2,
Zk = 2k − 3
k Xk − k − 3 k Xk−1,
Xk+1 = X − ( Yk − ks
2k + 4 f(Zk)
) .
For truncation order of this new method, we have the following theorem. The abovementioned procedure is presented in Appendix A.4 detailedly, as proof of Theorem 3. Theorem 3. If f has continuous second order derivative, the first and second derivative are bounded, and x(t) has continuous fourth derivative, then for fixed t, truncation error of (3.2) satisfies
L[x(tn);h] = O(h4).
The convergence of the new method and x(t) can be proved similar to Theorem 2.
3.2 Advantage of the new method
Since the new method has a truncation error of higher order than original Neaterov’s method, the iterations of the new method converge to the differential equation (1.2) when those of original
Nesterov’s method diverge. In another word, the new method is more stable for large step size. We present two numerical results in Figure 1 to confirm it.
Quadratic. F (x) = xTAx is a strongly convex function, in which x ∈ R2 and A is a 2× 2 matrix.
Linear regression. F (x) = n∑
i=1
(yi −wTi x)2, where n is the number of samples and (wi, yi) is the
ith sample.
In these examples, iterations of the new method converge to the minimum point, while those of original Nesterov’s method diverge, which confirms that the new method is more stable for large step size.
3.3 Absolute stability of Nesterov’s method and the new method
In this subsection, we explain the better stability of the new method with absolute stability theory. Firstly, recall the scheme of our new method
xn+1 = 10n2 + 9n+ 6
4n2 + 8n xn −
4n2 + 3
2n2 + 4n xn−1 + 2n− 1 4n+ 8 xn−2
− n 2n+ 4 sf
( 2n− 3
n xn − n− 3 n xn−1
) .
We use linear approximation
f ( xn +
n− 3 n
(xn − xn−1) ) = ∇F ( xn +
n− 3 n
(xn − xn−1) ) ≈ ∇2F · ( xn +
n− 3 n
(xn − xn−1) ) ,
and the characteristic equation of this finite scheme is approximately λ3− ( 10n2 + 9n+ 6
4n2 + 8n − s · ∇2F · 2n 2 − 3n 2n2 + 4n
) λ2+ ( 4n2 + 3
2n2 + 4n − s · ∇2F · n 2 − 3n 2n2 + 4n
) λ− 2n− 1
4n+ 8 = 0.
For large n, we can ignore the high order terms and the characteristic equation becomes λ3 − ( 5
2 − s · ∇2F ·
) λ2 + ( 2− s
2 · ∇2F
) λ− 1
2 = 0.
According to the absolute stability theory, the numerical stability of Nesterov’s scheme with respect to accumulated roundoff error is equivalent to this: all the roots of the characteristic equation lie in the unit circle (Leader, 2004). Noticing that the left hand of the equation can be factorized to(
λ− 1 2
)( λ2 − (2− s · ∇2F )λ+ 1 ) ,
the largest modulu of the roots is 1 when 0 ≤ s · ∇2F ≤ 4, and the absolutely stable region of the new method is s · ∇2F ∈ [0, 4]. When s·∇2F lies in the absoletely stable region, the related theory guarantees that the error caused by every iteration will not be magnified as the iteration number increases. To make the analysis more precise, we should consider the difference of the scheme between iterations caused by different n. We define the transfer matrix
Pn =  ( 10n2+9n+6 4n2+8n − s · ∇ 2F · 2n 2−3n 2n2+4n ) − ( 4n2+3 2n2+4n − s · ∇ 2F · n 2−3n 2n2+4n ) 2n−1 4n+8
1 0 0 0 1 0  and Qn = PnPn−1 · · ·P1. Error analysis shows that if the largest modulu of eigenvalues of Qn goes to zero, then error caused by iterations will be eliminated as the iteration number increases. Figure 2 presents the largest module of eigenvalues of Qn for different values of s · ∇2F . From the experiment we can see that the above condition is satisfied.
We then apply the same method to Nesterov’s method discussed in (Su et al., 2014) and conclude that the absolutely stable region of Nesterov’s method is [0, 43 ]. According to the above analysis, the absolutely stable region of the new method is four times as large as Nesterov’s method, so the new method is more stable, and we can choose larger step sizes to achieve faster convergence.

heading: 4 Application to matrix completion problem: modified FISTA
text: Our theory and numerical results show that the new method is more stable than original Nestrov’s method. So we can choose larger step size for new method and convergence to the optimal solution
can be faster, compared with original Nesterov’s method. In this section we apply the new method to matrix completion problem. We present a new algorithm which can be viewed as a modification of the well-konwn fast iterative shrinkage-thresholding algorithm (FISTA) (Beck & Teboulle, 2009). The performance of modified FISTA can also confirm the advantage of the new method. For matrix completion problem there exists a ‘true’ low rank matrix M . We are given some entries of M and asked to fill missing entries. There have been various algorithms to solve such problem (Candès & Recht, 2009; Keshavan et al., 2010). Besides, it is proposed that matrix completion can be transformed to the following unconstrained optimization problem (Mazumder et al., 2010)
minF (X) = 1
2 ∥Xobs −Mobs∥2 + λ∥X∥∗. (4.1)
Notice that F (X) is composed of a smooth term and a non-smooth term, so gradient-based algorithms cannot be used directly. Proximal gradient algorithms (Parikh & Boyd, 2014) are widely used in such composite optimization problems, and fast iterative shrinkage-thresholding algorithm (FISTA) is a successful algorithm. Moreover, FISTA has been extended to matrix completion case (Ji & Ye, 2009). For convenience, we set G(X) = 12∥Xobs − Mobs∥2, H(X) = λ∥X∥∗, and g(X) = ∇G(X). The idea of FISTA builds on Nesterov’s method. We also apply acclerated proximal gradient method (Parikh & Boyd, 2014) for our numerical experiment, which is composed of Nesterov’s method and proximal gradient descent. These two algorithms are presented in Appendix A.5. We find the performances of them are similar in our experiments. Our contribution is the third method (Algorithm 2), the new method (3.2) combined with proximal operator, which we call modified FISTA.
Algorithm 2 Modified FISTA Input: step size s Initial value: X2 = X1 = X0 ∈ M100. (k− 1)th iteration (k ≥ 2). Compute
Yk = 10k2 + 9k + 6
4k2 + 8k Xk −
4k2 + 3
2k2 + 4k Xk−1 + 2k − 1 4k + 8 Xk−2,
Zk = 2k − 3
k Xk − k − 3 k Xk−1,
Xk+1 = argmin X
{ 1
2 · 2k + 4 ks ∥∥∥∥X − (Yk − ks2k + 4g(Zk) )∥∥∥∥2 + λ∥X∥∗ } .
Notice that the minimizing problems in interations of above three algorithms can be solved directly by singular value decomposition (Cai & Candès, 2010). We take experiments on a simulated data set. We use fixed step sizes in the above three algorithms, and the performances are presented in Figure 3. We find empirically that for all methods, convergence is faster when step size is larger, so we choose the largest step sizes for all methods to compare their fastest convergence speed. Through experiments, we find the largest step size that makes modified FISTA convergent is 4.1 (accurate to one decimal place), while those for the first two algorithms are both 1.3. We also compare performances of the three methods with step sizes reduced from the largest in equal proportion. We find that when step sizes are chosen to be the largest or reduced from the largest in equal proportion (80%, 50%, 10%), convergence of modified FISTA is faster than the other two methods.
We also combine the three methods with backtracking (Beck & Teboulle, 2009) to choose step sizes automatically. We present modified FISTA with backtracking as Algorithm 3, and the other two algorithms are similar. Performances of the three algorithms with backtracking on abovementioned data set are presented in Figure 4. Convergence of modified FISTA is faster than the other two methods. Moreover, we find that the final step size of modified FISTA is larger.

heading: 5 Discussion
text: In this paper we prove that iterations of Nesterov’s accelerated method converge to solution of the derived differential equation as step sizes go to zero. We present a new accelerated method to make this convergence faster. We use numerical results to show that the new method is more stable, especially for large step sizes, and explan it using the order of truncation error. We then apply the new method to matrix completion problem and present a new algorithm which we call modified FISTA. Numerical experiments show that modified FISTA performs better than existing algorithms based on Nesterov’s acceleration because it can work with larger step sizes. We will also combine our new method with stochastic gradient-based algorithms and apply the new method to deep networks in the future.
Algorithm 3 Modified FISTA with backtracking Input: some β < 1 Initial value. X2 = X1 = X0 ∈ M100, step size s2. (k− 1)th iteration (k ≥ 2).
Yk = 10k2 + 9k + 6
4k2 + 8k Xk −
4k2 + 3
2k2 + 4k Xk−1 + 2k − 1 4k + 8 Xk−2,
Zk = 2k − 3
k Xk − k − 3 k Xk−1.
Find the smallest positive integer ik+1 such that with s = βik+1sk F (X̃) < F (Yk) + ⟨ X̃ − Yk, g(Zk) ⟩ + 1
2 · 2k + 4 ks ∥X̃ − Yk∥2,
where
X̃ = argmin X
{ 1
2 · 2k + 4 ks ∥∥∥∥X − (Yk − ks2k + 4g(Zk) )∥∥∥∥2 + λ∥X∥∗ } .
Set sk+1 = βik+1sk and compute
Xk+1 = X̃.
Our work shows that for an accelerated gradient method, the rate at which it converges to the derived differential equation is possibly related to its property as an optimization algorithm. We think this work suggests that more consideration should be given to the corresponding differential equations when studying optimization algorithms.

heading: A Appendix
text: A.1 Proof of Theorem 1
Theorem 1. Assume f satisfies L-Lipschitz condition, and solution x(t) of the derived differential equation (1.2) has a continuous third derivative. For fixed time t, the truncation error (2.2) satisfies
L[x(t);h] = O(h3). (A.1)
Proof. Notice that x(t− h) = x(t) +O(h).
Substitute this equation to the last term of L[x(t);h] to get
f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) = f ( x(t) +
t− 3h t
· O(h) ) .
Since f satisfies L-Lipschitz condition, we know
f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) =f(x(t)) +O(h)
=− ẍ(t)− 3 t ẋ(t) +O(h).
To get the second equality, we substitute the differential equation (1.2). Then we expend the first and third terms of L[x(t);h] to third order to get
x(t+ h) = x(t) + hx(1)(t) + h2
2 x(2)(t) +O(h3),
x(t− h) = x(t)− hx(1)(t) + h 2
2 x(2)(t) +O(h3).
Substitute these three equations to (2.2), we have L[x(t);h] = O(h3).
Remark 1. (A.1) can also be written as |L[x(t);h]| ≤ M1h3,
where M1 depends on sups≤t |x(1)(s)| and sups≤t |x(3)(s)|. Remark 2. Theorem 1 deals with the problem for fixed time t. To finish the proof of the convergence, we have to consider the situation that tn = nh, where n ≥ 1 is fixed. We set a fixed time t0 and assume that tn = nh < t0. Since x(t) has a continuous third derivative, x(t) and its first to third derivative are bounded in [0, t0]. We replace time t in the above proof by tn and expend the terms of (2.2). Now the term
−3h 3
2tn x(2)(tn)
obtained from the expansion of x(tn−1) cannot be viewed as O(h3), but there exists M2 > 0 such that ∣∣∣∣−3h32tn x(2)(tn) ∣∣∣∣ ≤ M2h2n . As a consequence, we have
|L[x(tn);h]| ≤ M1h3 +M2 h2
n , (A.2)
where M1 and M2 rely on t0.
A.2 Two lemmas for Theorem 2
For the proof of Theorem 2, we need the following two lemmas. Lemma 1. (Holte, 2009) For constant α, β > 0 and positive sequence {ηn}n≥0 satisfying
ηn ≤ β + α n−1∑ i=0 ηi, ∀n > 0,
the following inequality holds ηn ≤ eαn(β + αη0).
The above lemma is a classic result and refered to as discrete Gronwall inequality. Lemma 2. We define matrices Cn and Dn,l as
Cn =
( 2n−1 n+1 − n−2 n+1
1 0
) ,
Dn,l = CnCn−1 · · ·Cn−l+1,
where n ≥ 0 and 0 < l ≤ n+ 1. In addition, we set Dn,0 = I2. Then there exist positive constants M, M3 such that for all n, the following two inequalities hold, where the matrix norm is 2-norm.
sup 0≤l≤n+1
∥Dn,l∥ ≤ Mn,
Dn,n+1 ≤ M3. (A.3)
Proof. Since C2 = ( 1 0 1 0 ) ,
we notice that when n ≥ 2,
Dn,n−1 = ( 1 0 1 0 ) , Dn,n = ( 1 2
1 2
1 2
1 2
) , Dn,n+1 = ( 0 1 0 1 ) ,
having nothing to do with the value of n. So it is obvious that there exists M3 to make (A.3) holds and M4 > 0 such that for all n < 2 or n ≥ 2, l > n− 2 or l = 0,
∥Dn,l∥ ≤ M4n. (A.4)
Then we consider the condition when n ≥ 2, 0 < l ≤ n− 2. Notice that
Cn = ( 1 1 1 0 )( 1 1 0 n−2n+1 )( 1 1 1 0 )−1 .
For convenience, we write P = ( 1 1 1 0 ) .
Assume we have alreagy got Dn,l = P ( 1 an,l 0 bn,l ) P−1
satisfying 0 < an,l ≤ l, 0 < bn,l ≤ 1,
then since Dn,l+1 = Dn,lCn−l,
and 0 ≤ n−l−2n−l+1 < 1, Dn,l+1 has the same form
Dn,l+1 = P ( 1 an,l+1 0 bn,l+1 ) P−1,
satisfying 0 < an,l+1 ≤ l + 1, 0 < bn,l ≤ 1.
Then for fixed n, induce from l = 1, we get
Dn,l = PD̃n,lP −1 ≜ P ( 1 an,l 0 bn,l ) P−1,
satisfying 0 < an,l ≤ l ≤ n, 0 < bn,l ≤ 1, (A.5)
for all n ≥ 2, 0 < l ≤ n− 2. Then we can estimate ∥Dn,l∥. Notice that
D̃n,lD̃ T n,l = ( 1 + a2n,l an,lbn,l an,lbn,l a 2 n,l ) .
The eigenvalues of this matrix are
λ1,2 = 1 + a2n,l + b 2 n,l ±
√ (1 + a2n,l + b 2 n,l) 2 − 4b4
2 .
Combining this representation with (A.5), we get the estimation ∥D̃n,l∥ = √ max{|λ1|, |λ2|} ≤ √ 1 + a2n,l + b 2 n,l ≤ n+ 2.
So there exists M5 > 0, such that for all n ≥ 2, 0 < l ≤ n− 2, inequality ∥Dn,l∥ ≤ M5n (A.6)
holds. Combining (A.4) with (A.6), we finish the proof.
A.3 Proof of Theorem 2
Theorem 2. Under conditions in Theorem 1, for fixed time t, xt/h converges to x(t) as h goes to zero at a rate of O(h ln 1h ) if x0 = x(0) and x1 = x(h).
Proof. In this proof, we first calculate the error caused by a single iteration, which can be divided into an accumulation term and a truncation term. Then we use the estimation given by Theorem 1 and apply discrete Gronwall inequality to prove the convergence. Recall the recurrence relation
xn+1 = xn + n− 3 n
(xn − xn−1)− h2 · f ( xn +
n− 3 n
(xn − xn−1) )
and the definition of truncation error
x(tn+1) = x(tn) + n− 3 n
(x(tn)− x(tn−1))− h2f ( x(tn) +
n− 3 n
(x(tn)− x(tn−1)) ) + L[x(tn);h],
where tn = nh. Subtract the above two equations, and introduce overall error
en = x(tn)− xn,
we have en+1 =
2n− 3 n en − n− 3 n en−1 − h2bn−1 + L[x(tn);h],
which can also be written as
en+2 − 2n− 1 n+ 1 en+1 + n− 2 n+ 1 en = −h2bn + L[x(tn+1);h], (A.7)
where bn = f ( 2n− 1 n+ 1 xn+1 − n− 2 n+ 1 xn ) − f ( 2n− 1 n+ 1 x(tn+1)− n− 2 n+ 1 x(tn) ) . (A.8)
We will also use the notation
b∗n = − en+2 − 2n−1n+1 en+1 + n−2 n+1en
h2 .
Then we rewrite (A.7) into a form that is convenient for recurrence. We set
En = ( en+1 en ) , Cn = ( 2n−1 n+1 − n−2 n+1 1 0 ) , Bn = ( −h2b∗n 0 ) .
Then (A.7) can be written as En+1 = CnEn +Bn.
By recursive method, we have
En = Cn−1 · · ·C0E0 + n∑
l=1
Cn−1 · · ·Cn−l+1Bn−l.
With the notations introduced in Lemma 2, this equation can be written as
En = Dn−1,nE0 + n∑ l=1 Dn−1,l−1Bn−l. (A.9)
Now we need to estimate ∥Bn∥. Since f satisfies L-Lipschitz condition, from (A.8) we have |bn| ≤ L ( 2n− 1 n+ 1 |en+1|+ n− 2 n+ 1 |en| ) ≤ L (2|en+1|+ |en|) ≤ 3L∥En∥.
and ∥Bn∥ ≤ 3h2L∥En∥+ L[x(tn+1);h]. (A.10)
Take norm on both sides of (A.9) and substitute (A.10) and conclusion of Lemma 2, we have the following estimation
∥En∥ ≤ M3∥E0∥+M(n− 1) n−1∑ l=0 ( 3h2L∥El∥+ L[x(tl+1);h] ) ≤ M3∥E0∥+ 3Mnh2L
n−1∑ l=0 ∥El∥+Mn n−1∑ l=0 L[x(tl+1);h].
(A.11)
Now we deal with truncation errors. Recall (A.2) in remark of Theorem 1
|L[x(tl);h]| ≤ M1h3 +M2 h2
l .
Take sum to obtain n−1∑ l=0 |L[x(tl+1);h]| ≤ nM1h3 +M2h2 n−1∑ l=0 1 l + 1 . (A.12)
Notice the classic inequality n∑
i=1
1 i ≤ lnn+Me,
where Me refers to a positive constant. Substitute it to (A.12), we have n−1∑ l=0 |L[x(tl+1);h]| ≤ nM1h3 +M2h2(lnn+Me).
Substitute this inequality to (A.11), we get a control of ∥En∥
∥En∥ ≤ M3∥E0∥+ 3Mnh2L n−1∑ l=0 ∥El∥+MM1n2h3 +MM2Menh2 +MM2nh2 lnn
Using discrete Gronwall inequality, we have
∥En∥ ≤ e3Mn 2h2L ( M3∥E0∥+MM1n2h3 +MM2Menh2 +MM2nh2 lnn+ 3Mnh2L∥E0∥ ) .
Then for fixed t, we choose n = th to get
∥Et/h∥ ≤ e3Mt 2L ( (M3 + 3MthL)∥E0∥+ (MM1t2 +MM2Met)h+MM2th ln t
h
) .
Notice that lim h→0 h ln t h = 0,
so if E0 = 0, then the vector form of overall error Et/h satisfies
lim h→0
∥Et/h∥ = 0.
A.4 Proof of Theorem 3
Theorem 3. If f has continuous second order derivative, the first and second derivative are bounded, and x(t) has continuous fourth derivative, then for fixed t, truncation error of (3.2) satisfies
L[x(t);h] = O(h4).
Proof. Recall the proof of Throrem 1. Now we expand x(t− h) to first order
x(t− h) = x(t) + hx(1)(t) +O(h2).
Then we have
f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) =f ( x(t) + ( 1− 3h
t
) (hx(1)(t) +O(h2)) ) =f ( x(t) + hx(1)(t) +O(h2)
) =f ( x(t) + hx(1)(t) ) +O(h2).
We now expand f :
f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) = f(x(t)) + hx(1)(t)f (1)(x(t)) +O(h2).
To do this, we need f has continuous second derivative and the second derivative is bounded. Take derivetive on both sides of differential equation
ẍ+ 3
t ẋ+ f(x) = 0,
we have (f(x(t))) ′ = −x(3)(t)− 3
t x(2)(t) +
3 t2 x(1)(t).
So
f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) =− hx(3)(t)− ( 3h
t + 1
) x(2)(t)
+
( 3h
t2 − 3 t
) x(1)(t) +O(h2).
(A.13)
Expand x(t+ h), x(t− h), x(t− 2h) to the third order, we have( α1 + β1h
t +
γ1h 2
t2
) x(t+ h) = ( α1 + β1h
t +
γ1h 2
t2 ) [ x(t) + hx(1)(t) + h2
2 x(2)(t) +
h3
6 x(3)(t) +O(h4) ] ,(
α3 + β3h
t +
γ3h 2
t2
) x(t− h) = ( α3 + β3h
t +
γ3h 2
t2 ) [ x(t)− hx(1)(t) + h 2
2 x(2)(t)− h
3
6 x(3)(t) +O(h4) ] ,(
α4 + β4h
t +
γ4h 2
t2
) x(t− 2h) = ( α4 + β4h
t +
γ4h 2
t2 ) [ x(t)− 2hx(1)(t) + 2h2x(2)(t)− 4h 3
3 x(3)(t) +O(h4)
] .
Substitute these three equations and (A.13) to truncation error of recurrence relation (3.1)
L[x(t);h] = 4∑ i=1 ( αi + βih t + γih 2 t2 ) x(t+ (2− i)h)
+ h2f ( x(t) +
t− 3h t
(x(t)− x(t− h)) ) ,
then simple calculation shows that terms with order less than four will be eliminated if we choose coefficients according to the following equations
 α1 = 2 α2 = −5 α3 = 4
α4 = −1
,
 β1 = 9 2 − k β2 = −6 + 3k β3 = 3 2 − 3k
β4 = k
,
 γ1 = m1 γ2 = − 3m1 +m2 + 3 2 γ3 = m2
γ4 = m1 −m2 + 3
2
,
where k, m1, m2 can be chosen randomly. Notice that coefficients of recurrence relation (3.2) satisfy above equations.
A.5 Algorithms
Algorithm 4 FISTA Input: step size s Initial value: Y1 = X0 ∈ M100, t1 = 1. kth iteration (k ≥ 1). Compute
Xk = argmin X
{ 1
2s ∥X − (Yk − sg(Yk))∥2 + λ∥X∥∗
} ,
tk+1 = 1 + √ 1 + 4t2k 2 , Yk+1 = Xk + tk − 1 tk+1 (Xk −Xk−1).
Algorithm 5 Accelerated proximal gradient method Input: step size s Initial value: X1 = X0 ∈ M100. kth iteration (k ≥ 1). Compute
Yk = Xk−1 + k − 3 k (Xk−1 −Xk−2),
Xk+1 = argmin X
{ 1
2s ∥X − (Yk − sg(Yk))∥2 + λ∥X∥∗
} .
A.6 Details about Numerical Experiments in Section 4
Here we produce some details for our numerical experiments in Section 4. Our experiments are taken on a simulated data set. Firstly, we generated the ‘true’ low rank matrix M . To do this, we generate a random matrix M0. Entries of M0 are independent and uniformly distributed on (0, 20). Then we compute the singular value decomposition of M0, that is, M0 = UΣV T. After that, we set M = UΣ0V T, where Σ0 is a diagonal matrix with only three nonzero diagonal elements. It is not difficult to prove that M has rank 3.
Secondly, we generate the observation set. For every row of M , we choose randomly ten entrys to be observed. As a consequence, 10% entries are observed in total. After data generation step, we apply the abovementioned algorithms (accelerated proximal gradient method, FISTA and our modified FISTA) with fixed step sizes and backtracking to this data set. The parameter of the loss function (4.1) is λ = 0.005. For initial point, we simply choose the zero matrix (every entry equals to zero). For backtracking, we set the initial step size as 10 and the decay factor β = 0.1.

Review summaries:
This paper proposes an accelerated method that has a high-order truncation error $O(h^4)$ to the ordinary differential equation $\ddot{x} + \frac{3}{t}\dot{x} + f(x) = 0$ obtained from Nesterov's accelerated method by (Su et al., 2014), while Nesterov's method has $O(h^3)$ error. This implies that the iterates of the proposed method converge to the trajectory of the differential equation faster than those of Nesterov's method. The two toy numerical experiments illustrate such phenomenon for certain large step size. A matrix completion problem experiment is further included.

Review:  This paper refines the the truncation error analysis for discretizing the ODE to obtain accelerated optimization method.  The truncation results include higher order term. Built upon the analysis, the authors propose a new method which is claimed to be more stable for large step size and converges faster. Numerical evidence on matrix completion problem is provided.