Diverse Neural Network Learns True Target Functions

Introduction

Neural networks are a powerful class of nonlinear functions which have been successfully deployed in a variety of machine learning tasks. In the simplest form, neural networks with one hidden layer are linear combinations of nonlinear basis functions (units),

where $\sigma(w_{k}^{\top}x)$ is a basis function with weights $w_{k}$, and $v_{k}$ is the corresponding combination coefficient. Learning with neural networks involves adapting both the combination coefficients and the basis functions at the same time, usually by minimizing the empirical loss

with first-order methods such as (stochastic) gradient descent. It is believed that basis function adaptation is a crucial ingredient for neural networks to achieve more compact models and better performance [Barron, 1993, Yang et al., 2014].

However, the empirical loss minimization problem involved in neural network training is non-convex with potentially numerous local minima and saddle points. This makes formal analysis of training neural networks very challenging. Given the empirical success of neural networks, a sequence of important and urgent scientific questions need to be investigated: Can neural networks corresponding to stationary points of the empirical loss learn the true target function? If the answer is yes, then what are the key factors contributing to their nice optimization properties? Based on these understandings, can we design better regularization schemes and learning algorithms to improve the training of neural networks?

In this paper, we provide partial answers to these questions by analyzing one-hidden-layer neural networks with rectified linear units (ReLU) in a least-squares regression setting. We show that neural networks with diverse units have no spurious local minima. More specifically, we show that the training loss of neural networks decreases in proportion to $\left\|\partial L/\partial W\right\|^{2}/s_{m}^{2}(D)$ where $\partial L/\partial W$ is the gradient and $s_{m}(D)$ is the minimum singular value of the extended feature matrix $D$ (defined in Section 3.1). The minimum singular value is lower bounded by two terms, where the first term is related to the spectrum property of the kernel function associate with the activation $\sigma(\cdot)$, and the second term quantifies the diversity of the units, measured by the classical notion of geometric discrepancy of a set of vectors. Essentially, the slower the decay of the spectrum, the better the optimization landscape; the more diverse the unit weights, the more likely stationary points will result in small training loss and generalization error.

We bypass the hurdle of non-convexity by directly analyzing the first order optimality condition of the learning problem, which implies that there are no spurious local minima if the minimum singular value of the extended feature matrix is large enough. Bounding the singular value is challenging because it entangles the nonlinear activation function, the weights and data in a complicated way. Unlike most previous attempts, we directly analyze the effect of nonlinearity without assuming independence of the activation patterns from actual data; in fact, the dependence of the patterns on the data and the unit weights underlies the key connection to activation kernel spectrum and the diversity of the units.

We have constructed a novel proof, which makes use of techniques from geometric discrepancy and kernel methods, and have identified a new relation linking the smallest singular value to the diversity of the units and the spectrum of a kernel function associated with the unit. More specifically,

We identify and separate two factors in the minimum singular value: 1) an ideal spectrum that is related to the kernel of the activation function and an ideal configuration of diverse unit weights; 2) deviation from the ideal spectrum measured by how far away actual unit weights are from the diverse configuration. This new perspective reveals benign conditions in learning neural networks.

We characterize the deviation from the ideal diverse weight configuration using the concept of discrepancy, which has been extensively studied in the geometric discrepancy theory. This reveals an interesting connection between the discrepancy of the weights and the training loss of neural networks. Therefore, it serves as a clean tool to analyze and verify the learning and the generalization ability of the networks.

Our results also suggest a novel regularization scheme to promote unit diversity for potentially better generalization. In [MarSra15], it is shown that diversity of the neurons leads to smaller network size and better performance.

Whenever possible, we corroborate our theoretical analysis with numerical simulations. These numerical results include computing and verifying the relationship between the discrepancy of a learned neural network and the minimum singular value. Additionally, we measure the effects on the discrepancy with and without regularization. In all these examples, the experiments match with the theory nicely and they accord with the practice of using gradient descent to learn neural networks.

Related work

Kernel methods have many commonalities with one-hidden-layer neural networks. The random feature perspective [Rahimi and Recht, 2009, Cho and Saul, 2009] views kernels as linear combinations of nonlinear basis functions, similar to neural networks. The difference between the two is that the weights are random in kernels while in neural networks they are learned. Using learned weights leads to considerable smaller models as shown in [Barron, 1993]. However it is a non-convex problem and it is difficult to find the global optima. e.g., one-hidden-layer networks are NP-complete to learn in the worst case [Blum and Rivest, 1993]. We will make novel use of techniques from kernel methods to analyze learning in neural networks.

