Title: STOCHASTIC GRADIENT/MIRROR DESCENT: MINI-MAX OPTIMALITY AND IMPLICIT REGULARIZATION

Abstract: Stochastic descent methods (of the gradient and mirror varieties) have become increasingly popular in optimization. In fact, it is now widely recognized that the success of deep learning is not only due to the special deep architecture of the models, but also due to the behavior of the stochastic descent methods used, which play a key role in reaching "good" solutions that generalize well to unseen data. In an attempt to shed some light on why this is the case, we revisit some minimax properties of stochastic gradient descent (SGD) for the square loss of linear models-originally developed in the 1990's-and extend them to general stochastic mirror descent (SMD) algorithms for general loss functions and nonlinear models. In particular, we show that there is a fundamental identity which holds for SMD (and SGD) under very general conditions, and which implies the minimax optimality of SMD (and SGD) for sufficiently small step size, and for a general class of loss functions and general nonlinear models. We further show that this identity can be used to naturally establish other properties of SMD (and SGD), namely convergence and implicit regularization for over-parameterized linear models (in what is now being called the "interpolating regime"), some of which have been shown in certain cases in prior literature. We also argue how this identity can be used in the so-called "highly over-parameterized" nonlinear setting (where the number of parameters far exceeds the number of data points) to provide insights into why SMD (and SGD) may have similar convergence and implicit regularization properties for deep learning.

1. INTRODUCTION:

Deep learning has proven to be extremely successful in a wide variety of tasks (Krizhevsky et al., 2012; LeCun et al., 2015; Mnih et al., 2015; Silver et al., 2016; Wu et al., 2016) . Despite its tremendous success, the reasons behind the good generalization properties of these methods to unseen data is not fully understood (and, arguably, remains somewhat of a mystery to this day). Initially, this success was mostly attributed to the special deep architecture of these models. However, in the past few years, it has been widely noted that the architecture is only part of the story, and, in fact, the optimization algorithms used to train these models, typically stochastic gradient descent (SGD) and its variants, play a key role in learning parameters that generalize well.

In particular, it has been observed that since these deep models are highly over-parameterized, they have a lot of capacity, and can fit to virtually any (even random) set of data points (Zhang et al., 2016) . In other words, highly over-parameterized models can "interpolate" the data, so much so that this regime has been called the "interpolating regime" (Ma et al., 2018) . In fact, on a given dataset, the loss function often has (uncountably infinitely) many global minima, which can have drastically different generalization properties, and it is not hard to construct "trivial" global minima that do not generalize. Which minimum among all the possible minima we pick in practice is determined by the optimization algorithm that we use for training the model. Even though it may seem at first that, because of the non-convexity of the loss function, the stochastic descent algorithms may get stuck in local minima or saddle points, in practice they almost always achieve a global minimum (Kawaguchi, 2016; Zhang et al., 2016; Lee et al., 2016) , which perhaps can also be justified by the fact that these models are highly over-parameterized. What is even more interesting is that not only do these stochastic descent algorithms converge to global minima, but they converge to "special" ones that generalize well, even in the absence of any explicit regularization or early stopping (Zhang et al., 2016) . Furthermore, it has been observed that even among the common optimization algorithms, namely SGD or its variants (AdaGrad (Duchi et al., 2011) , RMSProp (Tieleman & Hinton, 2012) , Adam (Kingma & Ba, 2014) , etc.), there is a discrepancy in the solutions achieved by different algorithms and their generalization capabilities (Wilson et al., 2017) , which again highlights the important role of the optimization algorithm in generalization.

There have been many attempts in recent years to explain the behavior and properties of these stochastic optimization algorithms, and many interesting insights have been obtained (Achille & Soatto, 2017; Chaudhari & Soatto, 2018; Shwartz-Ziv & Tishby, 2017; Soltanolkotabi et al., 2017) . In particular, it has been argued that the optimization algorithms perform an implicit regularization (Neyshabur et al., 2017; Ma et al., 2017; Gunasekar et al., 2017; 2018a; Soudry et al., 2017; Gunasekar et al., 2018b) while optimizing the loss function, which is perhaps why the solution generalizes well. Despite this recent progress, most results explaining the behavior of the optimization algorithm, even for SGD, are limited to linear or very simplistic models. Therefore, a general characterization of the behavior of stochastic descent algorithms for more general models would be of great interest.

1.1. OUR CONTRIBUTION:

In this paper, we present an alternative explanation of the behavior of SGD, and more generally, the stochastic mirror descent (SMD) family of algorithms, which includes SGD as a special case. We do so by obtaining a fundamental identity for such algorithms (see Lemmas 2 and 5). Using these identities, we show that for general nonlinear models and general loss functions, when the step size is sufficiently small, SMD (and therefore also SGD) is the optimal solution of a certain minimax filtering (or online learning) problem. The minimax formulation is inspired by, and rooted, in H ∞ filtering theory, which was originally developed in the 1990's in the context of robust control theory (Hassibi et al., 1999; Simon, 2006; Hassibi et al., 1996) , and we generalize several results from this literature, e.g., (Hassibi et al., 1994; Kivinen et al., 2006) . Furthermore, we show that many properties recently proven in the learning/optimization literature, such as the implicit regularization of SMD in the over-parameterized linear case-when convergence happens- (Gunasekar et al., 2018a) , naturally follow from this theory. The theory also allows us to establish new results, such as the convergence (in a deterministic sense) of SMD in the over-parameterized linear case. We also use the theory developed in this paper to provide some speculative arguments into why SMD (and SGD) may have similar convergence and implicit regularization properties in the so-called "highly over-parameterized" nonlinear setting (where the number of parameters far exceeds the number of data points) common to deep learning.

