The Physical Systems Behind Optimization Algorithms

2016·Arxiv

Abstract

Abstract

1 Introduction

Many machine learning problems can be cast into an optimization problem of the following form:

where X ⊆ Rd and f : X → R is a continuously differentiable function. For simplicity, we assume that f is convex or approximately convex (more on this later). Perhaps, the earliest algorithm for solving (1.1) is the vanilla gradient descent (VGD) algorithm, which dates back to Euler and Lagrange. VGD is simple, intuitive, and easy to implement in practice. For large-scale problems, it is usually more scalable than more sophisticated algorithms (e.g. Newton).

Existing state-of-the-art analysis of VGD shows it achieves O(1/k) convergence for general convex functions and linear convergence rate for strongly convex functions, where k is the number of iterations (Nesterov, 2013). Recently, a class of so-called Nesterov’s accelerated gradient (NAG) algorithms has gained popularity in statistical signal processing and machine learning communities. These algorithms combine the vanilla gradient descent algorithm with an additional momentum term at each iteration. Such a modification, though simple, has a profound impact: the NAG algorithms attain faster convergence than VGD. Specifically, NAG achieves O(1/k2) convergence for general convex functions, and linear convergence with a better constant term for strongly convex functions (Nesterov, 2013).

Another closely related class of algorithms is randomized coordinate gradient descent (RCGD) algorithms. These algorithms conduct a gradient descent-type step in each iteration, but only with respect to a single coordinate. RCGD has similar convergence rates to VGD, but has a smaller overall computational complexity, since its computational cost per iteration of RCGD is much

smaller than VGD (Nesterov, 2012; Lu and Xiao, 2015). More recently, Lin et al. (2014); Fercoq

and Richt´arik (2015) applied Nesterov’s acceleration to RCGD, and proposed accelerated randomized coordinate gradient (ARCG) algorithms. Accordingly, they established similar accelerated convergence rates for ARCG.

Another line of research focuses on relaxing the convexity and strong convexity conditions for alternative regularity conditions, including restricted secant inequality, error bound, Polyak-Łojasiewicz, and quadratic growth conditions. These conditions have been shown to hold for many optimization problems in machine learning, and faster convergence rates have been estab-

