Title: A RANDOM MATRIX PERSPECTIVE ON MIXTURES OF NONLINEARITIES IN HIGH DIMENSIONS Anonymous

Abstract: One of the distinguishing characteristics of modern deep learning systems is that they typically employ neural network architectures that utilize enormous numbers of parameters, often in the millions and sometimes even in the billions. While this paradigm has inspired significant research on the properties of large networks, relatively little work has been devoted to the fact that these networks are often used to model large complex datasets, which may themselves contain millions or even billions of constraints. In this work, we focus on this high-dimensional regime in which both the dataset size and the number of features tend to infinity. We analyze the performance of a simple regression model trained on the random features F = f (W X + B) for a random weight matrix W and random bias vector B, obtaining an exact formula for the asymptotic training error on a noisy autoencoding task. The role of the bias can be understood as parameterizing a distribution over activation functions, and our analysis directly generalizes to such distributions, even those not expressible with a traditional additive bias. Intriguingly, we find that a mixture of nonlinearities can outperform the best single nonlinearity on the noisy autoecndoing task, suggesting that mixtures of nonlinearities might be useful for approximate kernel methods or neural network architecture design.

1. INTRODUCTION:

It is undeniable that in recent years deep learning systems have found widespread success in their applications to a diverse and ever-expanding set of domains. The foundational results on many tasks such as image recognition (Krizhevsky et al., 2012) , speech recognition (Hinton et al., 2012) , and machine translation (Wu et al., 2016) , have begun to make their way into higher-level products that people interact with and rely upon in their daily lives. Whether these products generate a medical diagnosis, a navigation decision, or some other important output, it is crucial to understand the inner-workings of the deep learning algorithms that generate them.

Unfortunately, our theoretical understanding of these deep learning algorithms continues to lag behind their impressive practical successes. One main challenge in building a fuller understanding stems from the fact that deep neural networks are complex nonlinear functions that employ millions or even billions (Shazeer et al., 2017) of parameters. Traditional wisdom would suggest that to this parameter complexity corresponds to an optimization difficulty. Recent work, however, suggests that as the width of a network's hidden layers becomes large, the loss function simplifies and a theoretical analysis becomes tractable (Jacot et al., 2018; Mei et al., 2018; Chizat & Bach, 2018a; Mei et al., 2019; Rotskoff et al., 2019; Rotskoff & Vanden-Eijnden, 2018) . In some scenarios, the simplification is such that throughout training the parameters of the model stay within an infinitesimal radius of their initial (random) values, implying that much about neural network training can be understood by studying the random initialization (Jacot et al., 2018; Chizat & Bach, 2018b; Lee et al., 2019) .

Another main challenge in building a rich understanding of deep learning systems stems from the fact that they are often trained on very large, complex datasets: even if the models themselves are very large, they may not be large in comparison to the number of constraints they are designed to satisfy. Indeed, many important phenomena may become apparent only by examining the high-dimensional regime in which the dataset size and width are both large and of the same order.

In this work, we focus on the high-dimensional regime and analyze the performance of a regression model trained on the random features F = f (W X + B) for a random weight matrix W and random bias vector B. We obtain an exact formula for the training error on a noisy autoencoding task in the limit that the width and dataset size both go to infinity. The result is determined by the resolvent of the kernel matrix F t F , whose properties we analyze via the resolvent method from random matrix theory. Our analysis also provides an exact formula for the eigenvalue density of the kernel matrix, which may be of independent interest since it provides a characterization for how spectral properties of the data covariance matrix propagate through neural network layers at initialization.

1.1. OUR CONTRIBUTIONS:

The main contribution of our work is an exact characterization of the training error of a ridgeregularized random feature regression model on a noisy autoencoder task in the high-dimensional regime. This is one of the first non-trivial models to be solved exactly in the joint limit of large data and large width and provides an interesting testing ground in which to analyze this regime. Some of our additional contributions include,

