One of the main recent trends within machine learning and data analytics using massive data sets is to leverage the inferential strength of the vast amounts of data by using relatively simple, but fast, optimization methods as algorithmic primitives [1]. Many of these optimization methods are modifications of the basic gradient descent (GD) method. Indeed, computationally more heavy approaches, such as interior point methods, are often infeasible for a given limited computational budget [2].

Moreover, the rise of deep learning has brought a significant boost for the interest in gradient methods. Indeed, a major insight within the theory of deep learning is that for typical high-dimensional models, e.g., those represented by deep neural networks, most of the local minima of the cost function (e.g., the empirical loss or training error) are reasonably close (in terms of objective value) to the global optimum [3]. These local minima can be found efficiently by gradient methods such as stochastic gradient descent (SGD), which is considered the de-facto standard algorithmic primitive for training deep neural networks [3].

This paper elaborates on the interpretation of some basic gradient methods such as GD and its variants as fixed-point iterations. These fixed-point iterations are obtained for operators associated with the convex objective function. Emphasizing the connection to fixed-point theory unleashes some powerful tools, e.g., on the acceleration of fixed-point iterations [4] or inexact fixed-point iterations [5, 6], for the analysis and construction of convex optimization methods.

In particular, we detail how the convergence of the basic GD iterations can be understood from the contraction properties of a specific operator which is associated naturally with a differentiable objective function. Moreover, we work out in some detail how the basic GD method can be accelerated by modifying the operator underlying GD in a way that preserves its fixed-points but decreases the contraction factor which implies faster convergence by the contraction mapping theorem.

A function f(·):ℝⁿ → ℝ is convex if

holds for any x, y ∈ ℝⁿ and α ∈ [0, 1] [2]. For a differentiable function f(·) with gradient ∇f(x), a necessary and sufficient condition for convexity is [7, p. 70]

which has to hold for any x, y ∈ ℝⁿ.

Our main object of interest in this paper is the optimization problem

Given a convex function f(x), we aim at finding a point x₀ with lowest function value f(x₀), i.e., f(x0)=minxf(x).

In order to motivate our interest in optimization problems like Equation (1), consider a machine learning problem based on training data X :={z(i)}i=1N consisting of N data points z⁽ⁱ⁾ = (d⁽ⁱ⁾, y⁽ⁱ⁾) with feature vector d⁽ⁱ⁾ ∈ ℝⁿ (which might represent the RGB pixel values of a webcam snapshot) and output or label y⁽ⁱ⁾ ∈ ℝ (which might represent the local temperature during the snapshot). We wish to predict the label y⁽ⁱ⁾ by a linear combination of the features, i.e.,

The choice for the weight vector x ∈ ℝⁿ is typically based on balancing the empirical risk incurred by the predictor (2), i.e.,

with some regularization term, e.g., measured by the squared norm ‖x‖². Thus, the learning problem amounts to solving the optimization problem

The learning problem (3) is precisely of the form (1) with the convex objective function

By choosing a large value for the regularization parameter λ, we de-emphasize the relevance of the training error and thus avoid overfitting. However, choosing λ too large induces a bias if the true underlying weight vector has a large norm [3, 8]. A principled approach to find a suitable value of λ is cross validation [3, 8].

Let us now show how one of the most basic methods for solving the problem (1), i.e., the GD method, can be obtained naturally as fixed-point iterations involving the gradient operator ∇f (cf. 5).

Our point of departure is the necessary and sufficient condition [7]

for a vector x0∈ℝn to be optimal for the problem (1) with a convex differentiable objective function f(·)∈SnL,U.

Lemma 1. We have ∇f(x) = 0 if and only if the vector x ∈ ℝⁿ is a fixed point of the operator

for an arbitrary but fixed non-zero α ∈ ℝ \ {0}. Thus,

Proof. Consider a vector x such that ∇f(x) = 0. Then,

Conversely, let x be a fixed point of T(α), i.e.,

Then,

     □

According to Lemma (1), the solution x₀ of the optimization problem (1) is obtained as the fixed point of the operator T(α) (cf. 14) with some non-zero α. As we will see shortly, the freedom in choosing different values for α can be exploited in order to compute the fixed points of T(α) more efficiently.

A straightforward approach to finding the fixed-points of an operator T(α) is via the fixed-point iteration

By tailoring a fundamental result of analysis (cf. [10], Theorem 9.23], we can characterize the convergence of the sequence x^(k) obtained from (16).

Lemma 2. Assume that for some q∈[0, 1), we have