In an attempt to make the paper easier to follow, we first describe the main ideas and results in a simpler setting, namely, SGD on the square loss of linear models, in Section 3, and mention the connections to H ∞ theory. The full results, for SMD on a general class of loss functions and for general nonlinear models, are presented in Section 4. We demonstrate some implications of this theory, such as deterministic convergence and implicit regularization, in Section 5, and we finally conclude with some remarks in Section 6. Most of the formal proofs are relegated to the appendix.

2. PRELIMINARIES:

Denote the training dataset by {(x i , y i ) : i = 1, . . . , n}, where x i ∈ R d are the inputs, and y i ∈ R are the labels. We assume that the data is generated through a (possibly nonlinear) model f i (w) = f (x i , w) with some parameter vector w ∈ R m , plus some noise v i , i.e., y i = f (x i , w) + v i for i = 1, . . . , n. The noise can be due to actual measurement error, or it can be due to modeling error (if the model f (x i , •) is not rich enough to fully represent the data), or it can be a combination of both. As a result, we do not make any assumptions on the noise (such as stationarity, whiteness, Gaussianity, etc.).

Since typical deep models have a lot of capacity and are highly over-parameterized, we are particularly interested in the over-parameterized (so-caled interpolating) regime, i.e., when m > n. In this case, there are many parameter vectors w (in fact, uncountably infinitely many) that are consistent with the observations. We denote the set of these parameter vectors by

W = {w ∈ R m | y i = f (x i , w), i = 1, . . . , n} .

(1) (Note the absence of the noise term, since in this regime we can fully interpolate the data.) The set W is typically an (m -n)-dimensional manifold and depends only on the training data {(x i , y i ) : i = 1, . . . , n} and nonlinear model f (•, •).

The total loss on the training set (empirical risk) can be denoted by L(w) = n i=1 L i (w), where L i (•) is the loss on the individual data point i. We assume that the loss L i (•) depends only on the residual, i.e., the difference between the prediction and the true label. In other words,

EQUATION

where l(•) can be any nonnegative differentiable function with l(0) = 0. Typical examples of l(•) include square (l 2 ) loss, Huber loss, etc. We remark that, in the interpolating regime, every parameter vector in the set W renders each individual loss zero, i.e., L i (w) = 0, for all w ∈ W.

3. WARM-UP: REVISITING SGD ON SQUARE LOSS OF LINEAR MODELS:

In this section, we describe the main ideas and results in a simple setting, i.e., stochastic gradient descent (SGD) for the square loss of a linear model, and we revisit some of the results from H ∞ theory (Hassibi et al., 1999; Simon, 2006) . In this case, the data model is y i = x T i w + v i , i = 1, . . . , n (where there is no assumption on v i ) and the loss function is L i (w) =1 2 (y i -x T i w)2 . Assuming the data is indexed randomly, the SGD updates are defined as w i = w i-1 -η∇L i (w i-1 ), where η > 0 is the step size or learning rate. 1 The update in this case can be expressed as

EQUATION

for i ≥ 1 (for i > n, we can either cycle through the data, or select them at random).

Remark. We should point out that, when the step size η is fixed, the SGD recursions have no hope of converging, unless there exists a weight vector w which perfectly interpolates the data {(x i , y i ) : i = 1, . . . , n}. The reason being that, if this is not the case, for any estimated weight vector in SGD there will exist at least one data point that has a nonzero instantaneous gradient and that will therefore move the estimate by a non-vanishing amount. 2 It is for this reason that the results on the convergence of SGD and SMD (Sections 3.3 and 5) pertain to the interpolating regime.

3.1. CONSERVATION OF UNCERTAINTY:

Prior to the i-th step of any optimization algorithm, we have two sources of uncertainty: our uncertainty about the unknown parameter vector w, which we can represent by w -w i-1 , and our uncertainty about the i-th data point (x i , y i ), which we can represent by the noise v i . After the i-th step, the uncertainty about w is transformed to w -w i . But what about the uncertainty in v i ? What is it transformed to? In fact, we will view any optimization algorithm as one which redistributes the uncertainties at time i -1 to new uncertainties at time i. The two uncertainties, or error terms, we will consider are e i and e p,i , defined as follows.

e i := y i -x T i w i-1 , and e p,i :=

EQUATION

)

e i is often referred to as the innvovations and is the error in predicting y i , given the input x i . e p,i is sometimes called the prediction error, since it is the error in predicting the noiseless output x T i w, i.e., in predicting what the best output of the model is. In the absence of noise, e i and e p,i coincide.

One can show that SGD transforms the uncertainties in the fashion specified by the following lemma, which was first noted in (Hassibi et al., 1996) . Lemma 1. For any parameter w and noise values {v i } that satisfy y i = x T i w + v i for i = 1, . . . , n, and for any step size η > 0, the following relation holds for the SGD iterates {w i } given in Eq. ( 3)

EQUATION