• An exact characterization of the spectral density of the randm feature matrix F = f (W X + B), extending prior results of Pennington & Worah (2017) to non-Gaussian data distributions and to non-zero bias distributions. • One interpretation of the random additive bias is that it induces a distribution of activation functions parameterized by B, i.e. f (Z; B) := f (Z + B). Our analysis trivially extends to any distribution of activation functions f (• ; B) parameterized by B. • We show that there exists a non-trivial distribution over activation functions that outperforms the best possible single activation function on a noisy autoencoding task. • Our method of proof introduces a surrogate "linearization" of F , F lin , that possesses the same spectral information as F . F lin and its properties are likely to be of further interest and utility in analyzing neural networks in high dimensions.

1.2. RELATED WORK:

Neural networks have been studied from the perspective of high-dimensional statistics in a number of recent works. Most prior work has focused on the bias-free case. Pennington & Worah (2017) studied the spectrum of the activation matrix f (W X) for iid Gaussian data and derived an analytic expression for the training error of a ridge-regularized random feature model trained on pure noise. It is natural to consider incorporating biases by appending a constant feature 1 to X. Unfortunately, this leads to biases that are the same order as the weights, and so the effect disappears in the large dataset limit. Moreover, this modification on the data violates the assumptions of Pennington & Worah (2017) . Hastie et al. (2019) study ridgeless interpolation in high-dimensional interpolation for linear features as well as nonlinear random features of iid Gaussian data. Louart et al. (2018) derived a deterministic equivalent for the resolvent of the kernel matrix F t F which allowed for a characterization of the asymptotic training and test performance of linear ridge regression of random feature models.

Other work has investigated learning dynamics and generalization in the high-dimensional regime (Liao & Couillet, 2018a; Lampinen & Ganguli, 2018; Advani & Ganguli, 2016; Advani & Saxe, 2017) as well as the spectra of more complicated objects such as the Hessian (Pennington & Bahri, 2017) and Fisher information matrix (Pennington & Worah, 2018) . From the mathematical perspective, random matrix theory provides natural tools (Silverstein & Bai, 1995a) for analyzing the behavior of neural networks in the high-dimensional regime. Liao & Couillet (2018b) examined spectra for data drawn from Gaussian mixture models; see also (El Karoui, 2010) on the spectra of random kernel matrices.

2. PRELIMINARIES:

Consider a dataset X ∈ R n0×m and the random feature matrix,

F = f (W X; B),

generated by a single hidden-layer network with iid Gaussian weights W ∈ R n1×n0 (W ak ∼ N (0, σ 2 W /n 0 )), activation function f , and biases B = b1 m T ∈ R n1×m (for b ∈ R n1 ). We regard the second argument of f as parametrizing (continuously or discretely) an ensemble of activation functions. We refer to B (or b) as the bias, in reference to the important special case f (W X + B) (additive bias). In general, however, we only assume the parameters b a are such that the measure1 n1 a δ ba converges in distribution to some limiting distribution µ B , effecting an arbitrary distribution over activation functions.

We assume that E |f (N ; b)| k for N ∼ N (0, σ) is finite for all 1 ≤ k ≤ 3, σ > 0, and b ∈ support(µ B ). When µ B is a single Dirac mass at location b 0 , the activation function can be written as f (W X; B) = f (W X; b 0 ) = g(W X) for some single-argument function g (for an additive bias, g(W X) = f (W X + b 0 )). When this is the case, we say the model has a single activation function, as opposed to a mixture or distribution of activation functions.

The quantity of interest for the investigations below is the kernel matrix 1 n1 F t F , and in particular its resolvent,

G(z) = 1 n1 F t F -zI -1

.

(1)

As we review in Sec. 4, the optimal regression coefficients of a linear model on the random features F is a simple function of this resolvent.

The high-dimensional regime that we study is the one in which the dataset size m, feature dimensionality n 0 , and hidden layer width n 1 all go to infinity at the same rate. In particular, as is standard in the random matrix literature, we assume that we can parameterize the limit in terms of the dataset size m in such a way that there exist two positive constants,

EQUATION

Note that the resolvent is a random matrix, but as m grows large, its normalized trace becomes a deterministic quantity. In the limit that m → ∞, this quantity is known as the Stieltjes transform 1 ,

EQUATION

Together with an auxiliary transform m(z), defined below, these deterministic quantities completely characterize the asymptotic training error of kernel ridge regression on a noisy autoencoder task in this high-dimensional regime.