for any x, y ∈ ℝⁿ. Then, the operator T(α) has a unique fixed point x₀ and the iterates x^(k) (cf. 16) satisfy

Proof. Let us first verify that the operator T(α) cannot have two different fixed points. Indeed, assume there would be two different fixed points x, y such that

This would imply, in turn,

However, since q < 1, this inequality can only be satisfied if ‖x − y‖ = 0, i.e., we must have x = y. Thus, we have shown that no two different fixed points can exist. The existence of one unique fixed point x₀ follows from Rudin [10, Theorem 9.23].

The estimate (18) can be obtained by induction and noting

Here, step (a) is valid since x₀ is a fixed point of T(α), i.e., x0=T(α)x0.     □

In order to apply Lemma (2) to (16), we have to ensure that the operator T(α) is a contraction, i.e., it satisfies (17) with some contraction coefficient q ∈ [0, 1). For the operator T(α) (cf. 14) associated with the function f(·)∈SnL,U this can be verified by standard results from vector analysis.

Lemma 3. Consider the operator T(α):x↦x-α∇f(x) with some convex function f(·)∈SnL,U. Then,

with contraction factor

Proof. First,

using z = ηx + (1 − η)y with some η ∈ [0, 1]. Here, we used in step (a) the mean value theorem of vector calculus [10, Theorem 5.10].

Combining (21) with the submultiplicativity of Euclidean and spectral norm [11, p. 55] yields

The matrix M^(α): = I−α∇²f(z) is symmetric (M^(α) = (M^(α))^T) with real-valued eigenvalues [11]

Since also

we obtain from (22)

     □

It will be handy to write out the straightforward combination of Lemma (2) and Lemma (3).

Lemma 4. Consider a convex function f(·)∈SnL,U with the unique minimizer x₀, i.e., f(x0)=minxf(x). We then construct the operator T(α):x↦x-α∇f(x) with a step size α such that

Then, starting from an arbitrary initial guess x⁽⁰⁾, the iterates x^(k) (cf. 16) satisfy

According to Lemma (4), and also illustrated in Figure 1, starting from an arbitrary initial guess x⁽⁰⁾, the sequence x^(k) generated by the fixed-point iteration (16) is guaranteed to converge to the unique solution x₀ of (1), i.e., limk→∞x(k)=x0. What is more, this convergence is quite fast, since the error ||x(k)-x0|| decays at least exponentially according to (25). Loosely speaking, this exponential decrease implies that the number of additional iterations required to have on more correct digit in x^(k) is constant.

Let us now work out the iterations (16) more explicitly by inserting the expression (14) for the operator T(α). We then obtain the following equivalent representation of (16):

This iteration is nothing but plain vanilla GD using a fixed step size α [3].

Since the GD iteration (26) is precisely the fixed-point iteration (16), we can use Lemma (4) to characterize the convergence (rate) of GD. In particular, convergence of GD is ensured by choosing the step size of GD (26) such that q(α) = max{|1−Uα|, |1−Lα|} < 1. Moreover, in order to make the convergence as fast as possible we need to chose the step size α = α^* which makes the contraction factor q(α) (cf. 20) as small as possible.

In Figure 2, we illustrate how the quantifies |1 − αL| and |1 − αU| evolve as the step size α (cf. 26) is varied. From Figure 2 we can easily read off the optimal choice

yielding the smallest possible contraction factor

We have arrived at the following characterization of GD for minimizing convex functions f(·)∈SnL,U.

Theorem 5. Consider the optimization problem (1) with objective function f(·)∈SnL,U, where the parameters L and U are fixed and known. Starting from an arbitrarily chosen initial guess x⁽⁰⁾, we construct a sequence by GD (26) using the optimal step size (27). Then,

In what follows, we will use the shorthand T :=T(α*) for the gradient operator T(α) (cf. 14) obtained for the optimal step size α = α^* (cf. 27).

Without a computational model taking into account a finite amount of resources, the study of the computational complexity inherent to (1) becomes meaningless. Consider having unlimited computational resources at our disposal. Then, we could build an “optimization device” which maps each function f(·)∈SnL,U to its unique minimum x₀. Obviously, this approach is infeasible since we cannot perfectly represent such a mapping, let alone its domain SnL,U, using a physical hardware which allows us only to handle finite sets instead of continuous spaces like SnL,U.

Let us further illustrate the usefulness of using a computational model in the context of machine learning from massive data sets (big data). In particular, as we have seen in the previous section, the regularized linear regression model (3) amounts to minimizing a convex quadratic function (8) with the particular choices (12). Even for this most simple machine learning model, it is typically infeasible to have access to a complete description of the objective function (8).

