Year: 2021

Abstract:
Deep neural networks (DNNs) have shown great success in many machine learning tasks. Their training is challenging since the loss surface of the network architecture is generally non-convex, or even non-smooth. How and under what assumptions is guaranteed convergence to a global minimum possible? We propose a reformulation of the minimization problem allowing for a new recursive algorithmic framework. By using bounded style assumptions, we prove convergence to an ε-(global) minimum using Õ(1/ε) gradient computations. Our theoretical foundation motivates further study, implementation, and optimization of the new algorithmic framework and further investigation of its non-standard bounded style assumptions. This new direction broadens our understanding of why and under what circumstances training of a DNN converges to a global minimum.

heading: 1 INTRODUCTION
text: In recent years, deep neural networks (DNNs) have shown a great success in many machine learning tasks. However, training these neural networks is challenging since the loss surface of network architecture is generally non-convex, or even non-smooth. Thus, there have been a long-standing question on how optimization algorithms may converge to a global minimum. Many previous work have investigated Gradient Descent algorithm and its stochastic version for over-parameterized setting (Arora et al., 2018; Soudry et al., 2018; Allen-Zhu et al., 2019; Du et al., 2019a; Zou & Gu, 2019). Although these works have shown promising convergence results under certain assumptions, there is still a lack of new efficient methods that can guarantee global convergence for machine learning optimization. In this paper, we address this problem using a different perspective. Instead of analyzing the traditional finite-sum formulation, we adopt a new composite formulation that exactly depicts the structure of machine learning where a data set is used to learn a common classifier.
Representation. Let { (x(i), y(i)) }n i=1
be a given training set with x(i) ∈ Rm, y(i) ∈ Rc, we investigate the following novel representation for deep learning tasks:
min w∈Rd
{ F (w) = 1
n n∑ i=1 φi(h(w; i))
} , (1)
where h(·; i) : Rd → Rc, i ∈ [n] = {1, . . . , n}, is the classifier for each input data x(i); and φi : Rc → R, i ∈ [n], is the loss function corresponding to each output data y(i). Our composite formulation (1) is a special case of the finite-sum problem minw∈Rd { F (w) = 1n ∑n i=1 f(w; i) } where each individual function f(·; i) is a composition of the loss function φi and the classifier h(·; i). This problem covers various important applications in machine learning, including logistic regression and neural networks. The most common approach for the finite-sum problem is using first-order methods such as (stochastic) gradient algorithms and making assumptions on the component functions f(·; i). As an alternative, we further investigate the structure of the loss function φi and narrow our assumption on the classifier h(·; i). For the purpose of this work, we first consider convex and Lipschitz-smooth loss functions while the classifiers can be non-convex. Using this representation, we propose a new framework followed by two algorithms that guarantee global convergence for the minimization problem.
Algorithmic Framework. Representation (1) admits a new perspective. Our key insight is to (A) define z(t)i = h(w (t); i), where t is an iteration count of the outer loop in our algorithmic framework.
Next (B), we want to approximate the change z(t+1)i − z (t) i in terms of a step size times the gradient ∇φi(z(t)i ) = (∂φi(z)/∂za)a∈[c] ∣∣ z=z
(t) i
,
and (C) we approximate the change h(w(t+1); i)− h(w(t); i) in terms of the first order derivative
H (t) i = (∂ha(w; i)/∂wb)a∈[c],b∈[d] ∣∣ w=w(t) .
Finally, we combine (A), (B), and (C) to equate the approximations of z(t+1)i − z (t) i and h(w(t+1); i) − h(w(t); i). This leads to a recurrence on w(t) of the form w(t+1) = w(t) − η(t)v(t), where η(t) is a step size and which involves computing v(t) by solving a convex quadratic subproblem, see the details in Section 4. We explain two methods for approximating a solution for the derived subproblem. We show how to approximate the subproblem by transforming it into a strongly convex problem by adding a regularizer which can be solved in closed form. And we show how to use Gradient Descent (GD) on the subproblem to find an approximation v(t) of its solution.
Convergence Analysis. Our analysis introduces non-standard bounded style assumptions. Intuitively, we assume that our convex and quadratic subproblem has a bounded solution. This allows us to prove a total complexity of Õ( 1ε3 ) to find an ε-(global) solution that satisfies F (ŵ)− F∗ ≤ ε, where F∗ is the global minimizer of F . Our analysis applies to a wide range of applications in machine learning: Our results hold for squared loss and softmax cross-entropy loss and applicable for a range of activation functions in DNN as we only assume that the h(·; i) are twice continuously differentiable and their Hessian matrices (second order derivatives) as well as their gradients (first order derivatives) are bounded.
Contributions and Outline. Our contributions in this paper can be summarized as follows.
• We propose a new representation (1) for analyzing the machine learning minimization problem. Our formulation utilizes the structure of machine learning tasks where a training data set of inputs and outputs is used to learn a common classifier. Related work in Section 2 shows how (1) is different from the classical finite-sum problem.
• Based on the new representation we propose a novel algorithm framework. The algorithmic framework approximates a solution to a subproblem for which we show two distinct approaches.
• For general DNNs and based on bounded style assumptions, we prove a total complexity of Õ( 1ε3 ) to find an ε-(global) solution that satisfies F (ŵ)−F∗ ≤ ε, where F∗ is the global minimizer of F .
We emphasize that our focus is on developing a new theoretical foundation and that a translation to a practical implementation with empirical results is for future work. Our theoretical foundation motivates further study, implementation, and optimization of the new algorithmic framework and further investigation of its non-standard bounded style assumptions. This new direction broadens our understanding of why and under what circumstances training of a DNN converges to a global minimum.
The rest of this paper is organized as follows. Section 2 discusses related work. Section 3 describes our setting and deep learning representation. Section 4 explains our key insight and derives our Framework 1. Section 5 presents our algorithms and their global convergence. All technical proofs are deferred to the Appendix.

heading: 2 RELATED WORK
text: Formulation for Machine Learning Problems. The finite-sum problem is one of the most important and fundamental problems in machine learning. Analyzing this model is the most popular approach in the machine learning literature and it has been studied intensively throughout the years (Bottou et al., 2018; Reddi et al., 2016; Duchi et al., 2011b). Our new formulation (1) is a special case of the finite-sum problem, however, it is much more complicated than the previous model since it involves the data index i both inside the classifiers h(·; i) and the loss functions φi. For a comparison, previous works only consider a common loss function l(ŷ, y) for the predicted value
ŷ and output data y (Zou et al., 2018; Soudry et al., 2018). Our modified version of loss function φi is a natural setting for machine learning. We note that when h(w; i) is the output produced by a model, our goal is to match this output with the corresponding target y(i). For that reason, the loss function for each output has a dependence on the output data y(i), and is denoted by φi. This fact reflects the natural setting of machine learning where the outputs are designed to fit different targets, and the optimization process depends on both outer function φi and inner functions h(·; i). This complication may potentially bring a challenge to theoretical analysis. However, with separate loss functions, we believe this model will help to exploit better the structure of machine learning problems and gain more insights on the neural network architecture.
Other related composite optimization models are also investigated thoroughly in (Lewis & Wright, 2016; Zhang & Xiao, 2019; Tran-Dinh et al., 2020). Our model is different from these works as it does not have a common function wrapping outside the finite-sum term, as in (Lewis & Wright, 2016). Note that a broad class of variance reduction algorithms (e.g. SAG (Le Roux et al., 2012), SAGA (Defazio et al., 2014), SVRG (Johnson & Zhang, 2013), SARAH (Nguyen et al., 2017)) is designed specifically for the finite-sum formulation and is known to have certain benefits over Gradient Descent. In addition, the multilevel composite problem considered in (Zhang & Xiao, 2021) also covers empirical risk minimization problem. However our formulation does not match their work since our inner function h(w; i) is not an independent expectation over some data distribution, but a specific function that depends on the current data.
Global Convergence for Neural Networks. A recent popular line of research is studying the dynamics of optimization methods on some specific neural network architectures. There are some early works that show the global convergence of Gradient Descent (GD) for simple linear network and two-layer network (Brutzkus et al., 2018; Soudry et al., 2018; Arora et al., 2019; Du et al., 2019b). Some further works extend these results to deep learning architectures (Allen-Zhu et al., 2019; Du et al., 2019a; Zou & Gu, 2019). These theoretical guarantees are generally proved for the case when the last output layer is fixed, which is not standard in practice. A recent work (Nguyen & Mondelli, 2020) prove the global convergence for GD when all layers are trained with some initial conditions. However, these results are for neural networks without bias neurons and it is unclear how these analyses can be extended to handle the bias terms of deep networks with different activations. Our novel framework and algorithms do not exclude learning bias layers as in (Nguyen & Mondelli, 2020).
Using a different algorithm, Brutzkus et al. (2018) investigate Stochastic Gradient Descent (SGD) for two-layer networks in a restricted linearly separable data setting. This line of research continues with the works from Allen-Zhu et al. (2019); Zou et al. (2018) and later with Zou & Gu (2019). They justify the global convergence of SGD for deep neural networks for some probability depending on the number of input data and the initialization process.
Over-Paramaterized Settings and other Assumptions for Machine Learning. Most of the modern learning architectures are over-parameterized, which means that the number of parameters are very large and often far more than the number of input data. Some recent works prove the global convergence of Gradient Descent when the number of neurons are extensively large, e.g. (Zou & Gu, 2019) requires Ω(n8) neurons for every hidden layer, and (Nguyen & Mondelli, 2020) improves this number to Ω(n3). If the initial point satisfies some special conditions, then they can show a better dependence of Ω(n). In Allen-Zhu et al. (2019), the authors initialize the weights using a random Gaussian distribution where the variance depends on the dimension of the problem. In non-convex setting, they prove the convergence of SGD using the assumption that the dimension depends inversely on the tolerance . We will discuss how these over-paramaterized settings might be a necessary condition to develop our theory.
Other standard assumptions for machine learning include the bounded gradient assumption (Nemirovski et al., 2009; Shalev-Shwartz et al., 2007; Reddi et al., 2016; Tran et al., 2021). It is also common to assume all the iterations of an algorithm stays in a bounded domain (Duchi et al., 2011a; Levy et al., 2018; Gürbüzbalaban et al., 2019; Reddi et al., 2018; Vaswani et al., 2021). Since we are analyzing a new composite formulation, it is understandable that our assumptions may also not be standard. However, we believe that there is a strong connection between our assumptions and the traditional setting of machine learning. We will discuss this point more clearly in Section 4.