The Stieltjes transform frequently arises in random matrix methods as a way to encode the spectra of matrices. In particular, if λ i are the eigenvalues of 1 n1 F t F and the empirical distribution of eigenvalues converges in distribution to some deterministic limiting density as m → ∞,

EQUATION

then (with appropriate technical assumptions), the limiting spectral density itself can be recovered from the Stieltjes transform m(z) via the inversion formula,

EQUATION

The Stieltjes transform then substitutes convergence in distribution for pointwise convergence for all z such that z > 0.

2.1. METHODS FOR COMPUTING THE STIELTJES TRANSFORM:

We briefly review two standard methods for computing the Stieltjes transform m(z), the resolvent method and the moments method.

The resolvent method is an approach for computing the Stieltes transform based on the application of the Schur complement formula to the resolvent itself (or to a closely-related block matrix). Intuitively, as the matrix size becomes large, the minors of the matrix are similar in distribution to the larger matrix, and, moreover, the Cauchy interlacing theorem guarantees that their Stieltjes transforms are We get excellent agreement of theory and simulation for spectral densities for any bias distribution. We set φ = 1.5, ψ = 0.8, σ X = σ W = 1, and f = ReLU. Simulations are performed on matrices of size m = 2 14 . Left Gaussian distribution over the biases with distribution N (0, 1). Right Bernoulli distribution over biases with distribution p = 0.5. close as well. This allows for the derivation of a self-consistent equation (SCE) in which the Stieltjes transform appears on the left-hand side as the trace of the resolvent, and on the right-hand side as the trace of one of its minors.

The moments method is more combinatorial in nature and involves expanding the resolvent for large z and computing the traces of each term,

EQUATION

The traces themselves are expended out as

EQUATION

where the sum runs over matrix indices a 1 , . . . , a k , α 1 , . . . α k . The essence of the moment method involves analyzing the asymptotic contribution of each term in the sum based on its combinatorial type and the details of F , and resumming the results to obtain m(z).

We refer the reader to (Erdos & Yau, 2017; Tao, 2012) for more details about these methods and additional background on random matrix theory.

3 RESULT FOR STIELTJES TRANSFORM

3.1. MAIN THEOREM:

We make the following assumptions on the data matrix X and bias vector b:

1. 1 n0 a X aα X aβ = δ αβ σ 2 X + O 1/ √ n 0 for all α and β 2. the empirical eigenvalue distribution of 1 n0 X t X converges in distribution to a measure µ X 3. 1 n1 n1 a=1 δ ba → µ B in distribution. Theorem 1. Define σ Z = σ W σ X and resolvent G(z) = 1 n1 F t F -zI -1

. Then under the above assumptions and for all z such that z > 0, the transforms

EQUATION

converge in probability to the unique solution, m(z) and m(z), of the Eqn. (9) that map C + to C + :

m(z) = E S∼µ X 1 C 0 (z) + SC 1 (z) and m(z) = E S∼µ X S C 0 (z) + SC 1 (z) (9)

where, 

EQUATION

D(B):

EQUATION

and where η(B) and ζ(B) are the Gaussian expectations

EQUATION

)

The proof is quite involved and is presented in the supplementary material. The basic idea is to derive a multivariate Gaussian random matrix model with the same correlation structure as F , then derive a self-consistent equation (SCE) using the resolvent method for this linearized version of F , F lin .

Remark 1. The self-consistent equations consist of two coupled equations involving the Stieltjes transform m(z) and an auxiliary object m(z), (cf. (Paul & Silverstein, 2009, Eq. (2))), which we will see in Cor. 1 essentially measures the autoencoding capacity of the network.

Remark 2. Note that the self-consistent equations contain an expectation over the limiting spectral density of the input data. While the assumptions on the data matrix X in Thm. 1 are quite general, they may not be optimal. See Sec. 3.3, where we show strong agreement with empirical data from MNIST and CIFAR-10 and for a range of synthetic distributions. This suggests that the theorem may hold for even more general data distributions.

3.2. ALTERNATE REPRESENTATION AND LIMITING RESULTS:

When the data distribution is iid Gaussian, the expectations in Eqn. ( 9) can be expressed in closed form, though one must be careful to choose the correct branch of the resulting function. For simplicity and future reference, we focus on the setting where 0 < φ ≤ ψ ≤ 1, in which case we have the coupled algebraic equations,

