Title: THEORETICAL PROPERTIES OF THE GLOBAL OPTIMIZER OF TWO-LAYER NEURAL NETWORK Anonymous authors Paper under double-blind review

Abstract: In this paper, we study the problem of optimizing a two-layer artificial neural network that best fits a training dataset. We look at this problem in the setting where the number of parameters is greater than the number of sampled points. We show that for a wide class of differentiable activation functions (this class involves most nonlinear functions and excludes piecewise linear functions), we have that arbitrary first-order optimal solutions satisfy global optimality provided the hidden layer is non-singular. We essentially show that these non-singular hidden layer matrix satisfy a "good" property for these big class of activation functions. Techniques involved in proving this result inspire us to look at a new algorithmic, where in between two gradient step of hidden layer, we add a stochastic gradient descent (SGD) step of the output layer. In this new algorithmic framework, we extend our earlier result and show that for all finite iterations the hidden layer satisfies the"good" property mentioned earlier therefore partially explaining success of noisy gradient methods and addressing the issue of data independency of our earlier result. Both of these results are easily extended to hidden layers given by a flat matrix from that of a square matrix. Results are applicable even if network has more than one hidden layer provided all inner hidden layers are arbitrary, satisfy non-singularity, all activations are from the given class of differentiable functions and optimization is only with respect to the outermost hidden layer. Separately, we also study the smoothness properties of the objective function and show that it is actually Lipschitz smooth, i.e., its gradients do not change sharply. We use smoothness properties to guarantee asymptotic convergence of O(1/number of iterations) to a first-order optimal solution.

1. INTRODUCTION:

Neural networks architecture has recently emerged as a powerful tool for a wide variety of applications. In fact, they have led to breakthrough performance in many problems such as visual object classification (Krizhevsky et al., 2012) , natural language processing (Collobert & Weston, 2008) and speech recognition (Mohamed et al., 2012) . Despite the wide variety of applications using neural networks with empirical success, mathematical understanding behind these methods remains a puzzle. Even though there is good understanding of the representation power of neural networks (Barron, 1994) , training these networks is hard. In fact, training neural networks was shown to be NP-complete for single hidden layer, two node and sgn(•) activation function (Blum & Rivest, 1988) . The main bottleneck in the optimization problem comes from non-convexity of the problem. Hence it is not clear how to train them to global optimality with provable guarantees. Neural networks have been around for decades now. A sudden resurgence in the use of these methods is because of the following: Despite the worst case result by Blum & Rivest (1988) , first-order methods such as gradient descent and stochastic gradient descent have been surprisingly successful in training these networks to global optimality. For example, Zhang et al. (2016) empirically showed that sufficiently over-parametrized networks can be trained to global optimality with stochastic gradient descent. Neural networks with zero hidden layers are relatively well understood in theory. In fact, several authors have shown that for such neural networks with monotone activations, gradient based methods will converge to the global optimum for different assumptions and settings (Mei et al., 2017; Hazan et al., 2015; Kakade et al., 2011; Kalai & Sastry, 2009) .