lished (e.g. Luo and Tseng (1993); Liu and Wright (2015); Necoara et al. (2015); Zhang and Yin

Although various theoretical results have been established, the algorithmic proof of convergence and regularity conditions in these analyses rely heavily on algebraic tricks that are sometimes arguably mysterious to be understood. To address this concern, differential equation approaches recently have attracted enormous interests on the analysis of optimization algorithms, because they provide a clear interpretation for the continuous approximation of the algorithmic systems (Su et al., 2014; Wibisono et al., 2016). Su et al. (2014) propose a framework for studying discrete algorithmic systems under the limit of infinitesimal time step. They show that Nesterov’s accelerated gradient (NAG) algorithm can be described by an ordinary differential equation (ODE) under the limit that time step tends to 0. Wibisono et al. (2016) study a more general family of ODE’s that essentially correspond to accelerated gradient algorithms. All these analyses, however, lack a link to a natural physical system behind the optimization algorithms. Therefore, they do not clearly explain why the momentum leads to acceleration. Meanwhile, these analyses only consider general convex conditions and gradient descent-type algorithms, and are NOT applicable to either the aforementioned relaxed conditions or coordinate-gradient-type algorithms (due to the randomized coordiante selection).

Our Contribution (I): We provide some new physics insights for the differential equation approaches. Particularly, we connect the differential equations to natural physical systems. Such a new connection allows us to establish a unified theory for understanding these optimization algorithms. Specifically, we consider the VGD, NAG, RCGD, and ARCG algorithms. All these algorithms are associated with damped oscillator systems with different particle mass and damping coefficients. For example, VGD corresponds to a massless particle system, NAG corresponds to a massive particle system. A damped oscillator system has a natural dissipation of its mechanical energy. The decay rate of the mechanical energy in the system essentially connects to the convergence rate of the algorithm. Our results match the convergence rates of all considered algorithms in existing literature. For a massless system, the convergence rate only depends on the gradient (force field) and smoothness of the function, whereas a massive particle system has an energy decay rate proportional to the ratio between the mass and damping coefficient. We further show that the optimal algorithm such as NAG correspond to an oscillator system near critical damping. Such a phenomenon is known in the physical literature that the critically damped system undergoes the fastest energy dissipation. Thus, this approach can potentially help us to design new optimization algorithms in a more intuitive way. As pointed out by the anonymous reviewers, although some of the intuitions provided in this paper are also presented in ?, we provide a more detailed analysis in this paper.

Our Contribution (II): We provide new analysis for more general optimization problems beyond general convexity and strong convexity, as well as more general algorithms. Specifically, we provide several concrete examples: (1) VGD achieves linear convergence under the Polyak-Łojasiewicz (PL) condition (possibly nonconvex), which matches the state-of-art result in Karimi et al. (2016); (2) NAG achieves accelerated linear convergence (with a better constant term) under both general convex and quadratic growth conditions, which matches the state-of-art result in Zhang (2016); (3) Coordinate-gradient-type algorithms share the same ODE approximation with gradient-type algorithms, and our analysis involves a more refined infinitesimal analysis; (4) Newton algorithm achieves linear convergence under the strongly convex and self-concordance conditions. See a summary in Table 1 as follows. Due to space limit, we present the extension to the nonsmooth composite optimization problem in Appendix.

Table 1: Our contribution compared with Su et al. (2014); Wibisono et al. (2016).

Recently, an independent paper of Wilson et al. (2016) studies a similar framework for analyzing first order optimization algorithms, and they focus on bridging the gap between discrete algorithmic analysis and continuous approximation. While we focus on understanding the physical systems behind the optimization. Both perspectives are essentially complementary to each other.

Before we proceed, we first introduce assumptions on the objective f .

Assumption 1.1 (L-smooth). There exits a constant L such that for any x, y ∈ Rd, we have ∥∇f (x)− ∇f (y)∥ ≤ L∥x − y∥.

Assumption 1.2 (µ-strongly convex). There exits a constant µ such that for any x, y ∈ Rd, we have f (x) ≥ f (y) + ⟨∇f (y),x − y⟩ + µ2∥x − y∥2.

Assumption 1.3 . (Lmax-coordinate-smooth) There exits a constant Lmax such that for any x, y ∈ Rd, we have |∇jf (x) − ∇jf (x\j,yj)| ≤ Lmax(xj − yj)2 for all j = 1,...,d.

The Lmax-coordinate-smooth condition has been shown to be satisfied by many machine learning problems such as Ridge Regression and Logistic Regression. For convenience, we define κ . Note that we also have Lmax ≤ L ≤ dLmax and κmax ≤ κ ≤ dκmax.

2 From Optimization Algorithms to ODE

We develop a unified representation for the continuous approximations of the aforementioned optimization algorithms. Our analysis is inspired by Su et al. (2014), where the NAG algorithm for general convex function is approximated by an ordinary differential equation under the limit of infinitesimal time step. We start with VGD and NAG, and later we will show that RCGD and ARCG can be approximated by the same ODE. For self-containedness, we present a brief review for popular optimization algorithms in Appendix A (VGD, NAG, RCGD, ARCG, and Newton).

2.1 A Unified Framework for Continuous Approximation Analysis

By considering an infinitesimal step size, we rewrite VGD and NAG in the following generic form:

convex. We then rewrite (2.1) as

When considering the continuous-time limit of the above equation, it is not immediately clear how the continuous-time is related to the step size k. We thus let h denote the time scaling factor

and study the possible choices of h later on. With this, we define a continuous time variable

where k is the iteration index, and X(t) from t is a trajectory characterizing the dynamics of the algorithm. Throughout the paper, we may omit (t) if it is clear from the context.

Note that our definition in (2.3) is very different from Su et al. (2014), where t is defined as t , i.e., fixing h . There are several advantages by using our new definition: (1) The new definition leads to a unified analysis for both VGD and NAG. Specifically, if we follow the same notion as Su et al. (2014), we need to redefine t for VGD, which is different from t NAG; (2) The new definition is more flexible, and leads to a unified analysis for both gradient-type (VGD and NAG) and coordinate-gradient-type algorithms (RCGD and ARCG), regardless of their different step sizes, e.g η = 1/L for VGD and NAG, and η for RCGD and ARCG; (3) The new definition is equivalent to Su et al. (2014) only when h . We will show later that, however, h ≍ √η is a natural requirement of a massive particle system rather than an artificial choice of h.

We then proceed to derive the differential equation for (2.2). By Taylor expansion

Taking the limit of h → 0, we rewrite (2.4) in a more convenient form,

Here (2.5) describes exactly a damped oscillator system in d dimensions with

Let us now consider how to choose h for different settings. The basic principle is that both m and c are finite under the limit h,η → 0. In other words, the physical system is valid. Taking VGD as an example, for which we have α under which, m → 0 and c → c0 for some constant c0. We call such a particle system massless. For NAG, it can also be verified that only h ) results in a valid physical system and it is massive (0 < m < ∞,0 ≤ c < ∞). Therefore, we provide a unified framework of choosing the correct time scaling factor h.