EQUATION

)

When µ B in Eqn. ( 10) is trivial, i.e. a single Dirac mass, then the result should reduce to the single activation function case with F = f (W X), which was studied in (Pennington & Worah, 2017) .

Indeed, writing η = E B [η(B)] and ζ = E B [ζ(b)

] for such a distribution, and eliminating m(z) from Eqns. ( 13) and ( 14), we find that m(z) satisfies the following quartic polynomial:

EQUATION

which agrees with the result in (Pennington & Worah, 2017) upon identifying m(z) = -(1φ/ψ)/z -φ/ψG(z).

3.3. SPECTRAL DENSITY ESTIMATES:

The self-consistent equations in Thm. 1 can be solved numerically by iterating Eqn. (9) until convergence, using numerical integration. By utilizing the Stieltjes inversion formula, Eqn.

(5), we can extract predictions for the limiting eigenvalue density of 1 n1 F t F . The results show close agreement with empirical spectral simulations from several interesting practical datasets and synthetic distributions, see Figs. 1 and 2 .

4. RIDGE-REGULARIZED NOISY AUTOENCODER:

We consider the problem of kernel ridge regression with random features given by F = f (W X; B) and noisy regression targets given by Y = AX + for some A ∈ R n2×n0 and independent Gaussian noise ∈ R n2×m such that ij ∼ N (0, σ 2 ). As is common in the literature of high-dimensional statistics, we assume an isotropic prior on A such that A t A → n2 n0 σ 2 A I as m, n 0 , n 2 → ∞. Note that when σ A = 0, we have the pure memorization setting studied by Pennington & Worah (2017) . Corollary 1. Let W * 2 be the minimizer for regularized training loss

EQUATION

)

with random features F = f (W X; B). Then the asymptotic training error converges in probability as

EQUATION

In particular, we see that the derivative of the Stieltjes transform m (-γ) measures the capacity to learn noisy labels, whereas m (-γ) measures pure autoencoding capacity. See Fig. 3 for a comparison between these theoretical predictions and simulation.

Proof. The optimal weights for the regularized loss are given by

W * 2 = 1 n1 Y G(γ)F t for G(γ) = ( 1 n1 F t F + γI) -1 , resulting in training error E train = 1 n2m ||Y -W * 2 F || 2 = γ 2 1 n2m tr Y t Y G(γ) 2 = γ 2 σ 2 A 1 n0m tr X t XG(γ) 2 + γ 2 σ 2 1 m tr G(γ) 2 = γ 2 σ 2 A d dγ ( m(-γ)) + γ 2 σ 2 d dγ (m(-γ)) = -γ 2 σ 2

A m (-γ) + σ2 m (-γ) . Remark 3. As in (Pennington & Worah, 2017) , there is a scaling homogeneity in the E train : an increase in the regularization constant γ can be compensated by a decrease in scale of W 2 , which, in turn, can be compensated by increasing the scale of F , which is equivalent to increasing η(b) and ζ(b). Owing to this homogeneity, we are free to choose a normalization of the activation function for which E b∼µ B [η(b)] = 1.

4.1. NONLINEAR MIXTURES CAN OUTPERFORM SINGLE NONLINEARITIES:

The bias term in our random feature model can be viewed as one way of defining a distribution over activation functions. The choice of distribution, in general, affects the performance of the model on a given task -as quantified by the expectations in Eqn. (10). A key benefit of our model and analytical approach is that it permits nontrivial distributions of nonlinear activation functions.

In this section, we build on our results for the training loss on noisy autoencoder tasks to examine the benefits of utilizing nonlinear mixtures. The goal here is not to identify "good" mixtures, since this will clearly be a dataset-and architecture-dependent question, nor is it to identify a large performance gap. Instead, we merely seek to demonstrate a proof-of-principle, namely that there exist non-trivial distributions over nonlinearities that can provably outperform the best possible single nonlinearity.

For this analysis, we consider the simplest possible nontrivial distribution over activation functions: a Bernoulli mixture of two different functions. To each of these functions we associate two constants, η and ζ, which derive from Eqn. ( 12) but have no B-dependence since each function is a single nonlinearity. Concretely, let