The empirical success of training neural networks with simple algorithms such as gradient descent has motivated researchers to explain their surprising effectiveness. In [Choromanska et al., 2015], the authors analyze the loss surface of a special random neural network through spin-glass theory and show that for many large-size networks, there is a band of exponentially many local optima, whose loss is small and close to that of a global optimum. The analyzed polynomial network is different from the actual neural network being used which typically contains ReLU nowadays. Moreover, the analysis does not lead to a generalization guarantee for the learned neural network.

A similar work shows that all local optima are also global optima in linear neural networks [Kawaguchi, 2016]. However their analysis for nonlinear neural networks hinges on independence of the activation patterns from the actual data, which is unrealistic. Some other works try to argue that gradient descent is not trapped in saddle points [Lee et al., 2016, Ge et al., 2015], as was suggested to be the major obstacle in optimization [Dauphin et al., 2014]. There is also a seminal work using tensor method to avoid the non-convex optimization problem in neural network [Janzamin et al., 2015]. However, the resulting algorithm is very different from typically used algorithms where only gradient information of the empirical loss $L$ is used.

[Soudry and Carmon, 2016] is the closest to our work, which shows that zero gradient implies zero loss for all weights except an exception set of measure zero. However, this is insufficient to guarantee low training loss since small gradient can still lead to large loss. Furthermore, their analysis does not characterize the exception set and it is unclear a priori whether the set of local minima fall into the exception set.

Problem setting and preliminaries

where $\sigma(\cdot)=\max\{0,\cdot\}$ is the rectified linear unit (ReLU) activation function, $\left\{w_{k}\right\}$ and $\left\{v_{k}\right\}$ are the unit weights and combination coefficients respectively, $n$ is the number of units, and $C_{W}$ is some constant. We restrict $v_{k}\in\left\{-1,1\right\}$ due to the positive homogeneity of ReLU,

That is, the magnitude of $v_{k}$ can always be scaled into the corresponding $w_{k}$. For convenience, let

be the column concatenation of the unit parameters; also let $W$ denote the set $\left\{w_{i}:1\leq i\leq k\right\}$. Let

denote the feasible set of $W$’s. A function $f\in\mathcal{F}$ will depend on $v$ and $W$, and it can be written as $f(x;v,W)$. But when clear from the context, it is shorten as $f(x)$.

Our primary goal is to identify conditions under which there are no spurious local minima. We need to identify a set $\mathcal{G}_{W}\subseteq\mathcal{F}_{W}$ such that when gradient descent outputs a solution $W\in\mathcal{G}_{W}$ with the gradient norm $\left\|\partial L/\partial W\right\|$ smaller than $\epsilon$, then the training and test errors can be bounded by $O(\epsilon^{2})$. Ideally, $\mathcal{G}_{W}$ should have clear characterization that can be easily verified, and should contain most $W$ in the parameter space (especially those solutions obtained in practice).

On notation, we will use $c$, $c^{\prime}$ or $C$, $C^{\prime}$ to denote constants and its value may change from line to line.

In this section, we will rewrite the set of first order conditions for minimizing the empirical loss $L$. This rewriting motivates the direction of our later analysis. More specifically, the gradient of the empirical loss w.r.t. $w_{k}$ is

for all $k=1,\ldots,n$. We will express this collection of gradient equations using matrix notation. Define the “extended feature matrix” as

However, in practice, we will not have the gradient being exactly zero because, e.g., we stop the algorithm in finite steps or because we use stochastic gradient descent (SGD). In other words, typically we only have $\left\|\partial L/\partial W\right\|\leq\epsilon$, and $D$ being full rank is insufficient since small gradient can still lead to large loss. More specifically, let $s_{m}(D)$ be the minimum singular value of $D$, we have

We can see that $s_{m}(D)$ needs to be large enough for the residual to be small. Thus it is important to identify conditions to lower bound $s_{m}(D)$ away from zero, which will be the focus of the paper.

2 Spectrum decay of activation kernel

We will later show that $s_{m}(D)$ is related to the decay rate of the kernel spectrum associated with the activation function. More specifically, for an activation function $\sigma(w^{\top}x)$, we can define the following kernel function

In particular, for ReLU, the kernel has a closed form

In fact, it is a dot-product kernel and its spectrum can be obtained through spherical harmonic decomposition:

where the eigenvalues are ordered $\gamma_{1}\geq\cdots\geq\gamma_{m}\geq\cdots\geq 0$ and the bases $\phi_{u}(x)$ are spherical harmonics. The $m$-th eigenvalue $\gamma_{m}$ will be related to $s_{m}(D)$.