As illustrated in Figure 1 , this means that each step of SGD can be thought of as a lossless transformation of the input uncertainties to the output uncertainties, with the specified coefficients.

Once one knows this result, proving it is straightforward. To see that, note that we can write

v i = y i -x T i w as v i = (y i -x T i w i-1 ) -(x T i w -x T i w i-1

). Multiplying both sides by √ η, we have

EQUATION

On the other hand, subtracting both sides of the update rule (3) from w yields

EQUATION

Squaring both sides of ( 6) and ( 7), and subtracting the results leads to Equation (5).

A nice property of Equation ( 5) is that, if we sum over all i = 1, . . . , T , the terms w -w i 2 and w -w i-1 2 on different sides cancel out telescopically, leading to the following important lemma.

Lemma 2. For any parameter w and noise values {v i } that satisfy y i = x T i w + v i for i = 1, . . . , n, any initialization w 0 , any step size η > 0, and any number of steps T ≥ 1, the following relation holds for the SGD iterates {w i } given in Eq. (3)

EQUATION

As we will show next, this identity captures most properties of SGD, and implies several important results in a very transparent fashion. For this reason, this relation can be viewed as a "fundamental identity" for SGD.

3.2. MINIMAX OPTIMALITY OF SGD:

For a given horizon T , consider the following minimax problem:

EQUATION

This minimax problem is motivated by the theory of H ∞ control and estimation (Francis, 1987; Hassibi et al., 1999; Bas ¸ar & Bernhard, 2008) . The denominator of the cost function can be interpreted as the energy of the uncertainties and consists of two terms, w -w 0 2 , the energy of our uncertainty of the unknown weight vector at the beginning of learning when we have not yet observed the data, and

T i=1 v 2

i , the energy of the uncertainty in the measurements. The numerator denotes the energy of the estimation errors in an online setting. The first term, w -w T 2 , is the energy of our uncertainty of the unknown weight vector after we have observed T data points, and the second term,

T i=1 e 2 p,i = T i=1 (x T i w -x T i w i-1 ) 2

, is the energy of the prediction error, i.e., how well we can predict the true uncorrupted output x T i w using measurements up to time i -1. The parameter η weighs the two energy terms relative to each other. In this minimax problem, nature has access to the unknown weight vector w and the noise sequence v i and would like to maximize the energy gain from the uncertainties to prediction errors (so that the estimator behaves poorly), whereas the estimator attempts to minimize the energy gain. Such an estimator is referred to as H ∞ -optimal and is robust because it safeguards against the worst-case noise. It is also conservative-for the exact same reason.3 Theorem 3. For any initialization w 0 , any step size 0 < η ≤ min i 1 xi 2 , and any number of steps T ≥ 1, the stochastic gradient descent iterates {w i } given in Eq. (3) are the optimal solution to the minimax problem (9). Furthermore, the optimal minimax value (achieved by SGD) is 1. This theorem explains the observed robustness and conservatism of SGD. Despite the conservativeness of safeguarding against the worst-case disturbance, this choice may actually be the rational thing to do in situations where we do not have much knowledge about the disturbances, which is the case in many machine learning tasks.

Theorem 3 holds for any horizon T ≥ 1. A variation of this result, i.e., when T → ∞ and without the w -w T 2 term in the numerator, was first shown in (Hassibi et al., 1994; 1996) . In that case, the ratio

η ∞ i=1 e 2 p,i w-w0 2 +η ∞ i=1 v 2 i

in the minimax problem is in fact the H ∞ norm of the transfer operator that maps the unknown disturbances (w -w 0 , { √ ηv i }) to the prediction errors { √ ηe p,i }.

We end this section with a stochastic interpretation of SGD (Hassibi et al., 1996) . Assume that the true weight vector has a normal distribution with mean w 0 and covariance matrix ηI, and that the noise v i are iid standard normal. Then SGD solves

EQUATION

and no exponent larger than 1 2 is possible, in the sense that no estimator can keep the expected cost finite. This means that, in the Gaussian setting, SGD minimizes the expected value of an exponential quadratic cost. The algorithm is thus very adverse to large estimation errors, as they are penalized exponentially larger than moderate ones.

3.3. CONVERGENCE AND IMPLICIT REGULARIZATION:

The over-parameterized (interpolating) linear regression regime is a simple but instructive setting, recently considered in some papers (Gunasekar et al., 2018a; Zhang et al., 2016) . In this setting, we can show that, for sufficiently small step, i.e. 0 < η ≤ min i 1 xi 2 , SGD always converges to a special solution among all the solutions W, in particular to the one with the smallest l 2 distance from w 0 . In other words, if, for example, initialized at zero, SGD implicitly regularizes the solution according to an l 2 norm. This result follows directly from Lemma 2.

To see that, note that in the interpolating case the v i are zero, and we have

EQUATION

for all w ∈ W. By dropping the w -w T 2 term and taking T → ∞, we have η

∞ i=1 2 -η x i 2 e 2 i ≤ w -w 0 2

, which implies that, for 0 < η < min i 2 xi 2 , we must have e i → 0 as i → ∞. When e i = y i -x T i w i-1 goes to zero, the updates in (3) vanish and we get convergence, i.e., w → w ∞ . Further, again because e i → 0, all the data points are being fit, which means w ∞ ∈ W. Moreover, it is again very straightforward to see from (11) that the solution converged to is the one with minimum Euclidean norm from the initial point. To see that, notice that the summation term in Eq. ( 11) is independent of w (it depends only on x i , y i and w 0 ). Therefore, by taking T → ∞ and minimizing both sides with respect to w ∈ W, we get w ∞ = arg min w∈W w -w 0 .