Despite the hardness of training the single hidden layer (or two-layer) problem, enough literature is available which tries to reduce the hardness by making different assumptions. E.g., Choromanska et al. (2014) made a few assumptions to show that every local minimum of the simplified objective is close to the global minimum. They also require some independent activations assumption which may not be satisfied in practice. For the same shallow networks with (leaky) ReLU activations, it was shown in Soudry & Carmon (2016) that all local minimum are global minimum of the modified loss function, instead of the original objective function. Under the same setting, Xie et al. (2016) showed that critical points with large "diversity" are near global optimal. But ensuring such conditions algorithmically is difficult. All the theoretical studies have been largely focussed on ReLU activation but other activations have been mostly ignored. In our understanding, this is the first time a theoretical result will be presented which shows that for almost all nonlinear activation functions including softplus, an arbitrary first-order optimal solution is also the global optimal provided certain "simple" properties of hidden layer. Moreover, we show that a stochastic gradient descent type algorithm will give us those required properties for free for all finite number of iterations hence even if the hidden layer variables are data dependent, we still get required properties. Our assumption on data distribution is very general and can be reasonable for practitioners. This comes at two costs: First is that the hidden layer of our network can not be wider than the dimension of the input data, say d. Since we also look at this problem in over-parametrized setting (where there is hope to achieve global optimality), this constraint on width puts a direct upper-bound of d 2 on the number of data points that can be trained. Even though this is a strong upper bound, recent results from margin bounds (Neyshabur et al., 2017) show that if optimal network is closer to origin then we can get an upper bound on number of samples independent of dimension of the problem which will ensure closeness of population objective and training objective. Second drawback of this general setting is that we can prove good properties of the optimization variables (hidden layer weights) for only finite iterations of the SGD type algorithm. But as it is commonly known, stochastic gradient descent converges to first order point asymptotically so ideally we would like to prove these properties for infinitely many iterations. We compare our results to some of the prior work of Xie et al. (2016) and Soudry & Carmon (2016) . Both of these papers use similar ideas to examine first order conditions but give quite different results from ours. They give results for ReLU or Leaky ReLU activations. We, on the other hand, give results for most other nonlinear activations, which can be more challenging. We discuss this in section 3 in more detail. We also formally show that even though the objective function for training neural networks is nonconvex, it is Lipschitz smooth meaning that gradient of the objective function does not change a lot with small changes in underlying variable. To the best of our knowledge, there is no such result formally stated in the literature. Soltanolkotabi et al. (2017) discuss similar results, but there constant itself depends locally on w max , a hidden layer matrix element, which is variable of the the optimization function. Moreover, there result is probabilistic. Our result is deterministic, global and computable. This allows us to show convergence results for the gradient descent algorithm, enabling us to establish an upper bound on the number of iterations for finding an ε-approximate first-order optimal solution ( ∇f () ≤ ε). Therefore our algorithm will generate an ε-approximate first-order optimal solution which satisfies aforementioned properties of the hidden layer. Note that this does not mean that the algorithm will reach the global optimal point asymptotically. As mentioned before, when number of iterations tend to infinity, we could not establish "good" properties. We discuss technical difficulties to prove such a conjecture in more detail in section 5 which details our convergence results. At this point we would also like to point that there is good amount of work happening on shallow neural networks. In this literature, we see variety of modelling assumptions, different objective functions and local convergence results. Li & Yuan (2017) focuses on a class of neural networks which have special structure called "Identity mapping". They show that if the input follows from Gaussian distribution then SGD will converge to global optimal for population objective of the "identity mapping" network. Brutzkus & Globerson (2017) show that for isotropic Gaussian inputs, with one hidden layer ReLU network and single non-overlapping convolutional filter, all local minimizers are global hence gradient descent will reach global optimal in polynomial time for the population objective. For the same problem, after relaxing the constraint of isotropic Gaussian inputs, they show that the problem is NP-complete via reduction from a variant of set splitting problem. In both of these studies, the objective function is a population objective which is significantly different from training objective in over parametrized domain. In over-parametrized regime, Soltanolkotabi et al. (2017) shows that for the training objective with data coming from isotropic Gaussian distribution, provided that we start close to the true solution and know maximum singular value of optimal hidden layer then corresponding gradient descent will converge to the optimal solution. This is one of its kind of result where local convergence properties of the neural network training objective function have studied in great detail. Our result differ from available current literature in variety of ways. First of all, we study the training problem in the over-parametrized regime. In that regime, the training objective can be significantly different from population objective. Moreover, we study the optimization problem for many general non-linear activation functions. Our result can be extended to deeper networks when considering the optimization problem with respect to outermost hidden layer. We also prove that stochastic noise helps in keeping the aforementioned properties of hidden layer. This result, in essence, provides justification for using stochastic gradient descent. Another line of study looks at the effect of over-parametrization in the training of neural networks (Haeffele & Vidal, 2015; Nguyen & Hein, 2017) . These result are not for the same problem as they require huge amount of over-parametrization. In essence, they require the width of the hidden layer to be greater than number of data points which is unreasonable in many settings. These result work for fairly general activations as do our results but we require a moderate over-parametrization, width × dimension ≥ number of data population, much more reasonable in practice as pointed before from margin bound results. They also work for deeper neural network as do our results when optimization is with respect to outermost hidden layer (and aforementioned technical properties are satisfied for all hidden layers).

2. NOTATION AND PROBLEM OF INTEREST:

We define set [q] := {1, . . . , q}. For any matrix A ∈ R a×b , we write vect(A) ∈ R ab×1 as vector form of the matrix A. For any vector z ∈ R k , we denote h(z)

:= h(z[1]), . . . , h(z[k])

T , where

z[i]

is the i-th element in vector z. B i (r) represents a l i -ball of radius r, centred at origin. We define component-wise product of two vectors with operator . We say that a collection of vectors,

{v i } N i=1 ∈ R d , is full rank if rank v 1 . . . v N = min{d, N }.

Similarly, we say that collection of matrices,

{M i } N i=1 ∈ R n×d , is full rank if rank [vect(M 1 ) . . . vect(M k )] = min{N, nd}.

A fully connected two-layer neural network has three parameters: hidden layer W , output layer θ and activation function h. For a given activation function, h, we define neural network function as

φ W,θ (u) := θ T h(W u).

In the above equation, W ∈ R n×d is hidden layer matrix, θ ∈ R n is the output layer. Finally h : R → R is an activation function.

The main problem of interest in this paper is the two-layer neural network problem given by min

W ∈R n×d θ∈R n f (W, θ) := 1 2N N i=1 (v i -φ W,θ (u i )) 2 . (2.1)

In this paper, we assume that

(u i , v i ) ∈ R d × R, i ∈ [N ]

are independently distributed data point and each u i is sampled from a d-dimensional Lebesgue measure.

3. THE BASIC IDEA AND THE ALGORITHM:

First-order optimality condition for the problem defined in (2.1), with respect to W [j, k] (j-th row, k-th column element of matrix

W ) ∀ j ∈ [n], ∀ k ∈ [d] is ∇ W f (W, θ)[j, k] = 1 N N i=1 {v i -θ T h(W u i )}h (W [j, :]u i )θ[j]u i [k] = 0. (3.1) Equation (3.1) is equivalent to N i=1 {v i -θ T h(W u i )} h (W u i ) θ u i T = 0. (3.2)

(3.1) can also be written in a matrix vector product form:

EQUATION

where

D :=    h (W [1, :]u 1 )θ[1]u 1 . . . h (W [1, :]u N )θ[1]u N . . . . . . . . . h (W [d, :]u 1 )θ[d]u 1 . . . h (W [d, :]u N )θ[d]u N    and s :=    v 1 -θ T h(W u 1 ) . . . v N -θ T h(W u N )    .

Notice that if matrix D ∈ R nd×N is of full column rank (which implies nd ≥ N , i.e., number of samples is less than number of parameters) then it immediately gives us that s = 0 which means such a stationary point is global optimal. This motivates us to investigate properties of h under which we can provably keep matrix D full column rank and develop algorithmic methods to help maintain such properties of matrix D.

Note that similar approach was explored in the works of Soudry & Carmon (2016) Moreover, we look at the algorithmic process of finding W differently. We show that SGD will achieve full rank property of matrix D with probability 1 for all finite iterations. So this is essentially a property of the algorithm and not of the model. Even if that is the case, to show global optimality, we need to prove that matrix D is full column rank in asymptotic sense.