For each spherical harmonic of order $t$, there are N(d,t)=\frac{2t+d-2}{t}\left(\begin{array}[]{c}t+d-3\\ t-1\end{array}\right) basis functions sharing the same eigenvalue. Therefore, the spectrum has a step like shape where each step is of length $N(d,t)$. Especially, for high dimensional input $x$, the number of such basis functions with large eigenvalues can be very large. Figure 2 illustrates the spectrum of the kernel for $d=1500$, and it is about $\Omega(m^{-1})$ for a large range of $m$. For more details about the decomposition, please refer to Appendix A.

Such step like shape also appears in the Gram matrix associated with the kernel. Figure 2 compares the spectra of the kernel of $d=15$ and the corresponding Gram matrix with $m=3000$. We can see the spectrum of the Gram matrix closely resembles that of the kernel. Such concentration phenomenon underlies the reason why the spectrum of $D^{\top}D$ is closely related to the corresponding kernel.

3 Weight discrepancy

Essentially, each $S_{xy}$ defines a slice-shaped area on the sphere which is carved out by the two half spaces $w^{\top}x\geq 0$ and $w^{\top}y\geq 0$.

Based on the collection ${\mathcal{S}}$, we can define two discrepancy measures relevant to ReLU units. $L_{\infty}$ discrepancy of $W$ w.r.t. ${\mathcal{S}}$ is defined as

where the expectation is taken over $x,y$ uniformly on the sphere. We use $L_{\infty}(W)$ and $L_{2}(W)$ as their shorthands. Both discrepancies measure how diverse the points $W$ are. The more diverse the points, the smaller the discrepancy.

For our analysis, we slightly generalize the discrepancy for $w_{k}$’s not on unit sphere, by setting

Main results

Our main result is a bound on the training loss and the generalization error, assuming sufficiently large $n,d$. To state the theorem, first recall that $\beta\in(0,1)$ is the decay exponent of the spectrum of the activation kernel in (15), that is, $\gamma_{m}$ is the $m$-th eigenvalue of the kernel and satisfies $\gamma_{m}=\Omega(m^{-\beta})$. Also recall that $\mathcal{F}_{W}$ denote the set of feasible values of $W$.

then there exists a set $\mathcal{G}_{W}\subseteq\mathcal{F}_{W}$ which takes up $1-\delta^{\prime}$ fraction of measure of $\mathcal{F}_{W}$, such that with probability at least $1-cm^{-\log m}-\delta$ the following holds. For any $W\in\mathcal{G}_{W}$ and any $v\in\left\{-1,1\right\}^{n}$, we have

Note that an immediate corollary is that any critical point in $\mathcal{G}_{W}$ is global optimum.

Furthermore, a randomly sampled set of weights $W$ are likely to fall into this set. This suggests a reason for the practical success of training with random initialization: after initialization, the parameters w.h.p. fall into the set, then would stay inside during training, and finally get to a point with small gradient, which by our analysis, has small error.

Remark 2.

An important feature about our result is that the set $\mathcal{G}_{W}$ has a simple explicit form:

where $c_{g}>0$ is a universal constant. Furthermore, $(L_{2}(W))^{2}$ has a simple closed form (See Theorem 7). Therefore, it is possible to directly check if a solution $W$ is in $\mathcal{G}_{W}$, or design regularization that make sure $W$ stays in the set $\mathcal{G}_{W}$.

Remark 3.

The above theorem is a special case of the following more general result.

and any $v\in\left\{-1,1\right\}^{n}$, we have

Analysis roadmap

Our key technical result is a lower bound on the smallest singular value of $D$ based on the spectrum of the activation kernel defined in (15) and the discrepancy of the weights defined in (20). Once the lower bound is obtained, we can use (13) to bound the training loss, and use Rademachar complexity to bound the generalization error.

It is interesting to compare the theorem to the results in [Soudry and Carmon, 2016], which shows that $D$ is full rank with probability one under small perturbations. However, full-rankness alone is not sufficient since its smallest singular value could be extremely small leading to possibly huge training loss. Instead, we directly bound the smallest singular value and relate it to the activation and the diversity of the weights.

Here we describe the high level intuition for bounding the minimum singular value. It is necessarily connected to the activation function and the diversity of the weights. For example, if $\sigma^{\prime}(t)$ is very small for all $t$, then the smallest singular value is expected to be very small. For the weights, if $d<m$ (the interesting case) and all $w_{k}$’s are the same, then $D$ cannot have rank $m$. If $w_{k}$’s are very similar to each other, then one would expect the smallest singular value to be very small or even zero. Therefore, some notion of diversity of the weights is needed.