Indeed, in order to fully specify the quadratic function in (8), we need to fully specify the matrix Q ∈ ℝ^n×n and the vector q ∈ ℝⁿ. For the (regularized) linear regression model (3) this would require to compute Q_LR (cf. 12) from the training data X={z(i)}i=1N. Computing the matrix Q_LR in a naive way, i.e., without exploiting any additional structure, amounts to a number of arithmetic operations on the order of N · n². This might be prohibitive in a typical big data application with N and n being on the order of billions and using distributed storage of the training data X [12].

There has emerged a widely accepted computational model for convex optimization which abstracts away the details of the computational (hard- and software) infrastructure. Within this computational model, an optimization method for solving (1) is not provided with a complete description of the objective function, but rather it can access the objective function only via an “oracle” [2, 13].

We might think of an oracle model as an application programming interface (API), which specifies the format of queries which can be issued by a convex optimization method executed on an application layer (cf. Figure 3). There are different types of oracle models but one of the most popular type (in particular for big data applications) is a first order oracle [13]. Given a query point x ∈ ℝⁿ, a first order oracle returns the gradient ∇f(x) of the objective function at this particular point.

A first order method (FOM) aims at solving (1) by sequentially querying a first order oracle, at the current iterate x^(k), to obtain the gradient ∇f(x^(k)) (cf. Figure 3). Using the current and past information obtained from the oracle, a FOM then constructs the new iterate x^{(k + 1)} such that eventually limk→∞x(k)=x0. For the sake of simplicity and without essential loss in generality, we will only consider FOMs whose iterates x^(k) satisfy [13]

According to Section (3), solving (1) can be accomplished by the simple GD iterations (26). The particular choice α^* (27) for the step size α in (26) ensures the convergence rate (κ−1κ+1)k with the condition number κ = U/L of the function class SnL,U. While this convergence is quite fast, i.e., the error decays exponentially with iteration number k, we would, of course, like to know how efficient this method is in general.

As detailed in Section (4), in order to study the computational complexity and efficiency of convex optimization methods, we have to define a computational model such as those underlying FOMs (cf. Figure 3). The next result provides a fundamental lower bound on the convergence rate of any FOM (cf. 29) for solving (1).

Theorem 6. Consider a particular FOM, which for a given convex function f(·)∈SnL,U generates iterates x^(k) satisfying (29). For fixed L, U there is a sequence of functions fn(·)∈SnL,U (indexed by dimension n) such that

with a sequence δ(n) such that limn→∞|δ(n)|=0.

Proof. see Section (8.1).     □

There is a considerable gap between the upper bound (28) on the error achieved by GD after k iterations and the lower bound (30) which applies to any FOM which is run for the same number iterations. In order to illustrate this gap, we have plotted in Figure 4 the upper and lower bound for the (quite moderate) condition number κ = 100.

Thus, there might exist a FOM which converges faster than the GD method (28) and comes more close to the lower bound (30). Indeed, in the next section, we will detail how to obtain an accelerated FOM by applying a fixed point preserving transformation to the operator T (cf. 16), which is underlying the GD method (26). This accelerated gradient method is known as the heavy balls (HB) method [14] and effectively achieves the lower bound (30), i.e., the HB method is already optimal among all FOM's for solving (1) with an objective function f(·)∈SnL,U.

Let us now show how to modify the basic GD method (26) in order to obtain an accelerated FOM, whose convergence rate essentially matches the lower bound (30) for the function class SnL,U with condition number κ = U/L>1 (cf. 7) and is therefore optimal among all FOMs.

Our derivation of this accelerated gradient method, which is inspired by the techniques used in Ghadimi et al. [15], starts from an equivalent formulation of GD as the fixed-point iteration

with the operator

As can be verified easily, the fixed-point iteration (31) starting from an arbitrary initial guess x¯(0)=(z(0)y(0)) is related to the GD iterate x^(k) (cf. 26), using initial guess z⁽⁰⁾, as

for all iterations k ≥ 1.

By the equivalence (33), Theorem (5) implies that for any initial guess x̄(0) the iterations (31) converge to the fixed point

with x₀ being the unique minimizer of f(·)∈SnL,U. Moreover, the convergence rate of the fixed-point iterations (31) is precisely the same as those of the GD method, i.e., governed by the decay of (κ−1κ+1)k, which is obtained for the optimal step size α = α^* (cf. 27).

We will now modify the operator T¯ in (32) to obtain a new operator M:ℝ2n→ℝ2n which has the same fixed points (34) but improved contraction behavior, i.e., the fixed point iteration