That question was partly answered in Xie et al. (2016) . They show that matrix D is full column rank by achieving a lower bound on smallest singular value of matrix D. But to get this, they need two facts. First, the activation function has to be ReLu so that they can find the spectrum of corresponding kernel matrix. Second, they require a bound on discrepancy of weights W . These conditions are strong in the sense that they restrict the analysis to a non-differentiable activation function and finding an algorithm satisfying discrepancy constraint on W can be a difficult task. On other hand, our results are proved for a simple SGD type algorithm which is easy to implement. But we do not get a lower bound on singular value of D in asymptotic sense. There are obvious pluses and minuses for both types of results.

For the rest of the discussion, we will assume that n = d (our results can be extended to case n ≤ d easily) and hence W is a square matrix. In this setting, we develop the following algorithm whose output is a provable first-order approximate solution. Here we present the algorithm and in next sections we will discuss conditions that are required to satisfy full rank property of matrix D as well as convergence properties of the algorithm.

In Algorithm 1, we use techniques inspired from alternating minimization to minimize with respect to θ and W . For minimization with respect to θ, we add gaussian noise to the gradient information. This will be useful to prove convergence of this algorithm. We use randomness in θ to ensure some "nice" properties of W which help us in proving that matrix D generated along the trajectory of Algorithm 1 is full column rank. More details will follow in next section.

The algorithm has two loops. An outer loop implements a single gradient step with respect to hidden layer, W . For each outer loop iteration, there is an inner loop which optimizes objective function with respect to θ using a stochastic gradient descent algorithm. In the stochastic gradient descent, we generate a noisy estimated of ∇ θ f (W, θ) as explained below.

Let ξ ∈ R d be a vector whose elements are i.i.d. Gaussian random variable with zero mean. Then for a given value of W we define stochastic gradient w.r.t. θ as follows:

G W (θ, ξ) = ∇ θ f (W, θ) + ξ. (3.4) Then we know that E[G W (θ, ξ)] = ∇ θ f (W, θ). We can choose a constant σ > 0 such that following holds E G W (θ, ξ) -∇ θ f (W, θ) 2 ≤ σ 2 . (3.5)

Moreover, in the algorithm we consider a case where θ ∈ R. Note that R can be kept equal to R d but that will make parameter selection complicated. In our convergence analysis, we will use

EQUATION

for some constant R, to make parameter selection simpler. We use prox-mapping P x : R d → R as follows:

P x (y) = argmin z∈R y, z -x + 1 2 z -x 2 . (3.7)

In case R is a ball centred at origin, solution of (3.7) is just projection of x -y on that ball. For case where R = R d then the solution is quantity x -y itself.

Algorithm 1 SGD-GD Algorithm procedure W 0 ← Random d × d matrix θ 0 ← Random d vector

Initialize N o to predefined iteration count for outer ietaration Initialize N i to predefined iteration count for inner iteration Begin outer iteration:

for k = 0, 1, 2, . . . , N o do θ 1 ← θ k Begin inner iteration: for i = 1, 2, . . . , N i do θ i+1 ← P θi (β i G W k (θ i , ξ k i )) θ av i+1 = i τ =1 β τ -1 i τ =1 β τ θ τ +1 end for θ k+1 ← θ av Ni+1 W k+1 ← W k -γ k ∇ W f (W k , θ k+1 ) end for return {W No+1 ; θ No+1 } end procedure

Notice that the problem of minimization with respect to θ is a convex minimization problem. So we can implement many procedures developed in the Stochastic optimization literature to get the convergence to optimal value (Nemirovski et al., 2009) . In the analysis, we note that one does not even need to implement complete inner iteration as we can skip the stochastic gradient descent suboptimally given that we improve the objective value with respect to where we started, i.e.,

f (W k , θ k+1 ) ≤ f (W k , θ k ). (3.8)

In essence, if evaluation of f for every iteration is not costly then one might break out of inner iterations before running N i iterations. If it is costly to evaluate function values then we can implement the whole SGD for convex problem with respect to θ as specified in inner iteration of the algorithm above. In each outer iteration, we take one gradient decent step with respect to variable W . We have total of N o outer iterations. So essentially we evaluate

∇ θ f (W, •) a total of N o N i times and ∇ W f (•, θ) total of N o times.

Overall, this algorithm is new form of alternate minimization, where one iteration can be potentially left suboptimally and other one is only one gradient step.

4. ARBITRARY POINT SATISFYING FIRST ORDER OPTIMALITY CONDITIONS:

We prove in this section that arbitrary first order optimal points are globally optimal. One does not expect to have an arbitrary first order optimal point because it has to depend on data. We still would like to put our analysis here because that inspires us to consider a new algorithmic framework in Section 3 providing similar results for all finite iterations of the algorithm. We say that h : R → R satisfy the condition "C1" if

• ∀ interval (a, b), {c 1 , c 2 , c 3 } ∈ R 3 s.t. {h (x) = c 1 , ∀x ∈ (a, b)} or {(x + c 2 )h (x) + h(x) = c 3 , ∀x ∈ (a, b)}.

One can easily notice that most activation functions used in practice e.g.,

• (Softplus) h(x) := ln(1 + e x ), • (Sigmoid) h(x) := 1 1+e -x , • (Sigmoid symmetric) h(x) := 1-e -x 1+e -x , • (Gaussian) h(x) := e -x 2 , • (Gaussian Symmetric) h(x) := 2e -x 2 -1, • (Elliot) h(x) := x 2(1+|x|) + 0.5, • (Elliot Symmetric) h(x) := x 1+|x| , • (Erf) h(x) := 2 √ π x 0 e -t 2 /2 dt, • (Hyperbolic tangent) h(x) := tanh(x),