where $\left\|G-G_{n}\right\|$ is the spectral norm of the difference.

For the first term in the lower bound, we observe that $G$ has a particular nice form: $G(i,j)=g(x_{i},x_{j})$, the kernel defined in (15). This allows us to apply the eigendecomposition of the kernel and positive definite matrix concentration inequality to bound $\lambda_{m}(G)$, which turns out to be around $m\gamma_{m}/2$.

For the second term, when $w_{k}$’s are indeed from the uniform distribution over the sphere, this can be bounded by concentration bounds. It turns out that when $w_{k}$’s are not too far away from that, it is still possible to do so. Therefore, we use the geometric discrepancy to measure the diversity of the weights, and show that when they are sufficiently diverse, $\left\|G-G_{n}\right\|$ is small. In particular, the entries in $G-G_{n}$ can be viewed as the kernel of some U-statistics, hence concentration bounds can be applied. The expected U-statistics turns out to be the $(L_{2}(W))^{2}$, which has a closed form and can be shown to be small.

Outline.

Theorem 3 is proved in Section 6, $L_{2}(W)$ and $\mathcal{G}_{W}$ are characterized in Section 6.2, and the proof sketch of Theorem 1 and Theorem 2 is provided in Section 7. We describe the proof sketch for the lemmas and provide the remaining proofs in the appendix.

Bounding the smallest singular value

Theorem 3 can be obtained from the following technical lemma.

With probability $\geq 1-m\exp\left(-m\gamma_{m}/8\right)-2m^{2}\exp\left(-4\log^{2}d\right)-\delta$, we have

Lemma 4 is meaningful only when $cn\rho(W)$ is small compared to $nm\gamma_{m}/2$. This requires $L_{2}(W)$ to be sufficiently small. In the following we will first provide the proof sketch of Lemma 4, and then bound that $L_{2}(W)$ in Section 6.2.

To prove Lemma 4, it is sufficient to bound the smallest eigenvalue of $G_{n}=D^{\top}D/n$. Note that $v_{k}\in\left\{-1,1\right\}$, so $v_{k}^{2}=1$, and thus the $(i,j)$-th entry of $G_{n}$ is

For ReLU, $\sigma^{\prime}(w^{\top}x)$ does not depend on the norm of $w$ so without loss of generality, we assume $\left\|w\right\|=1$. Consider a related matrix $G$ whose $(i,j)$-th entry is defined as

Note that $G(i,j)=g(x_{i},x_{j})$ where $g$ is the kernel defined in (15). This allows us to reason about the eigenspectrum of $G$, denoted as $\lambda_{1}(G)\geq\ldots\geq\lambda_{m}(G)$.

Therefore, our strategy is to first bound $\lambda_{m}(G)$ in Lemma 5 and then bound $|\lambda_{m}(G)-\lambda_{m}(G_{n})|$ in Lemma 6. Combining the two immediately leads to Lemma 4.

First, consider $\lambda_{m}(G)$. We consider a truncated version of spherical harmonic decomposition:

and the corresponding matrix $G^{[m]}$. On one hand, it is clear that $\lambda_{m}(G)\geq\lambda_{m}(G^{[m]})$. On the other hand, $G^{[m]}=AA^{\top}$ where $A$ is a random matrix whose rows are

Next, we bound $\lambda_{m}(G^{[m]})$ by matrix Chernoff bound [Tropp, 2012], and it is better than previous work [Braun, 2006]. This leads to the following lemma.

With probability at least $1-m\exp\left(-m\gamma_{m}/8\right)$,

Next, bound $|\lambda_{m}(G)-\lambda_{m}(G_{n})|$. By Weyl’s theorem, this is bounded by $\left\|G-G_{n}\right\|$. To simplify the notation, denote

Then $G(i,j)-G_{n}(i,j)=\left\langle x_{i},x_{j}\right\rangle E_{ij}$, and thus

where the last inequality holds with high probability since $x_{i}$’s are uniform over the unit sphere and thus we can apply sub-gaussian concentration bounds.

Note that $\sum_{i\neq j}E_{ij}^{2}/(m(m-1))$ is a U-statistics where the summands are dependent and typical concentration inequality for i.i.d. entries does not apply. Instead we use a Bernstein inequality for U-statistics [Peel et al., 2010] to show that with probability at least $1-\delta$, it is bounded by

The key observation is that the quantities in the above lemma are related to discrepancy:

Plugging (31)-(33) into (30) and (29), we have

The following inequality holds with probability at least $1-2m^{2}\exp\left(-\log^{2}d\right)-\delta$,

where $\rho(W)$ is as defined in Lemma 4.

2 Characterizing the discrepancy

In this subsection, we present a bound for $L_{2}(W)$ and show that the $\mathcal{G}_{W}$ defined in the following covers most $W$’s. Recall that

for $0<\delta^{\prime}<1$ and a proper constant $c_{g}>0$. The constant $c_{g}$ is the constant in Lemma 8. $\delta^{\prime}$ will be clear from the context where $\mathcal{G}_{W}$ is used.

First we provide a closed form for $L_{2}$ discrepancy of slices defined in (18). The proof is provided in the appendix.

The closed form is simple and intuitive. The kernel $k(w_{i},w_{j})$ measures how similar two units are. The discrepancy is the difference between the average pairwise similarity and the expected one over uniform distribution.

There exists a constant $c_{g}$, such that for any $0<\delta^{\prime}<1$, with probability at least $1-\delta^{\prime}$ over $W=\left\{w_{i}\right\}_{i=1}^{n}$ that are sampled from the unit sphere uniformly at random,

Final bound on generalization error

Here we provide the proof sketch of Theorem 2 and Theorem 1. More details of the proof are in Appendix F.

First, we prove Theorem 2. Suppose a solution $W$ satisfies the assumption and has small gradient $\left\|\partial L/\partial W\right\|$. Using (13), we have $\left\|r\right\|\leq{\left\|\partial L/\partial W\right\|}/{s_{m}(D)}.$ By Theorem 3 and the assumption in Theorem 2, with high probability $s_{m}^{2}(D)=\Omega(nm^{1-\beta})$. This implies the training loss is

The generalization error can then be derived using McDiamid’s inequality and Rademacher complexity. First, we need an upper bound on the difference of the loss for two data points for the McDiamid’s inequality. Since $\left\|x\right\|_{2}\leq 1$ and $\sum_{k}\left\|w_{k}\right\|_{2}\leq C_{W}$, we have $\left|f\right|\leq C_{W}$. Thus

where in the last inequality we use the fact that the true function $|y|\leq Y$. Next, we use the composition rules to compute the Rademacher complexity. Since the complexity of linear functions $\left\{w^{\top}x:\left\|w\right\|_{2}\leq b_{W},\left\|x\right\|_{2}\leq 1\right\}$ is $b_{W}/\sqrt{m}$ and $\sigma(\cdot)$ is 1-Lipschitz, and $\sum_{k}\left\|w_{k}\right\|_{2}\leq C_{W}$, the complexity $\mathcal{R}_{m}(\mathcal{F})\leq C_{W}/\sqrt{m}.$ Composing it with the loss function, and applying the bound in [Bartlett and Mendelson, 2002], we get the final generalization bound.

Discussions

In this section, we discuss and remark on further considerations and possible extensions of our current analysis.

2 Other activation functions

We can consider a family of activation functions of the form $\sigma(u)=\max\left\{u,0\right\}^{t}$, i.e., rectified polynomials [Cho and Saul, 2009, Krotov and Hopfield, 2016]. This requires two modifications to the analysis.

Examples for the first few $t$ are listed as follows.

Larger $t$ corresponds to more nonlinear activation functions and leads to slower decaying spectrum since there are more high frequency components.

We also need to change the definition of the discrepancy to accommodate the new kernels. Let

Therefore, the discrepancy is affected by how the kernels change due to change in activation functions.

The other modification is on the Rademacher complexity. Since the derivative $\sigma^{\prime}(u)=t\max\left\{u,0\right\}^{t-1}$, there is an additional factor of $t$ in front of the complexity. That is, larger $t$ leads to higher Rademacher complexity.

In summary, the best parameter $t$ depends on the balance between the two conflicting effects. On one hand, larger $t$ corresponds to slower decaying spectrum and makes the minimum singular value more likely to be larger. On the other hand, smaller $t$ leads to better generalization since the Rademacher complexity is smaller.

3 (Sub)gradient of the activation function

In summary, though for some $W\in\mathcal{G}_{W}$ the loss is not differentiable, one can define $\partial L/\partial W$ by using subgradients of ReLU $\sigma$ as follows:

for any $c\in$. Then under the conditions in our theorems, with high probability, for any $W\in\mathcal{G}_{W}$ and any definition of $\sigma^{\prime}$ in (41), the guarantees hold.