will converge faster than those obtained from T¯ in (31) (cf. Figure 5). In particular, this improved operator M is defined as

with

As can be verified easily, the fixed point (x0T,x0T)T of T¯ is also a fixed point of M.

Before we analyze the convergence rate of the fixed-point iteration (35), let us work out explicitly the FOM which is represented by the fixed-point iteration (35). To this end, we partition the kth iterate, for k ≥ 1, as

Inserting (38) into (35), we have for k ≥ 1

with the convention xHB(-1) :=0. The iteration (39) defines the HB method [14] for solving the optimization problem (1). As can be verified easily, like the GD method, the HB method is a FOM. However, contrary to the GD iteration (26), the HB iteration (39) also involves the penultimate iterate xHB(k-2) for determining the new iterate xHB(k).

We will now characterize the converge rate of the HB method (39) via its fixed-point equivalent (35). To this end, we restrict ourselves to the subclass of SnL,U given by quadratic functions of the form (8).

Theorem 7. Consider the optimization problem (1) with objective function f(·)∈SnL,U which is a quadratic (26). Starting from an arbitrarily chosen initial guess xHB(-1) and xHB(0), we construct a sequence xHB(k) via iterating (27). Then,

with

Proof. see Section (8.2).     □

The upper bound (40) differs from the lower bound (30) by the factor k. However, the discrepancy is rather decent as this linear factor in (40) grows much slower than the exponential (κ−1κ+1)k in (40) decays. In Figure 6, we depict the upper bound (40) on the error of the HB iterations (39) along with the upper bound (28) on the error of the GD iterations (26) and the lower bound (30) on the error of any FOM after k iterations.

We highlight that, strictly speaking, the bound (40) only applies to a subclass of smooth strongly convex functions f(·)∈SnL,U, i.e., it applies only to quadratic functions of the form (8). However, as discussed in Section (2), given a particular point x, we can approximate an arbitrary function f(·)∈SnL,U with a quadratic function f~(x) of the form (8). The approximation error ε(x) (cf. (11)) will be small for all points x sufficiently close to x₀. Making this reasoning more precise and using well-known results on fixed-point iterations with inexact updates [6], one can verify that the bound (40) essentially applies to any function f(·)∈SnL,U.

We have presented a fixed-point theory of some basic gradient methods for minimizing convex functions. The approach via fixed-point theory allows for a rather elegant analysis of the convergence properties of these gradient methods. In particular, their convergence rate is obtained as the contraction factor for an operator associated with the objective function.

The fixed-point approach is also appealing since it leads rather naturally to the acceleration of gradient methods via fixed-point preserving transformations of the underlying operator. We plan to further develop the fixed-point theory of gradient methods in order to accommodate stochastic variants of GD such as SGD. Furthermore, we can bring the popular class of proximal methods into the picture by replacing the gradient operator underlying GD with the proximal operator.

However, by contrast to FOMs (such as the GD method), proximal methods use a different oracle model (cf. Figure 3). In particular, proximal methods require an oracle which can evaluate the proximal mapping efficiently which is typically more expensive than gradient evaluations. Nonetheless, the popularity of proximal methods is due to the fact that for objective functions arising in many important machine learning applications, the proximal mapping can be evaluated efficiently.

In this section we present the (somewhat lengthy) proofs for the main results stated in Sections (5) and (6).

We collect some elementary results from linear algebra and analysis, which are required to prove our main results.

Lemma 8. Consider a matrix M = (ab10) ∈ ℝ^2×2 with spectral radius ρ(M). Then, there is an orthonormal matrix U ∈ ℂ^2×2 such that

where |λ₁|, |λ₂| ≤ ρ(M) and |d|(|a| + |b| + 1)ρ^k−1(M).

Proof. Consider an eigenvalue λ₁ of the matrix M with normalized eigenvector u=(u1,u2)H, i.e., Mu = λ₁u with ‖u‖ = 1. According to Golub and Van Loan [11, Lemma 7.1.2], we can find a normalized vector v=(v1,v2)H, orthogonal to u, such that

or equivalently

with some eigenvalue λ₂ of M. As can be read off (67), d = u₁(u₂a + v₂b) + v₁u₂ which implies (65) since |u₁|, |u₂|, |v₁|, |v₂| ≤ 1. Based on (66), we can verify (65) by induction.     □

Lemma 9. For any κ > 1 and k ∈ ℕ,

Proof. Let us introduce the shorthand z := exp(j2πθ) and further develop the LHS of (68) as

The denominator of the integrand in (69) can be factored as