(12)

Once again, this also implies that if SGD is initialized at the origin, i.e., w 0 = 0, then it converges to the minimum-l 2 -norm solution, among all the solutions.

4. MAIN RESULT: GENERAL CHARACTERIZATION OF STOCHASTIC MIRROR DESCENT:

Stochastic Mirror Descent (SMD) (Nemirovskii et al., 1983; Beck & Teboulle, 2003; Cesa-Bianchi et al., 2012; Zhou et al., 2017) is one of the most widely used families of algorithms for stochastic optimization, which includes SGD as a special case. In this section, we provide a characterization of the behavior of general SMD, on general loss functions and general nonlinear models, in terms of a fundamental identity and minimax optimality.

For any strictly convex and differentiable potential ψ(•), the corresponding SMD updates are defined as

EQUATION

where

EQUATION

) is the Bregman divergence with respect to the potential function ψ(•). Note that D ψ (•, •) is nonnegative, convex in its first argument, and that, due to strict convexity, D ψ (w, w ) = 0 iff w = w . Moreover, the updates can be equivalently written as

EQUATION

which are uniquely defined because of the invertibility of ∇ψ (again, implied by the strict convexity of ψ(•)). In other words, stochastic mirror descent can be thought of as transforming the variable w, with a mirror map ∇ψ(•), and performing the SGD update on the new variable. For this reason, ∇ψ(w) is often referred to as the dual variable, while w is the primal variable.

Different choices of the potential function ψ(•) yield different optimization algorithms, which, as we will see, result in different implicit regularizations. To name a few examples: For the potential function ψ(w) = 1 2 w 2 , the Bregman divergence is D ψ (w, w ) = 1 2 w -w 2 , and the update rule reduces to that of SGD. For ψ(w) = j w j log w j , the Bregman divergence becomes the unnormalized relative entropy (Kullback-Leibler divergence) D ψ (w, w ) = j w j log wj w j j w j + j w j , which corresponds to the exponentiated gradient descent (aka the exponential weights) algorithm. Other examples include ψ(w) = 1 2 w 2 Q = 1 2 w T Qw for a positive definite matrix Q, which yields D ψ (w, w ) = 1 2 (w-w ) T Q(w-w ), and the q-norm squared ψ(w) = 1 2 w 2 q , which with 1 p + 1 q = 1 yields the p-norm algorithms (Grove et al., 2001; Gentile, 2003) .

In order to derive an equivalent "conservation law" for SMD, similar to the identity (5), we first need to define a new measure for the difference between the parameter vectors w and w according to the loss function L i (•). To that end, let us define

EQUATION

which is defined in a similar way to a Bregman divergence for the loss function. 4 The difference though is that, unlike the potential function of the Bregman divergence, the loss function

L i (•) = (y i -f (x i , •))

need not be convex, even when (•) is, due to the nonlinearity of f (•, •). As a result, D Li (w, w ) is not necessarily non-negative. The following result, which is the general counterpart of Lemma 1, states the identity that characterizes SMD updates in the general setting. Lemma 4. For any (nonlinear) model f (•, •), any differentiable loss l(•), any parameter w and noise values {v i } that satisfy y i = f (x i , w) + v i for i = 1, . . . , n, and any step size η > 0, the following relation holds for the SMD iterates {w i } given in Eq. (15)

EQUATION

for all i ≥ 1, where

EQUATION

)

The proof is provided in Appendix A. Note that E i (w i , w i-1 ) is not a function of w. Furthermore, even though it does not have to be nonnegative in general, for η sufficiently small, it becomes nonnegative, because the Bregman divergence D ψ (., .) is nonnegative.

Summing Equation ( 17) over all i = 1, . . . , T leads to the following identity, which is the general counterpart of Lemma 2.

Lemma 5. For any (nonlinear) model f (•, •), any differentiable loss l(•), any parameter w and noise values {v i } that satisfy y i = f (x i , w) + v i for i = 1, . . . , n, any initialization w 0 , any step size η > 0, and any number of steps T ≥ 1, the following relation holds for the SMD iterates {w i } given in Eq. ( 15)

EQUATION

)

We should reiterate that Lemma 5 is a fundamental property of SMD, which allows one to prove many important results, in a direct way.

In particular, in this setting, we can show that SMD is minimax optimal in a manner that generalizes Theorem 3 of Section 3, in the following 3 ways: 1) General potential ψ(•), 2) General model f (•, •), and 3) General loss function l(•). The result is as follows.

Theorem 6. Consider any (nonlinear) model f (•, •), any non-negative differentiable loss l(•) with the property l(0) = l (0) = 0, and any initialization w 0 . For sufficiently small step size, i.e., for any η > 0 for which ψ(w) -ηL i (w) is convex for all i, and for any number of steps T ≥ 1, the SMD iterates {w i } given by Eq. (15), w.r.t. any strictly convex potential ψ(•), is the optimal solution to the following minimization problem

EQUATION

)

Furthermore, the optimal value (achieved by SMD) is 1.

The proof is provided in Appendix B. For the case of square loss and a linear model, the result reduces to the following form.