satisfy the condition C1. Note that h (x) also satisfy condition C1 for all of them. In fact, except for very small class of functions (which includes linear functions), none of the continuously differentiable functions satisfy condition C1. We first prove a lemma which establishes that columns of the matrix D (each column is a vector form of d×d matrix itself) are linearly independent when W = I d and h satisfies condition C1. We later generalise it to any full rank W using a simple corollary. The statement of following lemma is intuitive but its proof is technical. This means that if u i in the Problem (2.1) are coming from a Lebesgue measure then by corollary 4.2 we have h(W u i )u i T will be a full rank collection given that we have maintained full rank property of W . Now note that in the first-order condition, given in (3.3), row of matrix D are scaled by constant factors θ[j]'s, j ∈ [d] . Notice that we may assume θ[j] = 0 because otherwise there is no contribution of corresponding j-th row of W to the Problem (2.1) and we might as well drop it entirely from the optimization problem. Hence we can rescale rows of matrix D by factor 1

θ[j]

without changing the rank. In essence, corollary 4.2 implies that matrix D is full rank when W is full rank. So by our discussion in earlier section, we show that satisfying first-order optimality is enough to show global optimality under condition C1 for data independent W. Remark 4.3 Due to Lemma 4.1 and corollary 4.2 then, rank of collection h(u i )u i T is invariant under any rotation.

Remark 4.4 As a result of corollary above one can see that the collection of vectors h(W x i ) is full rank under the assumption that W is non-singular, x i ∈ R d are independently sampled from Lebesgue measure and h satisfies condition C1.

Remark 4.5 Since collection h(W u i ) is also full rank, we can say that z i := h(W 1 u i ) are independent and sampled from a Lebesgue measure for a non-singular matrix W 1 . Applying the Lemma to z i , we have collection of matrices g(W 2 z i )z i T are full rank with measure 1 for non-singular W 2 and g satisfying condition C1. So we see that for multiple hidden layers satisfying non-singularity, we can apply full rank property for collection of gradients with respect to outermost hidden layer.

Remark 4.6 If W ∈ R n×d is such that n ≤ d and W is full row rank, then we can extend its basis to create W and apply corollary 4.2 to get that h(W u i )u i T is full rank with measure 1. So this implies that h(W u i )u i T must have been full rank with probability 1 otherwise we will have contradiction.

Remark 4.7 We can extend corollary 4.2 to a general result that h(W u i )u i T has rank min{rank(W )d, N } with measure 1 by removing dependent rows and using remark 4.6.

5. CONVERGENCE RESULTS FOR STOCHASTIC GRADIENT DESCENT:

Even though we have proved that collection h(W u i )u i T N i=1 is full rank in the previous section, we need such W to be independent of data. In general, any algorithm will use data to find W and it appears that results in previous section are not meaningful in practice. However, the analysis of Lemma 4.1 motivates the idea that stochastic noise of θ might help in obtaining the required properties of W and D along the trajectory of Algorithm 1. In this section we first prove that by using random noise in stochastic gradient on θ gives a non-singular W k in every iteration. Then using this fact, we prove that matrix D generated along the algorithm is also full rank. The proof techniques are very similar to proof of Lemma 4.1. Later on, we will also show that overall algorithm will converge to approximate first-order optimal solution to Problem (2.1) by using smoothness properties. It should be noted however that this can not guarantee convergence to a global optimal solution. To prove such a result, one needs to analyze the smallest singular value of random matrix D, defined in (3.3). More specifically, we have to show that σ min (D) decreases at the rate slower than the first-order convergence rate of the algorithm so that the overall algorithm converges to the global optimal solution. Even if it is very difficult to prove such a result in theory, we think that such an assumption about σ min (D) is reasonable in practice. Now we analyze the algorithm. For the sake of simplicity of notation, let us define

ξ [k] := {ξ 1 [Ni] , . . . , ξ k [Ni] } (5.1) and ξ j [Ni] = {ξ j 1 . . . ξ j Ni }, (5.2)

where N i is the inner iteration count in Algorithm 1. Essentially ξ [k] contains the record of all random samples used until the k-th outer iteration in Algorithm 1 and ξ j [Ni] contains record of all random samples used in the inner iterations of j-th outer iteration.

Lemma 5.1 Pr{∃ v such that W k v = 0 ξ [k-1] } = 0, ∀ k ≥ 0,

where W k are matrices generated by Algorithm 1 and measure Pr{.

ξ [k-1] } is w.r.t. random variables ξ k [Ni] .

Now that we have proved that W k 's generated by the algorithm are full rank, we show that matrix D generated along the trajectory of the algorithm is full rank for any finite number of iterations. We use techniques inspired from Lemma 4.1 but this time we use Lebesgue measure over Θ rather than data. Over randomness of Θ, we can show that our algorithm will not produce any W such that corresponding matrix D is rank deficient. Since Θ is essentially designed to be independent of data so we will not produce rank deficient D throughout the process of randomized algorithm.

Lemma 5.2 Suppose W = W + D v Z where D v := diag(v[i], i ∈ [d]

) and v is a random vector with Lebesgue measure in R d . W , Z ∈ R d×d and Z = 0. Let h be a function which follows condition C1. Also assume that W is full rank with measure 1 over randomness of v. Then h(W u i )u i T is full rank with measure 1 over randomness of v.

Lemma 5.3 Collection of matrices h (W k+1 u i )u i T are full rank with measure 1, where the measure is over randomness of ξ k+1

[Ni]