Other activation functions such as rectified polynomials are differentiable and thus they do not have such issue.

4 Other input distribution

When the input distribution is not uniform, the spectrum of the kernel function defined in (14) will be different because the spherical harmonic bases are defined with respect to the input distribution. To ensure the spectrum decays slowly, we need to find a corresponding distribution of $W$ that “matches” the input distribution.

In Table 1, we compare the minimum eigenvalues with the two distributions. The uniform distribution on $F_{E}$ always leads to larger or the same minimum eigenvalues. However, as dimension increases, the difference becomes negligible. Note that the difference between the uniform distribution on the whole sphere and uniform on $F_{E}$ becomes exponentially small when the dimension $d$ increases, because the proportion of $E$ and $-E$ shrinks exponentially. This suggests that in high dimensions, uniform on the whole unit sphere is a reasonable distribution for $W$.

For a general input distribution $P(x)$, we can decompose it into small sets $dx$ and on every set, the distribution is uniform with measure $P(x)dx$. Then every small sets corresponds to a distribution of $W$. The final distribution of $W$ is the superposition of all such distributions, weighted by $P(x)dx$.

Numerical evaluation

In this section, we further investigate numerically the effects of gradient descent on the discrepancy and the effects of regularizing the weights using discrepancy measure.

One limitation of the analysis is that we have not analyzed how to obtain a solution $W\in\mathcal{G}_{W}$ with small gradient. The theoretical analysis of gradient descent is left for future work. Meanwhile we provide some numerical results supporting our claims.

Although the set $\mathcal{G}_{W}$ contains most $W$’s, it is still unclear whether the solutions given by gradient descent lie in the set. We design experiments to investigate this issue. The ground truth input data are of dimension $d=50$ and true function consists of $n=50$ units. We use networks of different $n$ to learn the true function and perform SGD with batch size 100 and learning rate 0.1 for 5000 iterations. Figure 4 shows how $\left(L_{2}(W)\right)^{2}$ changes as a function of $n$. It is slightly worse than $O(n^{-1})$ but scales better than $O(n^{-1/2})$, suggesting (stochastic) gradient descent outputs solutions with reasonable discrepancy.

2 Regularization

To reinforce solutions with small discrepancy, we propose a novel regularization term to minimize $L_{2}$ discrepancy:

It is essentially $L_{2}$ discrepancy without the constants.

To verify the effectiveness of the regularization term, we explore the relationship between the regularization and the minimum singular value. We first generate 20 random $W$’s, all with $n=100$ and $d=100$, and compute their discrepancy and singular values using $m=3000$. Then we optimize $R(W)$ and compare the quantities after optimization. The result is presented in Figure 4. We can see smaller regularization value corresponds to larger singular value.

We also conducts experiments to compare training and test errors with and without regularization. The ground truth data are of $d=100$ and $n=100$. We learn the true function by SGD with learning rate 0.1, momentum 0.9 and a total of 300,000 iterations. The regularization coefficients are chosen from $\left\{1,0.1,0.01,0.001\right\}$ and the best results are reported. We use neural networks of size $n\in\left\{100,150,200,300\right\}$ and for each $n$ we repeat five times with different random seeds. The result is summarized in Table 2. Regularization leads to lower training and test errors for most settings. Even in the case where the un-regularized one performs better, the errors are all small enough (within the same range as standard deviation), suggesting the noise begins to dominate.

We also compare the regularization effects on the MNIST dataset. The dataset contains 60,000 training and 10,000 test handwritten digits. To demonstrate the regularization effect, we train one hidden layer fully connected neural networks with $k=200,400,600,800$ units. The results are summarized in Table 3. Note that state-of-the-arts performance on MNIST are mostly obtained by convolutional neural networks. This experiment is not intended to achieve the state-of-the-arts but it tries to showcase the advantage of regularization on a real-world dataset.

From Table 3, we see regularization consistently leads to slightly better test error for all cases.

Conclusion

We have analyzed one-hidden-layer neural networks and identified novel conditions when local optima become global optima despite the non-convexity of the loss function. The key factors are the spectrum of the kernel associated with the activation function and the diversity of the units measured by discrepancy.

Although we focus on a least-square loss function and uniform input distribution, the analysis technique can be readily extended to other loss function and input distributions. At the moment, our analysis is still limited in the sense that it is independent of the actual algorithm. In the future work, we will explore the interplay between the discrepancy and gradient descent. In addition, we will further investigate the issue of designing an algorithm that guarantees good discrepancy thus small errors, possibly in a way similar to [Ge et al., 2016] in low-rank recovery problems.