Corollary 7. For L i (w) = 1 2 (y i -x T i w) 2 , for any initialization w 0 , any sufficiently small step size, i.e., 0 < η ≤ α xi 2 , and any number of steps T ≥ 1, the SMD iterates {w i } given by Eq. (15), w.r.t. any α-strongly convex potential ψ(•), is the optimal solution to

EQUATION

)

The optimal value (achieved by SMD) is 1.

We should remark that Theorem 6 and Corollary 7 generalize several known results in the literature.

In particular, as mentioned in Section 3, the result of (Hassibi et al., 1994 ) is a special case of Corollary 7 for ψ(w) = 1 2 w 2 . Furthermore, our result generalizes the result of (Kivinen et al., 2006) , which is the special case for the p-norm algorithms, again, with square loss and a linear model. Another interesting connection to the literature is that it was shown in (Hassibi & Kailath, 1995) that SGD is locally minimax optimal, with respect to the H ∞ norm. Strictly speaking, our result is not a generalization of that result; however, Theorem 6 can be interpreted as SGD/SMD being globally minimax optimal, but with respect to different metrics in the numerator and denominator. Namely, the uncertainty about the weight vector w is measured by the Bregman divergence of the potential, the uncertainty about the noise by the loss, and the prediction error by the "Bregman-divergencelike" expression of the loss.

5. CONVERGENCE AND IMPLICIT REGULARIZATION IN OVER-PARAMETERIZED MODELS:

In this section, we show some of the implications of the theory developed in the previous section.

In particular, we show convergence and implicit regularization, in the over-parameterized (so-called interpolating) regime5 , for general SMD algorithms. We first consider the linear interpolating case, which has been studied in the literature, and show that the known results follow naturally from our Lemma 5. Further, we shall obtain some new convergence results. Finally, we discuss the implications for nonlinear models, and argue that the same results hold qualitatively in highlyoverparameterized settings, which is the typical scenario in deep learning.

5.1. OVER-PARAMETERIZED LINEAR MODELS:

In this setting, the v i are zero, W = w | y i = x T i w, i = 1, . . . , n , and L i (w) = l(y i -x T i w), with any differentiable loss l(•). Therefore, Eq. ( 19) reduces to

EQUATION

for all w ∈ W, where

EQUATION

which is notably independent of w. As a result, we can easily minimize both sides of Eq. ( 22) with respect to w ∈ W, which for T → ∞ leads to the following result. Proposition 8. For any differentiable loss l(•), any initialization w 0 , and any step size η, consider the SMD iterates given in Eq. 

An equivalent form of Proposition 8 has been shown recently in, e.g., (Gunasekar et al., 2018a) . 6Other implicit regularization results have been shown in (Gunasekar et al., 2018b; Soudry et al., 2017) for classification problems, which are not discussed here. Note that the result of (Gunasekar et al., 2018a) does not say anything about whether the algorithm converges or not. However, our fundamental identity of SMD (Lemma 5) allows us to also establish convergence to the regularized point, for some common cases, which will be shown next.

What Proposition 8 says is that depending on the choice of the potential function ψ(•), the optimization algorithm can perform an implicit regularization without any explicit regularization term.

In other words, for any desired regularizer, if one chooses a potential function that approximates the regularizer, we can run the optimization without explicit regularization, and if it converges to a solution, the solution must be the one with the minimum potential.

In principle, one can choose the potential function in SMD for any desired convex regularization.

For example, we can find the maximum entropy solution by taking the potential to be the negative entropy. Another illustrative example follows.

Example [Compressed Sensing]::

In compressed sensing, one seeks the sparsest solution to an under-determined (over-parameterized) system of linear equations. The surrogate convex problem one solves is: min w 1 subject to y i = x T i w, i = 1, . . . n (28)

One cannot choose ψ(w) = w 1 , since it is neither differentiable nor strictly convex. However, ψ(w) = w 1+ , for any > 0, can be used. Figure 4 shows a compressed sensing example, with n = 50, m = 100, and sparsity k = 10. SMD was used with a step size of η = 0.001 and the potential function was ψ(•) = • 1.1 . SMD converged to the true sparse solution after around 10,000 iterations. On this example, it was an order of magnitude faster than standard l 1 optimization. Next we establish convergence to the regularized point for the convex case. Proposition 9. Consider the following two cases.

(i) l(•) is differentiable and convex and has a unique root at 0, ψ(•) is strictly convex, and η > 0 is such that ψ -ηL i is convex for all i.

(ii) l(•) is differentiable and quasi-convex, l (•) is zero only at zero, ψ(•) is α-strongly convex, and

0 < η ≤ min i α|yi-x T i wi-1| xi 2 |l (yi-x T i wi-1)| .

If either (i) or (ii) holds, then for any w 0 , the SMD iterates given in Eq. (15) converge to

EQUATION

)

The proof is provided in Appendix C.

5.2. DISCUSSION OF HIGHLY OVER-PARAMETERIZED NONLINEAR MODELS:

Let us consider the highly-overparameterized nonlinear model

EQUATION

where by highly-overparameterized we mean m n. Since the model is highly over-parameterized, it is assumed that we can perfectly interpolate the data points (x i , y i ) so that the noise v i is zero. In this case, the set of parameter vectors that interpolate the data is given by W = {w ∈ R m | y i = f (x i , w), i = 1, . . . , n}, and Eq. ( 19), again, reduces to