EQUATION

)

For the two functions themselves, we utilize (i) a "pure linear" activation function with η = 1, i.e. the identity function and (ii) a "pure nonlinear" (Hastie et al., 2019) activation function with η = 1 and ζ = 0. (The particular purely nonlinear function in (ii) is irrelevant, as our theory predicts and our experiments confirm; see Fig. 4 (a)). To be precise, we define for p ∼ Bernoulli(p),

EQUATION

)

where g ζ=0 is any function with η = 1 and ζ = 0 (see, e.g., the functions in Fig.

( 3) of (Pennington & Worah, 2017) ).

The task of computing E train for f p is cumbersome but purely algebraic. To see how to proceed, notice that the expectations in Eqn. ( 10) are simple for f p :

EQUATION

)

Plugging these equations into Eqns. ( 13) and ( 14), collecting terms and simplifying yields a set of coupled polynomial equations for m(z) and m(z). Taking the total derivative of these equations with respect to z yields two additional equations which can be solved to express m (z) and m (z) in terms of m(z) and m(z). Combining these results produces a polynomial system whose solution 2 encodes E train through Eqn. ( 17). The question is, does a non-trivial mixture actually outperform a single nonlinearity? First, we must understand the performance of the optimal single nonlinearity. We note that owing to the homogeneity of the the training loss in η, ζ, and γ, we can assume without loss of generality that η = 1. Therefore the entire effect of the nonlinearity should be encoded in the single constant ζ. In Fig. 4 (a), we plot our theoretical prediction for E train in solid lines and empirical simulations for finite networks as 1σ error bars. The activation function used for each simulation is chosen randomly, conditional on the value of ζ. So the good agreement in the left panel of Fig. 4 (a) demonstrates not just the correctness of our theoretical result but also the fact that E train depends on the activation function solely through the constant ζ. The blue stars in this figure indicate that the optimal single nonlinearity is neither purely linear (ζ = 1), nor purely nonlinear (ζ = 0), but rather something in between.

For this particular problem setup, the performance of the optimal single nonlinearity and the optimal Bernoulli mixture are rather close, as indicated by the top panel of Fig. 4 (b). However, owing to our precise analytical formulation, we can evaluate the training loss to high precision and observe that there is indeed a difference in performance between the two models, as shown in the bottom panel of Fig. 4(b) . This result establishes that there are some problems for which even the best single nonlinearity is outperformed by a mixture of nonlinearities.

5. CONCLUSIONS:

In this work we studied the feature matrix F = f (W X; B) where W is a random matrix with iid Gaussian entries. Under mild assumptions on X and B, we obtained an exact analytic formula, Eqn. ( 9), that characterizes the Stieltjes transform of the spectral density of F . The result allowed us to describe the exact training loss of a ridge-regularized noisy autoencoder in the high-dimensional, large-dataset limit, providing one of the first closed-form solutions to a non-trivial model in this limit. We found excellent agreement between the asymptotic predictions of Eqn. ( 9) and a variety of finite-dimensional empirical simulations.

We also advanced the interpretation of the bias B as one particular way of parameterizing a distribution of activation functions. Indeed, our derivations proceed completely unchanged whether this distribution is of the traditional additive form f (• + B) or the more general f (•; B). By examining the latter, we showed that there are configurations in which a non-trivial distribution over activation functions provably outperforms the best possible single activation function. This opens the door to future investigations regarding optimal methods for parameterizing distributions over activation functions for approximate kernel methods, and suggests the possibility that mixtures of nonlinearities could be a useful design consideration when constructing neural network architectures.

Remark 4. For simplicity, while writing f (i) we assume differentiability of f up to ith order. However, this assumption can be relaxed, as since the derivatives are only used to take a Gaussian expectation, they can be expressed using Stein's lemma. For example, ξ 1 (x) = EN f (σ X σ W N + x). Thus, similar results should hold for rectified linear units for example, regardless of their lack of differentiability at 0. Lemma 1. To leading order the correlation structure of F is

EQUATION

)

Proof of Lem. 1. Before we prove Lem. 1, we state a simple lemma that is easy to show using Taylor expansion.

Lemma 2. Suppose ε is a small constant and that N ∼ N (0, 1), then