Acknowledge

We thank Santosh Vempala, Lorenzo Rosasco and Jason Lee for valuable discussions. The research is supported in part by NSF/NIH BIGDATA 1R01GM108341, ONRN00014-15-1-2340, NSF IIS-1639792, NSF IIS-1218749, NSF CAREER IIS-1350983, Intel and NVIDIA, and by NSF grants CCF-0832797, CCF-1117309, CCF-1302518, DMS-1317308, Simons Investigator Award, and Simons Collaboration Grant.

References

Appendix A Spherical harmonic decomposition and kernel spectrum

Any function defined on the unit sphere has a spherical harmonic decomposition

For each order $t$, there are N(d,t)=\frac{2t+d-2}{t}\left(\begin{array}[]{c}t+d-3\\ t-1\end{array}\right) bases with the same coefficient. As a result, the spectrum $\gamma_{u}$ sorted by magnitude has the step like shape where each step is of length $N(d,t)$.

To compute the coefficients, we use the Legendre harmonics [Müller, 2012] with the following property

The spherical harmonics also form an orthonormal basis on the unit sphere:

Combining these properties, we can calculate the spectrum using

Appendix B Bounding λm​(G)\lambda_{m}(G) using matrix concentration bound: Proof of Lemma 5

For an integer $r>0$, define the truncated version of $g$ and the corresponding residue as

Let $c_{g}=\max_{x}g(x,x)$ then with probability at least $1-m\exp\left(-\frac{m\gamma_{m}}{8c_{g}}\right)$,

Therefore, matrix Chernoff bound (e.g., [Tropp, 2012]) gives

Choose $\epsilon=1/2$ and use the facts that $G^{[m]}=AA^{\top}$, $Y=A^{\top}A$ and $\lambda_{m}(G^{[m]})=\lambda_{m}(Y)$, we finish the proof. ∎

By Weyl’s theorem and the fact that $E^{m}$ is PSD,

Appendix C Bounding the difference between λm​(G)\lambda_{m}(G) and λm​(Gn)\lambda_{m}(G_{n}): Proof of Lemma 6

We are going to give an upper bound on $\left\|G_{n}-G\right\|$:

Our bound heavily relies on the inner products $\left|\left\langle x_{i},x_{j}\right\rangle\right|$ for all $i\neq j$ being small enough. In the next lemma, we provide such a result for uniformly distributed data.

Note that both $a$ and $b$ are sub-gaussian random variables with sub-gaussian norm $c/\sqrt{d}$ where $c$ is some constant [Vershynin, 2010].

The last inequality uses the independence of $a$ and $b$ and $\left\|\left\langle a,b\right\rangle\right\|_{\psi_{2}}\leq\left\|b\right\|_{2}\left\|a\right\|_{\psi_{2}}$ for a fixed $b$. ∎

Decomposing the sum into diagonal and off-diagonal terms gives us

Let $\mathcal{G}$ denote the event that for all $i\neq j\in[m]$, $|\left\langle x_{i},x_{j}\right\rangle|\leq O\left(\frac{\log d}{\sqrt{d}}\right)$, then by Lemma 10 and the union bound

Therefore, with probability at least $1-2m^{2}e^{-\log^{2}d}$, we have

Suppose $\left|E_{ij}\right|\leq B$, according to the concentration inequality (Theorem 2 in [Peel et al., 2010]), we have with probability at least $1-\delta$

Putting everything together, we have with probability at least $1-\delta-2m^{2}e^{-\log^{2}d}$

Appendix D Discrepancy of the weights

The $L_{\infty}$ discrepancy of $W$ with respect to ${\mathcal{S}}$ is

where the expectation is taken over $x,y$ uniformly on the sphere. We use $L_{\infty}(W)$ and $L_{2}(W)$ as their shorthands.

using the fact that $E_{ij}=\text{dsp}(W,S_{x_{i}x_{j}})$.

In the following subsections, we will discuss the discrepancies.

Consider the first term, which is equal to

where the third step is by invariance to $p$ and the fourth step is by Lemma 11. The theorem then follows. ∎

Theorem 7 lets us compute $L_{2}(W)$ for a fixed $W$. The next lemma gives a concrete bound for a special case where $W$ is uniformly distributed on the unit sphere.

Lemma 8 There exists a constant $c_{g}$, such that for any $0<\delta<1$, with probability at least $1-\delta$ over $W=\left\{w_{i}\right\}_{i=1}^{n}$ that are sampled from the unit sphere uniformly at random,