heading: 3 BACKGROUND
text: In this section, we discuss our formulation and notations in detail. Although this paper focuses on deep neural networks, our framework and theoretical analysis are general and applicable for other learning architectures.
Deep Learning Representation. Let {(x(i), y(i))}ni=1 be a training data set where x(i) ∈ Rm is a training input and y(i) ∈ Rc is a training output. We consider a fully-connected neural network with L layers, where the l-th layer, l ∈ {0, 1, . . . , L}, has nl neurons. We represent layer 0-th and L-th layer as input and output layers, respectively, that is, n0 = d and nL = c. For l ∈ {1, . . . , L}, let W (l) ∈ Rnl−1×nl and b(l) ∈ Rnl , where {(W (l), b(l))Ll=1} represent the parameters of the neural network. A classifier h(w; i) is formulated as
h(w; i) = W (L)>σL−1(W (L−1)>σL−2(. . . σ1(W (1)>x(i) + b(1)) . . . ) + b(L−1)) + b(L),
wherew = vec({W (1), b(1), . . . ,W (L), b(L)}) ∈ Rd is the vectorized weight and {σl}L−1l=1 are some activation functions. The most common choices for machine learning are ReLU, sigmoid, hyperbolic tangent and softplus. For j ∈ [c], hj(·; i) : Rd → R denotes the component function of the output h(·; i), for each data i ∈ [n] respectively. Moreover, we define h∗i = arg minz∈Rc φi(z), i ∈ [n]. Loss Functions. The well-known loss functions in neural networks for solving classification and regression problems are softmax cross-entropy loss and square loss, respectively:
(Softmax) Cross-Entropy Loss: F (w) = 1n ∑n i=1 f(w; i) with
f(w; i) = −y(i)> log(softmax(h(w; i))). (2)
Squared Loss: F (w) = 1n ∑n i=1 f(w; i) with
f(w; i) = 1
2 ‖h(w; i)− y(i)‖2. (3)
We provide some basic definitions in optimization theory to support our theory.
Definition 1 (L-smooth). Function φ : Rc → R is Lφ-smooth if there exists a constant Lφ > 0 such that, ∀x1, x2 ∈ Rc,
‖∇φ(x1)−∇φ(x2)‖ ≤ Lφ‖x1 − x2‖. (4)
Definition 2 (Convex). Function φ : Rc → R is convex if ∀x1, x2 ∈ Rc,
φ(x1)− φ(x2) ≥ 〈∇φ(x2), x1 − x2〉. (5)
The following corollary shows the properties of softmax cross-entropy loss (2) and squared loss (3).
Corollary 1. For softmax cross-entropy loss (2) and squared loss (3), there exist functions h(·; i) : Rd → Rc and φi : Rc → R such that, for i ∈ [n], φi(z) is convex and Lφ-smooth with Lφ = 1, and
f(w; i) = φi(h(w; i)) = φi(z) ∣∣ z=h(w;i) . (6)

heading: 4 NEW ALGORITHM FRAMEWORK

heading: 4.1 KEY INSIGHT
text: We assume f(w; i) = φi(h(w; i)) with φi convex and Lφ-smooth. Our goal is to utilize the convexity of the outer function φi. In order to simplify notation, we write ∇zφi(h(w(t); i)) instead of ∇zφi(z) ∣∣ z=h(w(t);i) and denote z(t)i = h(w (t); i). Starting from the current weight w(t), we would like to find the next point w(t+1) that satisfies the following approximation for all i ∈ [n]:
h(w(t+1); i) = z (t+1) i ≈ z (t) i − α (t) i ∇zφi(z (t) i ) = h(w (t); i)− α(t)i ∇zφi(h(w (t); i)). (7)
We can see that this approximation is a “noisy” version of a gradient descent update for every function φi, simultaneously for all i ∈ [n]. In order to do this, we use the following update
w(t+1) = w(t) − η(t)v(t), (8)
where η(t) > 0 is a learning rate and v(t) is a search direction that helps us approximate equation (7). If the update term η(t)v(t) is small enough, and if h(·; i) has some nice smooth properties, then from basic calculus we have the following approximation:
h(w(t+1); i) = h(w(t) − η(t)v(t); i) ≈ h(w(t); i)−H(t)i ( η(t)v(t) ) , (9)
where H(t)i is a matrix in Rc×d with first-order derivatives. Motivated by approximations (7) and (9), we consider the following optimization problem:
v (t) ∗ = arg min v∈Rd 1 2 1 n n∑ i=1 ‖H(t)i ( η(t)v ) − α(t)i ∇zφi(h(w (t); i))‖2. (10)
Hence, by solving for the solution v(t)∗ of problem (10) we are able to find a search direction for the key approximation (7). This yields our new algorithmic Framework 1, see below.
Framework 1 New Algorithm Framework
Initialization: Choose an initial point w(0) ∈ Rd; for t = 0, 1, · · · , T − 1 do
Solve for an approximation v(t) of the solution v(t)∗ of the problem in (10)
v (t) ∗ = arg min v∈Rd 1 2 1 n n∑ i=1 ‖η(t)H(t)i v − α (t) i ∇zφi(h(w (t); i))‖2
Update w(t+1) = w(t) − η(t)v(t) end for

heading: 4.2 TECHNICAL ASSUMPTIONS
text: Assumption 1. The loss function φi is convex and Lφ-smooth for i ∈ [n]. Moreover, we assume that it is lower bounded, i.e. infz∈Rc φi(z) > −∞ for i ∈ [n].
We have shown the convexity and smoothness of squared loss and softmax cross-entropy loss in Section 3. The bounded property of φi is required in any algorithm for the well-definedness of (1). Now, in order to use the Taylor series approximation, we need the following assumption on the neural network architecture h:
Assumption 2. We assume that h(·; i) is twice continuously differentiable for all i ∈ [n] (i.e. the second-order partial derivatives of all scalars hj(·; i) are continuous for all j ∈ [c] and i ∈ [n]), and that their Hessian matrices are bounded, that is, there exists a G > 0 such that for all w ∈ Rd, i ∈ [n] and j ∈ [c],
‖Mi,j(w)‖ = ‖Jw (∇whj(w; i))‖ ≤ G, (11)
where Jw denotes the Jacobian1. Remark 1 (Relation to second-order methods). Although our analysis requires an assumption on the Hessian matrices of h(w; i), our algorithms do not use any second order information or try to approximate this information. Our theoretical analysis focused on the approximation of the classifier and the gradient information, therefore is not related to the second order type algorithms. It is currently unclear how to apply second order methods into our problem, however, this is an interesting research question to expand the scope of this work.
1For a continuously differentiable function g(w) : Rd → Rc we define the Jacobian Jw(g(w)) as the matrix (∂ga(w)/∂wb)a∈[c],b∈[d].
Assumption 2 allows us to apply a Taylor approximation of each function hj(·; i) with which we prove the following Lemma that bounds the error in equation (9):
Lemma 1. Suppose that Assumption 2 holds for the classifier h. Then for all i ∈ [n] and 0 ≤ t < T ,
h(w(t+1); i) = h(w(t) − η(t)v(t); i) = h(w(t); i)− η(t)H(t)i v (t) + (t) i , (12)
where H
(t) i = Jw(h(w; i))|w=w(t) ∈ R c×d (13)
is defined as the Jacobian matrix of h(w; i) at w(t) and entries (t)i,j , j ∈ [c], of vector (t) i satisfy
| (t)i,j | ≤ 1
2 (η(t))2‖v(t)‖2G. (14)
In order to approximate (7) combined with (9), that is, to make sure the right hand sides of (7) and (9) are close to one another, we consider the optimization problem (10):
v (t) ∗ = arg min v∈Rd 1 2 1 n n∑ i=1 ‖η(t)H(t)i v − α (t) i ∇zφi(h(w (t); i))‖2.
The optimal value of problem (10) is equal to 0 if there exists a vector v(t)∗ satisfying η(t)H (t) i v (t) ∗ = α (t) i ∇zφi(h(w(t); i)) for every i ∈ [n]. Since the solution v (t) ∗ is in Rd and ∇zφi(h(w(t); i)) is in Rc, this condition is equivalent to a linear system with n · c constraints and d variables. In the overparameterized setting where dimension d is sufficiently large (d n · c) and there are no identical data, there exists almost surely a vector v(t)∗ that interpolates all the training set, see the Appendix for details.
Let us note that an approximation of v(t)∗ serves as the search direction for Framework 1. For this reason, the solution v(t)∗ of problem (10) plays a similar role as a gradient in the search direction of (stochastic) gradient descent method. It is standard to assume a bounded gradient in the machine learning literature (Nemirovski et al., 2009; Shalev-Shwartz et al., 2007; Reddi et al., 2016). Motivated by these facts, we assume the following Assumption 3, which implies the existence of a near-optimal bounded solution of (10):
Assumption 3. We consider an over-parameterized setting where dimension d is sufficiently large enough to interpolate all the data and the tolerance ε. We assume that there exists a bound V > 0 such that for ε > 0 and 0 ≤ t < T as in Framework 1, there exists a vector v̂(t)∗ε with ‖v̂(t)∗ε ‖2 ≤ V so that
1
2
1
n n∑ i=1 ‖η(t)H(t)i v̂ (t) ∗ε − α(t)i ∇zφi(h(w (t); i))‖2 ≤ ε2.
Our Assumption 3 requires a nice dependency on the tolerance ε for the gradient matrices H(t)i and ∇zφi(h(w(t); i)). We note that at the starting point t = 0, these matrices may depend on ε due to the initialization process and the dependence of d on ε. This setting is similar to previous works, e.g. Allen-Zhu et al. (2019).

heading: 5 NEW ALGORITHMS AND CONVERGENCE RESULTS