EQUATION

for all w ∈ W. Our proofs of convergence and implicit regularization for SGD and SMD in the linear case relied on two facts: (i) D Li (w, w i-1 ) was non-negative (this allowed us to show convergence), and (ii) D Li (w, w i-1 ) was independent of w (this allowed us to show implicit regularization). Unfortunately, neither of these hold in the nonlinear case.

However, they do hold in a local sense. In other words, (i) D Li (w, w i-1 ) ≥ 0 for w i-1 "close enough" to w (see Figure 3 ), and (ii) D Li (w, w i-1 ) is weakly dependent on w for w i-1 "close enough." (Both statements can be made precise.)

Figure 3 : Non-negativity of D Li (w, w i-1 ) for w i-1 "close enough" to w.

Now define:

EQUATION

)

Then one can show the following result. Theorem 10. There exists an > 0, such that if w * -w 0 < , then for sufficiently small step size η > 0:

1. SMD iterates converge to a point w ∞ ∈ W 2. w ∞ -w * = o( )

This shows that if the initial condition is close enough, then we have convergence to a point w ∞ that interpolates the data, and that w ∞ is an order of magnitude closer to w * (the implicitly regularized solution) than the initial w 0 was. At first glance, this result seems rather dissatisfying. It relies on w 0 being close to the manifold W which appears hard to guarantee. We would now like to argue that in deep learning w 0 being close to W is often the case.

In the highly-overparameterized regime, m n, and so the dimension of the manifold W is m -n, which is very large. Now if the x i are sufficiently random, then the tangent space to W at w * will be a randomly oriented affine subspace of dimension m -n. This means that any randomly chosen w 0 will whp have a very large component when projected onto W. In particular, it can be shown that w * -w 0

2 = O( n m ) • y -f (x, w) 2

, where y = vec(y i , i = 1, . . . , n) and f (x, w) = vec(f (x i , w), i = 1, . . . , n). Thus, we may expect that, when m n, the distance of any randomly chosen w 0 to W will be small and so SMD will converge to a point on W that approximately performs implicit regularization.

The gist of the argument is that (i) When m n, any random initial condition is "close" to the n -m dimensional solution manifold W, (ii) when w 0 is "close" to w * , then SMD converges to a point w ∞ ∈ W, (iii) w ∞ is "an order of magnitude closer" to w * than w 0 was, and (iv) thus, when highly overparamatrized, SMD converges to a point that exhibits implicit regularization.

Of course, this was a very heuristic argument that merits a much more careful analysis. But it is suggestive of the fact that SGD and SMD, when performed on highly-overparameterized nonlinear models, as occurs in deep learning, may exhibit implicit regularization.

6. CONCLUDING REMARKS:

We should remark that all the results stated throughout the paper extend to the case of time-varying step size η i , with minimal modification. In particular, it is easy to show that in this case, the identity (the counterpart of Eq. ( 19)) becomes

EQUATION

where E i (w i , w i-1 ) = D ψ (w i , w i-1 ) -η i D Li (w i , w i-1 ) + η i L i (w i ). As a consequence, our main result will be the same as in Theorem 6, with the only difference that the small-step-size condition in this case is the convexity of ψ(w) -η i L i (w) for all i, and the SMD with time-varying step size will be the optimal solution to the following minimax problem

EQUATION

)

Similarly, the convergence and implicit regularization results can be proven under the same conditions (See Appendix D for more details on the time-varying case).

This paper opens up a variety of important directions for future work. Most of the analysis developed here is general, in terms of the model, the loss function, and the potential function. Therefore, it would be interesting to study the implications of this theory for specific classes of models (such as different neural networks), specific losses, and specific mirror maps (which induce different regularization biases). Something for future work.

Clearly, the optimal solution(s) and the optimal value of this problem can, and will, be a function of {w i }. Similarly, we can also choose feasible points that depend on {w i }. Any choice of a feasible point ( ŵ, {v i }) gives a lower bound on the value of the problem. Before choosing a feasible point, let us first expand the D Li (w, w i-1 ) term in the numerator, according to its definition.

EQUATION

where we have used the fact that l(y i -f i (w)) = l(v i ) for all consistent w, in the first term. Now, we choose a feasible point as follows

EQUATION

)

where ŵ is the choice of w, as will be described soon. The reason for choosing this value for the noise is that it "fools" the estimator by making its loss on the corresponding data point zero. In other words, for this choice, we have

D Li (w, w i-1 ) = l(v i ) -l(0) + l (0)∇f (w i-1 ) T ( ŵ -w i-1 ) = l(v i )

because l(0) = l (0) = 0. It should be clear at this point that this choice makes the second terms in the numerator and the denominator equal, independent of the choice of ŵ. What remains to do, in order to show the 1 lower-bound, is to take care of the other two terms, i.e., D ψ (w, w T ) and D ψ (w, w 0 ). As we would like to make the ratio equal to one, we would like to have D ψ (w, w T ) = D ψ (w, w 0 ), which is equivalent to having

ψ(w) -ψ(w T ) -∇ψ(w T ) T (w -w T ) = ψ(w) -ψ(w 0 ) -∇ψ(w 0 ) T (w -w 0 )

which is, in turn, equivalent to