Proof. We know that

W k+1 = W k + γ k Θ k+1 N j=1 h (W k u j )u j T (v i -θ T h(W k u j )).

Now apply Lemma 5.2 to obtain the required result.

Note that Lemma 5.3 is very similar to the result in section 4. Some remarks are in order.

Remark 5.4 As a result of Lemma 5.3 above, one can see that collection of vectors h(W k u i ) is full rank for all finite iterations of Algorithm 1.

Remark 5.5 If we have a neural network with multiple hidden layer, we can assume that inner layers are initialized to arbitrary full rank matrices and we are optimizing w.r.t. outermost hidden layer. Corollary 4.2 and Remark 4.4 give us that input to outermost hidden layer are independent vectors sampled from some lebesgue measure. Then applying Algorithm 1 to optimize w.r.t. outermost hidden layer will give us similar results as mentioned in Lemma 5.3.

Hence we showed that algorithm will generate full rank matrix D for any finite iteration. Now to prove convergence of the algorithm, we need to analyze the function f (defined in (2.1)) itself. We show that f is a Lipschitz smooth function for any given instance of data {u i , v i } N i=1 . This will give us a handle to estimate convergence rates for the given algorithm.

Lemma 5.6 Assuming that h : R → R is such that its gradients, hessian as well as values are bounded and data {u i , v i } N i=1 is given then there exists a constant L such that

∇ W f (W 1 , θ) -∇ W f (W 2 , θ) F ≤ L W 1 -W 2 F . (5.3)

Moreover, a possible upper bound on L can be as follows:

L ≤ 1 N θ max L h N i=1 u i 2 2 |v i | + √ 2dL hh θ 2 N i=1 u i 2 2

Remark 5.7 Before staing the proof, we should stress that assumptions on h is satisfied by most activation functions e.g., sigmoid, sigmoid symmetric, gaussian, gaussian symmetric, elliot, elliot symmetric, tanh, Erf.

Remark 5.8 Note that one can easily estimate value of L given data and θ. Moreover, if we put constraints on θ 2 then L is constant in every iteration of the algorithm 1. As mentioned in section 3, this will provide an easier way to analyze the algorithm.

Lemma 5.9 Assuming that scalar function h is such that |h(•)| ≤ u then there exists L θ s.t.

∇ w f (W, θ 1 ) -∇ w f (W, θ 2 ) 2 ≤ L θ θ 1 -θ 2 2 (5.4)

Notice that Lemma 5.9 gives us value of L θ irrespective of value of W or data. Also observe that f (W, •) is convex function since hessian

∇ 2 θ f (W, θ) = 1 N N i=1 h(W u i )h(W u i ) T ,

which is the sum of positive semidefinite matrices. By Lemma 5.9, we know that f (W, •) is smooth as well. So we can use following convergence result provided by Lan (2012) for stochastic composite optimization. A simplified proof can be found in appendix.

Theorem 5.10 Assume that stepsizes β i satisfy 0 < β i ≤ 1/2L θ , ∀ i ≥ 1. Let {θ av i+1 } i≥1 be the sequence computed according to Algorithm 1. Then we have,

E[f (W k , θ av i+1 ) -f (W k , θ * W k )] ≤ K 0 (i), ∀ i ≥ 1, ∀ k ≥ 0, (5.5)

where K 0 (i) :

= i τ =1 β τ -1 θ 1 -θ * W k 2 2 + σ 2 i τ =1

β 2 i where θ 1 is the starting point for inner iteration and σ is defined in (3.5). Now we look at a possible strategy of selecting stepsize β i . Suppose we adopt a constant stepsize policy then we have

β i = β, ∀ i ∈ [N i ]. Then we have E[f (W k , θ av Ni+1 ) -f (W k , θ * W k )] ≤ θ 1 -θ * W k 2 N i β + σ 2 β. Now if we choose β = min 1 2L θ , 1 N i σ 2 ,

(5.6)

we get

E[f (W k , θ av Ni+1 ) -f (W k , θ * W k )] ≤ θ 1 -θ * W k 2 2L θ N i + σ √ N i + σ √ N i .

By Lemma 5.6, the objective function for neural networks is Lipschitz-smooth with respect to the hidden layer, i.e., it satisfies eq (5.

3). Notice that it is equivalent to saying

f (W 2 , w) -f (W 1 , w) -∇ W f (W 1 , w), W 2 -W 1 ≤ L 2 W 1 -W 2 2 F , ∀ W 1 , W 2 ∈ R d×d .

(5.7) Since we have a handle on the smoothness of objective function, we can provide a convergence result for the overall algorithm.

Theorem 5.11 Suppose γ k < 2 L then we have

EQUATION

) where R/2 is the radius of origin centred ball, R in algorithm, defined as

R := {r ∈ R d : r 2 ≤ R 2 }.

In view of Theorem 5.11, we can derive a possible way of choosing γ k , σ and N i to obtain a convergence result. More specifically, if

N i = N o , σ = 1 √ Ni , γ k = 1

L and β k τ is chosen according to (5.6) then we have

E min k=0,...,N ∇ W f (W k , θ k+1 ) 2 F ≤ 2L f (W 0 , θ 0 ) + R 2 (L θ + 1/2) + 1 N o

Note that in the algorithm 1, we have proved that having a stochastic noise helps keeping matrix D full rank for all finite iterations. Then in Theorem 5.11, we showed a methodical way of achieving approximate first order optimal. So essentially at the end of the finitiely many steps of algorithm 1, we have a point W which satisfies full rank property of D and is approximately first order optimal.