heading: 5.1 APPROXIMATING THE SOLUTION USING REGULARIZER
text: Since problem (10) is convex and quadratic, we consider the following regularized problem:
min v∈Rd
{ Ψ(v) = 1
2
1
n n∑ i=1 ‖η(t)H(t)i v − α (t) i ∇zφi(h(w (t); i))‖2 + ε 2 2 ‖v‖2
} , (15)
for some small ε > 0 and t ≥ 0. It is widely known that problem (15) is strongly convex, and has a unique minimizer v(t)∗ reg. The global minimizer satisfies∇vΨ(v(t)∗ reg) = 0. We have
∇vΨ(v) = 1
n n∑ i=1 [η(t)H (t) i >H (t) i η (t)v − α(t)i η (t)H (t) i >∇zφi(h(w(t); i))] + ε2 · v
=
( 1
n n∑ i=1 η(t)H (t) i >H (t) i η (t) + ε2I
) v − ( 1
n n∑ i=1 α (t) i η (t)H (t) i >∇zφi(h(w(t); i))
) .
Therefore,
v (t) ∗ reg =
( 1
n n∑ i=1 η(t)H (t) i >H (t) i η (t) + ε2I
)−1( 1
n n∑ i=1 α (t) i η (t)H (t) i >∇zφi(h(w(t); i))
) . (16)
If ε2 is small enough, then v(t)∗ reg is a close approximation of the solution v (t) ∗ for problem (10). Our first algorithm updates Framework 1 based on this approximation.
Algorithm 1 Solve for the exact solution of the regularized problem
Initialization: Choose an initial point w(0) ∈ Rd, tolerance ε > 0; for t = 0, 1, · · · , T − 1 do
Update the search direction v(t) as the solution v(t)∗ reg of problem in (15):
v(t) = v (t) ∗ reg =
( 1
n n∑ i=1 η(t)H (t) i >H (t) i η (t) + ε2I
)−1( 1
n n∑ i=1 α (t) i η (t)H (t) i >∇zφi(h(w(t); i))
)
Update w(t+1) = w(t) − η(t)v(t) end for
The following Lemma shows the relation between the regularized solution v(t)∗ reg and the optimal solution of the original convex problem v̂(t)∗ε . Lemma 2. For given ε > 0, suppose that Assumption 3 holds for bound V > 0. Then, for iteration 0 ≤ t < T , the optimal solution v(t)∗ reg of problem (15) satisfies ‖v(t)∗ reg‖2 ≤ 2 + V and
1
2
1
n n∑ i=1 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2 ≤ (1 + V 2 )ε2. (17)
Based on Lemma 2, we guarantee the global convergence of Algorithm 1 and prove our first theorem. Since it is currently expensive to solve for the exact solution of problem (15), our algorithm serves as a theoretical method to obtain the global convergence for the finite-sum minimization.
Theorem 1. Let w(t) be generated by Algorithm 1 where we use the closed form solution for the search direction. We execute Algorithm 1 for T = βε outer loops for some constant β > 0. We assume Assumption 1 holds. Suppose that Assumption 2 holds for G > 0 and Assumption 3 holds for V > 0. We set the step size equal to η(t) = D √ ε for some D > 0 and choose a learning rate α (t) i = (1 + ε)α (t−1) i = (1 + ε) tα (0) i . Based on β, we define α (0) i = α eβLφ
with α ∈ (0, 13 ). Let F∗ be the global minimizer of F , and h∗i = arg minz∈Rc φi(z), i ∈ [n]. Then
1
T T−1∑ t=0 [F (w(t))− F∗] ≤ eβLφ(1 + ε) 2(1− 3α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ(3ε+ 2) 8α(1− 3α) [ c(4 + (V + 2)GD2)2 + 8 + 4V ] · ε. (18)
We note that β is a constant for the purpose of choosing the number of iterations T . The analysis can be simplified by choosing β = 1 with T = 1ε . Notice that the common convergence criteria for
finding a stationary point for non-convex problems is 1T ∑T t=1 ||∇F (wt)||2 ≤ O(ε). This criteria has been widely used in the existing literature for non-convex optimization problems. Our convergence criteria 1T ∑T t=1[F (wt) − F∗] ≤ O(ε) is slightly different, in order to find a global solution for non-convex problems.
Our proof for Theorem 1 is novel and insightful. It is originally motivated by the Gradient Descent update (7) and the convexity of the loss functions φi. For this reason it may not be a surprise that Algorithm 1 can find an ε-global solution after O ( 1 ε ) iterations. However, computing the exact solution in every iteration might be extremely challenging, especially when the number of samples n is large. Therefore, we present a different approach to this problem in the following section.

heading: 5.2 APPROXIMATION USING GRADIENT DESCENT
text: In this section, we use Gradient Descent (GD) algorithm to solve the strongly convex problem (15). It is well-known that if ψ(x) − µ2 ‖x‖
2 is convex for ∀x ∈ Rc, then ψ(x) is µ-strongly convex (see e.g. Nesterov (2004)). Hence Ψ(·) is ε2-strongly convex. For each iteration t, we use GD to find a search direction v(t) which is sufficiently close to the optimal solution v(t)∗ reg in that
‖v(t) − v(t)∗ reg‖ ≤ ε. (19)
Our Algorithm 2 is described as follows.
Algorithm 2 Solve the regularized problem using Gradient Descent
Initialization: Choose an initial point w(0) ∈ Rd, tolerance ε > 0; for t = 0, 1, · · · , T − 1 do
Use Gradient Descent algorithm to solve Problem (15) and find a solution v(t) that satisfies
‖v(t) − v(t)∗ reg‖ ≤ ε
Update w(t+1) = w(t) − η(t)v(t) end for
Since Algorithm 2 can only approximate a solution within some ε-preciseness, we need a supplemental assumption for the analysis of our next Theorem 2:
Assumption 4. Let H(t)i be the Jacobian matrix defined in Lemma 1. We assume that there exists some constant H > 0 such that, for i ∈ [n], ε > 0, and 0 ≤ t < T as in Algorithm 2,
‖H(t)i ‖ ≤ H√ ε . (20)
Assumption 4 requires a mild condition on the bounded Jacobian of h(w; i), and the upper bound may depend on ε. This flexibility allows us to accommodate a good dependence of ε for the theoretical analysis. We are now ready to present our convergence theorem for Algorithm 2.
Theorem 2. Let w(t) be generated by Algorithm 2 where v(t) satisfies (19). We execute Algorithm 2 for T = βε outer loops for some constant β > 0. We assume Assumption 1 holds. Suppose that Assumption 2 holds for G > 0, Assumption 3 holds for V > 0 and Assumption 4 holds for H > 0. We set the step size equal to η(t) = D √ ε for some D > 0 and choose a learning rate α (t) i = (1 + ε)α (t−1) i = (1 + ε) tα (0) i . Based on β, we define α (0) i = α eβLφ
with α ∈ (0, 14 ). Let F∗ be the global minimizer of F , and h∗i = arg minz∈Rc φi(z), i ∈ [n]. Then
1
T T−1∑ t=0 [F (w(t))− F∗] ≤ eβLφ(1 + ε) 2(1− 4α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ(4ε+ 3) 2α(1− 4α) [ D2H2 + c(2 + (V + ε2 + 2)GD2)2 + 2 + V ] · ε.
Theorem 2 implies Corollary 2 which provides the computational complexity for Algorithm 2. Note that for (Stochastic) Gradient Descent, we derive the complexity in terms of component gradient calculations for the finite-sum problem (1). As an alternative, for Algorithm 2 we compare the number of component gradients in problem (15). Such individual gradient has the following form:
∇vψi(v) = η(t)H(t)i >H (t) i η (t)v − α(t)i η (t)H (t) i >∇zφi(h(w(t); i)).
In machine learning applications, the gradient of f(·; i) is calculated using automatic differentiation (i.e. backpropagation). Since f(·; i) is the composition of the network structure h(·; i) and loss function φi(·), this process also computes the Jacobian matrix H(t)i and the gradient∇zφi(h(w(t); i)) at a specific weight w(t). Since matrix-vector multiplication computation is not expensive, the cost for computing the component gradient of problem (15) is similar to problem (1). Corollary 2. Suppose that the conditions in Theorem 2 hold with η(t) = D √ ε̂√ N
for some D > 0 and 0 < ε̂ ≤ N (that is, we set ε = ε̂/N ), where
N = eβLφ
∑n i=1 ‖h(w (0);i)−h∗i ‖ 2
n(1−4α)αβ + 7eβLφ[D2H2+c(2+(V+3)GD2)2+2+V ] 2α(1−4α) .
Then, the total complexity to guarantee min0≤t≤T−1[F (w(t))−F∗] ≤ 1T ∑T−1 t=0 [F (w
(t))−F∗] ≤ ε̂ is O ( nN
3β ε̂3 (D 2H2 + (ε̂2/N)) log(Nε̂ ) ) .
Remark 2. Corollary 2 shows that O (1/ε̂) outer loop iterations are needed in order to reach an ε̂-global solution, and it proves that each iteration needs the equivalent of O ( n ε̂2 log( 1 ε̂ ) )
gradient computations for computing an approximate solution. In total, Algorithm 2 has total complexity O ( n ε̂3 log( 1 ε̂ ) ) for finding an ε̂-global solution.
For a comparison, Stochastic Gradient Descent uses a total of O( 1ε2 ) gradient computations to find a stationary point satisfying E[‖∇F (ŵ)‖2] ≤ ε for non-convex problems (Ghadimi & Lan, 2013). Gradient Descent has a better complexity in terms of ε, i.e. O(nε ) such that ‖∇F (ŵ)‖
2 ≤ ε (Nesterov, 2004). However, both methods may not be able to reach a global solution of (1). In order to guarantee global convergence for nonconvex settings, one may resort to use Polyak-Lojasiewicz (PL) inequality (Karimi et al., 2016; Gower et al., 2021). This assumption is widely known to be strong, which implies that every stationary point is also a global minimizer.

heading: 6 FURTHER DISCUSSION AND CONCLUSIONS
text: This paper presents an alternative composite formulation for solving the finite-sum optimization problem. Our formulation allows a new way of exploiting the structure of machine learning problems and the convexity of squared loss and softmax cross entropy loss, and leads to a novel algorithmic framework that guarantees global convergence (when the outer loss functions are convex and Lipschitz-smooth). Our analysis is general and can be applied to various different learning architectures, in particular, our analysis and assumptions match practical neural networks; in recent years, there has been a great interest in the structure of deep learning architectures for over-parameterized settings (Arora et al., 2018; Allen-Zhu et al., 2019; Nguyen & Mondelli, 2020). Algorithm 2 demonstrates a gradient method to solve the regularized problem, however, other methods can be applied to our framework (e.g. conjugate gradient descent).
Our theoretical foundation motivates further study, implementation, and optimization of the new algorithmic framework and further investigation of its non-standard bounded style assumptions. Possible research directions include more practical algorithm designs based on our Framework 1, and different related methods to solve the regularized problem and approximate the solution. This potentially leads to a new class of efficient algorithms for machine learning problems. This paper presents a new perspective to the research community.

heading: ETHICS STATEMENT
text: This paper does not contain ethics concerns.

heading: APPENDIX

heading: A TABLE OF NOTATIONS
text: Notation Meaning
F∗ Global minimization function of F in (1)
F∗ = minw∈Rd F (w)
h∗i h ∗ i = arg minz∈Rc φi(z), i ∈ [n]
v (t) ∗ Solution of the convex problem in (10)
minv∈Rd 1 2 1 n ∑n i=1 ‖η(t)H (t) i v − α (t) i ∇zφi(h(w(t); i))‖2
v(t) An approximation of v(t)∗ which is used as the search direction in Framework 1
v̂ (t) ∗ε A vector that satisfies
1 2 1 n ∑n i=1 ‖η(t)H (t) i v − α (t) i ∇zφi(h(w(t); i))‖2 ≤ ε2
for some ε > 0 and ‖v̂(t)∗ε ‖2 ≤ V , for some V > 0. v (t) ∗ reg Solution of the strongly convex problem in (15)
minv∈Rd { 1 2 1 n ∑n i=1 ‖η(t)H (t) i v − α (t) i ∇zφi(h(w(t); i))‖2 + ε 2 2 ‖v‖ 2 }

heading: B USEFUL RESULTS
text: The following lemmas provide key tools for our results. Lemma 3 (Squared loss). Let b ∈ Rc and define φ(z) = 12‖z − b‖
2 for z ∈ Rc. Then φ is convex and Lφ-smooth with Lφ = 1. Lemma 4 (Softmax cross-entropy loss). Let index a ∈ [c] and define
φ(z) = log
[ c∑
k=1
exp(zk − za) ] = log [ c∑
k=1
exp(w>k z)
] ,
for z = (z1, . . . , zc)> ∈ Rc, wherewk = ek−ea with ei representing the i-th unit vector (containing 1 at the i-th position and 0 elsewhere). Then φ is convex and Lφ-smooth with Lφ = 1.
The following lemma is a standard result in (Nesterov, 2004). Lemma 5 ((Nesterov, 2004)). If φ is Lφ-smooth and convex, then for ∀z ∈ Rc,
‖∇φ(z)‖2 ≤ 2Lφ(φ(z)− φ(z∗)), (21) where z∗ = arg minz φ(z).
The following useful derivations could be used later in our theoretical analysis. Since φi is convex, by Definition 2 we have
φi(h(w; i)) ≥ φi(h(w′; i)) + 〈 ∇zφi(z) ∣∣∣ z=h(w′;i) , h(w; i)− h(w′; i) 〉 . (22)
If φi is convex and Lφ-smooth, then by Lemma 5∥∥∥∥∇zφi(z)∣∣∣ z=h(w;i) ∥∥∥∥2 ≤ 2Lφ [φi(h(w; i))− φi(h∗i )] , (23) where h∗i = arg minz∈Rc φi(z).
We compute gradients of f(w; i) in term of φi(h(w; i)).
• Gradient of softmax cross-entropy loss:
∇φi(z) ∣∣ z=h(w;i) =
( ∂φi(z)
∂z1
∣∣∣ z=h(w;i) , . . . , ∂φi(z)
∂zc
∣∣∣ z=h(w;i) )> ,
where for j ∈ [c],
∂φi(z)
∂zj
∣∣∣ z=h(w;i) =  exp ( [h(w;i)]j−[h(w;i)]I(y(i)) ) ∑c k=1 exp ( [h(w;i)]k−[h(w;i)]I(y(i)) ) , j 6= I(y(i)) − ∑ k 6=I(y(i)) exp ( [h(w;i)]k−[h(w;i)]I(y(i))
) ∑c k=1 exp ( [h(w;i)]k−[h(w;i)]I(y(i)) ) , j = I(y(i)) . (24)
• Gradient of squared loss: ∇φi(z) ∣∣ z=h(w;i) = h(w; i)− y(i). (25)

heading: C ADDITIONAL DISCUSSION

heading: C.1 ABOUT ASSUMPTION 2
text: We make a formal assumption for the case h(·; i) is closely approximated by k(·; i). Assumption 5. We assume that for all i ∈ [n] there exists some approximations k(w; i) : Rd → Rc such that
|kj(w; i)− hj(w; i)| ≤ ε, ∀w ∈ Rd, i ∈ [n] and j ∈ [c], (26)
where k(·; i) are twice continuously differentiable (i.e. the second-order partial derivatives of all scalars kj(·; i) are continuous for all i ∈ [n]), and that their Hessian matrices are bounded:
‖Mi,j(w)‖ = ‖Jw (∇wkj(w; i))‖ ≤ G, ∀w ∈ Rd, i ∈ [n] and j ∈ [c]. (27)
Assumption 5 allows us to prove the following Lemma that bound the error in equation (9): Lemma 6. Suppose that Assumption 5 holds for the classifier h. Then for all i ∈ [n] and 0 ≤ t < T , we have:
h(w(t+1); i) = h(w(t) − η(t)v(t); i) = h(w(t); i)− η(t)H(t)i v (t) + (t) i , (28)
where H(t)i is defined to be the Jacobian matrix of the approximation k(w; i) at w (t):
H (t) i := Jwk(w; i)|w=w(t) =  ∂k1(w;i) ∂w1 . . . ∂k1(w;i)∂wd . . . . . . . . .
∂kc(w;i) ∂w1 . . . ∂kc(w;i)∂wd
 ∣∣∣∣∣ w=w(t) ∈ Rc×d. (29)
Additionally we have,
| (t)i,j | ≤ 1
2 (η(t))2‖v(t)‖2G+ 2ε, j ∈ [c]. (30)
Note that these result recover the case when h(·; i) is itself smooth. Hence we analyze our algorithms using the result of Lemma 6, which generalizes the result from Lemma 1.

heading: C.2 ABOUT ASSUMPTION 3
text: In this section, we justify the existence of the search direction in Assumption 3 (almost surely). We argue that there exists a vector v̂(t)∗ε satisfying
1
2
1
n n∑ i=1 ‖η(t)H(t)i v̂ (t) ∗ε − α(t)i ∇zφi(h(w (t); i))‖2 ≤ ε2.
It is sufficient to find a vector v satisfying that
η(t)H (t) i v = α (t) i ∇zφi(h(w (t); i)) for every i ∈ [n].
Since the solution v is in Rd and ∇zφi(h(w(t); i)) is in Rc, this condition is equivalent to a linear system with n ·c constraints and d variables. LetA and b be the following stacked matrix and vector:
A =  H (t) 1 η (t) . . .
H (t) n η(t)
 ∈ Rn·c×d, and b =  α (t) 1 ∇zφ1(h(w(t); i)) . . .
α (t) n ∇zφn(h(w(t); i))  ∈ Rn·c, then the problem reduce to finding the solution of the equation Av = b. In the over-parameterized setting where dimension d is sufficiently large (d n · c), then rank A = n · c almost surely and there exists almost surely a vector v that interpolates all the training set.
To demonstrate this fact easier, we consider a simple neural network where the classifier h(w; i) is formulated as
h(w; i) = W (2)>σ(W (1)>x(i)),
where c = 1, W (1) ∈ Rm×l and W (2) ∈ Rl×1, w = vec({W (1),W (2)}) ∈ Rd is the vectorized weight where d = l(m+ 1) and σ is sigmoid activation function.
H (t) i is defined to be the Jacobian matrix of h(w; i) at w (t):
H (t) i := Jwh(w; i)|w=w(t) = [ ∂h(w;i) ∂w1 . . . ∂h(w;i)∂wd ] ∣∣∣∣∣ w=w(t) ∈ R1×d,
then
A = η(t)  H (t) 1 . . .
H (t) n
 = η(t)  ∂h(w;1) ∂w1 . . . ∂h(w;1)∂wd . . . . . . . . .
∂h(w;n) ∂w1 . . . ∂h(w;n)∂wd  ∈ Rn×d. We want to show that A has full rank, almost surely. We consider the over-parameterized setting where the last layer has at least n neuron (i.e. l = n and the simple version when c = 1. We argue that rank of matrix A is greater than or equal to rank of the submatrix B created by the weights of the last layer W (2) ∈ Rn:
B =  ∂h(w;1) ∂W (2) 1 . . . ∂h(w;1) ∂W (2) n
. . . . . . . . . ∂h(w;n)
∂W (2) 1
. . . ∂h1(w;n) ∂W (2) n  ∈ Rn×n. Note that h(·, i) is a linear function of the last weight layers (in this simple case W (2) ∈ Rn and σ(W (1)>x(i)) ∈ Rn), we can compute the partial derivatives as follows:
∂h(w; i)
∂W (2) = σ(W (1)>x(i)); i ∈ [n].
Hence
B =  σ(W (1)>x(1)) . . .
σ(W (1)>x(n))  ∈ Rn×n. Assuming that there are no identical data, and σ is the sigmoid activation, the set of weights W (1) that make matrix B degenerate has measure zero. Hence B has full rank almost surely, and we have the same conclusion for A. Therefore we are able to prove the almost surely existence of a solution v of the linear equation Av = b for simple two layers network. Using the same argument, this result can be generalized for larger neural networks where the dimension d is sufficiently large (d nc).
C.3 INITIALIZATION EXAMPLE
Our Assumption 3 requires a nice dependency on the tolerance ε for the gradient matrices H(0)i and ∇zφi(h(w(0); i)). We note that at the starting point t = 0, these matrices may depend on ε due to the initialization process and the dependence of d on ε. In order to accommodate the choice of learning rate η(0) = D √ ε in our theorems, in this section we describe a network initialization that
satisfies ‖H(0)i ‖ = Θ ( 1√ ε ) where the gradient norm ‖∇zφi(h(w(0); i))‖ is at most constant order with respect to ε. To simplify the problem, we only consider small-dimension data and networks without activation.
About the target vector: We choose φi to be the softmax cross-entropy loss. By Lemma 7 (see below), we have that the gradient norm is upper bounded by a constant c, where c is the output dimension of the problem and is not dependent on ε. Note that when we stack all gradients for n data points, then the size of new vector is still not dependent on ε.
About the network architecture: For simplicity, we consider the following classification problem where
• The input data is in R2. There are only two data points {x(1), x(2)}. Input data is bounded and non-degenerate (we will clarify this property later).
• The output data is (categorical) in R2: {y(1) = (1, 0), y(2) = (0, 1)}.
We want to have an over-parameterized setting where the dimension of weight vector is at least nc = 4. We consider a simple network with two layers, no biases and no activation functions. Let the number of neurons in the hidden layer bem. The flow of this network is (in) R2 → Rm → R2 (out). First, we consider the case where m = 1.
• The first layer has 2 parameters (w1, w2) and only 1 neuron that outputs z(i) = w1x (i) 1 +
w2x (i) 2 (the subscript is for the coordinate of input data x (i)). • The second layer has 2 parameters (w3, w4). The final output is
h(w, i) = [w3(w1x (i) 1 + w2x (i) 2 ), w4(w1x (i) 1 + w2x (i) 2 )] > ∈ R2,
with w = [w1, w2, w3, w4]> ∈ R4. This network satisfies that the Hessian matrices of h(w; i) are bounded. Let Q and b be the following stacked matrix and vector:
Q =
[ H
(0) 1 H (0) 2
] ∈ R4×4, and b = [ ∇zφ1(h(w(0); 1)) ∇zφ2(h(w(0); 2)) ] ∈ R4,
Then we have the following:
Q = Q(w) =
[ H
(0) 1 H (0) 2
] =  ∇w[w3(w1x(1)1 + w2x (1) 2 )] ∇w[w4(w1x(1)1 + w2x (1) 2 )] ∇w[w3(w1x(2)1 + w2x (2) 2 )]
∇w[w4(w1x(2)1 + w2x (2) 2 )]

=  w3x (1) 1 w3x (1) 2 w1x (1) 1 + w2x (1) 2 0 w4x (1) 1 w4x (1) 2 0 w1x (1) 1 + w2x (1) 2 w3x (2) 1 w3x (2) 2 w1x (2) 1 + w2x (2) 2 0
w4x (2) 1 w4x (2) 2 0 w1x (2) 1 + w2x (2) 2
 .
The determinant of this matrix is a polynomial of the weight w and the input data. Under some mild non-degenerate condition of the input data, we can choose some base point w′ that made this matrix invertible (note that if this condition is not satisfied, we can rescale/add a very small noise to the data - which is the common procedure in machine learning).
Hence the system Qu = b always has a solution. Now we consider the following two initializations:
1. We choose to initialize the starting point at w(0) = 1√ ε w′ and note that Q(w) is a linear function of w and Q(w′) is independent of ε. Then the norm of matrix Q(w(0)) has the same scale with 1√ ε .
2. Instead of choosing m = 1, we consider an over-parameterized network where m = 1ε (recall that m is the number of neurons in the hidden layer). The hidden layer in this case is:
z =  z (i) 1 = w (1) 1,1x (i) 1 + w (1) 2,1x (i) 2 . . .
z (i) m = w (1) 1,mx (i) 1 + w (1) 2,mx (i) 2
.
The output layer is:{ y (i) 1 = z (i) 1 w (2) 1,1 + · · ·+ z (i) m w (2) m,1 = (w (1) 1,1x (i) 1 + w (1) 2,1x (i) 2 )w (2) 1,1 + · · ·+ (w (1) 1,mx (i) 1 + w (1) 2,mx (i) 2 )w (2) m,1
y (i) 2 = z (i) 1 w (2) 1,2 + · · ·+ z (i) m w (2) m,2 = (w (1) 1,1x (i) 1 + w (1) 2,1x (i) 2 )w (2) 1,2 + · · ·+ (w (1) 1,mx (i) 1 + w (1) 2,mx (i) 2 )w (2) m,2
with w = [w(1)1,1, . . . , w (1) 1,m, w (1) 2,1, . . . , w (1) 2,m, w (2) 1,1, w (2) 1,2, . . . , w (2) m,1, w (2) m,2] > ∈ R4m. Hence,
Q(w) =  w (2) 1,1x (1) 1 . . . w (2) m,1x (1) 1 w (2) 1,1x (1) 2 . . . w (2) m,1x (1) 2 z (1) 1 0 . . . z (1) m 0 w (2) 1,2x (1) 1 . . . w (2) m,2x (1) 1 w (2) 1,2x (1) 2 . . . w (2) m,2x (1) 2 0 z (1) 1 . . . 0 z (1) m w (2) 1,1x (2) 1 . . . w (2) m,1x (2) 1 w (2) 1,1x (2) 2 . . . w (2) m,1x (2) 2 z (2) 1 0 . . . z (2) m 0
w (2) 1,2x (2) 1 . . . w (2) m,2x (2) 1 w (2) 1,2x (2) 2 . . . w (2) m,2x (2) 2 0 z (2) 1 . . . 0 z (2) m
 .
Hence, the number of (possibly) non-zero elements in each row is 3m = 3ε . For matrix A of rank r, we have ‖A‖2 ≤ ‖A‖F ≤ √ r‖A‖2. Since the rank of Q(w) is at most 4 (nc = 4, independent of ε), we only need to find the Frobenius norm of Q(w). We have
‖Q(w)‖F = √√√√ 4∑ i=1 4m∑ j=1 |qij |2.
Let qmin and qmax be the element with smallest/largest magnitude of Q(w). Suppose that x(i) 6= (0, 0) and choose w 6= 0 such that z 6= 0, qmin > 0 and independent of ε. Hence, √ 8√ ε |qmin| ≤ ‖Q(w)‖F ≤ √ 12√ ε |qmax|.
Hence, ‖Q(w)‖ = Θ (
1√ ε
) . Therefore this simple network initialization supports the dependence
on ε for our Assumption 3. We note that a similar setting is found in (Allen-Zhu et al., 2019), where the authors initialize the weights using a random Gaussian distribution with a variance depending on the dimension of the problem. In non-convex setting, they prove the convergence of SGD using the assumption that the number of neurons m depends inversely on the tolerance ε.
Lemma 7. For softmax cross-entropy loss, and x = h(w; i) ∈ Rc, for ∀w ∈ Rd and i ∈ [n], we have ∥∥∥∥∇zφi(x)∣∣∣
x=h(w;i) ∥∥∥∥2 ≤ c. (31) Proof. By (24), we have for i = 1, . . . , n,
• For j 6= I(y(i)):( ∂φi(x)
∂xj
∣∣∣ x=h(w;i) )2 = ( exp ( [h(w; i)]j − [h(w; i)]I(y(i)) )∑c k=1 exp ( [h(w; i)]k − [h(w; i)]I(y(i)) ))2
=
( exp ( [h(w; i)]j − [h(w; i)]I(y(i)) ) 1 + ∑ k 6=I(y(i)) exp ( [h(w; i)]k − [h(w; i)]I(y(i)) ))2 ≤ 1. • For j = I(y(i)):(
∂φi(x)
∂xj
∣∣∣ x=h(w;i) )2 = (∑ k 6=I(y(i)) exp ( [h(w; i)]k − [h(w; i)]I(y(i)) )∑c k=1 exp ( [h(w; i)]k − [h(w; i)]I(y(i)) ) )2
=
( ∑ k 6=I(y(i)) exp ( [h(w; i)]k − [h(w; i)]I(y(i)) ) 1 + ∑ k 6=I(y(i)) exp ( [h(w; i)]k − [h(w; i)]I(y(i)) ))2 ≤ 1 Hence, for i = 1, . . . , n,∥∥∥∥∇zφi(x)∣∣∣
x=h(w;i) ∥∥∥∥2 = c∑ j=1 ( ∂φi(x) ∂xj ∣∣∣ x=h(w;i) )2 ≤ c.
This completes the proof.

heading: D PROOFS OF LEMMAS AND COROLLARY 1

heading: PROOF OF LEMMA 1
text: Proof. Since h(·; i) are twice continuously differentiable for all i ∈ [n], we have the following Taylor approximation for each component outputs hj(·; i) where j ∈ [c] and i ∈ [n]:
hj(w (t+1); i) = hj(w (t) − η(t)v(t); i)
= hj(w (t); i)− Jwhj(w; i)|w=w(t)η(t)v(t) +
1 2 (η(t)v(t))>Mi,j(w̃ (t))(η(t)v(t)),
(32)
where Mi,j(w̃(t)) is the Hessian matrices of hj(·; i)at w̃(t) and w̃(t) = αw(t) + (1 − α)w(t+1) for some α ∈ [0, 1]. This leads to our desired statement:
h(w(t+1); i) = h(w(t) − η(t)v(t); i) = h(w(t); i)− η(t)H(t)i v (t) + (t) i ,
where
(t) i,j =
1 2 (η(t)v(t))>Mi,j(w̃ (t))(η(t)v(t)), j ∈ [c],
Hence we get the final bound:
| (t)i,j | ≤ 1
2 ∣∣∣(η(t)v(t))>Mi,j(w̃(t))(η(t)v(t))∣∣∣ ≤ 1
2 (η(t))2‖v(t)‖2 · ‖Mi,j(w̃(t))‖
(11)
≤ 1 2 (η(t))2‖v(t)‖2G, j ∈ [c].

heading: PROOF OF LEMMA 2
text: Proof. From Assumption 3, we know that there exists v̂(t)∗ε so that
1
2
1
n n∑ i=1 ‖η(t)H(t)i v̂ (t) ∗ε − α(t)i ∇zφi(h(w (t); i))‖2 ≤ ε2,
and ‖v̂(t)∗ε ‖2 ≤ V , for some V > 0. Hence,
1
2
1
n n∑ i=1 ‖η(t)H(t)i v̂ (t) ∗ε − α(t)i ∇zφi(h(w (t); i))‖2 + ε 2 2 ‖v̂(t)∗ε ‖2 ≤ ε2 + ε2 2 V = (1 + V 2 )ε2.
Since v(t)∗ reg is the optimal solution of the problem in (15) for 0 ≤ t < T , we have
1
2
1
n n∑ i=1 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2 + ε 2 2 ‖v(t)∗ reg‖2 ≤ (1 + V 2 )ε2.
Therefore, we have (17) and ‖v(t)∗ reg‖2 ≤ 2 + V for 0 ≤ t < T .

heading: PROOF OF LEMMA 3
text: Proof. 1. We want to show that for any α ∈ [0, 1]
φ(αz1 + (1− α)z2) ≤ αφ(z1) + (1− α)φ(z2), ∀z1, z2 ∈ Rc, (33)
in order to have the convexity of φ with respect to z (see (Nesterov, 2004)).
For any α ∈ [0, 1], we have for ∀z1, z2 ∈ Rc,
α‖z1 − b‖2 + (1− α)‖z2 − b‖2 − ‖α(z1 − b) + (1− α)(z2 − b)‖2
= α‖z1 − b‖2 + (1− α)‖z2 − b‖2 − α2‖z1 − b‖2 − (1− α)2‖z2 − b‖2
− 2α(1− α)〈z1 − b, z2 − b〉 ≥ α(1− α)‖z1 − b‖2 + (1− α)α‖z2 − b‖2 − 2α(1− α)‖z1 − b‖ · ‖z2 − b‖ = α(1− α) (‖z1 − b‖ − ‖z2 − b‖)2 ≥ 0,
where the first inequality follows according to Cauchy-Schwarz inequality 〈a, b〉 ≤ ‖a‖·‖b‖. Hence,
1 2 ‖αz1 + (1− α)z2 − b‖2 ≤ α 2 ‖z1 − b‖2 + (1− α) 2 ‖z2 − b‖2.
Therefore, (33) implies the convexity of φ with respect to z.
2. We want to show that ∃Lφ > 0 such that
‖∇φ(z1)−∇φ(z2)‖ ≤ Lφ‖z1 − z2‖, ∀z1, z2 ∈ Rc. (34)
Notice that∇φ(z) = z − b, then clearly ∀z1, z2 ∈ Rc,
‖∇φ(z1)−∇φ(z2)‖ = ‖z1 − z2‖.
Therefore, (34) implies the Lφ-smoothness of φ with respect to z with Lφ = 1.

heading: PROOF OF LEMMA 4
text: Proof. 1. For ∀z1, z2 ∈ Rc and 1 ≤ k ≤ c, denote uk,1 = exp(w>k z1) and uk,2 = exp(w>k z2) and using Holder inequality
c∑ k=1 ak · bk ≤
( c∑
k=1
|ak|p ) 1 p (
c∑ k=1
|bk|q ) 1 q , where 1
p +
1 q = 1, (35)
we have
φ(αz1 + (1− α)z2) = log
[ c∑
k=1
exp(w>k (αz1 + (1− α)z2))
] = log [ c∑
k=1
uαk,1 · u (1−α) k,2 ] (35)
≤ log ( c∑ k=1 u α· 1α k,1 )α( c∑ k=1 u (1−α)· 1 (1−α) k,2 )1−α = α log [ c∑
k=1
exp(w>k z1)
] + (1− α) log [ c∑
k=1
exp(w>k z2) ] = αφ(z1) + (1− α)φ(z2),
where the first inequality since log(x) is an increasing function for ∀x > 0 and exp(v) > 0 for ∀v ∈ R. Therefore, (33) implies the convexity of φ with respect to z. 2. Note that ‖∇2φ(z)‖ ≤ Lφ if and only if φ(z) is Lφ-smooth (see (Nesterov, 2004)). First, we compute gradient of φ(z):
• For i 6= a:
∂φ(z)
∂zi = exp(zi − za)∑c k=1 exp(zk − za) .
• For i = a:
∂φ(z) ∂zi = − ∑ k 6=a exp(zk − za)∑c k=1 exp(zk − za) = − ∑c k=1 exp(zk − za) + 1∑c k=1 exp(zk − za)
= −1 + 1∑c k=1 exp(zk − za) = −1 + exp(zi − za)∑c k=1 exp(zk − za) .
We then calculate ∂ 2φ(z) ∂zj∂zi = ∂∂zj ( ∂φ(z) ∂zi ) • For i = j:
∂2φ(z) ∂zj∂zi =
exp(zi − za)[ ∑c k=1 exp(zk − za)]− exp(zi − za) exp(zi − za)
[ ∑c k=1 exp(zk − za)]2
= exp(zi − za)[
∑c k=1 exp(zk − za)− exp(zi − za)]
[ ∑c k=1 exp(zk − za)]2
.
• For i 6= j:
∂2φ(z) ∂zj∂zi = − exp(zj − za) exp(zi − za) [ ∑c k=1 exp(zk − za)]2 .
Denote that yi = exp(zi − za) ≥ 0, i ∈ [c], we have:
• For i = j: ∣∣∣∣∂2φ(z)∂zj∂zi ∣∣∣∣ = ∣∣∣∣yi(∑ck=1 yk − yi)(∑ck=1 yk)2 ∣∣∣∣ . • For i 6= j: ∣∣∣∣∂2φ(z)∂zj∂zi ∣∣∣∣ = |yiyj |(∑ck=1 yk)2 .
Recall that for matrix A = (aij) ∈ Rc×c: ‖A‖2 ≤ ‖A‖2F = ∑c i=1 ∑c j=1 |aij |2. We have:
c∑ j=1 ∣∣∣∣∂2φ(z)∂zj∂zi ∣∣∣∣2 ≤ 1(∑ck=1 yk)4 y2i ( c∑ k=1 yk − yi)2 + ∑ j 6=i (yiyj) 2  = 1
( ∑c k=1 yk) 4 y2i ( c∑ k=1 yk) 2 − 2y2i c∑ k=1 yk.yi + y 4 i + ∑ j 6=i (yiyj) 2  = 1
( ∑c k=1 yk) 4
[ y2i ( c∑ k=1 yk) 2 − 2y3i c∑ k=1 yk + y 2 i c∑ k=1 y2k ] Therefore,
‖∇2φ(z)‖2 ≤ c∑ i=1 c∑ j=1 ∣∣∣∣∂2φ(z)∂zj∂zi ∣∣∣∣2
≤ 1 ( ∑c k=1 yk) 4
[ (
c∑ i=1 y2i )( c∑ k=1 yk) 2 − 2( c∑ i=1 y3i )( c∑ k=1 yk) + ( c∑ i=1 y2i )( c∑ k=1 y2k)
]
≤ ( ∑c i=1 y 2 i )( ∑c k=1 yk) 2
( ∑c k=1 yk) 4 ≤
( ∑c k=1 yk) 4
( ∑c k=1 yk) 4 = 1,
where the last inequality holds since
( c∑ i=1 y2i )( c∑ k=1 y2k) ≤ ( c∑ i=1 y3i )( c∑ k=1 yk)⇔ ( c∑ k=1 y2k) ≤ √√√√( c∑ i=1 y3i )( c∑ k=1 yk),
which follows by the application of Holder inequality (35) with p = 2, q = 2, ak = y 3/2 k , and bk = y 1/2 k (Note that yk ≥ 0, k ∈ [c]). Hence, ‖∇2φ(z)‖ ≤ Lφ with Lφ = 1 which is equivalent to Lφ-smoothness of φ.

heading: PROOF OF LEMMA 6
text: Proof. Since k(·; i) are twice continuously differentiable for all i ∈ [n], we have the following Taylor approximation for each component outputs kj(·; i) where j ∈ [c] and i ∈ [n]:
kj(w (t+1); i) = kj(w (t) − η(t)v(t); i)
= kj(w (t); i)− Jwkj(w; i)|w=w(t)η(t)v(t) +
1 2 (η(t)v(t))>Mi,j(w̃ (t))(η(t)v(t)),
(36)
where Mi,j(w̃(t)) is the Hessian matrices of kj(·; i)at w̃(t) and w̃(t) = αw(t) + (1 − α)w(t+1) for some α ∈ [0, 1]. Shifting this back to the original function hj(·; i) we have:
hj(w (t+1); i) = kj(w (t+1); i) + (hj(w (t+1); i)− kj(w(t+1); i))
(36) = kj(w (t); i)− Jwkj(w; i)|w=w(t)η(t)v(t) + 1
2 (η(t)v(t))>Mi,j(w̃ (t))(η(t)v(t))
+ (hj(w (t+1); i)− kj(w(t+1); i)),
= hj(w (t); i)− Jwkj(w; i)|w=w(t)η(t)v(t) +
1 2 (η(t)v(t))>Mi,j(w̃ (t))(η(t)v(t))
+ (hj(w (t+1); i)− kj(w(t+1); i)) + (kj(w(t); i)− hj(w(t); i)),
which leads to our desired statement:
h(w(t+1); i) = h(w(t) − η(t)v(t); i) = h(w(t); i)− η(t)H(t)i v (t) + (t) i ,
where
(t) i,j =
1 2 (η(t)v(t))>Mi,j(w̃ (t))(η(t)v(t))
+ (hj(w (t+1); i)− kj(w(t+1); i)) + (kj(w(t); i)− hj(w(t); i)), j ∈ [c],
Hence we get the final bound:
| (t)i,j | ≤ 1
2 ∣∣∣(η(t)v(t))>Mi,j(w̃(t))(η(t)v(t))∣∣∣ + |hj(w(t+1); i)− kj(w(t+1); i)|+ |kj(w(t); i)− hj(w(t); i)|
(26)
≤ 1 2 ∣∣∣(η(t)v(t))>Mi,j(w̃(t))(η(t)v(t))∣∣∣+ 2ε, ≤ 1
2 (η(t))2‖v(t)‖2 · ‖Mi,j(w̃(t))‖+ 2ε
(11)
≤ 1 2 (η(t))2‖v(t)‖2G+ 2ε, j ∈ [c].

heading: PROOF OF COROLLARY 1
text: Proof. The proof of this corollary follows directly by the applications of Lemmas 3 and 4.

heading: E TECHNICAL PROOFS FOR THEOREM 1
text: Lemma 8. Suppose that Assumption 2 holds for G > 0 and Assumption 3 holds for V > 0, and v(t) = v (t) ∗ reg. Consider η(t) = D √ ε for some D > 0 and ε > 0. For i ∈ [n] and 0 ≤ t < T , we have
‖ (t)i ‖ 2 ≤ 1
4 c(4 + (V + 2)GD2)2ε2. (37)
Proof. From (14), for i ∈ [n], j ∈ [c], and for 0 ≤ t < T , by Lemma 1 and Lemma 6 we have
| (t)i,j | ≤ 1 2 (η(t))2‖v(t)‖2G+ 2ε ≤ 1 2 (V + 2)GD2ε+ 2ε = 1 2 ε(4 + (V + 2)GD2),
where the last inequality follows by the fact ‖v(t)‖2 = ‖v(t)∗ reg‖2 ≤ 2 + V of Lemma 2 and η(t) = D √ ε. Hence,
‖ (t)i ‖ 2 = c∑ j=1 | (t)i,j | 2 ≤ 1 4 c(4 + (V + 2)GD2)2ε2.
Lemma 9. Let w(t) be generated by Algorithm 1 where we use the closed form solution for the search direction. We execute Algorithm 1 for T = βε outer loops for some constant β > 0. We assume Assumption 1 holds. Suppose that Assumption 2 holds for G > 0 and Assumption 3 holds for V > 0. We set the step size equal to η(t) = D √ ε for some D > 0 and choose a learning rate α (t) i ≤ αLφ , for some α ∈ (0, 1 3 ). For i ∈ [n] and 0 ≤ t < T , we have
‖h(w(t+1); i)− h∗i ‖2 ≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 3α)α (t) i [φi(h(w (t); i))− φi(h∗i )]
+ (3ε+ 2)
4 c(4 + (V + 2)GD2)2 · ε
+ 3ε+ 2
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2 (38)
Proof. Note that we have the optimal solution v(t)∗ reg for the optimization problem (15) for 0 ≤ t < T . From (12), we have, for i ∈ [n],
h(w(t+1); i) = h(w(t) − η(t)v(t)∗ reg; i)
= h(w(t); i)− η(t)H(t)i v (t) ∗ reg + (t) i
= h(w(t); i)− α(t)i ∇zφi(h(w (t); i)) + (t) i − [η (t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))].
Hence, we have
‖h(w(t+1); i)− h∗i ‖2
= ‖h(w(t); i)− h∗i − α (t) i ∇zφi(h(w (t); i)) + (t) i − [η (t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))]‖2
= ‖h(w(t); i)− h∗i ‖2 + (α (t) i ) 2‖∇zφi(h(w(t); i))‖2
+ ‖ (t)i ‖ 2 + ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
− 2 · 〈h(w(t); i)− h∗i , α (t) i ∇zφi(h(w (t); i))〉
+ 2 · 〈h(w(t); i)− h∗i , (t) i 〉
− 2 · 〈h(w(t); i)− h∗i , η(t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉
− 2 · 〈α(t)i ∇zφi(h(w (t); i)), (t) i 〉
+ 2 · 〈α(t)i ∇zφi(h(w (t); i)), η(t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉
− 2 · 〈 (t)i , η (t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉,
where we expand the square term. Now applying Young’s inequalities: 2|〈u, v〉| ≤ ‖u‖ 2
ε/2 +(ε/2)‖v‖ 2
for ε > 0 and 2|〈u, v〉| ≤ ‖u‖2 + ‖v‖2 we have:
‖h(w(t+1); i)− h∗i ‖2
= ‖h(w(t); i)− h∗i ‖2 + (α (t) i ) 2‖∇zφi(h(w(t); i))‖2
+ ‖ (t)i ‖ 2 + ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
− 2α(t)i 〈h(w (t); i)− h∗i ,∇zφi(h(w(t); i))〉
+ ε
2 ‖h(w(t); i)− h∗i ‖2 +
2 ε ‖ (t)i ‖ 2
+ ε
2 ‖h(w(t); i)− h∗i ‖2 +
2 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
+ 2(α (t) i ) 2‖∇zφi(h(w(t); i))‖2 + 2‖ (t)i ‖ 2
+ 2‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
(22)
≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 + 3(α (t) i ) 2‖∇zφi(h(w(t); i))‖2
+ ( 3 + 2
ε
) ‖ (t)i ‖ 2 + ( 3 + 2
ε
) ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
− 2α(t)i [φi(h(w (t); i))− φi(h∗i )]
Note that from (23) we get that ‖∇zφi(h(w(t); i))‖2 ≤ 2Lφ[φi(h(w(t); i))−φi(h∗i )]. Applying this and using the fact that α(t)i ≤ αLφ , for some α ∈ (0, 1 3 ), we are able to derive:
‖h(w(t+1); i)− h∗i ‖2
≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 3α)α (t) i [φi(h(w (t); i))− φi(h∗i )]
+ 3ε+ 2
ε ‖ (t)i ‖
2 + 3ε+ 2
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 3α)α (t) i [φi(h(w (t); i))− φi(h∗i )]
+ 3ε+ 2
ε
1 4 c(4 + (V + 2)GD2)2ε2 + 3ε+ 2 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
where the last inequality follows by Lemma 8.
Lemma 10. Let w(t) be generated by Algorithm 1 where we use the closed form solution for the search direction. We execute Algorithm 1 for T = βε outer loops for some constant β > 0. We assume Assumption 1 holds. Suppose that Assumption 2 holds for G > 0 and Assumption 3 holds for V > 0. We set the step size equal to η(t) = D √ ε for some D > 0 and choose a learning rate α (t) i = (1 + ε)α (t−1) i = (1 + ε) tα (0) i . Based on β, we define α (0) i = α eβLφ with α ∈ (0, 13 ).
We have
1
T T−1∑ t=0 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )] ≤ eβLφ(1 + ε) 2(1− 3α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ
8α(1− 3α) (3ε+ 2)
[ c(4 + (V + 2)GD2)2 + 8 + 4V ] · ε.
(39)
Proof. Rearranging the terms in Lemma 9, we have
φi(h(w (t); i))− φi(h∗i ) ≤
1
2(1− 3α)
( (1 + ε)
α (t) i
‖h(w(t); i)− h∗i ‖2 − 1
α (t) i
‖h(w(t+1); i)− h∗i ‖2 )
+ 1 8(1− 3α) · 1 α (t) i · ε(3ε+ 2)c(4 + (V + 2)GD2)2 + 1
2(1− 3α) · 1 α (t) i · 3ε+ 2 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
≤ 1 2(1− 3α)
( (1 + ε)
α (t) i
‖h(w(t); i)− h∗i ‖2 − (1 + ε)
α (t+1) i
‖h(w(t+1); i)− h∗i ‖2 )
+ eβLφ
8α(1− 3α) · ε(3ε+ 2)c(4 + (V + 2)GD2)2
+ eβLφ 2α(1− 3α) · 3ε+ 2 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2.
The last inequality follows because the learning rate satisfies α(0)i = α eβLφ ≤ αLφ and for t = 1, . . . , T = βε for some β > 0
α (t) i = (1 + ε)α (t−1) i = (1 + ε) tα (0) i ≤ (1 + ε) Tα (0) i = (1 + ε)
β/ε α eβLφ ≤ α Lφ ,
since (1 + x)1/x ≤ e, x > 0. Moreover, we have 1 α
(t) i
≤ 1 α
(0) i
= eβLφ α , t = 0, . . . , T − 1.
Taking the average sum from t = 0, . . . , T − 1, we have
1
T T−1∑ t=0 [φi(h(w (t); i))− φi(h∗i )] ≤
1 2(1− 3α)T · (1 + ε) α (0) i ‖h(w(0); i)− h∗i ‖2
+ eβLφ
8α(1− 3α) · ε(3ε+ 2)c(4 + (V + 2)GD2)2
+ eβLφ 2α(1− 3α) · 3ε+ 2 ε 1 T T−1∑ t=0 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
= eβLφ(1 + ε)
2(1− 3α)αβ ε · ‖h(w(0); i)− h∗i ‖2
+ eβLφ
8α(1− 3α) · ε(3ε+ 2)c(4 + (V + 2)GD2)2
+ eβLφ 2α(1− 3α) · 3ε+ 2 ε 1 T T−1∑ t=0 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2.
Taking the average sum from i = 1, . . . , n, we have
1
T T−1∑ t=0 1 n n∑ i=1 [φi(h(w (t); i))− φi(h∗i )]
≤ e βLφ(1 + ε) 2(1− 3α)αβ ε · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2
+ eβLφ
8α(1− 3α) · ε(3ε+ 2)c(4 + (V + 2)GD2)2
+ eβLφ 2α(1− 3α) · 3ε+ 2 ε 1 T T−1∑ t=0 1 n n∑ i=1 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
(17) ≤ e βLφ(1 + ε) 2(1− 3α)αβ ε · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2
+ eβLφ
8α(1− 3α) · ε(3ε+ 2)c(4 + (V + 2)GD2)2
+ eβLφ 2α(1− 3α) · 3ε+ 2 ε (2 + V )ε2. (40)
Note that
1
T T−1∑ t=0 1 n n∑ i=1 [φi(h(w (t); i))− φi(h∗i )] = 1 T T−1∑ t=0 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]. (41)
Therefore, applying (41) to (40), we have
1
T T−1∑ t=0 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]
≤ e βLφ(1 + ε) 2(1− 3α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ
8α(1− 3α) (3ε+ 2)
[ c(4 + (V + 2)GD2)2 + 8 + 4V ] · ε.
which is our desired result.

heading: PROOF OF THEOREM 1
text: Proof. We have
F∗ = min w∈Rd F (w) = min w∈Rd
( 1
n n∑ i=1 fi(w)
) = 1
n min w∈Rd ( n∑ i=1 fi(w) )
≥ 1 n n∑ i=1 min w∈Rd (fi(w)) = 1 n n∑ i=1 f∗i ≥ 1 n n∑ i=1 φi(h ∗ i ). (42)
Hence F∗ − 1n ∑n i=1 φi(h ∗ i ) ≥ 0. Therefore
1
T T∑ t=1 [F (w(t))− F∗] = 1 T T∑ t=1
( 1
n n∑ i=1 [f(w(t); i)− φi(h∗i )]−
[ F∗ − 1
n n∑ i=1 φi(h ∗ i )
])
≤ 1 T T∑ t=1 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]
(39) ≤ e βLφ(1 + ε) 2(1− 3α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ(3ε+ 2) 8α(1− 3α) [ c(4 + (V + 2)GD2)2 + 8 + 4V ] · ε.

heading: F TECHNICAL PROOFS FOR THEOREM 2
text: Lemma 11. For 0 ≤ t < T , suppose that Assumption 3 holds for V ≥ 0 and v(t) satisfies (19). Then
‖v(t)‖2 ≤ 2(ε2 + V + 2).
Proof. From ‖v(t) − v(t)∗ reg‖ ≤ ε. Using ‖a‖2 ≤ 2‖a− b‖2 + 2‖b‖2, we have
‖v(t)‖2 ≤ 2‖v(t) − v(t)∗ reg‖2 + 2‖v(t)∗ reg‖2 (19) ≤ 2ε2 + 4 + 2V.
where the last inequality follows since ‖v(t)∗ reg‖2 ≤ 2 + V for some V > 0 in Lemma 2.
Lemma 12. Suppose that Assumption 2 holds for G > 0 and Assumption 3 holds for V > 0. Consider η(t) = D √ ε for some D > 0 and ε > 0. For i ∈ [n] and 0 ≤ t < T , we have
‖ (t)i ‖ 2 ≤ c(2 + (V + ε2 + 2)GD2)2ε2. (43)
Proof. From (14), for i ∈ [n], j ∈ [c], and for 0 ≤ t < T , by Lemma 1 and Lemma 6 we have
| (t)i,j | ≤ 1 2 (η(t))2‖v(t)‖2G+ 2ε ≤ 1 2 2(ε2 + V + 2)GD2ε+ 2ε = ε(2 + (V + ε2 + 2)GD2),
where the last inequality follows by the application of Lemma 11 and η(t) = D √ ε. Hence,
‖ (t)i ‖ 2 = c∑ j=1 | (t)i,j | 2 ≤ c(2 + (V + ε2 + 2)GD2)2ε2.
Lemma 13. Let w(t) be generated by Algorithm 2 where v(t) satisfies (19). We execute Algorithm 2 for T = βε outer loops for some constant β > 0. We assume Assumption 1 holds. Suppose that Assumption 2 holds for G > 0, Assumption 3 holds for V > 0 and Assumption 4 holds for H > 0. We set the step size equal to η(t) = D
√ ε for some D > 0 and choose a learning rate α(t)i ≤ αLφ ,
for some α ∈ (0, 14 ). For i ∈ [n] and 0 ≤ t < T , we have
‖h(w(t+1); i)− h∗i ‖2 ≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 4α)α (t) i [φi(h(w (t); i))− φi(h∗i )] + ε(4ε+ 3) [ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] + 4ε+ 3
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2 (44)
Proof. Note that v(t) is obtained from the optimization problem (15) for 0 ≤ t < T . From (9), we have, for i ∈ [n],
h(w(t+1); i) = h(w(t) − η(t)v(t); i)
= h(w(t); i)− η(t)H(t)i v (t) + (t) i
= h(w(t); i)− η(t)H(t)i (v (t) − v(t)∗ reg)− α(t)i ∇zφi(h(w (t); i)) + (t) i
− [η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))].
Hence, we have
‖h(w(t+1); i)− h∗i ‖2
= ‖h(w(t); i)− h∗i − η(t)H (t) i (v (t) − v(t)∗ reg)− α(t)i ∇zφi(h(w (t); i))
+ (t) i − [η (t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))]‖2
= ‖h(w(t); i)− h∗i ‖2 + ‖η(t)H (t) i (v (t) − v(t)∗ reg)‖2 + (α(t)i ) 2‖∇zφi(h(w(t); i))‖2
+ ‖ (t)i ‖ 2 + ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
− 2 · 〈h(w(t); i)− h∗i , η(t)H (t) i (v (t) − v(t)∗ reg)〉
− 2 · 〈h(w(t); i)− h∗i , α (t) i ∇zφi(h(w (t); i))〉
+ 2 · 〈h(w(t); i)− h∗i , (t) i 〉
− 2 · 〈h(w(t); i)− h∗i , η(t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉
+ 2 · 〈η(t)H(t)i (v (t) − v(t)∗ reg), α(t)i ∇zφi(h(w (t); i))〉
− 2 · 〈η(t)H(t)i (v (t) − v(t)∗ reg), (t)i 〉
+ 2 · 〈η(t)H(t)i (v (t) − v(t)∗ reg), η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉
− 2 · 〈α(t)i ∇zφi(h(w (t); i)), (t) i 〉
+ 2 · 〈α(t)i ∇zφi(h(w (t); i)), η(t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉
− 2 · 〈 (t)i , η (t)H (t) i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))〉,
where we expand the square term. Now applying Young’s inequalities: 2|〈u, v〉| ≤ ‖u‖ 2
ε/3 +(ε/3)‖v‖ 2
for ε > 0 and 2|〈u, v〉| ≤ ‖u‖2 + ‖v‖2 we have:
‖h(w(t+1); i)− h∗i ‖2
= ‖h(w(t); i)− h∗i ‖2 + ‖η(t)H (t) i (v (t) − v(t)∗ reg)‖2 + (α(t)i ) 2‖∇zφi(h(w(t); i))‖2
+ ‖ (t)i ‖ 2 + ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
+ ε
3 ‖h(w(t); i)− h∗i ‖2 +
3 ε ‖η(t)H(t)i (v (t) − v(t)∗ reg)‖2
− 2α(t)i 〈h(w (t); i)− h∗i ,∇zφi(h(w(t); i))〉
+ ε
3 ‖h(w(t); i)− h∗i ‖2 +
3 ε ‖ (t)i ‖ 2
+ ε
3 ‖h(w(t); i)− h∗i ‖2 +
3 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
+ 3(η(t))2‖H(t)i (v (t) − v(t)∗ reg)‖2 + 3(α(t)i ) 2‖∇zφi(h(w(t); i))‖2 + 3‖ (t)i ‖ 2
+ 3‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
(22)
≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 + 4(α (t) i ) 2‖∇zφi(h(w(t); i))‖2
+ ( 4 + 3
ε
) ‖η(t)H(t)i (v (t) − v(t)∗ reg)‖2 + ( 4 + 3
ε
) ‖ (t)i ‖ 2
+ ( 4 + 3
ε
) ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
− 2α(t)i [φi(h(w (t); i))− φi(h∗i )]
Note that from (23) we get that ‖∇zφi(h(w(t); i))‖2 ≤ 2Lφ[φi(h(w(t); i))−φi(h∗i )]. Applying this and using the fact that α(t)i ≤ αLφ , for some α ∈ (0, 1 4 ), we are able to derive:
‖h(w(t+1); i)− h∗i ‖2
≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 4α)α (t) i [φi(h(w (t); i))− φi(h∗i )]
+ 4ε+ 3
ε ‖η(t)H(t)i (v
(t) − v(t)∗ reg)‖2 + 4ε+ 3
ε ‖ (t)i ‖ 2
+ 4ε+ 3
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
(a) ≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 4α)α (t) i [φi(h(w (t); i))− φi(h∗i )]
+ 4ε+ 3
ε D2ε
H2
ε ‖v(t) − v(t)∗ reg‖2 +
4ε+ 3
ε ‖ (t)i ‖ 2
+ 4ε+ 3
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
(b) ≤ (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 4α)α (t) i [φi(h(w (t); i))− φi(h∗i )]
+ 4ε+ 3 ε D2H2 · ε2 + 4ε+ 3 ε · c(2 + (V + ε2 + 2)GD2)2ε2 + 4ε+ 3
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
= (1 + ε)‖h(w(t); i)− h∗i ‖2 − 2(1− 4α)α (t) i [φi(h(w (t); i))− φi(h∗i )] + ε(4ε+ 3) [ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] + 4ε+ 3
ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
where (a) follows by using matrix vector inequality ‖Hv‖ ≤ ‖H‖‖v‖, where H ∈ Rc×d and v ∈ Rd and Assumption 4 in (20) and η(t) = D √ ε for some D > 0 and ε > 0; (b) follows by the fact that ‖v(t) − v(t)∗ reg‖2 ≤ ε2 in (19) and Lemma 12.
Lemma 14. Let w(t) be generated by Algorithm 2 where v(t) satisfies (19). We execute Algorithm 2 for T = βε outer loops for some constant β > 0. We assume Assumption 1 holds. Suppose that Assumption 2 holds for G > 0, Assumption 3 holds for V > 0 and Assumption 4 holds for H > 0.We set the step size equal to η(t) = D √ ε for some D > 0 and choose a learning rate α (t) i = (1 + ε)α (t−1) i = (1 + ε) tα (0) i . Based on β, we define α (0) i = α eβLφ with α ∈ (0, 14 ).
We have
1
T T−1∑ t=0 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]
≤ e βLφ(1 + ε) 2(1− 4α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ(4ε+ 3) 2α(1− 4α) [ D2H2 + c(2 + (V + ε2 + 2)GD2)2 + 2 + V ] · ε. (45)
Proof. Rearranging the terms in Lemma 13, we have
φi(h(w (t); i))− φi(h∗i ) ≤
1
2(1− 4α)
( (1 + ε)
α (t) i
‖h(w(t); i)− h∗i ‖2 − 1
α (t) i
‖h(w(t+1); i)− h∗i ‖2 )
+ 1 2(1− 4α) · 1 α (t) i
· ε(4ε+ 3) [ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] + 1
2(1− 4α) · 1 α (t) i · 4ε+ 3 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
≤ 1 2(1− 4α)
( (1 + ε)
α (t) i
‖h(w(t); i)− h∗i ‖2 − (1 + ε)
α (t+1) i
‖h(w(t+1); i)− h∗i ‖2 )
+ eβLφ
2α(1− 4α) · ε(4ε+ 3)
[ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] +
eβLφ 2α(1− 4α) · 4ε+ 3 ε ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2.
(46)
The last inequality follows because the learning rate satisfies α(0)i = α eβLφ ≤ αLφ and for t = 1, . . . , T = βε for some β > 0
α (t) i = (1 + ε)α (t−1) i = (1 + ε) tα (0) i ≤ (1 + ε) Tα (0) i = (1 + ε)
β/ε α eβLφ ≤ α Lφ ,
since (1 + x)1/x ≤ e, x > 0. Moreover, we have 1 α
(t) i
≤ 1 α
(0) i
= eβLφ α , t = 0, . . . , T − 1.
Taking the average sum from t = 0, . . . , T − 1, we have
1
T T−1∑ t=0 [φi(h(w (t); i))− φi(h∗i )] ≤
1 2(1− 4α)T · (1 + ε) α (0) i ‖h(w(0); i)− h∗i ‖2
+ eβLφ
2α(1− 4α) · ε(4ε+ 3)
[ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] +
eβLφ 2α(1− 4α) · 4ε+ 3 ε 1 T T−1∑ t=0 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
= eβLφ(1 + ε)
2(1− 4α)αβ ε · ‖h(w(0); i)− h∗i ‖2
+ eβLφ
2α(1− 4α) · ε(4ε+ 3)
[ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] +
eβLφ 2α(1− 4α) · 4ε+ 3 ε 1 T T−1∑ t=0 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2.
Taking the average sum from i = 1, . . . , n, we have
1
T T−1∑ t=0 1 n n∑ i=1 [φi(h(w (t); i))− φi(h∗i )]
≤ e βLφ(1 + ε) 2(1− 4α)αβ ε · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2
+ eβLφ
2α(1− 4α) · ε(4ε+ 3)
[ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] +
eβLφ 2α(1− 4α) · 4ε+ 3 ε 1 T T−1∑ t=0 1 n n∑ i=1 ‖η(t)H(t)i v (t) ∗ reg − α(t)i ∇zφi(h(w (t); i))‖2
(17) ≤ e βLφ(1 + ε) 2(1− 4α)αβ ε · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2
+ eβLφ
2α(1− 4α) · ε(4ε+ 3)
[ D2H2 + c(2 + (V + ε2 + 2)GD2)2 ] +
eβLφ 2α(1− 4α) · 4ε+ 3 ε (2 + V )ε2. (47)
Note that
1
T T−1∑ t=0 1 n n∑ i=1 [φi(h(w (t); i))− φi(h∗i )] = 1 T T−1∑ t=0 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]. (48)
Therefore, applying (48) to (47), we have
1
T T−1∑ t=0 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]
≤ e βLφ(1 + ε) 2(1− 4α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ
2α(1− 4α) (4ε+ 3)
[ D2H2 + c(2 + (V + ε2 + 2)GD2)2 + 2 + V ] · ε.

heading: PROOF OF THEOREM 2
text: Proof. From (42) we have F∗ − 1n ∑n i=1 φi(h ∗ i ) ≥ 0. This leads to
1
T T∑ t=1 [F (w(t))− F∗] = 1 T T∑ t=1
( 1
n n∑ i=1 [f(w(t); i)− φi(h∗i )]−
[ F∗ − 1
n n∑ i=1 φi(h ∗ i ) ]) (42)
≤ 1 T T∑ t=1 1 n n∑ i=1 [f(w(t); i)− φi(h∗i )]
(45) ≤ e βLφ(1 + ε) 2(1− 4α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ eβLφ(4ε+ 3) 2α(1− 4α) [ D2H2 + c(2 + (V + ε2 + 2)GD2)2 + 2 + V ] · ε.
(49)

heading: PROOF OF COROLLARY 2
text: Proof. For each iteration 0 ≤ t < T , we need to find v(t) satisfying the following criteria:
‖v(t) − v(t)∗ reg‖2 ≤ ε2,
for some ε > 0. Using Gradient Descent we need O(nLµ log( 1 ε2 )) = O(2n L µ log( 1 ε )) number of gradient evaluations (Nesterov, 2004), where L and µ = ε2 are the smooth and strongly convex constants, respectively, of Ψ. Let
ψi(v) = 1
2 ‖η(t)H(t)i v − α (t) i ∇zφi(h(w (t); i))‖2, i ∈ [n]. (50)
Then, for any v ∈ Rc
∇vψi(v) = η(t)H(t)i >[η(t)H (t) i v − α (t) i ∇zφi(h(w (t); i))], i ∈ [n]. (51)
Consider η(t) = D √ ε for some D > 0 and ε > 0, we have for i ∈ [n] and 0 ≤ t < T
‖∇2vψi(v)‖ = (η(t))2‖H (t) i >H (t) i ‖ ≤ (η (t))2‖H(t)i ‖ · ‖H (t) i ‖
(20)
≤ D2H2.
Hence, ‖∇2vΦ(v)‖ ≤ ‖∇2vψi(v)‖ + ε2 for any v ∈ Rc which implies that L = D2H2 + ε2 (Nesterov (2004)) and Lµ = D2H2+ε2 ε2 . Therefore, the complexity to find v (t) for each iteration t is O(2nD 2H2+ε2
ε2 log( 1 ε )).
Let us choose 0 < ε ≤ 1. From (49), we have
1
T T−1∑ t=0 [F (w(t))− F∗] ≤ eβLφ (1− 4α)αβ · 1 n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 · ε
+ 7eβLφ 2α(1− 4α) [ D2H2 + c(2 + (V + 3)GD2)2 + 2 + V ] · ε = Nε,
where
N = eβLφ (1−4α)αβ 1
n n∑ i=1 ‖h(w(0); i)− h∗i ‖2 + 7eβLφ 2α(1−4α) [ D2H2 + c(2 + (V + 3)GD2)2 + 2 + V ] .
Let ε̂ = Nε with 0 < ε̂ ≤ N . Then, we need T = Nβε̂ for some β > 0 to guarantee min0≤t≤T−1[F (w (t)) − F∗] ≤ 1T ∑T−1 t=0 [F (w
(t)) − F∗] ≤ ε̂. Hence, the total complexity is O ( nN
3β ε̂3 (D 2H2 + (ε̂2/N)) log(Nε̂ ) ) .

Review summaries:
This paper presents a new optimization method for finding global minima of nonconvex finite sum problems. In particular, the summands are functions of the form $\phi_{i}\circ h$ where $\phi_{i}$ is convex and Lipschitz smooth, while $h$ is nonconvex. Each iteration of the method consists of solving an auxiliary regularized least squares (RLS) problem, followed by a gradient step. Additional analysis is given for when the RLS problem is solved inexactly. Finally, a claimed $\tilde{O}(\varepsilon^{-3})$ complexity is established under strong boundedness assumptions on various solution sets.

The paper provides a new gradient-based algorithm. The algorithm is based on the observation that a loss function for a single sample can be written as composition of two functions (the logits and the actual loss function). It computes the direction by means of solving a quadratic MSE problem. They provide a convergent analysis of the algorithm (there is a version where the quadratic problem is solved explicitly through a closed form expression and a version where gradient descent is applied to solve the problem approximately). There are no computational experiments.   The main contributions are in the algorithms themselves and the accompanying analyses. It's unclear how novel are the proof techniques since everything resembles second order algorithms. The authors also claim the actual reformulation to be a novel contribution however such a formulation is straightforward and used in many contexts (some of my lecture slides from years ago show the formulation used by authors as a possible formulation for the overall loss function). Despite of this, the design of the algorithm should get credit.