EQUATION

)

Since ∇ψ is an invertible function, ∇ψ(w T ) -∇ψ(w 0 ) = 0, if w T = w 0 . Therefore, the above equation has a solution for w, if w T = w 0 . As a result, choosing ŵ to be a solution to (46) makes D ψ ( ŵ, w T ) = D ψ ( ŵ, w 0 ), if w T = w 0 . For the case when w T = w 0 , it is trivial that D ψ ( ŵ, w T ) = D ψ ( ŵ, w 0 ) for any choice of ŵ. In this case, we only need to choose ŵ to be different from w 0 , to avoid making the ratio 0 0 . Hence, we have the following choice ŵ = a solution of (46) for w T = w 0 w 0 + δw for some δw = 0 for w T = w 0 (47)

Choosing the feasible point ŵ, {v i } according to ( 47) and ( 45) leads to

max w,{vi} D ψ (w, w T ) + η T i=1 D Li (w, w i-1 ) D ψ (w, w 0 ) + η T i=1 l(v i ) ≥ D ψ ( ŵ, w T ) + η T i=1 l(f i (w i-1 ) -f i ( ŵ)) D ψ ( ŵ, w 0 ) + η T i=1 l(f i (w i-1 ) -f i ( ŵ)) . (48)

Taking the minimum of both sides with respect to {w i }, we have

min {wi} max w,{vi} D ψ (w, w T ) + η T i=1 D Li (w, w i-1 ) D ψ (w, w 0 ) + η T i=1 l(v i ) ≥ min {wi} D ψ ( ŵ, w T ) + η T i=1 l(f i (w i-1 ) -f i ( ŵ)) D ψ ( ŵ, w 0 ) + η T i=1 l(f i (w i-1 ) -f i ( ŵ)) = 1. (49)

The equality to 1 comes from the fact the that the optimal solution of the minimization either has w * T = w 0 or w * T = w 0 , and in both cases the ratio is equal to 1.

2. Now we prove that, under the small step size condition (convexity of ψ(w) -ηL i (w) for all i), SMD makes the minimax value at most 1, which means that it is indeed an optimal solution. Recall from Lemma 5 that

D ψ (w, w 0 ) + η T i=1 l(v i ) = D ψ (w, w T ) + T i=1 E i (w i , w i-1 ) + η T i=1 D Li (w, w i-1 )

,

where E i (w i , w i-1 ) = D ψ (w i , w i-1 ) -ηD Li (w i , w i-1 ) + ηL i (w i ).

It is easy to check that when ψ(w)

-ηL i (w) is convex, D ψ (w i , w i-1 ) -ηD Li (w i , w i-1 )

is in fact a Bregman divergence (i.e. the Bregman divergence with respect to the potential ψ(w) -ηL i (w)), and therefore it is nonnegative for any w i and w i-1 . Furthermore, we know that the loss L i (w i ) is also nonnegative for all w i . It follows that E i (w i , w i-1 ) is nonnegative for all values of w i , w i-1 and i. As a result, we have the following bound.

EQUATION

Since the Bregman divergence D ψ (w, w 0 ) and the loss l(v i ) are nonnegative, the left-hand side expression is nonnegative, and it follows that

EQUATION

In fact, this means that independent of the choice of the maximizer (i.e. for all {v i } and w), as long as the step size condition is met, SMD makes the ratio less than or equal to 1.

Combining the results of 1 and 2 above concludes the proof.

B.1 PROOF OF THEOREM 3:

Proof. This result is a special case of Theorem 6, which was proven above. In this case, ψ(w)

= 1 2 w 2 , f (x i , w) = x T i w, and l(z) = 1 2 z 2 . Therefore, D ψ (w, w T ) = 1 2 w -w T 2 , D ψ (w, w 0 ) = 1 2 w -w 0 2 , D Li (w, w i-1 ) = 1 2 (x T i w -x T i w i-1

) 2 , and l(v i ) = 1 2 v 2 i , which leads to the result.

C PROOF OF PROPOSITION 9:

Proof. To prove convergence, we appeal again to Equation ( 22), i.e.

D ψ (w, w 0 ) = D ψ (w, w T ) + T i=1 (E i (w i , w i-1 ) + ηD Li (w, w i-1 )) ,

for all w ∈ W. We prove the two cases separately.

1. The proof of case (i) is straightforward. When l(•) is differentiable and convex, L i is also convex, and therefore D Li (w, w i-1 ) is nonnegative. Moreover, when ψ -ηL i is convex, E i (w i , w i-1 ) is also nonnegative. Therefore, the entire summand in Eq. ( 52) is nonnegative, and has to go to zero for i → ∞. That is because as T → ∞, the sum should remain bounded, i.e., ∞ i=1 (E i (w i , w i-1 ) + ηD Li (w, w i-1 )) ≤ D ψ (w, w 0 ). As a result of the non-negativity of both terms in the sum, we have both E i (w i , w i-1 ) → 0 and D Li (w, w i-1 ) → 0 as i → ∞, the latter of which implies L i (w i-1 ) → 0. This implies that the updates in (15) vanish and we get convergence, i.e., w → w ∞ . Further, again because L i (w i-1 ) → 0, and 0 is the unique root of l(•), all the data point are being fit, which means w ∞ ∈ W.