EQUATION

EQUATION

) and

EQUATION

)

Recall ξ 0 (x) = 0 by assumption. Using the observation in Eqn. ( 21), we repeatedly rewrite the expectations in F t F as we can reduce the expectation over W to an expectation over at most two correlated Gaussian random variables. Let N and N denote independent, standard Gaussian random variables. Using Lem. 2, we see

EQUATION

Again, using Lem. 2, we get

EQUATION

For the covariance calculation, we define ρ αβ := c αβ / √ c αα and recall that ρ αβ is order 1/ √ m. Thus,

EQUATION

)

For the first term, we get

EQUATION

)

For the second term, we have

EQUATION

Thus, to second order we have

EQUATION

A.1 LINEARIZED MODEL Define F lin := CΘ 1 Σ + V 2 -C 2 1/2 Θ 2 , (40) where: Θ 1 and Θ 2 are n 1 × m matrices that have iid Gaussian entries with mean 0 and variance 1/n 1 ; 2) V and C are n 1 × n 1 diagonal matrices with entries

EQUATION

for a bias vector b; 3) Σ is a matrix square root of 1 n0 X t X. A simple calculation shows that the entries of the matrix F lin match the first and second mixed moments of F given in Lem. 1. However, unlike F , F lin is linear in the random matrices Θ 1 and Θ 2 ; in this sense, it is a linearization of F .

B DERIVATION OF SELF-CONSISTENT EQUATION FOR F LIN:

Theorem 2. Let m m (z) be the Stieltjes transform of 1 n1 F lin t F lin , i.e. 1 m trG(z). Then with probability 1, as n 1 → ∞, for all z such that z > 0, m n1 (z) → m(z), where m(z) is the solution of the coupled equations

EQUATION

EQUATION

We want to study the eigenvalues of 1 n1 F lin t F lin . Note without loss of generality it is sufficient to consdier diagonal Σ, since we can diagonalize Σ using some orthogonal matrices O and O as it is the square root of a positive definite matrix. Moreover, these orthogonal matrices, when applied to either Θ 1 or Θ 2 do not change their distributions. Thus, diagonalizing Σ, so that Σ αα = λ X α , where λ X α are the eigenvalues of X t X/n 0 , results in an equivalent matrix ensemble in distribution (see Silverstein & Bai (1995b) for more detail). With this simplification, F lin has independent entries given by

EQUATION

To make the derivation easier, we can partly linearize the problem by studying the matrix

EQUATION

)

By the Schur complement formula, one easily finds

EQUATION

and z mm (z) :

EQUATION

)

where G is the inverse of H. Again by the Schur complement formula

EQUATION

and

EQUATION

for α ∈ {1, . . . , m} and a ∈ {m + 1, . . . , m + n 1 } and G (a) is the inverse of the minor H (a) . Since G(α) is independent of Θ 1 1α , . . . , Θ 1 n1α and Θ 2 1α , . . . , Θ 2 n1α , we see by taking the expectation over these variables that

EQUATION

Moreover, standard concentration inequalities and the Ward identity (see Erdos & Yau (2017) ) show n1 a,b=1 We may also replace G (α) with G at the expense of another small error that can be bounded using the Cauchy interlacing theorem:

EQUATION

EQUATION

G (α)

m+a,m+a -G m+a,m+a ≤ O(1/m). Using this control over these sums, we see In particular α = -1 is a linear function, α = 0 is regular ReLU, and α = 1 is the a scaled absolute value function. We normalize all functions so that E b [η(b)] = 1. We get excellent agreement with theory from only a single sample.

1 G αα = -z -

:

The distinction between the Stieltjes transform m(z) and the dataset size m should be clear from context.

Special care must be taken in selecting the correct root of this equation, in accordance with the condition that m(z) ∼ -1/z for large |z|.

Figures and Tables:

Figure1: We get excellent agreement of theory and simulation for spectral densities for any bias distribution. We set φ = 1.5, ψ = 0.8, σ X = σ W = 1, and f = ReLU. Simulations are performed on matrices of size m = 2 14 . Left Gaussian distribution over the biases with distribution N (0, 1). Right Bernoulli distribution over biases with distribution p = 0.5.

