{'fraction_non_alphanumeric': 0.04664127741580251, 'fraction_numerical': 0.01681655501514178, 'mean_word_length': 4.110681205299566, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Over-parametrized neural networks as under-determined linear systems

We draw connections between simple neural networks and under-determined linear systems to comprehensively explore several interesting theoretical questions in the study of neural networks. First, we emphatically show that it is unsurprising such networks can achieve zero training loss. More specifically, we provide lower bounds on the width of a single hidden layer neural network such that only training the last linear layer suffices to reach zero training loss. Our lower bounds grow more slowly with data set size than existing work that trains the hidden layer weights. Second, we show that kernels typically associated with the ReLU activation function have fundamental flaws -- there are simple data sets where it is impossible for widely studied bias-free models to achieve zero training loss irrespective of how the parameters are chosen or trained. Lastly, our analysis of gradient descent clearly illustrates how spectral properties of certain matrices impact both the early iteration and long-term training behavior. We propose new activation functions that avoid the pitfalls of ReLU in that they admit zero training loss solutions for any set of distinct data points and experimentally exhibit favorable spectral properties.

1. Introduction. Neural networks are among the predominant mathematical models in machine learning, though theoreticians still struggle to understand their efficacy fully. While the origins of neural networks go back decades, they have recently benefited from modern computing architectures, the growth of labeled data sets, and their breadth of applicability. In this work, we use a simple mathematical model to understand the properties of "over-parameterized" neural networks (roughly, those with more model parameters than are seemingly necessary) through the lens of underdetermined systems of equations in numerical analysis.
Motivated by empirical observations [58], there has been significant work on showing that neural networks can achieve zero training loss when the weights are trained using a simple optimization algorithm such as (stochastic) gradient descent [1,16,17,36]. However, such work depends delicately on analyzing networks in what is essentially a "linearized regime" [11,19,28]. In this regime, the training dynamics mostly follow a kernel method [26] called the neural tangent kernel (NTK). Surprisingly, the property of achieving zero training error is effectively independent of training the weights of the network and the theoretical analysis relies on the fact that the hidden layer weights hardly change from their random initialization.
In this paper, we show that just training the last layer is sufficient to achieve zero training loss in most cases. The lower bounds we find on the network width to acquire this property are significantly smaller than those previously developed for training the hidden layer weights. More specifically, for a simple two-layer neural network, we provide conditions for the existence of zero training loss solutions and gradient descent's convergence to it when only the last layer weights are trained. Notably, our results hold for randomly selected weights and are distinct from results on finite sample expressiveness [58,Sect. 4], where hidden layer weights are initialized in a data-dependent manner. When a network's weights are chosen at random, training the last layer is equivalent to solving an under-determined linear system, which is a well-understood problem in numerical analysis. The linear system corresponds to the so-called random feature regime [39,40,41], except with feature choices motivated by the structure of a single layer fully connected neural network (see Section 2). The under-determined linear systems viewpoint turns out to be reasonably flexible and readily applicable to the use of "pre-trained" models.
To achieve our theoretical results, we use matrix concentration inequalities to connect the under-determined linear system to a square system in the infinite width limit. The study of the corresponding infinite width limit is used in classical analysis of neural networks via ridge function approximation theory [6,14,29,37] and is closely related to kernel methods [34,54]. We can show that in the infinite width limit, one can achieve zero training loss by showing that the kernel induced by certain activation functions and random weights is strictly positive definite.
In Section 3, we rigorously characterize when an activation function corresponds to a positive definite kernel. A critical insight is that the commonly used rectified linear unit (ReLU) activation function is not strictly positive definite for a bias-free network. We provide a set of eight data points that demonstrates this. Furthermore, The same eight data points explicitly show that the NTK is also only positive semidefinite when using ReLU. Importantly, our results hold for any finite width and the infinite width limit. In the finite width setting, our results are independent of the choice of hidden layer weights and conclude that there are benign datasets for which a bias-free single hidden layer neural network with ReLU cannot generically achieve zero training loss.
Motivated by these observations, in Section 3 we introduce a new set of activation functions from the radial basis function literature that yield positive definite kernels in the infinite limit and experimentally produce well-conditioned finite width models. Our analysis is closely related to recent connections between random features models and kernel methods [4,5,15], spectral properties of neural networks [10,38,43], and generalization performance of kernel methods [9]. Our new activation functions complement these results by experimentally exhibiting favorable spectral properties. Furthermore, we make concrete connections between the spectral properties of finite systems and the training process in Section 4.
More specifically, in Section 4, we show that achieving zero training error via gradient descent applied to the last layer is relatively simple to accomplish. Thus, achieving zero training loss with over-parameterized networks and simple optimization methods is not surprising. If gradient descent is initialized appropriately, the minimal norm solution is found. This fact alone is not particularly useful as it is easy to characterize the minimal norm solution of a consistent under-determined linear system. More interesting are the connections we develop between training the last layer and the so-called Landweber iteration [27]. Our analysis shows the connection between spectral properties of the kernel and optimization performance, thereby motivating the choice of activation functions that correspond to positive definite kernels and have favorable spectral properties. The suggested activation functions have such properties, and simple numerical experiments show their efficacy both in early iteration performance and ability to achieve zero training error via gradient descent.
Our results show that there is effectively no gap between what is achievable by random features models and fully trained neural networks in terms of training loss. This does not mean the two models are equivalent more generally -fully trained neural networks are more complex models and nominally more capable. In fact, recent research explores the theoretical distinctions between the random features and NTK regimes [19,20]. Further theoretical insight on the random features regime and the models properties can be drawn from its connections to "ridgeless" regression and interpolation [7,8,23]. However, many of these results are asymptotic (in data dimension, network width, and data set size) and make statistical assumptions on the data. In this work, we explicitly focus on finite problem sizes and deliberately make minimal assumptions on the input data. Lastly, while there is extensive work on the generalization properties of neural networks, even in the random features regime [31,44]), we explicitly omit such a discussion here.
2. An over-parameterized two-layer neural network. Suppose that we are given n distinct data points x 1 , . . . , x n ∈ R d that are normalized so that x i 2 = 1 for 1 ≤ i ≤ n, 1 where each data point is assigned a label y 1 , . . . , y n ∈ R. We let X ∈ R d×n be the data matrix, where X(:, i) = x i and y = y 1 · · · y n is the vector of labels.
We consider the task of training a two-layer neural network, G : R d → R, with a single fully-connected hidden layer of width m and linear last layer such that G(x i ) ≈ y i for 1 ≤ i ≤ n. More specifically, for a fixed continuous scalar-valued activation function γ : R → R and integer m, G takes the form 2 where W ∈ R d×m is the weight matrix, α ∈ R m is the last layer, and γ(x T W ) is the 1 × m vector obtained by applying γ to x T W ∈ R 1×m entrywise. Throughout this paper, we assume the activation function γ is Lipschitz continuous with Lipschitz constant L, and its magnitude is bounded by C on [−1, 1], i.e., |γ(z)| ≤ C for z ∈ [−1, 1]. For many common activation functions it suffices to take C = 1.
One usually fits G by jointly learning W and α via the least-squares problem It is standard to try and solve (2.1) by gradient descent or stochastic gradient descent. There are many popular choices for the activation function γ, and we will analyze specific choices later. The choice of γ affects how computationally expensive it is to solve (2.1) and the quality of the final discovered G. When m > n, we consider the model to be over parametrized since there are more hidden nodes than training data points. When there are more parameters than data points in the neural network, it is quite plausible to find an interpolating G so that G(x i ) = y i for 1 ≤ i ≤ n. In other words, we may expect the existence of W and α such that Surprisingly, we find that some of the most popular choices of γ do not allow one to find an interpolating G; regardless of how over-parameterized the network is. Concretely, we show this for the commonly used ReLU activation function, which is defined as γ(z) = (z) + , where (z) + = z if z > 0 and 0 for z ≤ 0. Theorem 2.2 shows that for this model ReLU does not admit interpolating solutions when using the well-separated data points in Definition 2.1. where 0 ∈ R d−3 is the zero vector and s = 1/ √ 3.
Therefore, E T v = 0 and E is of rank at most 7.
The implications of Theorem 2.2 are quite significant, particularly when we generalize the result to properties of the induced kernels in Section 3, and are summarized in Corollary 2.3. In particular, there are "nice" data sets X containing distinct and well-separated points for which (2.1) may not exactly fit all possible labels y when the ReLU activation function is used. More precisely, it shows that for ReLU, γ(X T R W ) is not full row-rank regardless of the width m and choice of weights W. Consequently, there are sets of labels y outside the range of γ(X T R W ), which precludes the possibility of generically achieving zero training loss regardless of how W and α are trained. Practically, Remark 1 shows how more commonly used models may avoid this theoretical pitfall. Nevertheless, the existence of such data sets may help shed light on what makes problems difficult for practical single hidden layer neural networks.
Importantly, planting (2.3) into any other data set makes solving (2.1) challenging as γ(X T W ) remains rank-deficient. An additional implication of having these "singular" point sets is that for any n ≥ 8, we can make the minimal singular value of γ(X T W ) arbitrarily small while enforcing that X has no parallel points. This observation has consequences on the width of networks shown to achieve zero training loss in prior work [16,17].
and (X, y) be any dataset that contains the points in Definition 2.1. Then, there exists a non-zero vector q ∈ R n such that if y T q = 0 no bias-free single hidden layer neural network using ReLU activation, i.e., γ(z) = (z) + in (2.1) can achieve zero error for (2.2).
Proof. Without a loss of generality, assume that X is ordered such that X = X R X C and y = y R y C , where y R ∈ R 8 are the labels associated with X R from T as in the proof of Theorem 2.2.
We will show that the fixed vector q = v 0 T satisfies the criteria of the Corollary; assume that y T q = 0, which implies that y T R v = 0. Now, suppose that there are weights W , and α that achieved zero error for (2.2). This implies that m −1/2 γ(X T R W )α = y R and, therefore, m −1/2 v T γ(X T R W )α = v T y R . By Theorem 2.2, we find that v T γ(X T R W ) = 0, so y T R v = 0 leads to a contradiction. Remark 1. While we formally allow for the inclusion of bias terms through augmentation of the data, this technique formally constrains the set of allowable X. These constrains may suffice to ensure that (X T W ) + is full row-rank for some W provided the columns of X are distinct. Consequently, the bias terms are theoretically quite important for a ReLU network to be generically able to achieve zero training loss -this agrees with the universal approximation theorem [29] and finite width expressiveness results [58,Sect. 4].
2.1. The random features regime. In light of (2.2), when is it reasonable to expect an interpolating solution to exist? And, if it does, then can gradient descent find it? To answer this question, we study (2.1) in the so-called random features regime [39], where W is arbitrary and kept fixed while one optimizes over α. In this regime, provided that m ≥ n, (2.1) reduces to an under-determined least-squares problem 3 given by Characterizing solutions to (2.4) requires understanding the properties of the matrix γ(X T W ) for a fixed W . For example, if γ(X T W ) has full row-rank, then the minimal norm solution to (2.4) is given by α = 1 √ m γ(W T X)z, where z solves the non-singular linear system This means that it is relatively easy to obtain zero training loss provided that γ(X T W ) has full row-rank or, equivalently, the matrix in (2.5) is non-singular. Mirroring standard initialization techniques, we study (2.4) when W is a random matrix with each column is drawn independently from the uniform distribution over the sphere, i.e., w i ∼ U(S d−1 ). Denoting the matrix in (2.5) as H, it is given by We would like to ensure that H is invertible with high probability when m is sufficiently large and accomplish this by connecting H to a matrix that characterizes the infinite width limit of our simple neural network.
2.1.1. Infinite width limit. For a given γ, the entries of the population-level matrix H ∈ R n×n are given by Since H is a Gram matrix, we know that H is symmetric and H 0 for all continuous γ. Notably, the matrix H appears when analyzing the infinite width limit of neural networks and their relation to Gaussian processes [34,54]. Throughout this paper, we often write E w instead of E w∼U (S d−1 ) to keep the notation compact. As we will show, if H 0 it follows that γ(X T W ) is full row-rank for sufficiently large m. This implies that (2.4) has a unique minimum norm solution with a zero objective value. The specific choice of γ can impact the properties of H and how they relate to training (see Section 4). In Section 3, we chose a specific γ motivated by bounding λ min (H) away from zero.
Remark 2. While it is tempting to try to find a connection between the matrix H in (2.7) and a kernel interpolation matrix, i.e., is a positive definite reproducing kernel, it is typically not possible. To see this, note that the uniform distribution prescribes the inner product used to construct entries of H, and it is often not in agreement with the inner product needed to produce a proper kernel matrix [5,25]. The more appropriate analogy is to that of Random Fourier Features [39] and Random Kitchen Sinks [41]. First, we control λ min ( H) with high probability by using a matrix Chernoff bound.
are independent and identically distributed positive semi-definite matrices. Furthermore, they satisfy E w γ(X T w)γ(w T X) = H, and we have that where C is the upper bound on |γ| from Section 2. Based on these observations, the statement follows from a matrix Chernoff bound [49, Thm. 5.1.1] combined with a union bound to simultaneously control the probability bounds on both the minimal and maximal eigenvalue.
Using a matrix Bernstein inequality, we can understand how well H approximates H with high probability in the spectral norm.
Here, κ(H) denotes the condition number of H.
Proof. First, note that Since the weights w i ∈ S d−1 are independently drawn from U(S d−1 ), we may invoke a matrix Bernstein inequality (see [ for any > 0. The result follows when m satisfies (2.8) by setting = λ min (H)/2 and using the fact that H 2 /λ min (H) = κ(H).
It is useful to inspect the lower bounds on m in Lemmas 2.4 and 2.5. First, since λ min (H) ≤ e T i He i ≤ C 2 , the lower bounds on m grow like (n log n)/λ min (H). For fixed data points, this means that optimizing the lower bounds on m is purely a question of choosing γ such that λ min (H) is not too small, and H is well-conditioned. Moreover, the appearance of C 2 /λ min (H) is natural as the definition of H includes two copies of γ. If we scale γ by some constant p, so that C becomes pC, then we also scale λ min (H) by p 2 . Therefore, our bounds are properly invariant with respect to scalings of γ.
Remark 3. While we have restricted our discussion here to w ∼ U(S d−1 ), analogous results are possible to derive for w ∼ N (0, σ 2 I) by using sub-Gaussian versions of the matrix Bernstein inequality (see Appendix A.1).
3. Activation functions and their properties. Lemmas 2.4 and 2.5 show that for sufficiently large m the stochastic estimator H is a good approximation to the population matrix H, and many of the properties of H are reflected in the simple model (2.2). This makes the choice of γ crucial in both the random features regime and the lazy training regime [11]. First, we show that the popular ReLU activation function can exhibit catastrophically bad behavior even for relatively benign data sets and we extend that analysis to the NTK. This result also extends to the swish activation function (see Appendix A.2). Second, we characterize the activation functions that led to a strictly positive definite H for all data sets.
3.1. Failure of the ReLU activation function. In Section 2, our simple model shows that the ReLU activation function has some severe shortcomings when used in (2.1) (see Theorem 2.2). We now show that the 8 × 8 matrix H in (2.7) is also singular for the data points given in Definition 2.1. 4 Proof. First, we note that Now, let x 1 , . . . , x 8 be the point set given in (2.3). Let X ∈ R 8×d be the matrix such that the kth row of X is x k for 1 ≤ k ≤ 8. The matrix of inner-products is given by T is an eigenvector corresponding to an eigenvalue of 0. Note that This means that Hv = 0 and the matrix H is singular.
3.2. Failure of the ReLU Neural Tangent Kernel. In situations where the weights are randomly drawn and then trained, the NTK has gained popularity [26] and it is particularly relevant in regimes where training leads to small changes in the weights [11,17,28]. We find that the NTK for ReLU has the same severe shortcomings as the activation function itself. When the weights are initialized on S d−1 , the matrix of interest for the NTK is 1 z>0 denotes the indicator function for z > 0 (i.e., 1 z>0 = 1 if z > 0 and 0 otherwise). Our key observation is that the point set in Definition 2.1 also causes G to be singular.
Since the point set in Definition 2.1 causes H, H, G, and G to be singular and Theorem 3.1 holds for any W , one cannot overcome this catastrophic behavior by jointly training W and α in (2.1). Moreover, this shows that without placing additional assumptions on the data or modifying the simple model it is impossible to generically assume H and G are non-singular. 7 Moreover, prior results [32,33] adapted to this setting suggests there are many point sets that share the properties of those in Definition 2.1 when the kernel function does not have both an even and odd part that is non-polynomial. Lastly, For H and G these results are related to observations in [43], where the relative effectiveness of learning distinct components of an underlying function is studied.
Remark 4. While we have followed the closely related literature (particularly that on achieving zero training loss) by omitting bias terms explicitly, as alluded to in Remark 1 the addition of bias terms effectively constrain x i to certain parts of the sphere. It is possible the minimal eigenvalues of H and G are bounded away from zero when these additional restrictions are placed on x i , making the inclusion of bias terms vital theoretically. Notably, the specific kernels we analyze are extensively studied (see, e.g., [2] and [28, Appendix C]) and therefore these results have important consequences for studying neural networks via approximation properties of the kernel methods induced by their infinite width limit.
3.3. General activation functions and relations to kernels. While Subsections 3.1 and 3.2 demonstrate fundamental problems with ReLU, we would like to find sufficient conditions on γ to ensure that H in (2.7) is strictly positive definite for all data sets containing distinct points. Since the distribution in the expectation in (2.7) is U(S d−1 ), which is a rotationally invariant distribution, we know that for any continuous γ there exists a function φ : We can use the reproducing kernel literature on S d−1 to give a sufficient condition [45,25].
Theorem 3.3. For d ≥ 3, the matrix H is strictly positive definite for any distinct data points x 1 , . . . , such that Z γ contains infinitely many odd integers as well as infinitely many even integers. Here, C (τ ) k is the ultraspherical polynomial of degree k with parameter τ > 0.
Proof. To ensure that H is strictly positive definite for any (t), then letting Γ(·) denote the Gamma function, we have the following absolutely convergent expansion for φ: From [57], we find that φ(x, y) = φ(x T y) is a strictly positive kernel on the sphere if and only ifZ φ = {k ∈ Z + : c k > 0} contains infinitely many odd integers and infinitely many even integers. This completes the proof asZ φ contains infinitely many odd integers and infinitely many even integers if Z γ = {k ∈ Z + : a k = 0} does.

Wendland kernels.
Motivated by our observations of the ReLU activation function and analysis in Subsections 2.1 and 2.2, we seek a choice of activation function γ that "optimizes" two quantities: (1) The minimal eigenvalue of H should be as large as possible, and (2) The condition number of H should be as small as possible. The first point has immediate implications for the existence of solutions with zero training loss, while the second impacts both the actual training process and allows us to invoke classical results for robustness and stability. Towards this end, we advocate for the non-linearities derived from so-called Wendland kernels, which are compactly supported radial basis functions [52,53].

The Kernels.
While there are infinitely many classes of Wendland kernels, we restrict our attention to the least smooth variants [53, Tab. 9.1]. Written for data on the sphere the two we consider are given by and the corresponding kernels are Φ d,k (x, y) = φ d,k 2 − 2x T y . Notably, Φ d,k is a strictly positive definite kernel on S d−1 [53]. To modulate the width of the kernels we introduce the parameter ζ > 0 and define Φ d,k,ζ (x, y) = φ d,k ( 2 − 2x T y/ζ). While all positive definite kernels must have eigenvalues that decay to zero as the number of data points grows, Wendland kernels have significantly slower decay than other common choices. This property is the impetus for our choice as we may expect that use of a feature map motivated by the Wendland kernels will lead to a well-conditioned γ(X T W ) for m not much larger than n.
Remark 5. A rather interesting feature of (3.6) is that the form of the kernel depends on the dimension d. This dependence is necessary to ensure the desired eigenvalue decay and maintain representational power as d grows. Moreover, this is not a typical feature of activations functions used in the literature.

The activation functions.
We build an activation function from the Wendland kernel in the manner suggested by Reproducing Kernel Hilbert Spaces (RKHS). Specifically, we consider As before, it is important to note that using φ d,k,ζ as a non-linearity does not imply that the population level matrix H is the associated kernel matrix. In fact, in general, we have as the inner-product over the sphere induced by the uniform distribution is not the same inner-product associated with the RKHS constructed from φ d,k,ζ . Nevertheless, Theorem 3.4 does show that H is strictly positive definite and, therefore, when this activation function is used with distinct data points λ min (H) > 0.
There is often a rather stark contrast between the conditioning of γ(X T W ) when either a ReLU or Wendland activation function is used (see Subsection 3.5). In general, the use of a Wendland activation function leads to a better conditioned matrix γ(X T W ) once m is a small multiple of n. While formally ζ should depend on properties of the data set X, we find that ζ = √ 2 is a natural starting point as it gives the Wendland based activation function the same support as ReLU, i.e., it is non-zero if x T i x j > 0.

Numerical experiments.
To complement our theoretical results and explore properties of (2.1) with various activation functions, we provide several numerical experiments. We consider four data sets and four activation functions (see Table 1).
1. Synthetic-1. We take the data points to be 500 i.i.d. points sampled from a uniform distribution on S 9 . The label for a point is assigned to be x i is x T i e, where e is the vector of all ones. 8 2. Synthetic-2. We take the data points to be 92 i.i.d. points sampled from a uniform distribution on S 2 , combined with the bad point set in Definition 2.1, for a total of 100 points: 50 points are labeled 1 (including all of the points in the bad point set), and the other 50 points are labeled −1. 3. Prostate-Cancer. We consider the 97 data points from 8 clinical measures (features) used to predict prostate-specific antigen levels relevant to prostate cancer [46]. Each feature is transformed into its standard score, and then the data points are projected onto S 7 . The labels are the log of the antigen levels. 4. Fashion-MNIST. We generate the data points from 28 × 28 images of fashion items [56]. We first sample 250 images labeled as sneakers and 250 images labeled as t-shirts. We project this subset onto its first 10 principal components and then project each embedding point onto S 9 . The labels are +1 and −1 corresponding to sneakers and t-shirts, respectively. We explore the condition number of H for varying network widths and choices of the activation function. As a point of comparison, we provide information from H either via an analytical expression or high accuracy quadrature. For quadrature, we use the Funk-Hecke formula paired with a numerically computed truncated ultraspherical polynomial expansion of γ. There are many reasons to consider the condition number when studying (2.1). First, we observe that a more ill-conditioned population matrix H means that a larger width may be necessary for H to approximate H well (in the sense of Lemma 2.5). Second, the conditioning of H governs the convergence of gradient descent (see Section 4). Third, since the condition number of H is the square of the condition number of the matrix in the least-squares problem in (2.4), which bounds the sensitivity of the solution to perturbations in the labels [24,Chapter 21].
In our experiments, we vary the network width m and measure the condition number κ( H), where H comes from the activation functions in Table 1. In all of our experiments, we set ζ = √ 2 in the Wendland kernels. The swish activation is from a parametric family of activations designed to be smooth approximations of ReLU [42]; here, we use the default implemented in TensorFlow [48]. For each width and activation function, we repeat the experiment 100 times with weights sampled i.i.d. from U(S d−1 ).
As the network width becomes sufficiently large, the Wendland kernels lead to well-conditioned systems (see Figure 1). When using ReLU or swish, the conditioning shows little dependence on the network width, and the condition number of H converges to that of H when m ≥ 10n. Although the swish activation is essentially a smooth ReLU approximation [42], we can see that the smoothness produces poorly conditioned matrices. This aligns with the kernel perspective: generally, smooth kernels generate ill-conditioned kernel matrices. The results on Synthetic-2 verify Theorem 2.2, as H is numerically singular for all widths. 9 We also find that when the swish activation function is used γ(X T W ) does not have full row-rank, which aligns with the results in Appendix A.2.
4. Training loss. In practice, how a network is trained is as much a part of the model as the network architecture itself. Modern deep networks are often vastly over parametrized, and the optimization algorithms used during training can have numerous hyper-parameters and nebulous stopping criteria. While some recent work [17,16] proves that zero training error is achievable for (2.1) via gradient descent, the results are for unrealistically wide (large m) regimes and contain opaque consideration of the dependence on the minimal eigenvalue of H (or, more appropriately, G). 10 4.1. Zero training error in the random features regime. One can characterize the set of simple neural networks (2.1) with width m ≥ n that can achieve zero training loss for arbitrary y. This is not a new observation, as it is simply the theory of under-determined linear systems.
for any y ∈ R n .
Proof. For any A with full row-rank, In Proposition 4.1, we have effectively fixed the training data X. Furthermore, for fixed W it is easy to write down an exact characterization of the minimal 2-norm solution α of (2.4) (using, e.g., the SVD). However, practically, there are still numerous questions to address.
4.1.1. Ensuring full row-rank. For any γ and X such that λ min (H) > 0 we can easily ensure γ(X T W ) is full row-rank by using Lemma 2.4; if m ≥ 10nC 2 λ min (H) log(2n/δ) then γ(X T W ) is full row-rank with probability 1 − δ. 11 In Section 3, the assumption that λ min (H) > 0 does not hold for ReLU. Nevertheless, for any fixed X it is possible that λ min (H) > 0 and we frame our results in terms of this quantity. 9 As expected based on Theorem 3.1, H is numerically singular when its entries are computed explicitly. 10 As we saw previously, for finite n and d ≥ 3 this quantity can be made arbitrarily small, requiring arbitrarily wide neural networks to reach zero training error. 11 There is a broader class of γ for which one can concoct W based on the training data such that γ(X T W ) is full row-rank even with m = n. We are uninterested in such schemes.

4.2.
Conditioning and the convergence of gradient descent. Once we have a set of weights W for which γ(X T W ) is full row-rank, optimizing over the finallayer weights α reduces to the under-determined linear least-squares problem (2.4). To simplify the exposition in this section we let Z ≡ 1 √ m γ(X T W ) be a fixed matrix. While there are standard direct methods for solving such linear least-squares problems [21], we consider what happens with gradient-based methods common in machine learning.
When m ≥ n and Z has full row rank, gradient descent converges to a minimizer that achieves zero training error for a small enough fixed step size η. 12 First, we find it illustrative to show the convergence of gradient descent in Lemma 4.2 explicitly as it both clearly shows the dependence of η on the conditioning of Z and that convergence occurs at the expected rate. Theorem 4.3 then formalizes a statement about gradient descent, achieving zero training loss in our setting. Notably, in contrast to prior work, we have convergence to a specific solution (the least-norm one) that achieves zero training loss, not just convergence of the loss to zero.
Lemma 4.2. Assume Z ∈ R n×m has full row-rank and consider (2.4). For any α (0) , the sequence of gradient descent iterates given by Proof. First, we relate the objective function at step k + 1 to step k as This immediately implies that All that remains to ensure convergence, i.e., Zα ( ) − y 2 → 0 as → ∞, is to pick an η sufficiently small so that (I − ηZZ T ) 2 < 1. To accomplish this, we minimize where σ min (Z) and σ max (Z) are the smallest and largest singular values of Z. Concretely, the minimizer occurs when 1 − ησ 2 min (Z) = ησ 2 max (Z) − 1 so that η = 2/(σ 2 min (Z) + σ 2 max (Z)). Using this choice of η, we find that and the result follows.
When m > n, there are an infinite number of global minimizers for (2.4), and initializing with α (0) = 0 is a sufficient condition for converge to the minimizer α ln of (2.4) with the minimum 2-norm. Similar results hold for certain stochastic gradient methods [30,47].
Building on Lemma 4.2, Theorem 4.3 provides a clear statement of when gradient descent applied to the single hidden layer model (2.1) can achieve zero training loss.
. . , n, and x i = x j for i = j; • γ and H as in (2.7) with γ ensuring that λ min (H) > 0; and • a failure probability 0 < δ < 1. Let m ≥ 10nC 2 log(2n/δ)/λ min (H) and W ∈ R d×m be a weight matrix with columns drawn i.i.d. from U(S d−1 ). Then, there is a unique least-norm solution to (2.4) that achieves zero objective value with probability at least 1 − δ. Furthermore, when there is a unique least-norm solution the sequence of iterates {α (k) } ∞ k=0 generated by gradient descent on the function with α (0) = 0 and step size η satisfying 0 < η < 2/λ max ( H) converges to the leastnorm solution of (2.4) as k → ∞, where H is given in (2.6).
Proof. First, by Lemma 2.4 we have that H 0 and therefore γ(X T W ) is full row-rank. This implies that (2.4) has a unique solution. Now, we let where U ∈ R n×n and V ∈ R m×n have orthonormal columns and Σ ∈ R n×n is diagonal with entries σ j > 0 for 1 ≤ j ≤ n. Starting with α (0) = 0, by Lemma A.6, we have that where f (z) = 1 − 1 − ηz 2 . Using the fact that the unique minimal norm solution to (2.4) is α (ln) = V Σ −1 U T y, we have that Since σ 2 j = λ j ( H) for 1 ≤ j ≤ n choosing 0 < η < 2/λ max ( H) ensures that max 1≤j≤n |1 − ησ 2 j | < 1 and the result follows. Remark 6. The optimal step size in Theorem 4.3 (with respect to an upper bound on convergence to the least-norm solution) is η = 2/ λ min ( H) + λ max ( H) . In this case we have that , n δ , where G is the NTK matrix for the chosen activation function. Moreover, these results require smaller step sizes to ensure convergence. Irrespective of the known theoretical results, one may expect that Theorem 4.3 could be improved if we simultaneously optimize over W and α since γ(W T X) may become better conditioned as the weights are "learned". In Appendix A.3, we discuss this possibility and, in practice, observe the opposite effect. The condition number of γ(W T X) is observed to grow during training. This may suggest that zero training loss is achieved more quickly by simply fixing the weights at their initialization.

Numerical experiments.
Using the datasets from Subsection 3.5, we explore the implications of Lemma 4.2 and Theorem 4.3. Specifically, we demonstrate that gradient descent can eventually converge to zero training error and illustrate the dependence of that behavior on the conditioning of γ(X T W ).
First, we using the Prostate-Cancer dataset with m = 200n. Figure 2 shows the dependence between the condition number and the convergence rate of gradient descent. We observe that the condition number of κ( H) is smaller for the Wendland activation functionss (see Figure 1). Indeed, gradient descent with the fixed step size given in Lemma 4.2 shows that the Wendland activation functions lead to a faster convergence rate than ReLU and swish.
At the same time, the number of gradient descent steps needed for convergence is quite large -tens of thousands of steps for the Wendland activation and over a million steps for ReLU. With the swish activation function, gradient descent does not even converge after 10 million steps. This is far larger than the number of steps one might use in practice [58], and the decrease in error stagnates after just 20 iterations (see Figure 2, right).

Early iteration behavior as Landweber iteration.
To explain the differences in early iterations of gradient descent during training, we investigate the Landweber iteration [22]. As before, denote 1 √ m γ(X T W ) by Z. Let Z = U ΣV T be the "reduced" SVD of Z, where U ∈ R n×n , Σ ∈ R n×n , and V ∈ R m×n . Lemma A.6 shows that the iterates of gradient descent for (2.4) with fixed step size η are given by Notably, taking η as in Lemma 4.2 we have that Moreover, we also have that The key observation is that if f k (σ j ) ≈ 1 for all the singular values then α (k) ≈ V Σ −1 U T y, which is the least-norm solution of (2.4). This is accomplished for small k if the eigenvalues of H are relatively flat such that 2λ j ( H) ≈ λ min ( H) + λ max ( H) for all j. We empirically observe that the desired slow eigenvalue decay from the Wendland activation functions (see Figure 3 (left)). Figure 3 (right) shows that many Landweber filter values f k (σ j ) ≈ 1 by k = 20. The faster eigenvalue decay seen when using ReLU or swish explain the relative benefits of the Wendland activation function during the early iterations.
A more nuanced view of this discussion is that gradient descent quickly "resolves" certain components of the solution that are associated with the components of y in the direction of left singular vectors corresponding to singular values close to 1. This is further illustrated by (4.4) and Lemma A.6, where we explicitly see the interplay between singular values, the step size, and the decomposition of y in the basis of left singular vectors.

Conclusions.
Our theory and experiments highlight the limits of certain common analysis techniques and metrics, advocate for careful consideration of activation functions, and provide insight on properties that aid in understanding the training of certain simple neural networks. By considering the random features regime we were able to extensively leverage existing theory for under-determined linear system and, from this perspective, it is entirely unsurprising that simple optimization methods are able to find solutions that achieve zero-training loss. Thus, careful understanding of the properties of models learned during training beyond their training loss is essential to discern the relative differences between models. For simple over-parameterized models, cleanly characterizing sets of solutions that achieve small training error and have preferable properties (e.g., favorable robustness or generalization behavior) provides a major opportunity.
We can again answer how well H approximates H. To do so, we first slightly expand our assumption on the activation function γ, we now require that For example, (A.3) is satisfied with C = 1 for the ReLU and swish activation functions. Next, the theoretical setup will depend on sub-Gaussian random variables.
Definition A.1 (Sub-Gaussian random vector). A random vector y in R d is sub-Gaussian with parameter J if The h k in (A.2) are sub-Gaussian random vectors. Proof. For any z ∈ S n−1 , where the second inequality follows from (A.3). Since W has i.i.d. N (0, 1) entries and each data point in X has unit 2-norm, y r = X T w r ∼ N (0, XX T ) and Cy T r z ∼ N (0, C 2 z T XX T z). Thus, The first inequality follows from the fact that a univariate normal random variable is sub-Gaussian with sub-Gaussian norm equal to an absolute constant times its variance [51, Example 2.5.8].
We use the following bound on the second moments of sub-Gaussian random variables.
Proposition A.3. Let Z 1 , . . . , Z m be i.i.d. random vectors in R n that have sub-Gaussian distributions with parameter J and let m ≥ n. Then, with probability at least 1 − 2e 2n−c m 2 /J 2 , where c is an absolute constant.
Proof. See the proof of Proposition 2.1 in [50].
An immediate corollary is a lower bound on the width that makes H close to H. In this result, σ 1 (X) 2 /λ min (H) is the analog of κ(H) in Lemma 2.5. We can use this fact to develop an analogous result for zero training error with weights initialized as i.i.d. Gaussian entries.
Theorem A.5 (Analog of Theorem 4.3). Consider • a dataset {(x i , y i )} n i=1 with x i ∈ R d , y i ∈ R, x i 2 = 1 for all i = 1, . . . , n, and x i = x j for i = j; • γ and H as in (A.3) and (A.1) with γ ensuring that the λ min (H) > 0; and • a failure probability 0 < δ < 1. Let m satisfy the bound in (A.5) and W ∈ R d×m be a weight matrix with i.i.d. random entries W ij ∼ N (0, 1). Then, there is unique least-norm solution of (2.4) that achieves zero training loss with with probability at least 1 − δ. Furthermore, the sequence of iterates {α (k) } ∞ k=0 generated by gradient descent on the function A.2. Failure of the swish activation function. Here, we show that the swish activation function has the same issues as ReLU (see Theorem 3.1) with the same point set (see Definition 2.1). The swish feature map is γ(t) = t/(1 + e −t ). Note that (γ(t) + γ(−t))/2 = z/2 so that γ(t) = γ even (t) + t/2, i.e., Therefore, we cannot conclude that H is a strictly positive definite matrix. In fact, we find that the point set in (2.3) also causes the matrix H jk = φ(x T j x k ) to be singular. To see this, one can repeat the argument in Theorem 3.1. In the same way we need to show that φ(1) − 3φ(s 2 ) + 3φ(−s 2 ) − φ(−1) = 0 when s = 1/ √ 3, which holds since φ even (1) − 3φ even (s 2 ) + 3φ even (−s 2 ) − φ even (−1) = 0 and hence, A.3. Joint weight training. In Subsection 4.2, we observe that ReLU can lead to ill-conditioned systems. However, the hidden layer weights were fixed, and one might suspect that jointly training all of the weights -along with a more sophisticated gradient method -might lead to better conditioning. To test this, we jointly learned all weights on the Synthetic-1 dataset with width m = 1.5n, using a stochastic gradient method with batch size equal to 10, momentum equal to 0.9, weight decay equal to 1e-5, 32-bit floats instead of 64-bit floats, and a learning rate decay of 0.99 after each epoch. We found that the condition number of γ(W T X) actually increases during training (see Figure 4). where Z is an n × m matrix, n < m, with full row rank. We consider computing the minimizer of (A.6) using gradient descent with a fixed step size η from the starting point α (0) = 0. The gradient descent iterates are (A.7) α (k+1) = α (k) − η(Z T Zα (k) − Z T y) = k j=0 (I − ηZ T Z) j · ηZ T y.
This is known as the Landweber iteration [22, Chapter 6]. Landweber iteration is typically analyzed for overdetermined least-squares problems. Below, Lemma A.6 works out an expression for the iterates in the underdetermined setting.
Lemma A.6. Consider Z ∈ R n×m and y ∈ R m with m ≥ n. Assume Z has full row-rank, 0 < η < Z 2 2 , and let Z = U ΣV T be the reduced SVD of Z. The sequence of gradient descent iterates {α (k) } ∞ k=1 for