We think this kind of result can be extended to variety of different first order methods developed for Lipschitz-smooth non-convex optimization problems. More specifically, accelerated gradient method such as unified accelerated method proposed by Ghadimi et al. (2015) or accelerated gradient method by Ghadimi & Lan (2016) can be applied in outer iteration. We can also use stochastic gradient descent method for outer iteration. For this, we need a stochastic algorithm that is designed for non-convex and Lipschitz smooth function optimization. Randomized stochastic gradient method, proposed by Ghadimi & Lan (2013) , Stochastic variance reduction gradient method (SVRG) by Reddi et al. or Simplified SVRG by Allen-Zhu & Hazan can be employed in outer iteration. Convergence of these new algorithms will follow immediately from the convergence results of respective studies. Some work may be needed to prove that they hold matrix D full rank. We leave that for the future work. We also leave the problem of proving a bound on singular value for future. This will close the gap between empirical results and theory. Value of Lipschitz constant, L, puts a significant impact on the running time of the algorithm. Notice that if L increases then correspondingly N o and N i increase linearly with L. So we need methods by which we can reduce the value of the estimate of L. One possible idea would be to use l 1 -ball for feasible region of θ. More specifically, if R = B 1 (R/2) then we can possibly enforce sparsity on θ which will allow us to put better bound on L.

A PROOFS OF AUXILIARY RESULTS:

In this appendix, we provide proofs for auxiliary results.

A.1 PROOF OF LEMMA 4.1

The result is trivially true for d =1, we will show this using induction on d. Define

v i := vect(h(x i )x i T ), i ∈ [N ]. Note that it suffices to prove independence of vector v i , i ∈ [N ] for N ≤ d 2 .

Now for sake of contradiction assume that v i , i ∈ [N ], are linearly dependent with positive joint measure on x i , i ∈ [N ] which is equivalent to positive measure on individual x i , ∀ i ∈ [N ] due to independence of vectors x i . Since x i 's are sampled from Lebesgue measure so positive measure on x i , ∀i ∈ [N ], implies there exists a d-dimensional volume for each x i such that corresponding v i are linearly dependent. We can assume volume to be d-dimensional hyper-cuboid

Z i := {x ∈ R d : a i < x < b i }, ∀ i ∈ [N ]

(otherwise we can inscribe a hyper-cuboid in that volume). Notice that since Z i is a d-dimensional hyper-cuboid so

a i [k] < b i [k], ∀ i ∈ [N ], ∀ k ∈ [d].

Moreover, for any collection satisfying x i ∈ Z i , corresponding collection of vector v i are linearly dependent, i.e.,

v 1 = µ 2 v 2 + • • • + µ N v N , such that ∀i ∈ [N ], x i ∈ Z i . (A.1)

Noticing the definition of Z 1 , we can choose > 0 s.t. x 1 := x 1 + e 1 ∈ Z 1 . Since we ensure that x 1 ∈ Z 1 then by (A.1) we have

v 1 := vect(h( x 1 ) x 1 T ) = µ 2 v 2 + • • • + µ N v N . (A.2)

So using (A.1) and (A.2) we get

v 1 -v 1 = λ 2 v 2 + • • • + λ N v N . (A.3) Since h(x i )x i T [j, k] = h(x i [j])x i [k], we have h(x 1 )x 1 T [j, k] = h( x 1 ) x 1 T [j, k], ∀j ∈ {2, .

. . , d}, k ∈ {2, . . . , d}. So we have (d -1) 2 components of v 1 -v 1 are zero. Let us define:

w 1 =                                        (x 1 [1] + )h(x 1 [1] + ) -x 1 [1]h(x 1 [1]) h(x 1 [2]) . . . h(x 1 [d]) (2d -1) x 1 [2](h(x 1 [1] + ) -h(x 1 [1]))

. . .

x 1 [d](h(x 1 [1] + ) -h(x 1 [1])) , z 1 =        0 . . . (d -1) 2 0 , and notice that v 1 -v 1 = w 1 z 1 . Since > 0, w 1 = 0 with measure 1. Let y i := x i [2 : d] then last (d -1) 2 equations in (A.3) gives us λ 2 h(y 2 )y 2 T + • • • + λ N h(y N )y N T = z 1 = 0 (A.4)

By definition we have y i ∈ R d-1 are independently sampled from (d -1)-dimensional Lebesgue measure. So by inductive hypothesis, rank of collection of matrices h(y i )y

i T , i ∈ {2, . . . , N } = min{(d -1) 2 , N -1}. So if N -1 ≤ (d -1) 2 then λ 2 = • • • = λ N = 0

with measure 1, then by (A.3) we have w 1 = 0 with measure 1, which is contradiction to the fact that w 1 = 0 with measure 1. This gives us N > (d -1) 2 + 1 (A.5)

Notice that (A.4) in its matrix form can be written as linear system vect(h(y 2 )y 2 T ) . . . vect(h(y

N )y N T )    λ 2 . . . λ N    = 0 (A.6)

By (A.6), we have that vector of λ's lies in the null space of the matrix. Finally by inductive hypothesis and (A.5) we conclude that the dimension of that space is

N -1 -(d -1) 2 {> 0}. Let u 1 , . . . , u N -1-(d-1) 2 ∈ R N -1 be the basis of that null space i.e. vect(h(y 2 )y 2 T ) . . . vect(h(y N )y N T ) u j = 0, ∀ j ∈ {1, N -1 -(d -1) 2 }

Define t i ∈ R 2d-1 as:

t i :=          x i [1]h(x i [1]) . . . x i [1]h(x i [d]) x i [2]h(x i [1]) . . . x i [d]h(x i [1])         

then we can rewrite (A.3) as