2.2 A Physical System: Damped Harmonic Oscillator

In classic mechanics, the harmonic oscillator is one of the first mechanic systems, which admit an exact solution. This system consists of a massive particle and restoring force. A typical example is a massive particle connecting to a massless spring.

The spring always tends to stay at the equilibrium position. When it is stretched or compressed, there will be a force acting on the object that stretches or compresses it. The force is

always pointing toward the equilibrium position. The energy stored in the spring is

where X denotes the displacement of the spring, and K is the Hooke’s constant of the spring. Here V (x) is called the potential energy in existing literature on physics.

Figure 1: An illustration of the

harmonic oscillators: A mas-

sive particle connects to a mass-

less spring. (Top) Undamped

harmonic oscillator; (Bottom)

Damped harmonic oscillator.

The harmonic oscillator is closely related to optimization algorithms. As we will show later, all our aforementioned optimization algorithms simply simulate a system, where a particle is falling inside a given potential. From a perspective of optimization, the equilibrium is essentially the minimizer of the quadratic potential function V (x. The desired property of the system is to stop the particle at the minimizer. However, a simple harmonic oscillator would not be sufficient and does not correspond to a convergent algorithm, since the system never stops: the particle at the equilibrium has the largest kinetic energy, and the inertia of the massive particle would drive it away from the equilibrium.

One natural way to stop the particle at the equilibrium is adding damping to the system, which dissipates the mechanic energy, just like the real-world mechanics. A simple damping is a force proportional to the negative velocity of the particle (e.g. submerge the system in some viscous

fluid) defined as

where c is the viscous damping coefficient. Suppose the potential energy of the system is f (x), then

the differential equation of the system is,

for under damped or nearly critical damped system (e.g. c2 ≲ 4mK).

For an over damped system (i.e. c2 > 4mK), the energy decay is

For extremely over damping cases, i.e., c2 ≫ 4mK, we have cm −�

not depend on the particle mass. The system exhibits a behavior as if the particle has no mass. In the language of optimization, the corresponding algorithm has linear convergence. Note that the convergence rate does only depend on the ratio c/m and does not depend on K when the system is under damped or critically damped. The fastest convergence rate is obtained, when the system is critically damped, c2 = 4mK.

2.3 Sufficient Conditions for Convergence

For notational simplicity, we assume that x∗ potential energy of the particle system is simply defined as V (t) := V (X(t)) := f (X(t)). If an algorithm converges to optimal, a sufficient condition is that the corresponding potential energy V decreases over time. The decreasing rate determines the convergence rate of the corresponding algorithm.

Theorem 2.1. Let γ(t) > 0 be a nondecreasing function of t and Γ(t) ≥ 0 be a nonnegative function.

Suppose that γ(t) and Γ(t) satisfy

Then the convergence rate of the algorithm is characterized by 1γ(t).

This further implies f (X) ≤ V (t) + Γ(t) ≤ γ(0+)(f (X(0+))+Γ(0+))γ(t) .

In words, γ(t)[V (t)+Γ(t)] serves as a Lyapunov function of system. We say that an algorithm is (1/γ)-convergent, if the potential energy decay rate is O(1/γ). For example, γ(tcorresponds to linear convergence, and γ = at corresponds to sublinear convergence, where a is a constant and independent of t. In the following section, we apply Theorem 2.1 to different problems by choosing different γ’s and Γ’s.

3 Convergence Rate in Continuous Time

We derive the convergence rates of different algorithms for different families of objective functions. Given our proposed framework, we only need to find γ and Γ to characterize the energy decay.

3.1 Convergence Analysis of VGD

We study the convergence of VGD for two classes of functions: (1) General convex function — Nesterov (2013) has shown that VGD achieves O(L/k) convergence for general convex functions; (2) A class of functions satisfying the Polyak-Łojasiewicz (PŁ) condition, which is defined as fol-

Assumption 3.1 . We say that f satisfies the µ-PŁ condition, if there exists a constant µ such that for any x ∈ Rd, we have 0 < f (x)∥∇f (x)∥2 ≤ 12µ.

Karimi et al. (2016) has shown that the PŁ condition is the weakest condition among the following conditions: strong convexity (SC), essential strong convexity (ESC), weak strong convexity (WSC), restricted secant inequality (RSI) and error bound (EB). Thus, the convergence analysis for the PŁ condition naturally extends to all the above conditions. Please refer to Karimi et al. (2016) for more detailed definitions and analyses as well as various examples satisfying such a condition in machine learning.

3.1.1 Sublinear Convergence for General Convex Function

By choosing Γ(t

where the last inequality follows from the convexity of f . Thus, Theorem 2.1 implies

Plugging t , we match the convergence rate in Nesterov

3.1.2 Linear Convergence Under the Polyak-Łojasiewicz Condition

By Theorem 2.1, for some constant C depending on x0, we obtain

which matches the behavior of an extremely over damped harmonic oscillator. Plugging t = kh and c , we match the convergence rate in Karimi et al. (2016):

for some constant C depending on x(0).

3.2 Convergence Analysis of NAG

We study the convergence of NAG for a class of convex functions satisfying the Polyak-Łojasiewicz (PŁ) condition. The convergence of NAG has been studied for general convex functions in Su et al. (2014), and therefore is omitted. Nesterov (2013) has shown that NAG achieves a linear convergence for strongly convex functions. Our analysis shows that the strong convexity can be relaxed as it does in VGD. However, in contrast to VGD, NAG requires f to be convex.

For a L-smooth convex function satisfying µ-PŁ condition, we have the particle mass and damping coefficient as m Karimi et al. (2016), under convex-

ity, PŁ is equivalent to quadratic growth (QG). Formally, we assume that f satisfies the following

condition.

Assumption 3.2 . We say that f satisfies the µ-QG condition, if there exists a constant µ such that for any x ∈ Rd, we have f (x) − f (x∗) ≥ µ2 ∥x − x∗∥2.

We then proceed with the proof for NAG. We first define two parameters, λ and σ. Let

Given properly chosen λ and σ, we show that the required condition in Theorem 2.1 is sat-isfied. Recall that our proposed physical system has kinetic energy m2 ∥ ˙X(t)∥2. In contrast to an un-damped system, NAG takes an effective velocity ˙X + σcX in the viscous fluid. By simple ma-

nipulation,

We then observe

Since c2 , we argue that if positive σ and λ satisfy

then we guarantee d(γ(t)(V (t)+Γ(t)))dt ≤ 0. Indeed, we obtain

By convexity of f , we have λc�1+ mσ2c2µ �f (X)−σc⟨X,∇f (X)⟩ ≤ σcf (X)−σc⟨X,∇f (X)⟩ ≤ 0. To make

(3.5) hold, it is sufficient to set σ . By Theorem 2.1, we obtain

for some constant C′′ depending on x(0). Plugging t

we have that

Comparing with VGD, NAG improves the constant term on the convergence rate for convex functions satisfying PŁ condition from L/µ to�L/µ. This matches with the algorithmic proof of Nes- terov (2013) for strongly convex functions, and Zhang (2016) for convex functions satisfying the QG condition.

3.3 Convergence Analysis of RCGD and ARCG

Our proposed framework also justifies the convergence analysis of the RCGD and ARCG algorithms. We will show that the trajectory of the RCGD algorithm converges weakly to the VGD algorithm, and thus our analysis for VGD directly applies. Conditioning on x(k), the updating

formula for RCGD is

where η is the step size and i is randomly selected from {1,2,...,d} with equal probabilities. Fixing

a coordinate i, we compute its expectation and variance as

We define the infinitesimal time scaling factor h ≤ η as it does in Section 2.1 and denote �Xh(t) := x(⌊t/h⌋). We prove that for each i ∈ [d], �Xhi (t) converges weakly to a deterministic function Xi(t) as

η → 0. Specifically, we rewrite (3.8) as,

Taking the limit of η → 0 at a fix time t, we have

Since ∥∇f (�Xh(t))∥2 is bounded at the time t, we have 1η Var��Xh(t + h) − �Xh(t)����Xh(t)�= O(h). Using an infinitesimal generator argument in Ethier and Kurtz (2009), we conclude that �Xh(t) converges to X(t) weakly as h → 0, where X(t) satisfies, ˙X(t) + 1d ∇f (X(t

by (3.4), we have

for some constant C1 depending on x(0). The analysis for general convex functions follows similarly. One can easily match the convergence rate as it does in (3.2), f (x(k)) ≤ c∥x0∥22kh

Repeating the above argument for ARCG, we obtain that the trajectory �Xh(t) converges weakly

to X(t), where X(t) satisfies

For general convex function, we have m . By the analysis of Su et al. (2014), we have f (xk) ≤ C2dk2 , for some constant C2 depending on x(0) and Lmax.

For convex functions satisfying µ-QG condition, m

f (xk) ≤ C3 exp�− 25d

3.4 Convergence Analysis for Newton

Newton’s algorithm is a second-order algorithm. Although it is different from both VGD and NAG, we can fit it into our proposed framework by choosing η and the gradient as L�∇2f (X)�−1 ∇f (X). We consider only the case f is µ-strongly convex, L-smooth and ν-self-concordant. By (2.5), if h/η

is not vanishing under the limit of h → 0, we achieve a similar equation,

where viscosity tensor of the system. In such a system, the function f not only determines the gradient field, but also determines a viscosity tensor field. The particle system is as if submerged in an anisotropic fluid that exhibits different viscosity along different directions. We release the particle at point x0 that is sufficiently close to the minimizer 0, i.e. ∥x0 − 0∥ ≤ ζ for some parameter ζ determined by ν, µ, and L. Now we consider the decay of the potential energy

V (X) := f (X). By Theorem 2.1 with γ(t

By simple calculus, we have ∇f (X. By the self-concordance condition, we have

where ∥v∥X 2. By integration and the convexity

of f , we have

Note that our proposed ODE framework only proves a local linear convergence for Newton method under the strongly convex, smooth and self concordant conditions. The convergence rate contains an absolute constant, which does not depend on µ and L. This partially justifies the superior local convergence performance of the Newton’s algorithm for ill-conditioned problems with very small µ and very large L. Existing literature, however, has proved the local quadratic convergence of the Newton’s algorithm, which is better than our ODE-type analysis. This is mainly because the discrete algorithmic analysis takes the advantage of “large” step sizes, but the ODE only characterizes “small” step sizes, and therefore fails to achieve quadratic convergence.

4 Numerical Illustration and Discussions

Due to the space limit, we present an numerical illustration in Figure 2. See more details on numerical results in Appendix C. We then give a more detailed interpretation of our proposed system from a perspective of physics:

Consequence of Particle Mass — As shown in Section 2, a massless particle system (mass m = 0) describes the simple gradient descent algorithm.

Figure 2: The algorithmic iter-

Damping and Convergence Rate — For a quadratic potential V (x, the system has a exponential energy decay, where the exponent factor depends on mass m, damping coefficient c, and the property of the function (e.g. PŁ-conefficient). As discussed in Section 2, the decay rate is the fastest when the system is critically damped, i.e, c2 . For either under or over damped system, the decay rate is slower. For a potential function f satisfying convexity and µ-PŁ condition, NAG corresponds to a nearly critically damped system, whereas VGD corresponds to an extremely over damped system, i.e., c2 ≫ 4mµ. Moreover, we can achieve different acceleration rate by choosing different m/c ratio for NAG, i.e., α for some absolute constant s > 0. However s = 1/2 achieves the largest convergence rate since it is exactly the critical damping: c2 = 4mµ.

Connecting PŁ Condition to Hooke’s law — The µ-PŁ and convex conditions together naturally mimic the property of a quadratic potential V , i.e., a damped harmonic oscillator. Specifically, the

µ-PŁ condition

guarantees that the force field is strong enough, since the left hand side of the above equation is exactly the potential energy of a spring based on Hooke’s law. Moreover, the convexity condition V (x) ≤ ⟨∇V (x),X⟩ guarantees that the force field has a large component pointing at the equilibrium point (acting as a restoration force). As indicated in Karimi et al. (2016), PŁ is a much weaker condition than the strong convexity. Some functions that satisfy local PŁ condition do not even satisfy convexity, e.g., matrix factorization. The connection between the PŁ condition and the Hookes law indicates that strong convexity is not the fundamental characterization of linear convergence. If there is another condition that employs a form of the Hookes law, it should employ linear convergence as well.

References

Ethier, S. N. and Kurtz, T. G. (2009). Markov processes: characterization and convergence, vol. 282. John Wiley & Sons.

Fercoq, O. and Richt´arik, P. (2015). Accelerated, parallel, and proximal coordinate descent. SIAM Journal on Optimization 25 1997–2023.

Gong, P. and Ye, J. (2014). Linear convergence of variance-reduced stochastic gradient without strong convexity. arXiv preprint arXiv:1406.1102 .

Karimi, H., Nutini, J. and Schmidt, M. (2016). Linear convergence of gradient and proximalgradient methods under the polyak-łojasiewicz condition. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases. Springer.

Lin, Q., Lu, Z. and Xiao, L. (2014). An accelerated proximal coordinate gradient method. In Advances in Neural Information Processing Systems.

Liu, J. and Wright, S. J. (2015). Asynchronous stochastic coordinate descent: Parallelism and convergence properties. SIAM Journal on Optimization 25 351–376.

Lu, Z. and Xiao, L. (2015). On the complexity analysis of randomized block-coordinate descent methods. Mathematical Programming 152 615–642.

Luo, Z.-Q. and Tseng, P. (1993). Error bounds and convergence analysis of feasible descent methods: a general approach. Annals of Operations Research 46 157–178.

Necoara, I., Nesterov, Y. and Glineur, F. (2015). Linear convergence of first order methods for non-strongly convex optimization. arXiv preprint arXiv:1504.06298 .

Nesterov, Y. (2012). Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization 22 341–362.

Nesterov, Y. (2013). Introductory lectures on convex optimization: A basic course, vol. 87. Springer Science & Business Media.

Nocedal, J. and Wright, S. (2006). Numerical optimization. Springer Science & Business Media.

Polyak, B. T. (1963). Gradient methods for minimizing functionals. Zhurnal Vychislitel’noi Matematiki i Matematicheskoi Fiziki 3 643–653.

Rockafellar, R. T. (2015). Convex analysis. Princeton university press.

Su, W., Boyd, S. and Candes, E. (2014). A differential equation for modeling nesterov’s accelerated gradient method: theory and insights. In Advances in Neural Information Processing Systems.

Wibisono, A., Wilson, A. C. and Jordan, M. I. (2016). A variational perspective on accelerated methods in optimization. arXiv preprint arXiv:1603.04245 .

Wilson, A. C., Recht, B. and Jordan, M. I. (2016). A lyapunov analysis of momentum methods in optimization. arXiv preprint arXiv:1611.02635 .

Zhang, H. (2016). New analysis of linear convergence of gradient-type methods via unifying error bound conditions. arXiv preprint arXiv:1606.00269 .

Zhang, H. and Yin, W. (2013). Gradient methods for convex minimization: better rates under weaker conditions. arXiv preprint arXiv:1303.4645 .

A A Brief Review of Popular Optimization Algorithms

A.1 Vanilla Gradient Descent Algorithm

A vanilla gradient descent (VGD) algorithm starts from an arbitrary initial solution x(0). At the

k-th iteration (k > 0), VGD takes

where η is a properly chosen step size. Since VGD only needs to calculate a gradient of f in each iteration, the computational cost per iteration is usually linearly dependent on d. For a L-smooth f , we can choose a constant step size such that η ≤ 1L to guarantee convergence.

VGD has been extensively studied in existing literature. Nesterov (2013) show that:

A.2 Nesterov’s Accelerated Gradient Algorithms

The Nesterov’s accelerated gradient (NAG) algorithms combines the vanilla gradient descent algorithm with an additional momentum at each iteration. Such a modification, though simple, enables NAG to attain better convergence rate than VGD. Specifically, NAG starts from an arbitrary initial solution x(0) along with an auxiliary solution y(0) -th iteration, NAG

takes

where α for general convex f and α for strongly convex f . Intuitively speaking,

NAG takes an affine combination of the current and previous solutions to compute the update

for the two subsequent iterations. This can be viewed as the momentum of a particle during its movement. Similar to VGD, NAG only needs to calculate a gradient of f in each iteration. Similar to VGD, we can choose η ≤ 1L for a L-smooth f to guarantee convergence.

NAG has also been extensively studied in existing literature. Nesterov (2013) show that:

A.3 Randomized Coordinate Gradient Descent Algorithm

A randomized coordinate gradient descent (RCGD) algorithm is closely related to VGD. RCGD starts from an arbitrary initial solution x(0). Different from VGD, RCGD takes a gradient descent step only over a coordinate. Specifically, at the k-th iteration (k > 0), RCGD randomly selects a

coordinate j from 1,...,d, and takes

where η is a properly chosen step size. Since RCGD only needs to calculate a coordinate gradient of f in each iteration, the computational cost per iteration usually does not scale with d. For a Lmax-coordinate-smooth f , we can choose a constant step size such that η ≤ 1Lmax to guarantee convergence.

RCGD has been extensively studied in existing literature. Nesterov (2012); Lu and Xiao (2015) show that:

A.4 Accelerated Randomized Coordinate Gradient Algorithms

Similar to NAG, the accelerated randomized coordinate gradient (ARCG) algorithms combine the randomized coordinate gradient descent algorithm with an additional momentum at each iteration. Such a modification also enables ARCG to attain better convergence rate than RCGD. Specifically, ARCG starts from an arbitrary initial solution x(0) along with an auxiliary solution y(0) -th iteration (k > 0), ARCG randomly selects a coordinate j from 1,...,d, and

takes

Here α is strongly convex, and α is general convex. Similar to RCGD, we can choose η ≤ 1Lmax for a Lmax-coordinate-smooth f to guarantee convergence.

ARCG has been studied in existing literature. Lin et al. (2014); Fercoq and Richt´arik (2015) show that:

A.5 Newton’s Algorithm

The Newton’s (Newton) algorithm requires f to be twice differentiable. It starts with an arbitrary

initial x(0). At the k-th iteration (k > 0), Newton takes

The inverse of the Hessian matrix adjusts the descent direction by the landscape at x(k−1). Therefore, Newton often leads to a steeper descent than VGD and NAG in each iteration, espcially for highly ill-conditioned problems.

Newton has been extensively studied in existing literature with an additional self-concordant assumption as follows:

Assumption A.1 . Suppose that f is smooth and convex. We define g(t) = f (x + tv). We say that f is self-concordant, if for any x ∈ Rd, v ∈ Rd, and t ∈ R, there exists a constant ν, which is

independent on f such that we have

Nocedal and Wright (2006) show that for a L-smooth, µ-strongly convex and ν-self-concordant, f , Newton attains a local quadratic convergence in conjunction. Specifically, given a suitable initial solution x(0) satisfying ∥x(0) − x∗∥2 ≤ ζ, where ζ < 1 is a constant depending on on L, µ, and

ν, there exists a constant ξ depending only on ν such that we have

Note that (A.9) is also referred as an iteration complexity of �O(loglog(1/ϵ)), where �O hides the constant term depending on L, µ, and ν. Since Newton needs to calculate the inverse of the Hessian matrix, its per iteration computation cost is at least O(d3). Thus, it outperforms VGD and NAG when we need a highly accurate solution, i.e., ϵ is very small.

B Extension to Nonsmooth Composite Optimization

Our framework can also be extended to nonsmooth composite optimization in a similar manner to Su et al. (2014). Let g be an L-smooth function, and h be a general convex function (not necessarily

smooth). For x ∈ Rd, the composite optimization problem solves

Analogously to Su et al. (2014), we define the force field as the directional subgradient G(x,p) of function f , where G : Rd × Rd → Rd is defined as G(x,p) ∈ ∂f (x) and ⟨G(x,p),p⟩ where ∂f (x) denotes the sub-differential of f . The existence of G(x,p) is guaranteed by Rockafellar

(2015). Accordingly, a new ODE describing the dynamics of the system is

Under the assumption that the solution to the ODE exists and is unique, we illustrate the analysis by VGD (the mass m

are straightforward. Specifically, a convex function f satisfies µ-proximal-PŁ if

where x∗ = 0 is the global minimum point of f . Slightly different from the definition of the proximal-PŁ condition in Karimi et al. (2016) involving a step size parameter, (B.1) does not involve any additional parameter. This is actually a more intuitive definition by choosing an appro-

priate subgradient. Let γ(t0, we study

By Taylor expansions and the local Lipschitz property of convex function f , we have

Combining the above two expansions, we obtain

This further implies

analysis follows exactly the same as it does in Section 3.1.2.

C Numerical Illustration

We present an illustration of our theoretical analysis in Figure 2. We consider a strongly convex

quadratic program

Obviously, f (x) is strongly convex and x∗ is the minimizer. We choose η and NAG, and η for RCGD and ARCG. The trajectories of VGD and NAG are obtained by the default method for solving ODE in MATLAB.

designed for accessibility and to further open science