Then let $\mathcal{G}$ denote the event that $|x|=|u^{\top}v|\leq c\sqrt{\log d/d}$ for a sufficient large constant $c>0$, so that by Lemma 10, $\Pr[\neg\mathcal{G}]\leq O(1/d^{4})$. Then

Then by Berstein’s inequality, we have with probability at least $1-\delta$ over the $W$ uniformly on the sphere,

A similar argument holds for $T_{2}=\frac{1}{n^{2}}\sum_{i,j=1}^{n}\left(\frac{1}{2}-d(w_{i},w_{j})\right)$. Note that

We have that with probability at least $1-\delta$ over the $W$ uniform from the sphere,

Below are some technical lemmas that are used in the analysis.

The first two are straightforward. The third is implicit in the proof of Theorem 1.21 in [Bilyk and Lacey, 2015]. ∎

Appendix E The spectrum of γm\gamma_{m}

The spectrum of the kernel matrix $g(x,y)=\left(\frac{1}{2}-\frac{\arccos\left\langle x,y\right\rangle}{2\pi}\right)\left\langle x,y\right\rangle$ is determined by the spherical decomposition coefficients.

We need $\gamma_{m}$ to decrease slower than $O(1/\sqrt{m})$ within a reasonable range, such as $m\leq 1000000$.

Although the kernel associated with ReLU decreases faster than the desired rate, we can choose from a family of such arccos kernels such that the spectrum decays slower than $1/\sqrt{m}$.

Larger $n$ corresponds to more nonlinear activation functions and leads to slower decaying spectrum.

Higher orders of $J_{n}(\theta)$ seems to be extremely complicated.

Although there is no analytical solution to the spectrum, we can compute them numerically.

Figure 2 illustrates the spectra of several arccos kernels compared to $O(1/m)$ and $O(1/\sqrt{m})$.

Appendix F Rademacher complexity and final error bounds: Proof of Theorem 2 and Theorem 1

We apply the argument in [Bartlett and Mendelson, 2002] to our setting to get Lemma 12. Combining it with Theorem 3 leads to Theorem 2. Further combining it with Lemma 8 leads to Theorem 1.

Suppose the data are bounded: $\left|y\right|\leq Y$ and $\left\|x\right\|_{2}\leq 1$. Let

Then with probability $\geq 1-\delta$, for any $f\in\mathcal{F}$,

where $\mathcal{L}$ is the set of loss functions

Let $S$ and $S^{\prime}$ be two datasets that differ by exactly one data point $(x_{i},y_{i})$ and $(x^{\prime}_{i},y^{\prime}_{i})$. Then we have a bound on the difference of loss functions. Since $\left\|x\right\|_{2}\leq 1$ and $\sum_{k}\left\|w_{k}\right\|_{2}\leq C_{W}$, we have $\left|f\right|\leq C_{W}$. Thus

Similarly, we can get the other side of the inequality and have $\left|\Phi(S)-\Phi(S^{\prime})\right|\leq\frac{Y^{2}+C_{W}^{2}}{m}$.

From McDiamids’ inequality, with probability at least $1-\delta$ we get

The first term on the right-hand side can be bounded by Rademacher complexity as shown in the book Foundations of Machine Learning (3.13). In the end, we have the bound

where $\mathcal{R}_{m}(\mathcal{L})$ is the Rademacher complexity of the function class $\mathcal{L}$.

We can find the Rademacher complexity by using composition rules. The Rademacher complexity of linear functions $\left\{w^{\top}x:\left\|w\right\|_{2}\leq b_{W},\left\|x\right\|_{2}\leq 1\right\}$ is $b_{W}/\sqrt{m}$, where $m$ is the number of data points. If a function $\phi$ is $L$-Lipschitz, then for any function class $\mathcal{H}$, we have $\mathcal{R}(\phi\circ\mathcal{H})\leq L\mathcal{R}(\mathcal{H})$. In addition, we also have $\mathcal{R}(cH)=\left|c\right|\mathcal{R}(H)$ and $\mathcal{R}(\sum_{k}F_{k})\leq\sum_{k}\mathcal{R}(F_{k})$.

So for the function class $\mathcal{F}$ that describes a neural network, we have

It is derived by using the fact that $\sigma^{\prime}(\cdot)$ is 1-Lipschitz and $\sum_{k}\left\|w_{k}\right\|_{2}\leq C_{W}$.

Finally composing on the loss function we get

using the fact that the ground truth in the loss should be bounded by $Y$ and the function bounded by $C_{W}$, thus the Lipschitz constant of the loss function is bounded by $Y+C_{W}$. ∎