w 1 z 1 = t 2 . . . t N vect(h(y 2 )y 2 T ) . . . vect(h(y N )y N T ) u 1 . . . u N -1-(d-1) 2    λ 2 . . . λ N -(d-1) 2    (A.7) which implies that w 1 = λ 2 v 2 + • • • + λ N -(d-1) 2 v N -(d-1) 2 (A.8)

where

v i = t 2 . . . t N u i-1 , i = 2, . . . , N -(d -1) 2 and z 1 part of the equation is already satisfied due to selection of null space. Since N ≤ d 2 ⇒ N -1 -(d -1) 2 ≤ 2d -2 then 2d -1 equations specified in (A.8) are consistent in ≤ (2d -2) variables.

Hence we get linearly dependent equations ∀x 1 1 ∈ (a 1 1 , b 1 1 ) and small enough. Since x 2 , . . . , x N are kept constant, v 2 , . . . , v N are constant. So t 2 , . . . , t N are constants and we can choose the same basis of null space u 1 , . . . , u N -1-(d-1) 2

. Hence we have v 2 , . . . , v (N -(d-1) 2 ) are constant. Let us define the set S to be the index set of linearly independent rows of matrix [ v 2 . . . v N -(d-1) 2 ] and every other row is a linear combination of rows in S. Since (A.8) is consistent so the same combination must be valid for the rows of w 1 . Now if N ≤ d 2 -1 then number of variables in (A.8) is ≤ 2d -3 but number of equations is 2d -1, therefore at least two equations are linearly dependent on other equation. This implies last (2d -2) equations then function must be dependent on each other:

d j=2 α j h(x (1) [j]) + h(x (1) [1] + ) -h(x (1) [1]) d j=2 β j x (1) [j] = 0

for some fixed combination α j , β j . If we divide above equation by and take the limit as → 0 then we see that h satisfies following differential equation on interval (a 1 1 , b 1 1 ): h (x) = c 1 which is a contradiction to the condition C1! Clearly this leaves only one case i.e. N = d 2 and (2d -1) equations must satisfy dependency of the following form for all x

(1) 1 ∈ (a

(1) 1 , b (1) 1 ): (x (1) [1] + )h(x (1) [1] + ) -x (1) [1]h(x (1) [1]) = d j=2 α j h(x (1) [j]) + h(x (1) [1] + ) -h(x (1) [1]) d j=2 β j x (1) [j]

Again by similar arguments, the combination is fixed. Let H(x) = xh(x) then dividing above equation by and taking the limit as → 0, we can see that h satisfies following differential equation:

H (x) = c 1 + c 2 h (x) ⇒ (x -c 2 )h (x) + h(x) = c 1 (A.9

) which is again a contradiction to the condition C1 So we conclude that for N ≤ d 2 there does not exist hyper-cuboids Z i such that vol(Z i ) > 0 and for all x i ∈ Z i , corresponding v i are linearly dependent. So we get rank of collection {v i } N i=1 = min{N, d 2 } with measure 1.

A.2 PROOF OF COROLLARY 4.2 Let us define x := W u be another random variable. Since W is full rank and u has Lebesgue measure ⇒ x has Lebesgue measure. Now we claim that the collection h(W u i )u i T is full rank iff the collection h(x i )x i T is full rank. This can observed as follows:

N i=1 λ i h(x i )x i T = 0 ⇔ N i=1 λ i h(W u i )u i T W T = 0 ⇔ N i=1 λ i h(W u i )u i T = 0

Here the second statement follows from the fact W is a non-singular matrix. Now by Lemma 4.1 we have that collection h(x i )x i T is linearly independent with measure 1. So h(W u i )u i T is linearly independent with measure 1. Since any rotation is U is a full rank matrix so we have the result.

A.3 PROOF OF LEMMA 5.1 This is true for k = 0 trivially since we are randomly sampling matrix W 0 . We now show this by induction on k.

Recall that gradient of f (W, θ) with respect to W can be written as

N i=1 {v i -θ T h(W u i )} h (W u i ) θ u i T

. Notice that in effect, we are multiplying i-th row of the rank one matrix h (W u i )u i T by i-th element of vector θ. So this can be rewritten as a matrix product

N i=1 {v i -θ T h(W u i )}Θh (W u i )u i T ,

where Θ := diag{θ[i], i = 1, . . . , d}. So iterative update of the algorithm can be given as

W k+1 = W k -γ k Θ k+1 ∇ W f (W k , θ k+1 ), ∀ k ≥ 0.

Notice that given ξ [k] , vector θ k+1 and corresponding diagonal matrix Θ k+1 are found by SGD in the inner loop so θ k+1 is a random vector. More specifically, since

{ξ k+1 i } Ni i=1 is sequence of (A.22)

It then follows from (A.20),(A.21) and (A.22) that

EQUATION

Finally, using (A.18), (A.19) and (A.22), we get a possible finite upper bound on the value of L:

L ≤ 1 N θ max L h N i=1 u i 2 2 |v i | + √ 2dL hh θ 2 N i=1 u i 2 2

Also note that this bound is valid even if W is not a square matrix.

A.6 PROOF OF LEMMA 5.9

Noting that

-∇ θ f (W, θ) = 1 N N i=1 {v i -θ T h(W u i )}h(W u i ),

we have