Figure 2: Empirical spectral densities agree with our predictions for varied data distributions and shape parameters. Top left: One class from CIFAR (airplane), mean subtracted. Top right: Classes {0, 8} from MNIST, mean subtracted. Bottom left: Gaussian input data, varying the NN parameter settings and activation function. Bottom right: Input data with a bimodal spectrum. All plots used f = ReLU, φ = 1.5, ψ = 0.8 and σ W = σ B = 1, except for the indicated modified parameter. Empirical densities were smoothed using a Gaussian KDE.

Figure 3: Comparisons of simulated ridge regression error and our theoretical prediction. We use ReLU with σ X = σ W = σ B = 1 for all plots and vary the shape parameters φ and ψ. For simulations we use m = 2 13 throughout. Note we also normalize the activation function so that E b [η(b)] = 1. (a): Our predictions for ridge regression with random labels are solid lines. Simulated losses are the red stars. (b) Autoencoder error. (c) Noisy autoencoder error with σ ε = 0.2.

Fig. 4(a)  shows the result of this calculation in solid lines for various values of γ, while the 1σ error bars show empirical simulations with finite networks. The red stars in the figure show that for many values of γ, the optimal mixture percentage is intermediate, i.e. 0 < p < 1.

Figure 4: Performance on ridge-regularized noisy autoencoder with σ = 1, φ = 1/2, and ψ = 1/2. (a) Theoretical predictions for training error (solid lines) and 1σ error bars for empirical simulations of finite networks (n 0 = 192, n 1 = 384, m = 384) for various values of ridge regularization constant γ as the activation function varies. In the left panel, a single activation function f is used. In the right panel, the non-linearity is f p , a Bernoulli(p)-mixture of a purely linear (ζ = 1) and purely non-linear (ζ = 0) function. Each simulation uses a randomly-chosen non-linearity having the specified values of ζ, demonstrating that E train depends on the non-linearity solely through this constant. Red and blue stars denote minima. (b) Training error as a function of γ for the optimal f and f p , as determined in (a). The bottom panel shows the difference in training error, demonstrating that the optimal Bernoulli(p)-mixture of non-linearities has smaller training error than the best single non-linearity.

Figure 5: Comparison of right singular vectors of X and F for MNIST (left) and CIFAR10 (right).Entry ij shows x i • f j . Although our theoretical results do not give predictions for how the singular vectors change, we found interesting behavior, with very little change to the largest singular vectors (which are nearly isolated in the spectrum), but more mixing of singular vectors in the dense part of the distribution.

Figure6: Comparisons of simulated ridge regression error and our theoretical prediction. Here we vary the activation function by changing the slope α in leaky ReLU. In particular α = -1 is a linear function, α = 0 is regular ReLU, and α = 1 is the a scaled absolute value function. We normalize all functions so that E b [η(b)] = 1. We get excellent agreement with theory from only a single sample.

Peer review (iclr2020):
This paper investigates the asymptotic spectral density of a random feature model F(Wx + B). This is an extension of existing result that analyzed a model without the bias term, i.e., F(WX). This extension requires a modification of the proof technique. In addition to that, it analyzed a mixture of linear and non-linear activation functions, and show that mixture is better than single nonlinearity in terms of expected training error for ridge regression estimators. Pros: - This paper investigates an interesting problem and it successfully extends the existing work. The theoretical curve well matches the simulated curve. - The finding that mixture of nonlinearities gives better expected training error is interesting. Cons: - The extension to the model with bias seems a bit incremental. In practice, we may consider an input with additional constant feature, X <- [X,1], to deal with both models in a unified manner. There should be more discussion about why this kind of trivial argument cannot be applied in the analysis. - The effect of mixture of activation functions is investigated in the *training error,* but I don*t see much significance on investigating the training error thoroughly. Instead, people are interested in the test error. I guess there does not appear such a trade-off for the test error and the linear activation function would be always better because the true function is the linear model. Hence, more expositions about why the training error is investigated should be provided. More minor comment: - I guess the definition of Etrain (Eq.(17)) requires an expectation with respect to the training data. - Assumptions of the activation function f should be provided; is it just assumed to be differentiable?, ReLU is included? - The definition of G(gamma) in page 6 had better to be consistent to that in previous pages.