Lemma 11. For any (nonlinear) model f (•, •), any differentiable loss l(•), any parameter w and noise values {v i } that satisfy y i = f (x i , w) + v i for i = 1, . . . , n, any initialization w 0 , any step size sequence {η i }, and any number of steps T ≥ 1, the following relation holds for the SMD iterates {w i } given in Eq. (64) D ψ (w, w 0 ) + T i=1 η i l(v i ) = D ψ (w, w T ) + T i=1 (E i (w i , w i-1 ) + η i D Li (w, w i-1 )) ,

Proof. The proof is straightforward by summing the following equation for all i = 1, . . . , T D ψ (w, w i-1 ) + η i l(v i ) = D ψ (w, w i ) + E i (w i , w i-1 ) + η i D Li (w, w i-1 ),

which can be easily shown in the same way as in the proof of Lemma 4 in Appendix A.

Theorem 12. Consider any general model f (•, •), and any differentiable loss function l(•) with property l(0) = l (0) = 0. For sufficiently small step size, i.e., for any sequence {η i } for which ψ(w) -η i L i (w) is convex for all i, the SMD iterates {w i } given by Eq. ( 64) are the optimal solution to the following minimization problem (67)

Furthermore, the optimal value (achieved by SMD) is 1.

Proof. The proof is similar to that of Theorem 6, as presented in Appendix B. The argument for the upper-bound of 1 is exactly the same. For the second part of the proof, we use the previous Lemma. It follows from the convexity of ψ(w) -η i L i (w) that E i (w i , w i-1 ) ≥ 0, and as a result we have

EQUATION

for SMD updates, which concludes the proof.

The convergence and implicit regularization results hold similarly, and can be formally stated as follows.

Proposition 13. Consider the following two cases.

(i) l(•) is differentiable and convex and has a unique root at 0, ψ(•) is strictly convex, and the positive sequence {η i } is such that ψ -η i L i is convex for all i.

(ii) l(•) is differentiable and quasi-convex and has zero derivative only at 0, ψ(•) is α-strongly convex, and 0 < η i ≤ α|yi-x T i wi-1| xi 2 |l (yi-x T i wi-1)| for all i. If either (i) or (ii) holds, then for any initialization w 0 , the SMD iterates given in Eq. 

Proof. The proof is similar to that of Proposition 9, as provided in Appendix C.

:

For the sake of simplicity of presentation, we present the results for constant step size. We show in the appendix that all the results extend to the case of time-varying step-size.

Of course, one may get convergence by having a vanishing step size ηi → 0. However, in this case, convergence is not surprising-since, effectively, after a while the weights are no longer being updated-and the more interesting question is "what" the recursion converges to.

The setting described is somewhat similar to the setting of online learning, where one considers the relative performance of an online learner who needs to predict, compared to a clairvoyant one who has access to the entire data set(Shalev-Shwartz, 2012;Hazan, 2016). In online learning, the relative performance is described as a difference, rather than as a ratio in H ∞ theory, and is referred to as regret.

It is easy to verify that for linear models and quadratic loss we obtain DL i (w, w ) = (x T i w -x T i w ) 2 .

In the classical under-parameterized (online streaming) case with white noise, the same theory can be used to establish convergence to the true parameter under the so-called Robbins-Monro conditions ( ∞ i=1 ηi = ∞, ∞ i=1 η 2 i < ∞) in a very direct and simple way (seeAzizan & Hassibi (2019)).

To be precise, the authors in(Gunasekar et al., 2018a) assume convergence to a global minimizer of the loss function L(w) = n i=1 l(yi -x T i w), which with their assumption of the loss function l(•) having a unique finite root is equivalent to assuming convergence to a point w∞ ∈ W.

Figures and Tables:

Figure 1: Illustration of Lemma 1. Each step of SGD can be viewed as a transformation of the uncertainties with the right coefficients.

(15) with respect to any strictly convex potential ψ(•). If the iterates converge to a solution w ∞ ∈ W, then w ∞ = arg min w∈W D ψ (w, w 0 ). (26) Remark. In particular, for the initialization w 0 = arg min w∈R m ψ(w), if the iterates converge to a solution w ∞ ∈ W, then w ∞ = arg min w∈W ψ(w).

Figure 2: The training loss and actual error of stochastic mirror descent for compressed sensing. SMD recovers the actual sparse signal.

D ψ (w, w T ) + T i=1 η i D Li (w, w i-1 ) D ψ (w, w 0 ) + T i=1 η i l(v i ).

(64) converge to w ∞ = arg min w∈W D ψ (w, w 0 ).

Peer review (iclr2019):
 The authors look at SGD, and SMD updates applied to various models and loss functions. They derive a fundamental identity lemma 2 for the case of linear model and squared loss + SGD and in general for non-linear models+ SMD + non squared loss functions. The main results shown are 1. SGD is optimal in a certain sense for squared loss and linear model. 2. SGD always converges to a solution closest to the starting point. 3. SMD when it converges, converges to a point closest to the starting point in the bregman divergence. The convergence of SMD iterates is shown for certain learning scenarios. Pros: Shows implicit regularization properties for models beyond linear case. Cons: 1. The notion of optimality is w.r.t. a metric that is pretty non-standard and it was not clear to me as to why the metric is important to study (the ratio metric in eq 9). 2. The result is not very surprising since SMD is pretty much a gradient descent w.r.t a different distance metric.