∇ θ f (W, θ 1 ) -∇ w f (W, θ 2 ) 2 = 1 N N i=1 {v i -θ T 1 h(W u i )}h(W u i ) -{v i -θ T 2 h(W u i )}h(W u i ) 2 = 1 N N i=1 {-h(W u i )h(W u i ) T θ 1 + h(W u i )h(W u i ) T θ 2 } 2 = 1 N N i=1 h(W u i )h(W u i ) T (θ 2 -θ 1 ) 2 ≤ 1 N N i=1 h(W u i )h(W u i ) T 2 . θ 1 -θ 2 2 = 1 N N i=1 h(W u i )h(W u i ) T 2 . θ 1 -θ 2 2 = 1 N λ max N i=1 h(W u i )h(W u i ) T . θ 1 -θ 2 2 ≤ 1 N N i=1 λ max h(W u i )h(W u i ) T . θ 1 -θ 2 2 ∵ Weyl's Inequality = 1 N N i=1 h(W u i ) 2 2 . θ 1 -θ 2 2 ≤ u 2 d θ 1 -θ 2 2

where u 1 and u 2 are upper bounds on scalar functions |h(•)| and |h (•)| respectively and d is rowdimension of W .

A.7 PROOF OF THEOREM 5.11

We know by Lemma 5.6 that f (•, θ) is a Lipschitz smooth function. So using (5.7) we have . Now taking expectation with respect to ξ [No] (which is defined in (5.1)), we have

f (W k+1 , θ k+1 ) ≤ f (W k , θ k+1 ) + L 2 vect(W k+1 -W k ) 2 + vect(∇ W f (W k , θ k+1 )), vect(W 1 k+1 -W 1 k ) = f (W k , θ k+1 ) -γ k - L 2 γ 2 k vect(∇ W f (W k , θ k+1 )) 2 ≤ f (W k , θ k ) + Ni τ =1 β k τ -1 1 2 θ k -θ * W k 2 2 + Ni τ =1 β k τ ξ k τ , θ * W k -θ k τ + Ni τ =1 β k 2 τ ξ k τ 2 2(1 -L θ β k τ ) -(γ k - L 2 γ 2 k ) vect(∇ W G(W k , θ k+1 ))

E ξ k τ , θ * W k -θ k τ ξ [k-1] ∪ ξ k [τ -1] = 0, which implies E ξ [No ] ξ k τ , θ * W k -θ k τ

= 0. We also have E ξ by simply writing first order necessary condition for problem (A.25). Also we note that objective function in (A.25) is strongly convex with parameter 1. Then we have

β τ G τ , d τ + 1 2 d τ 2 2 + 1 2 u -θ τ +1 2 2 ≤ β τ G τ , u -θ τ + 1 2 u -θ τ 2 2 .

(A.26)

We will use eq (A.26) along with smoothness and convexity of the function f to get the final convergence result. Notice that due to smoothness, we have

β τ f (θ τ +1 ) ≤ β τ [f (θ τ ) + g(θ τ ), d τ + L θ 2 d τ 2 ] = β τ [f (θ τ ) + g(θ τ ), d τ ] + 1 2 d τ 2 - (1 -L θ β τ ) 2 d τ 2 .

then due to (A.24) we have,

β τ f (θ τ +1 ) ≤ β τ [f (θ τ ) + G τ , d τ ] -β τ ξ k τ , d τ + 1 2 d τ 2 - (1 -L θ β τ ) 2 d τ 2 ≤ β τ [f (θ τ ) + G τ , d τ ] + 1 2 d τ 2 - (1 -L θ β τ ) 2 d τ 2 + β τ ξ k τ . d τ ≤ β τ f (θ τ ) + β τ G τ , d τ + 1 2 d τ 2 + β 2 τ ξ k τ 2

2(1 -L θ β τ ) .

:



Figures and Tables:

Suppose x i ∈ R d are independently chosen vectors from any d-dimensional Lebesgue measure and let h : R → R be any function that satisfies condition C1 then collection of matrices h(x i )x i T , i ∈ [N ] are full rank with measure 1.Now Lemma 4.1 gives us a simple corollary:Corollary 4.2 If W is a nonsingular square matrix and u i ∈ R d is independently sampled from a Lebesgue measure then the collection of matrices h(W u i )u i T N i=1 is full rank with measure 1.

where the last inequality follows from equation (A.28) and (A.29). From (3.6), we have θ ≤ R/2 so L is constant for each outer iteration. Summing (A.23) from k = 0 to N o and dividing both side by k ) vect(∇ W1 f (W k , θ k+1 ))

[No ] ξ k τ 2 ≤ σ 2 , and hence E min k=0,...,N ∇ W f (W k , θ k+1 ) simplicity of notation, we define f (•) := f (W k , •), g(•) := ∇f (•) = ∇ θ f (W k , •) and G W k (θ τ , ξ k τ ) := G τ .Then from (3.4) and Algorithm 1 we get G τ = g(θ τ ) + ξ k τ . (A.24) Also define d τ := θ τ +1 -θ τ . Notice that θ τ +1 is optimal solution to the problem min u∈R d β τ G τ , u -θ τ +

and Xie et al. (2016). To get the full rank property for matrix D,Soudry & Carmon (2016) use leaky ReLu function. Basically this leaky activation function adds noise to entries of matrix D which allows them to show matrix D is full rank and hence all local minimums are global minimums. So this is essentially a change of model. We, on the other hand, do not change model of the problem.

Table content: <table/>

Peer review (iclr2018):
 I only got access to the paper after the review deadline; and did not have a chance to read it until now. Hence the lateness and brevity. The paper tackles an important theoretical question; and it offers results that are complementary to existing results (e.g., Soudry et al). However, the paper does not properly relate their results, assumptions in the context of the existing literature. Much explanation is needed in the author reply in order to clear these questions. The work should not be evaluated from a practical perspective as it is of a theoretical nature. I agree with most of the criticism raised by other reviewers. However, I also believe the authors managed to clear essentially of the criticism in they reply. The paper lacks in clarity as currently written. The results are interesting, but more explanation is needed for the main message to be conveyed more clearly. I suggest 7, but the paper has a potential to become 8 in my eyes in a future resubmission.