Title: LEARNING BACKPROPAGATION-FREE DEEP ARCHIT-ECTURES WITH KERNELS

Abstract: One can substitute each neuron in any neural network with a kernel machine and obtain a counterpart powered by kernel machines. The new network inherits the expressive power and architecture of the original but works in a more intuitive way since each node enjoys the simple interpretation as a hyperplane (in a reproducing kernel Hilbert space). Further, using the kernel multilayer perceptron as an example, we prove that in classification and under certain losses, an optimal representation that minimizes the risk of the network can be characterized for each hidden layer. This result removes the need of backpropagation in learning the model and can be generalized to any feedforward kernel network. Moreover, unlike backpropagation, which turns models into black boxes, the optimal hidden representation enjoys an intuitive geometric interpretation, making the dynamics of learning in a deep kernel network transparent. Empirical results are provided to complement our theory.

1. INTRODUCTION:

Any neural network (NN) can be turned into a kernel network (KN) by replacing each artificial neuron (McCulloch & Pitts, 1943) , i.e., learning machine of the form f (x) = σ(w x + b), with a kernel machine, i.e., learning machine of the form f (x) = w, φ(x) + b with kernel function k(x, y) = φ(x), φ(y) . This combination of connectionism and kernel method enables the learning of hierarchical, distributed representations with kernels.

In terms of training, similar to NN, KN can be trained with backpropagation (BP) (Rumelhart et al., 1986) . In the context of supervised learning, the need for BP in learning a deep architecture is caused by the fact that there is no explicit target information to tune the hidden layers (Rumelhart et al., 1986) . Moreover, BP is usually computationally intensive and can suffer from vanishing gradient. And most importantly, BP results in hidden representations that are notoriously difficult to interpret or assess, turning deep architectures into "black boxes".

The main theoretical contribution of this paper is the following: Employing the simplest feedforward, fully-connected KN as an example, we prove that in classification and under certain losses, the optimal representation for each hidden layer that minimizes the risk of the network can be explicitly characterized. This result removes the need for BP and makes it possible to train the network in a feedforward, layer-wise fashion. And the same idea can be generalized to other feedforward KNs.

The layer-wise learning algorithm gives the same optimality guarantee as BP in the sense that it minimizes the risk. But the former is much faster and evidently less susceptible to vanishing gradient. Moreover, the quality of learning in the hidden layers can be directly assessed during or after training, providing more information about the model to the user. For practitioners, this enables completely new model selection paradigms. For example, the bad performance of the network can now be traced to a certain layer, allowing the user to debug the layers individually. Most importantly, the optimal representation for each hidden layer enjoys an intuitive geometric interpretation, making the learning dynamics in a deep KN more transparent than that in a deep NN. A simple acceleration method that utilizes the "sparse" nature of the optimal hidden representations is proposed to further reduce computational complexity.

Empirical results on several computer vision benchmarks are provided to demonstrate the competence of the model and the effectiveness of the greedy learning algorithm. 2) • F (1) on the right is better than that on the left and hence the F (2) on the right corresponds to the network-wise optimality of the second layer. But the F (2) on the left is clearly a better estimate to the layer-wise optimality T 2 if the quality of estimation is measured by the supremum distance.

2. FROM NEURAL TO KERNEL NETWORKS:

In this section, we discuss how to build a KN using a given NN. First, the generic approach is described in Fig. 1a . Note that KN inherits the expressive power of the original NN since a kernel machine is a universal function approximator under mild conditions (Park & Sandberg, 1991; Micchelli et al., 2006) and the two models share the same architecture. However, KN works in a more intuitive way since each node is a simple linear model in a reproducing kernel Hilbert space (RKHS).

We now concretely define the KN equivalent of an l-layer Multilayer Perceptron (MLP), which we shall refer to as the kernel MLP (kMLP). 1 Given a random sample (x n , y n ) N n=1 , where (x n , y n ) ∈ X 1 × Y ⊂ R d0 × R, denote (x n ) N n=1 as S and (y n ) N n=1 as Y S for convenience. For i = 1, 2, . . . , l, consider kernel k (i) : X i × X i → R, X i ⊂ R di-1 (for i > 1, d i-1 is determined by the width of the i -1 th layer). k (i) (x, y) = φ (i) (x), φ (i) (y) Hi , where φ (i) is a mapping into RKHS H i .

For i ≥ 1, the i th layer in a kMLP, denoted F (i) , is an array of d i kernel machines:

F (i) : X i → R di , F (i) = (f (i) 1 , f (i) 2 , . . . , f (i) di ), a d i -tuple. Let F (0) be the identity map on R d0 , each f (i) j : X i → R is a hyperplane in H i : f (i) j (x) = w (i) j , φ (i) F (i-1) •• • ••F (0) (x) Hi +b (i) j , w (i) j ∈ H i , b (i) j ∈ R. In practice, f (i) j (x) is usually implemented as N n=1 α (i) nj k (i) F (i-1) • • • • • F (0) (x), F (i-1) • • • • • F (0) (x n ) + b (i) j , where the α (i) nj , b (i) j

∈ R are the learnable parameters. 2 The set of mappings

{F ( ) • • • • • F (1) : α (i) nj , b (i) j ∈ R

for all admissible n, j, i} defines an l-layer kMLP. In the rest of this paper, we shall restrict our discussions to this kMLP.

3. ASSUMPTIONS AND NOTATIONS:

We now specify the assumptions that we impose on all kernels considered in this paper. First, we consider real, continuous, symmetric kernels only and we call a kernel positive semidefinite (PSD) or positive definite (PD) if for any S, the kernel matrix defined as (G) mn = k(x m , x n ) is PSD or PD, respectively. We shall always assume that any kernel considered is at least PSD and that

k (i) (x, x) = c < +∞ for all x ∈ X i and inf x,y∈Xi k (i) (x, y) = a > -∞.

It is straightforward to check using Cauchy-Schwarz inequality that the first condition implies max x,y∈Xi k (i) (x, y) = c. For each fixed x ∈ X i , we assume that k (i) (x, y), as a function of y, is L (i)

x -Lipschitz with respect to the Euclidean metric on X i . Let sup x∈Xi L (i)

x = L (i) , which we assume to be finite.

The following notations will be used whenever convenient: We use the shorthand

F (1) (S) for (F (1) (x n )) N n=1 . For i = 2, 3, . . . , l, F (i) (S) := F (i) • F (i-1) • • • • • F (1) (S)

and the same is true with S substituted by any x. Throughout this paper, notations such as F (i) can either be used to denote a set of functions or a specific function in some set depending on the context. Also, when there is no confusion, we shall suppress the dependency of any loss function on the example for brevity, i.e., for a loss function , instead of writing (f (x), y), we shall write (f ).

4. A LAYER-WISE LEARNING ALGORITHM:

To simplify discussion, we shall restrict ourselves to binary classification (Y = {+1, -1}) and directly give the result on classification with more than two classes in the end. A generalization to regression is left as future work. Again, we only focus on kMLP although the idea can be directly generalized to all feedforward KNs. We now discuss the layer-wise learning algorithm, beginning by addressing the difficulties with training a deep architecture layer-by-layer.

4.1. FUNDAMENTAL DIFFICULTIES:

There are two fundamental difficulties with learning a deep architecture layer-wise. First, the hidden layers do not have supervision (labels) to learn from. And it depends on BP to propagate supervision from the output backward (Rumelhart et al., 1986) . We shall prove that for kMLP, one can characterize the optimal target representation for each hidden layer, which induces a risk for that layer. The target is optimal in the sense that it minimizes the risk of the subsequent layer and eventually that of the network if all layers succeed in learning their optimal representations. This optimal representation defines what we call "layer-wise optimality".

The other difficulty with layer-wise learning is that for any given hidden layer, when the upstream layers create error, layer-wise optimality may not coincide with "network-wise optimality", i.e., the solution of this layer that eventually leads to the composed solution that minimizes the risk of the network in this suboptimal case. Indeed, when a hidden layer creates error, the objective of any layer after it becomes learning a solution that is a compromise between one that is close to the layer-wise optimality and one that prevents the error from the "bad" layer before it from "getting through" easily. And the best compromise is the network-wise optimality. The two solutions may not coincide, as shown in the toy example in Fig. 1b . Clearly, we would like to always learn the network-wise optimality at each layer, but the learner is blind to it if it is only allowed to work on one layer at a time. By decomposing the overall error of the network into error at each layer, we prove that in fact, network-wise optimality is learnable for each hidden layer even in a purely layer-wise fashion and that the proposed layer-wise algorithm learns network-wise optimality at each layer.

4.2. THE OPTIMAL HIDDEN REPRESENTATIONS:

We now address the first difficulty in layer-wise learning. The basic idea is first described in Section 4.2.1. Then we provide technical results in Section 4.2.2 and Section 4.2.3 to fill in the details.

4.2.1. BASIC IDEA:



Given a deep architecture:

F := F (l) • • • • • F (1)

and a loss function l defined for this network which induces a risk that we wish to minimize: R l = E l (F). BP views this problem in the following way: R l is a function of F. The learner tries to find an F that minimizes R l using the random sample S with labels Y S according to some learning paradigm such as Empirical Risk Minimization (ERM) or Structural Risk Minimization (SRM) (Vapnik, 2000; Shalev-Shwartz & Ben-David, 2014) . S is considered as fixed in the sense that it cannot be adjusted by the learner.

Alternatively, one can view R l as a function of F (l) and the learner tries to find an F (l) minimizing R l using random sample S l-1 with labels Y S according to some learning paradigm, where S l-1 := F (l-1) •• • ••F (1) (S)3 . The advantage is that the learner has the freedom to learn both the function F (l) and the random sample S l-1 . And since S l-1 determines the decision of the learning paradigm, which then determines R l , R l is now essentially a function of both F (l) and S l-1 : R l = E l (F (l) , S l-1 ).

The key result is that independently of the actual learning of F (l) , one can characterize the sufficient condition on S l-1 under which R l , as a function of S l-1 , is minimized, as we shall prove. In other words, the "global minimum" of R l w.r.t. S l-1 can be explicitly identified prior to any training. This gives the optimal S l-1 , which we denote as S l-1 .

Moreover, the characterization of S l-1 gives rise to a new loss function l-1 and thus also a new risk R l-1 that is a function of 1) . Consequently, the same reasoning would allow us to deduce S l-2 before the learner learns F (l-1) . And this analysis can be applied to each layer, eventually leading to a greedy learning algorithm that sequentially learns

F (l-1) • • • • • F (

F (1) , F (2) •F (1) * , . . . , F (l) • F (l-1) * •• • ••F (1)

* , in that order, where the asterisk on the superscript indicates that the corresponding layer has been learned and frozen.

The layer-wise learning algorithm provides a framework that enjoys great flexibility. To be specific, one could stop the above analysis at any layer i, then learn layers i + 1, . . . , l in a greedy fashion but still learn layers 1, . . . , i together with BP. Thus, it is easy to see that BP can be brought under this framework as a special case. Nevertheless, in later text, we shall stay on the one end of the spectrum where each layer is learned individually for clarity.

We now present the formal results that give the optimal hidden representations. By the reasoning above, the analysis starts from the last hidden layer (layer l -1) and proceeds backward.

4.2.2. FORMAL RESULTS: S l-1:

To begin with, we need to approximate the true classification error R l since it is not computable. To this end, we first review a well-known complexity measure. Definition 4.1 (Gaussian complexity (Bartlett & Mendelson, 2002) ). Let P be a probability distribution on a metric space X and suppose x 1 , . . . , x N are independent random elements distributed as P. Let F be a set of functions mapping from

X into R. Define ĜN (F) = E sup f ∈F 2 N N n=1 g n f (x i ) x 1 , . . . , x N ,

where g 1 , . . . , g N are independent standard normal random variables. The Gaussian complexity of

F is G N (F) = E ĜN (F)

Intuitively, Gaussian complexity quantifies how well elements in a given function class can be correlated with a noise sequence of length N , i.e., the g n (Bartlett & Mendelson, 2002) . Based on this complexity measure, we have the following bound on the expected classification error.

Theorem 4.2. (Bartlett & Mendelson, 2002) For each F mapping S into X l-1 , let F l,A = f : x → w, φ (l) F (x) H l + b w H l ≤ A, b ∈ R . Fix A,

γ > 0, with probability at least 1 -δ and for any N ∈ N, every function f (l) in F l,A satisfies

P (yf (l) (x) ≤ 0) ≤ Rl (f (l) ) + 2G N (F l,A ) + 8 γ + 1 log(4/δ) 2N ,

where Rl (f

(l) ) = 1 N N n=1 max(0, 1 -y n f (l) (x n )/γ

), the empirical hinge loss. Given the assumptions on k (l) , for any F , we have

G N (F l,A ) ≤ 2A c N .

Without loss of generality, we shall set hyperparameter γ = 1. We now characterize S l-1 . Note that for a given f (l) , A = w f (l) H l is the smallest nonnegative real number such that f (l) ∈ F l,A and it is immediate that this gives the tightest bound in Theorem 4.2. Let κ = 1 N N n=1 1 {yn=+} . Lemma 4.3 (optimal S l-1 ). Given a learning paradigm minimizing Rl (f (l) ) + τ w f (l) H l using representation S l-1 = F (S), where τ is any positive constant satisfying τ < 2(c -a) min(κ, 1κ). Denote as S l-1 any representation satisfying

k (l) (F (x + ), F (x -)) = a and k (l) (F (x), F (x )) = c (1)

for all pairs of x + , x -from distinct classes in S and all pairs of x, x from the same class. Suppose the learning paradigm returns f (l) under this representation. Let S • l-1 be another representation under which the learning paradigm returns f (l)• . If f (l)• achieves zero hinge loss on at least one example from each class, then for any

N ∈ N, Rl (f (l) ) + τ w f (l) H l ≤ Rl (f (l)• ) + τ w f (l)• H l .

The optimal representation S l-1 , characterized by Eq. 1, enjoys a straightforward geometric interpretation: Examples from distinct classes are as distant as possible in the RKHS whereas examples from the same class are as concentrated as possible (see proof (C) of Lemma 4.3 for a rigorous justification). Intuitively, it is easy to see that such a representation is the "easiest" for the classifier.

The conditions in Eq. 1 can be concisely summarized in an ideal kernel matrix G defined as

(G ) mn = a, if y m = y n ; (G ) mn = c, if y m = y n .

And to have the l -1 th layer learn S l-1 , it suffices to train it to minimize some dissimilarity measure between G and the kernel matrix computed from k (l) and F (l-1) (S), which we denote G l-1 . Empirical alignment (Cristianini et al., 2002) , L 1 and L 2 distances between matrices can all serve as the dissimilarity measure. To simplify discussion, we let the dissimilarity measure be the

L 1 distance l-1 (F (l-1) , (x m , y m ), (x n , y n )) = |(G ) mn -(G l-1 ) mn |.

This specifies Rl-1 (F (l-1) ) as the sample mean of ( l-1 (F (l-1) , (x m , y m ), (x n , y n ))) N m,n=1 and R l-1 as the expectation of l-1 over (X 1 , Y )×(X 1 , Y ). Note that due to the boundedness assumption on k (l) , l-1 ≤ 2 max(|c|, |a|).

4.2.3. FORMAL RESULTS: S:

l-2 , . . . , S 1 Similar to Section 4.2.2, we first need to approximate R l-1 . Lemma 4.4. For j = 1, 2, . . . , d l-1 , let f (l-1) j ∈ F l-1 , where F l-1 is a given hypothesis class.

There exists an absolute constant C > 0 such that for any N ∈ N, with probability at least

1 -δ, R l-1 (F (l-1) ) ≤ Rl-1 (F (l-1) ) + 4L (l) Cd l-1 max(|c|, |a|) G N (F l-1 ) + 8 log(2/δ) N .

We are now in a position to characterize S l-2 . For the following lemma only, we further assume that k (l-1) (x, y), as a function of (x, y), depends only on and strictly decreases in x -y 2 for all x, y ∈ X l-1 with k (l-1) (x, y) > a, and that the infimum

inf x,y∈X l-1 k (l-1) (x, y) = a is attained in X l-1 at all x, y with x -y 2 ≥ η. Also assume that inf x,y∈X l-1 ; x-y 2<η ∂k (l-1) (x, y)/∂ x - y 2 = ι (l-1

) is defined and is positive.

Consider an F mapping S into X l-1 , let F l-1,A = f : x → w, φ (l-1) F (x) H l-1 + b w H l-1 ≤ A, b ∈ R . For a given F (l-1) = f (l-1) 1 , . . . , f (l-1) d l-1

, it is immediate that

A = max j w f (l-1) j H l-1

is the smallest nonnegative real number such that f (l-1) j ∈ F l-1,A for all j, giving the tightest bound in Lemma 4.4 (recall the bound on

G N (F l-1,A ) in Theorem 4.2). Let ψ = N m,n=1 1 {ym =yn} /N 2 .

Lemma 4.5 (optimal S l-2 ). Given a learning paradigm minimizing Rl-1 (F (l-1) ) +

τ max j w f (l-1) j H l-1 using representation S l-2 = F (S)

, where τ is any positive constant satisfying τ < 2d l-1 (c -a)ψι (l-1) . Denote as S l-2 any representation satisfying

k (l-1) (F (x + ), F (x -)) = a and k (l-1) (F (x), F (x )) = c

for all pairs of x + , x -from distinct classes in S and all pairs of x, x from the same class. Suppose the learning paradigm returns

F (l-1) = f (l-1) 1 , . . . , f (l-1) d l-1

under this representation. Let S • l-2 be another representation under which the learning paradigm returns F (l-1)• . If F (l-1)• achieves zero loss on at least one pair of examples from distinct classes, then for any

N ∈ N, Rl-1 (F (l-1) ) + τ max j w f (l-1) j H l-1 ≤ Rl-1 (F (l-1)• ) + τ max j w f (l-1)• j H l-1 .

Applying this analysis to the rest of the hidden layers, it is evident that the i th layer, i = 1, 2, ..., l -1, should be trained to minimize the difference between G and the kernel matrix computed with k (i+1) and F (i) (S), denoted G i . Generalizing to classification with more than two classes requires no change to the algorithm since the definition of G is agnostic to the number of classes involved in the classification task. Also note that the sufficiency of expanding the kernel machines of each layer on the training sample (see Section 2) for the learning objectives in Lemma 4.3 and Lemma 4.5 is trivially justified since the generalized representer theorem directly applies (Schölkopf et al., 2001) . Now since the optimal representation is consistent across layers, the dynamics of layer-wise learning in a kMLP is clear: The network maps the random sample sequentially through layers, with each layer trying to map examples from distinct classes as far as possible in the RKHS while keeping examples from the same class in a cluster as concentrated as possible. In other words, each layer learns a more separable representation of the sample. Eventually, the output layer works as a classifier on the final representation and since the representation would be "simple" after the mappings of the lower layers, the learned decision boundary would generalize better to unseen data, as suggested by the bounds above.

4.3. LEARNING NETWORK-WISE OPTIMALITY:

We now discuss how to design a layer-wise learning algorithm that learns network-wise optimality at each layer. A rigorous description of the problem of layer-wise optimality drifting away from network-wise optimality and the search for a solution begins with the following bound on the total error of any two consecutive layers in a kMLP. Lemma 4.6. For any i = 2, . . . , l, let the target function and the approximation function be T i , F (i) * :

X i → R di , respectively. Let i = F (i) * -T i s := sup x∈Xi F (i) * (x) -T i (x) 2 , we have F (i) * • F (i-1) * -T i • T i-1 s ≤ i + √ i-1 2L (i) di j=1 w f (i) * j 2 Hi .

(2) By applying the above bound sequentially from the input layer to the output, we can decompose the error of an arbitrary kMLP into the error of the layers. This result gives a formal description of the problem: The hypothesis with the minimal norm minimizes the propagated error from upstream, but evidently, this hypothesis is not necessarily close to the layer-wise optimality T i .

Moreover, this bound provides the insight needed for learning network-wise optimality individually at each layer: For the i th layer, i ≥ 2, searching for network-wise optimality amounts to minimizing the r.h.s. of Eq. 2. Lemma 4.3 and Lemma 4.5 characterized T i for i < l and learning objectives that bound i were provided earlier in the text accordingly. Based on those results, the solution that minimizes the new learning objective Ri (F (i) ) + τ max j w f (i) j

Hi:

, where τ > 0 is a hyperparameter, provides a good approximate to the minimizer of the r.h.s. of Eq. 2 if, of course, τ is chosen well. Thus, taking this as the learning objective of the i th layer produces a layer-wise algorithm that learns network-wise optimality at this layer. Note that for BP, one usually also needs to heuristically tune the regularization coefficient for weights as a hyperparameter.

4.4. ACCELERATING THE UPPER LAYERS:

There is a natural method to accelerate the upper layers (all but the input layer): The optimal representation F (S) is sparse in the sense that φ(F

(x m )) = φ(F (x n )) if y m = y n and φ(F (x m )) = φ(F (x n )) if y m = y n (see the proof (C) of Lemma 4.3).

Since a kernel machine built on this representation of the given sample is a function in the RKHS that is contained in the span of the image of the sample, retaining only one example from each class would result in exactly the same hypothesis class because trivially, we have

{ N n=1 α n φ(F (x n ))|α n ∈ R} = {α + φ(F (x + )) + α -φ(F (x -))|α + , α -∈ R} for arbitrary x + ,

x -in S when the representation F (S) is optimal. Thus, after training a given layer, depending on how close the actual kernel matrix is to the ideal one, one can (even randomly) discard a large portion of centers for kernel machines of the next layer to speed up the training of it without sacrificing performance. As we will later show in the experiments, randomly keeping a fraction of the training sample as centers for upper layers produces performance comparable to or better than that obtained with using the entire training set.

5. RELATED WORKS:

The idea of combining connectionism with kernel method was initiated by Cho & Saul (2009) . In their work, an "arc cosine" kernel was so defined as to imitate the computations performed by a one-layer MLP. Zhuang et al. (2011) extended the idea to arbitrary kernels with a focus on MKL, using an architecture similar to a two-layer kMLP. As a further generalization, Zhang et al. (2017) independently proposed kMLP and the KN equivalent of CNN. However, they did not extend the idea to any arbitrary NN. Scardapane et al. (2017) proposed to reparameterize each nonlinearity in an NN with a kernel expansion, resulting in a network similar to KN but is trained with BP. There are other works aiming at building "deep" kernels using approaches that are different in spirit from those above. Wilson et al. (2016) proposed to learn the covariance matrix of a Gaussian process using an NN in order to make the kernel "adaptive". This idea also underlies the now standard approach of combining a deep NN with an SVM for classification, which was first explored by Huang & LeCun (2006) and Tang (2013) . Such an interpretation can be given to KNs as well, as we point out in Appendix B.5. Mairal et al. (2014) proposed to learn hierarchical representations by learning mappings of kernels that are invariant to irrelevant variations in images.

Much works have been done to improve or substitute BP in learning a deep architecture. Most aim at improving the classical method, working as add-ons for BP. The most notable ones are perhaps the unsupervised greedy pre-training techniques proposed by Hinton et al. (2006) and Bengio et al. (2007) . Among works that try to completely substitute BP, none provided a comparable optimality guarantee in theory as that given by BP. Fahlman & Lebiere (1990) pioneered the idea of greedily learn the architecture of an NN. In their work, each new node is added to maximize the correlation between its output and the residual error signal. Several authors explored the idea of approximating error signals propagated by BP locally at each layer or each node (Bengio, 2014; Carreira-Perpinan & Wang, 2014; Lee et al., 2015; Balduzzi et al., 2015; Jaderberg et al., 2016) . Kulkarni & Karande (2017) proposed to train NN layer-wise using an ideal kernel matrix that is a special case of that in our work. No theoretical results were provided to justify its optimality for NN. Zhou & Feng (2017) proposed a BP-free deep architecture based on decision trees, but the idea is very different from ours. Raghu et al. (2017) attempted to quantify the quality of hidden representations toward learning more interpretable deep architectures, sharing a motivation similar to ours.

6. EXPERIMENTS:

We compared kMLP learned using the proposed greedy algorithm with other popular deep architectures including MLP, Deep Belief Network (DBN) (Hinton & Salakhutdinov, 2006) and Stacked Autoencoder (SAE) (Vincent et al., 2010) , with the last two trained using a combination of unsupervised greedy pre-training and standard BP (Hinton et al., 2006; Bengio et al., 2007) . Note that we only focused on comparing with the standard, generic architectures because kMLP, as the KN equivalent of MLP, does not have a specialized architecture or features designed for specific application domains. Several optimization and training techniques were applied to the MLPs to boost performance. These include Adam (Kingma & Ba, 2014), RMSProp (Tieleman & Hinton, 2012) , dropout (Srivastava et al., 2014) and batch normalization (BN) (Ioffe & Szegedy, 2015) . kMLP accelerated using the proposed method (kMLP FAST ) was also compared. For these models, we randomly retained a subset of the centers of each upper layer before its training. As for the benchmarks used, rectangles, rectangles-image and convex are binary classification datasets, mnist (10k) and mnist (10k) rotated are variants of MNIST (Larochelle et al., 2007; LeCun et al., 2010) .

And fashion-mnist is the Fashion-MNIST dataset (Xiao et al., 2017) . To further test the proposed layer-wise learning algorithm and the acceleration method, we compared greedily-trained kMLP with MLP and kMLP trained using BP (Zhang et al., 2017) using the standard MNIST (LeCun et al., 2010) . Two popular acceleration methods for kernel machines were also compared on the same benchmark, including using a parametric representation (i.e., for each node in a kMLP, f (x) = k(w, x), w learnable) (kMLP PARAM ) and using random Fourier features (kMLP RFF ) (Rahimi & Recht, 2008) . More details for the experiments can be found in Appendix A4 .

From Table 1a , we see that the performance of kMLP is on par with some of the most popular and most mature deep architectures. In particular, the greedily-trained kMLPs compared favorably with their direct NN equivalents, i.e., the MLPs, even though neither batch normalization nor dropout was used for the former. These results also validate our earlier theoretical results on the layer-wise learning algorithm, showing that it indeed has the potential to be a substitute for BP with an equivalent optimality guarantee. Results in Table 1b further demonstrate the effectiveness of the greedy learning scheme. For both the single-hidden-layer and the two-hidden-layer kMLPs, the layer-wise algorithm consistently outperformed BP. It is worth noting that the proposed acceleration trick, despite being extremely simple, is clearly very effective and even produced models outperforming the original ones. This shows that kMLP together with the greedy learning scheme can be of practical interest even when dealing with the massive data sets in today's machine learning.

Last but not least, we argue that it is the practical aspects that makes the greedy learning framework promising. Namely, this framework of learning makes deep architectures more transparent and intuitive, which can serve as a tentative step toward more interpretable, easy-to-understand models with strong expressive power. Also, new design paradigms are now possible under the layer-wise framework. For example, each layer can now be "debugged" individually. Moreover, since learning becomes increasingly simple for the upper layers as the representations become more and more well-behaved, these layers are usually very easy to set up and also converge very fast during training.

each hidden layer, respectively. For DBN-3 and SAE-3, the sizes of the three hidden layers varied in intervals [500, 3000], [500, 4000] and [1000, 6000], respectively. DBN-1 used a much larger hidden layer than DBN-3 to obtain comparable performance. A simple calculation shows that the total numbers of parameters in the kMLPs were fewer than those in the corresponding DBNs and SAEs by orders of magnitude in all experiments. Like in the training for the kMLPs, the validation set were also reserved for early-stopping in final training. The DBNs and SAEs had been pre-trained unsupervisedly before the supervised training phase, following the algorithms described in (Hinton et al., 2006; Bengio et al., 2007) . More detailed settings for these models were reported in (Larochelle et al., 2007) .

APPENDIX B FURTHER ANALYSIS:

In this section, we provide some further analysis on kMLP and the layer-wise learning algorithm. Namely, in Appendix B.1, we give a bound on the Gaussian complexity of an l-layer kMLP, which describes the intrinsic model complexity of kMLP. In particular, the bound describes the relationship between the depth/width of the model and the complexity of its hypothesis class, providing useful information for model selection. In Appendix B.2, we give a constructive result stating that the dissimilarity measure being optimized at each hidden layer will not increase as training proceeds from the input layer to the output. This also implies that a deeper kMLP performs at least as well as its shallower counterparts in minimizing any loss function they are trained on. In Appendix B.3, a result similar to Lemma 4.3 is provided, stating that the characterization for the optimal representation can be made much simpler if one uses a more restricted learning paradigm. In fact, in contrast to Lemma 4.3, both necessary and sufficient conditions can be determined under the more restricted setting. In Appendix B.4, we provide a standard, generic method to estimate the Lipschitz constant of a continuously differentiable kernel, as this quantity has been repeatedly used in many of our results in this paper. In Appendix B.5, we state some advantages of kMLP over classical kernel machines. In particular, empirical results are provided in Appendix B.5.1, in which a two-layer kMLP consistently outperforms the classical Support Vector Machine (SVM) (Cortes & Vapnik, 1995) as well as several SVMs enhanced by Multiple Kernel Learning (MKL) algorithms (Bach et al., 2004; Gönen & Alpaydın, 2011) .

B.1 GAUSSIAN COMPLEXITY OF KMLP:

We first give a result on the Gaussian complexity of a two-layer kMLP.

Lemma B.1. Given kernel k : X 2 × X 2 → R, where X 2 ⊂ R d1 . Let F 2 = {f : X 2 → R, f (x) = m ν=1 α ν k(x ν , x) + b | α = (α 1 , . . . , α m ) ∈ R m , α 1 ≤ A, b ∈ R}, where the x ν are arbitrary examples from X 2 . Consider F 1 = {(f 1 , . . . , f d1 ) : X 1 → X 2 | f j ∈ Ω}, where Ω is a given hypothesis class that is closed under negation, i.e., if f ∈ Ω, then -f ∈ Ω. Define F 2 • F 1 = {h : x → m ν=1 α ν k(F (x ν ), F (x)) + b | α 1 ≤ A, b ∈ R, F ∈ F 1 }.

If the range of some element in Ω contains 0, we have

G N (F 2 • F 1 ) ≤ ALd 1 G N (Ω).

The above result can be easily generalized to kMLP with an arbitrary number of layers.

Lemma B.2. Given an l-layer kMLP, for each f

(i) j (x) = m ν=1 α (i) νj k (i) (F (i-1) (x ν ), F (i-1) (x))+b, α (i) j := (α (i) 1j , . . . , α (i) mj ) ∈ R m , b ∈ R, assume α (i) j 1 ≤ A i and let d l = 1.

Denote the class of functions implemented by this kMLP as F l , we have

G N (F l ) ≤ d 1 l i=2 A i L (i) d i G N (Ω).

Proof. It is trivial to check that the hypothesis class of each layer is closed under negation and that there exists a function in each of these hypothesis classes whose range contains 0. Then the result follows from repeatedly applying Lemma B.1.

B.2 NONINCREASING LOSS ACROSS LAYERS:

Lemma B.3. For i ≥ 2, assume k (i) is PD and fix layers 1, 2, . . . , i -1 at arbitrary states F (1) , F (2) , . . . , F (i-1) . Let the loss functions i , i-1 be the same up to their domains, and denote both i and i-1 as . Suppose layer i is trained with a gradient-based algorithm to minimize the loss (F (i) ). Denote the state of layer i after training by F (i) * . If d i-1 = d i , there is an initialization for F (i) such that (F (i) * ) ≤ (F (i-1) ).

(3)

Calculation for this initialization is specified in the proof.

For i = 1, under the further assumption that X 1 ⊂ R d0 and d 0 = d 1 , Eq. 3 becomes (F (1) * ) ≤ (F (0) ), where F (0) is the identity map on X 1 .

Remark B.3.1. For the greedily-trained kMLP, Lemma B.3 applies to the hidden layers and implicitly requires that k (i+1) = k (i) since the loss function for layer i, when viewed as a function of F (i) , takes the form (k (i+1) (F (i) )) and can be rewritten as

k (i+1) (F (i) ). Similarly, (k (i) (F (i-1) )) is k (i) (F (i-1)

). Since Lemma B.3 assumes to be the same across layers (otherwise it does not make sense to compare between layers), this forces

k (i+1) = k (i) . Then it follows that k (i+1) = k (i) .

Further, if k (i+1) and k (i) have the property that k(x, y) = k(x, ȳ), where x, ȳ denote the images of x, y under an embedding of R p into R q (p ≤ q) defined by the identity map onto a p-dimensional subspace of R q , then the condition

d i-1 = d i can be relaxed to d i-1 ≤ d i .

This lemma states that for a given kMLP, when it has been trained upto and including the i th hidden layer, the i + 1 th hidden layer can be initialized in such a way that the value of its loss function will be lower than or equal to that of the i th hidden layer after training. In particular, the actual hidden representation "converges" to the optimal represetation as training proceeds across layers. On the other hand, when comparing two kMLPs, this result implies that the deeper kMLP will not perform worse in minimizing the loss function than its shallower counterpart.

In deep learning literature, results analogous to Lemma B.3 generally state that in the hypothesis class of a NN with more layers, there exists a hypothesis that approximates the target function nontrivially better than any hypothesis in that of another shallower network (Sun et al., 2016) . Such an existence result for kMLP can be easily deduced from the earlier bound on its Gaussian complexity (see Lemma B.2). However, these proofs of existence do not guarantee that such a hypothesis can always be found through learning in practice, whereas Lemma B.3 is constructive in this regard. Nevertheless, one should note that Lemma B.3 does not address the risk R = E . Instead, it serves as a handy result that guarantees fast convergence of upper layers during training in practice.

B.3 SIMPLER OPTIMAL REPRESENTATION UNDER MORE RESTRICTED LEARNING PARADIGMS:

The following lemma states that if we are willing to settle with a more restricted learning paradigm, the necessary and sufficient condition that guarantees the optimality of a representation can be characterized and is simpler than that described in Lemma 4.3. The setup for this lemma is the same as that of Lemma 4.3 except that the assumption that the numbers of examples from the two classes are equal is not needed.

Lemma B.4. Consider a learning paradigm that minimizes Rl (f (l) ) + τ G N (F l,A ) using representation S l-1 = F (S) under the constraint that it returns an f (l) ∈ F l,A with Rl (f (l) ) = 0, where τ is any positive constant. For any N ∈ N, Rl (f (l) ) + τ G N (F l,A ) is minimized over all linearly separable representations if and only if the representation F (S) satisfies

k (l) (F (x + ), F (x -)) = a

for all pairs of x + , x -from distinct classes in S.

B.4 LIPSCHITZ CONSTANT FOR CONTINUOUSLY DIFFERENTIABLE KERNELS:

In general, for a continuously differentiable function f : R → R with derivative f and any a, b ∈ R, a < b, we have

|f (b) -f (a)| = b a f (x)dx ≤ b a |f (x)|dx ≤ max a≤x≤b |f (x)|(b -a).

This simple result can be used to bound the Lipschitz constant of a continuously differentiable kernel. For example, for Gaussian kernel k : X × X → R, X ⊂ R, k(x, y) = e -(x-y) 2 /σ 2 , we have ∂k(x, y)/∂y = 2(x -y)k(x, y)/σ 2 . Hence for each fixed x ∈ X, k(x, y) is Lipschitz in y with Lipschitz constant bounded by sup y∈X 2(x -y)k(x, y)/σ 2 . In practice, X is always compact and can be a rather small subspace of some Euclidean space after normalization of data, hence this would provide a reasonable approximation to the Lipschitz constant of Gaussian kernel.

B.5 COMPARING KMLP WITH CLASSICAL KERNEL MACHINES:

There are mainly two issues with classical kernel machines besides their usually high computational complexity. First, despite the fact that under mild conditions, they are capable of universal function approximation (Park & Sandberg, 1991; Micchelli et al., 2006) and that they enjoy a very solid mathematical foundation (Aronszajn, 1950) , kernel machines are unable to learn multiple levels of distributed representations (Bengio et al., 2013) , yet learning representations of this nature is considered to be crucial for complicated artificial intelligence (AI) tasks such as computer vision, natural language processing, etc. (Bengio, 2009; LeCun et al., 2015) . Second, in practice, performance of a kernel machine is usually highly dependent on the choice of kernel since it governs the quality of the accessible hypothesis class. But few rules or good heuristics exist for this topic due to its extremely task-dependent nature. Existing solutions such as MKL (Bach et al., 2004; Gönen & Alpaydın, 2011) view the task of learning an ideal kernel for the given problem to be separate from the problem itself, necessitating either designing an ad hoc kernel or fitting an extra trainable model on a set of generic base kernels, complicating training. kMLP learns distributed, hierarchical representations because it inherits the architecture of MLP. To be specific, first, we see easily that the hidden activation of each layer, i.e., F (i) (x) ⊂ R di , is a distributed representation (Hinton, 1984; Bengio et al., 2013) . Indeed, just like in an MLP, each layer of a kMLP consists of an array of identical computing units (kernel machines) that can be activated independently. Further, since each layer in a kMLP is built on top of the previous layer in exactly the same way as how the layers are composed in an MLP, the hidden representations are hierarchical (Bengio et al., 2013) .

Second, kMLP naturally combines the problem of learning an ideal kernel for a given task and the problem of learning the parameters of its kernel machines to accomplish that task. To be specific, kMLP performs nonparametric kernel learning alongside learning to perform the given task. Indeed, for kMLP, to build the network one only needs generic kernels, but each layer F (i) can be viewed as a part of a kernel of the form k (i+1) (F (i) (x), F (i) (y)). The fact that each F (i) is learnable makes this kernel "adaptive", mitigating to some extent any limitation of the fixed generic kernel k (i+1) . The training of layer i makes this adaptive kernel optimal as a constituent part of layer i + 1 for the task the network was trained for. And it is always a valid kernel if the generic kernel k (i+1) is. Note that this interpretation has been given in a different context by Huang & LeCun (2006) and Bengio et al. (2013) , we include it here only for completeness.

B.5.1 EMPIRICAL RESULTS:

We now compare a single-hidden-layer kMLP using simple, generic kernels with SVMs enhanced by MKL algorithms that used significantly more kernels to demonstrate the ability of kMLP to automatically learn task-specific kernels out of standard ones. The standard SVM and seven other SVMs enhanced by popular MKL methods were compared (Zhuang et al., 2011) , including the classical convex MKL (Lanckriet et al., 2004) with kernels learned using the extended level method proposed in (Xu et al., 2009 ) (MKL LEVEL ); MKL with L p norm regularization over kernel weights (Kloft et al., 2011) (L p MKL), for which the cutting plane algorithm with second order Taylor approximation of Table 2 : Average test error (%) and standard deviation (%) from 20 runs. Results with overlapping 95% confidence intervals (not shown) are considered equivalent. Best results are marked in bold. The average ranks (calculated using average test error) are provided in the bottom row. When computing confidence intervals, due to the limited sizes of the data sets, we pooled the 20 random samples. For kMLP, all results were achieved using a greedily-trained, one-hidden-layer model with the number of kernel machines ranging from 3 to 10 on the first layer for different data sets. The second layer was a single kernel machine. All kernel machines within one layer used the same Gaussian kernel, and the two kernels on the two layers differed only in kernel width σ. All hyperparameters were chosen via 5-fold cross-validation. As for the other models compared, for each data set, SVM used a Gaussian kernel. For the MKL algorithms, the base kernels contained Gaussian kernels with 10 different widths on all features and on each single feature and polynomial kernels of degree 1 to 3 on all features and on each single feature. For 2LMKL INF , one Gaussian kernel was added to the base kernels at each iteration. Each base kernel matrix was normalized to unit trace. For L p MKL, p was selected from {2, 3, 4}. For MKM, the degree parameter was chosen from {0, 1, 2}. All hyperparameters were selected via 5-fold cross-validation. From Table 2 , kMLP compares favorably with other models, which validates our claim that kMLP learns its own kernels nonparametrically hence can work well even without excessive kernel parameterization. Performance difference among models can be small for some data sets, which is expected since they are all rather small in size and not too challenging. Nevertheless, it is worth noting that only 2 Gaussian kernels were used for kMLP, whereas all other models except for SVM used significantly more kernels.

APPENDIX C PROOFS:

Proof of Lemma 4.3. Throughout this proof we shall drop the layer index l for brevity. Given that the representation satisfies Eq. 1, the idea is to first collect enough information about the returned f = (w , b ) such that we can compute R(f ) + τ w H and then show that for any other F (S) satisfying the condition in the lemma, suppose the learning paradigm returns f = (w , b ) ∈ F A , then R(f ) + τ w H ≥ R(f ) + τ w H . We now start the formal proof.

First, note that in the optimal representation, i.e., an F (S) such that Eq. 1 holds, it is easy to see that φ(F (x -)) -φ(F (x + )) H is maximized over all representations for all x -, x + .

Moreover, note that given the representation is optimal, we have φ(F (x)) = φ(F (x )) if y = y and φ(F (x)) = φ(F (x )) if y = y : Indeed, by Cauchy-Schwarz inequality, for all x, x ∈ S, k(F (x), F (x )) = φ(F (x)), φ(F (x )) H ≤ φ(F (x)) H φ(F (x )) H and the equality holds if where α ∈ A, the (x n , x n ) are a sample of size N from X 1 × X 2 and g 1 , . . . , g N , g 1,1 , . . . , g 2N,d are i.i.d. standard normal random variables.

Let arbitrary α, β ∈ A be given, define

X α -X β 2 2 = E(X α -X β ) 2

, where the expectation is taken over the g n . Define Y α -Y β 2 2 similarly and we have

X α -X β 2 2 = N n=1 ω(f α,1 (x n ), . . . , f α,d (x n ), y n ) -ω(f β,1 (x n ), . . . , f β,d (x n ), y n ) 2 ≤ L 2 N n=1 d i=1 (f α,i (x n ) -f β,i (x n )) 2 + (f α,i (x n ) -f β,i (x n )) 2 = Y α -Y β 2 2

. By Slepian's lemma (Pisier, 1999) and since the F i are closed under negation,

N 2 ĜN (ω • F) = E gn sup α∈A X α ≤ 2E gn,i,g N +n,i sup α∈A Y α ≤ N 2 2L d i=1 ( Ĝ1 N (F i ) + Ĝ2 N (F i )).

Taking the expectation of the x n and x n on both sides, we have

N 2 G N (ω • F) ≤ N 2 2L d i=1 (G 1 N (F i ) + G 2 N (F i )).

Proof of Lemma 4.4. Normalize l-1 to [0, 1] by dividing 2 max(|c|, |a|). Then the loss function becomes l-1 F (l-1) , (x m , y m ), (x n , y n ) = 1 2 max(|c|, |a|) k (l) F (l-1) (x m ), F (l-1) (x n ) -(G ) mn .

For each fixed (G ) mn , l-1 F (l-1) , (x m , y m ), (x n , y n )l-1 F (l-1) , (x m , y m ), (x n , y n )

≤ 1 2 max(|c|, |a|) k (l) F (l-1) (x m ), F (l-1) (x n ) -k (l) F (l-1) (x m ), F (l-1) (x n ) ≤ 1 2 max(|c|, |a|) k (l) F (l-1) (x m ), F (l-1) (x n ) -k (l) F (l-1) (x m ), F (l-1) (x n ) + k (l) F (l-1) (x m ), F (l-1) (x n ) -k (l) F (l-1) (x m ), F (l-1) (x n ) ≤ L (l) 2 max(|c|, |a|) F (l-1) (x n ) -F (l-1) (x n ) 2 + F (l-1) (x m ) -F (l-1) (x m ) 2 ≤ L (l) max(|c|, |a|) F (l-1) (x n ) -F (l-1) (x n ), F (l-1) (x m ) -F (l-1) (x m ) 2 . Hence l-1 is L (l) / max(|c|, |a|)-Lipschitz in (F (l-1) (x m ), F (l-1) (x n ))

with respect to the Euclidean metric on R 2d l-1 for each (G ) mn .

The result follows from Lemma C.1 with m = 2, d = d l-1 and Corollary 15 in (Bartlett & Mendelson, 2002) .

Proof of Lemma 4.5. This proof uses essentially the same idea as that of Lemma 4.3. Due to the complete dependence of k (l) (x, y) on x -y 2 , we can rewrite k (l) F (l-1) (F (x m )), F (l-1) (F (x n ))

as h (l) F (l-1) (F (x m )) -F (l-1) (F (x n )) 2 for some h (l) . Define µ

(l-1) mn = φ (l-1) (F (x m )) - φ (l-1) (F (x n )) H l-1 , we have Rl-1 (F (l-1) ) = 1 N 2 N m,n=1 h (l) F (l-1) (F (x m )) -F (l-1) (F (x n )) 2 -(G ) mn = 1 N 2 N m,n=1 h (l) d l-1 j=1 f (l-1) j (F (x m )) -f (l-1) j (F (x n )) 2 -(G ) mn = 1 N 2 N m,n=1 h (l) d l-1 j=1 µ (l-1) mn w f (l-1) j H l-1 cos θ F,w f (l-1) j 2 -(G ) mn .

Given that the representation is optimal, we have

Rl-1 (F (l-1) ) + τ max 1≤j≤d l-1 w f (l-1) j H l-1 = 1 N 2 ym =yn h (l) d l-1 j=1 µ (l-1) mn w f (l-1) j H l-1 cos θ F,w f (l-1) j 2 -a + τ max 1≤j≤d l-1 w f (l-1) j H l-1 .

The returned F (l-1) must satisfy

d l-1 j=1 µ (l-1) mn w f (l-1) j H l-1 cos θ F,w f (l-1) j 2 ≤ η; w f (l-1) 1 H l-1 = • • • = w f (l-1) d l-1 H l-1 =: w H l-1 ; cos θ F,w = 1.

The first observation is trivial since if otherwise, one can always reduce the largest w f (l-1)

j H l-1

to obtain equality to η, this gives the same Rl-1 (F (l-1) ) with a smaller τ max j w f (l-1) j H l-1 . Note that if during shrinking the largest w f (l-1) j H l-1 , this element ceases to be the largest among all j, we shall continue the process with the new largest instead. To see the rest of the claim, note that for the largest of the w f (l-1) j H l-1 , we must have θ F,w f (l-1) j = 0 since if not, one could shift and rotate the hyperplane and again obtain the same Rl-1 (F (l-1) ) with a smaller τ max j w f (l-1) j H l-1 . Reducing the largest w f (l-1) j H l-1 and increasing the second largest by scaling, one would get a smaller risk. It is immediate that the minimal risk (w.r.t. only the first and second largest w f (l-1) j H l-1 ) is attained when the largest and the second largest w f (l-1) j H l-1 are equal. Then a similar argument as before gives θ F,w f (l-1) j = 0 for both of them. The rest of the claim follows via repeatedly applying this argument to all the w f (l-1)

j H l-1 . Define |f (F (x + )) -f (F (x -))| = φ (l-1) (F (x + )) -φ (l-1) (F (x -)) H l-1 w H l-1 .

Then we have:

Rl-1 (F (l-1) ) + τ max

1≤j≤d l-1 w f (l-1) j H l-1 = ψ h (l) d l-1 |f (F (x + )) -f (F (x -))| -a + τ |f (F (x + )) -f (F (x -))| φ (l-1) (F (x + )) -φ (l-1) (F (x -)) H l-1 = ψ h (l) d l-1 |f (F (x + )) -f (F (x -))| -a + τ |f (F (x + )) -f (F (x -))| 2(c -a) .

As we have shown,

d l-1 |f (F (x + )) -f (F (x -))| ∈ [0, η]. Let λ = |f (F (x + )) -f (F (x -))|

and differentiate the r.h.s. of the above equation w.r.t. λ and using the assumption on τ , we have

d l-1 ψ dh (l) (t) dt + τ 2(c -a) < 0.

Hence the overall risk decreases in λ over [0, η/ d l-1 ], which implies that the returned F (l-1) must have λ = η/ d l-1 and

Rl-1 (F (l-1) ) + τ max

1≤j≤d l-1 w f (l-1) j H l-1 = τ η 2d l-1 (c -a) .

Now for any other representation S l-2 = F (S), suppose the learning paradigm returns F (l-1) .

Assume without loss of generality that the largest f

(l-1) j (F (x + )) -f (l-1) j (F (x -)

) over all j, x + and x -is attained at j = 1 and write w

= w f (l-1) 1 , f = f (l-1) 1 for convenience. Let x w + , x w - be the pair with the largest f (F (x + )) -f (F (x -)). Note that the assumption on F implies that f (F (x + )) -f (F (x -)) ≥ η/ d l-1 . Then we have Rl-1 (F (l-1) ) + τ max 1≤j≤d l-1 w f (l-1) j H l-1 ≥ τ w H l-1 = τ |f (F (x + )) -f (F (x -))| φ (l-1) (F (x + )) -φ (l-1) (F (x -)) H l-1 | cos θ F ,w | ≥ τ η 2d l-1 (c -a) = Rl-1 (F (l-1) ) + τ max 1≤j≤d l-1 w f (l-1) j H l-1 ,

proving the lemma.

Proof of Lemma 4.6. First, it is trivial that the so-defined s metric is indeed a metric. In particular, it satisfies the triangle inequality. For i = 2, . . . , l,

F (i) * • F (i-1) * -T i • T i-1 s ≤ F (i) * • F (i-1) * -F (i) * • T i-1 s + F (i) * • T i-1 -T i • T i-1 s ≤ sup x∈Xi di j=1 f (i) * j • F (i-1) * (x) -f (i) * j • T i-1 (x) 2 + i = sup x∈Xi di j=1 f (i) * j

(λ), k (i) (λ, F (i-1) * (x)) -k (i) (λ, T i-1 (x)) 

F (i-1) * (x) + L (i) Ti-1(x)

F (i-1) * (x) -T i-1 (x) Taking expectation with respect to the x n finishes the proof.

Proof of Lemma B.3. For i ≥ 1, given the assumption that d i-1 = d i , we first show that one can find an initialization for layer i, denoted F (i)• , such that F (i)• (S) = F (i-1) (S).

Without loss of generality, assume i = 1, then it amounts to showing that a set of parameters can be found for the kernel machines f

(1) 1 , . . . , f

(1) d1 of F (1) such that when these kernel machines are initialized with this certain set of parameters, this layer agrees with the identity function when restricted to S. Recall that f nm , G is the N × N kernel matrix defined as (G) nm = k (1) (x n , x m ), D is a d 0 × N matrix whose n th column is the n th example x n = (x n1 , x n2 , . . . , x nd0 ) for n = 1, 2, . . . , N . It is straightforward that the set of parameters A initializes F (1) to the desired F (1)• . Moreover, A can always be solved in closed-form by inverting G since k (1) being PD guarantees that G is invertible.

Combining the above result and the defining property of gradient of the loss function with respect to the parameters, i.e., that gradient points at the direction of the greatest rate of increase of the loss function, we then have (F (0) ) = (F (1)• ) ≥ (F (1) * ).

Proof of Lemma B.4. In this proof we shall again drop the layer index l for brevity.

Identify f ∈ F A using w, b upto R-scalar multiplication. Fix an F such that φ(F (S)) is linearly separable (otherwise the learning paradigm is not applicable).

Let x + , x -be the closest pair of examples from distinct classes in H. For a given f = (w, b), let x w -, x w + be a pair of examples from different classes that are closest to f in H in the sense that f (F (x w + )), f (F (x w -)) are the smallest in absolute value among all x + , x -, respectively. Denote the canonical form of f with respect to φ(F (S)) with equality attained by at least one of x w + , x w -as f F = (w F , b F ), i.e., f F satisfies yf F (x) ≥ 1 for all examples in (x n , y n ) N n=1 and for at least one of x w -, x w + , the equality is attained. Define A f F := inf b F ∈R w F H . Suppose the algorithm returns f = (w , b ) and recall that A is the smallest such that f ∈ F A , then for all f = (w, b) ∈ F A such that R(f ) = 0,

A ≥ w H ≥ w F H ≥ A f F .

And it is easy to see that the last equality holds if and only if b F is such that y w -f F (x w -) = y w + f F (x w + ). Then f F ∈ F A . Since R(f F ) = 0 and we have the freedom to choose b F such that the last equality holds, by construction of the learning paradigm and the minimality of A, we have

A = min f ∈F A R(f )=0 A f F = A f F .

In canonical form, for any f F ∈ F A , choose b F such that w F H = A f F , then

2 A f F = 1 A f F φ(F (x w + )) -φ(F (x w -)), w F H ≤ 1 A f F φ(F (x + )) -φ(F (x -)), w F H = φ(F (x + )) -φ(F (x -)) H cos θ F,w .

By construction of the learning paradigm, f must make the equality hold and satisfy θ F,w = 0. Hence we have

A = A f F = 2 φ(F (x + )) -φ(F (x -)) H ,

which proves that, as a function of F , A achieves its minimum if and only if F maximizes φ(F (x + )) -φ(F (x -)) H . Since where we have used the assumption on k, namely, that k(x, x) = φ(x), φ(x) H = φ(x) 2 H = c, for all x. It immediately follows that any minimizer F of A must minimize k(F (x + ), F (x -)) for all pairs of examples from opposite classes. This proves the desired result.

:

A PyTorch-based(Paszke et al.,

2017) library for implementing kMLP and the proposed layer-wise learning algorithm is available at: anonymized URL.2 The optimality of this expansion will be later justified in Section 4 under the layer-wise setting using representer theorem(Schölkopf et al., 2001).

This is in fact a set of random samples asF (l-1) • • • • • F (1) is a set of functions.

Scripts for all experiments performed in this paper are available at: anonymized URL.

Figures and Tables:

Figure 1: (a) Any NN (left, presented in the usual weight-nonlinearity abstraction) can be abstracted as a "graph" (right) with each node representing a neuron and each edge the input-output relationship between neurons. If a node receives multiple inputs, we view its input as a vector in some Euclidean space, as indicated by the colored rectangles. Under this abstraction, each neuron (f (x) = σ(w x + b)) can be directly replaced by a kernel machine (f (x) = w, φ(x) + b with kernel k(x, y) = φ(x), φ(y) ) mapping from the same Euclidean space into the real line without altering the architecture and functionality of the model. (b) Illustration for layer-wise optimality drifting away from network-optimality. Consider a two-layer network and let T 1 , T 2 be the target function of the first and second layer, respectively. If the first layer creates error, which is illustrated by F (1) (x) being far away from T 1 (x), the composed solution F(2) • F (1) on the right is better than that on the left and hence the F (2) on the right corresponds to the network-wise optimality of the second layer. But the F (2) on the left is clearly a better estimate to the layer-wise optimality T 2 if the quality of estimation is measured by the supremum distance.

) (λ, F (i-1) * (x)) -k (i) (λ, T i-1 (x)) ) (λ, F (i-1) * (x)) -k (i) (λ, T i-1 (x))

Proof of Lemma B.1. Since Ω and F 2 are both closed under negation, we haveĜN (F 2 • F 1 ) = E sup (x ν ), F (x n ))g n (x ν ), F (x n ))g n (x ν ), F (x n )) -k(F (x ν ), 0) g n (x n ) 1 g n ≤ ALd 1 ĜN (Ω).

jn k (i) (x n , x) and {α (i) jn ∈ R | n = 1, 2, . . . , N } is the set of trainable parameters for this kernel machine.Solve equation AG = D for A, where A is the d 0 ×N matrix of parameters defined as (A) nm = α(1)

-2 φ(F (x + )), φ(F (x -)) H = arg min F k F (x + ), F (x -) ,

Test errors (%) and 95% confidence intervals (%). When two results have overlapping confidence intervals, they are considered equivalent. Best results are marked in bold. The numbers following the model names indicate the number of hidden layers used. For kMLP FAST , we also include (in parentheses) the portion of centers retained, i.e., the number of training examples randomly chosen as centers for the given layer divided by the size of the training set.(a) Comparing kMLPs (trained layer-wise) with other popular deep architectures trained with BP and BP enhanced by unsupervised greedy pre-training. Further testing the proposed layer-wise learning algorithm and acceleration method using standard MNIST.

Table content: <table><tr><td/><td/><td/><td colspan="2">RECTANGLES</td><td>RECTANGLES-IMAGE</td><td>CONVEX</td><td colspan="2">MNIST (10K)</td><td colspan="3">MNIST (10K) ROTATED</td><td>FASHION-MNIST</td></tr><tr><td colspan="2">MLP-1 (SGD)</td><td/><td colspan="2">7.16± 0.23</td><td>33.20± 0.41</td><td>32.25± 0.41</td><td colspan="2">4.69± 0.19</td><td colspan="2">18.11± 0.34</td><td/><td>15.47± 0.71</td></tr><tr><td colspan="2">MLP-1 (ADAM)</td><td/><td colspan="2">5.37± 0.20</td><td>28.82± 0.40</td><td>30.07± 0.40</td><td colspan="2">4.71± 0.19</td><td colspan="2">18.64± 0.34</td><td/><td>12.98± 0.66</td></tr><tr><td colspan="2">MLP-1 (RMSPROP+BN)</td><td/><td colspan="2">5.37± 0.20</td><td>23.81± 0.37</td><td>28.60± 0.40</td><td colspan="2">4.57± 0.18</td><td colspan="2">18.75± 0.34</td><td/><td>14.55± 0.69</td></tr><tr><td colspan="3">MLP-1 (RMSPROP+DROPOUT)</td><td colspan="2">5.50± 0.20</td><td>23.67± 0.37</td><td>36.28± 0.42</td><td colspan="2">4.31± 0.18</td><td colspan="2">14.96± 0.31</td><td/><td>12.86± 0.66</td></tr><tr><td colspan="2">MLP-2 (SGD)</td><td/><td colspan="2">5.05± 0.19</td><td>22.77± 0.37</td><td>25.93± 0.38</td><td colspan="2">5.17± 0.19</td><td colspan="2">18.08± 0.34</td><td/><td>12.94± 0.66</td></tr><tr><td colspan="2">MLP-2 (ADAM)</td><td/><td colspan="2">4.36± 0.18</td><td>25.69± 0.38</td><td>25.68± 0.38</td><td colspan="2">4.42± 0.18</td><td colspan="2">17.22± 0.33</td><td/><td>11.48± 0.62</td></tr><tr><td colspan="2">MLP-2 (RMSPROP+BN)</td><td/><td colspan="2">4.22± 0.18</td><td>23.12± 0.37</td><td>23.28± 0.37</td><td colspan="2">3.57± 0.16</td><td colspan="2">13.73± 0.30</td><td/><td>11.51± 0.63</td></tr><tr><td colspan="3">MLP-2 (RMSPROP+DROPOUT)</td><td colspan="2">4.75± 0.19</td><td>23.24± 0.37</td><td>34.73± 0.42</td><td colspan="2">3.95± 0.17</td><td colspan="2">13.57± 0.30</td><td/><td>11.05± 0.61</td></tr><tr><td>DBN-1</td><td/><td/><td colspan="2">4.71± 0.19</td><td>23.69± 0.37</td><td>19.92± 0.35</td><td colspan="2">3.94± 0.17</td><td colspan="2">14.69± 0.31</td><td/><td>N/A</td></tr><tr><td>DBN-3</td><td/><td/><td colspan="2">2.60± 0.14</td><td>22.50± 0.37</td><td>18.63± 0.34</td><td colspan="2">3.11± 0.15</td><td colspan="2">10.30± 0.27</td><td/><td>N/A</td></tr><tr><td>SAE-3</td><td/><td/><td colspan="2">2.41± 0.13</td><td>24.05± 0.37</td><td>18.41± 0.34</td><td colspan="2">3.46± 0.16</td><td colspan="2">10.30± 0.27</td><td/><td>N/A</td></tr><tr><td>KMLP-1</td><td/><td/><td colspan="2">2.24± 0.13</td><td>23.29± 0.37</td><td>19.15± 0.34</td><td colspan="2">3.10± 0.15</td><td colspan="2">11.09± 0.28</td><td/><td>11.72± 0.63</td></tr><tr><td colspan="2">KMLP-1 FAST</td><td/><td colspan="2">2.36± 0.13 (0.05)</td><td>23.86± 0.37 (0.01)</td><td>20.34± 0.35 (0.17)</td><td colspan="2">2.95± 0.15 (0.1)</td><td colspan="2">12.61± 0.29 (0.1)</td><td/><td>11.45± 0.62 (0.28)</td></tr><tr><td>KMLP-2</td><td/><td/><td colspan="2">2.24± 0.13</td><td>23.30± 0.37</td><td>18.53± 0.34</td><td colspan="2">3.16± 0.15</td><td colspan="2">10.53± 0.27</td><td/><td>11.25± 0.62</td></tr><tr><td colspan="2">KMLP-2 FAST</td><td/><td colspan="10">2.21± 0.13 (0.3/0.3) 23.24± 0.37 (0.01/0.3) 19.32± 0.35 (0.005/0.03) 3.18± 0.15 (0.3/0.3) 10.94± 0.27 (0.1/0.7) 11.15± 0.62 (1/0.28)</td></tr><tr><td colspan="2">(b) MLP-1</td><td/><td colspan="2">MLP-1</td><td>MLP-2</td><td>MLP-2</td><td colspan="2">KMLP-1</td><td colspan="2">KMLP-1</td><td/><td>KMLP-1</td><td>RFF</td></tr><tr><td colspan="2">(RMSPROP+BN)</td><td colspan="3">(RMSPROP+DROPOUT)</td><td>(RMSPROP+BN)</td><td>(RMSPROP+DROPOUT)</td><td>(BP)</td><td/><td colspan="2">(GREEDY)</td><td/><td>(BP)</td></tr><tr><td colspan="2">2.05± 0.28</td><td/><td colspan="2">1.77± 0.26</td><td>1.58± 0.24</td><td>1.67± 0.25</td><td colspan="2">3.44± 0.36</td><td colspan="2">1.77± 0.26</td><td/><td>2.01± 0.28</td></tr><tr><td>KMLP-1</td><td>PARAM</td><td/><td>KMLP-1</td><td>FAST</td><td>KMLP-2</td><td>KMLP-2</td><td>KMLP-2</td><td>RFF</td><td>KMLP-2</td><td>PARAM</td><td/><td>KMLP-2</td><td>FAST</td></tr><tr><td colspan="2">(BP)</td><td/><td colspan="2">(GREEDY)</td><td>(BP)</td><td>(GREEDY)</td><td>(BP)</td><td/><td colspan="2">(BP)</td><td/><td>(GREEDY)</td></tr><tr><td colspan="2">1.88± 0.27</td><td colspan="3">1.75± 0.26 (0.54)</td><td>3.66± 0.37</td><td>1.56± 0.24</td><td colspan="2">1.92± 0.27</td><td colspan="2">2.45± 0.30</td><td colspan="2">1.47± 0.24 (1/0.19)</td></tr></table>

L p was adopted; Generalized MKL in(Varma & Babu, 2009) (GMKL), for which the target kernel class was the Hadamard product of single Gaussian kernel defined on each dimension; Infinite Kernel Learning in(Gehler & Nowozin, 2008) (IKL) with MKL LEVEL as the embedded optimizer for kernel weights; 2-layer Multilayer Kernel Machine in(Cho & Saul, 2009) (MKM); 2-Layer MKL (2LMKL) and Infinite 2-Layer MKL in(Zhuang et al., 2011) (2LMKL INF ).Eleven binary classification data sets that have been widely used in MKL literature were split evenly for training and test and were all normalized to zero mean and unit variance prior to training. 20 runs with identical settings but random weight initializations were repeated for each model. For each repetition, a new training-test split was selected randomly.

Table content: <table><tr><td/><td>Size/Dimension</td><td>SVM</td><td>MKL LEVEL</td><td>L p MKL</td><td>GMKL</td><td>IKL</td><td>MKM</td><td>2LMKL</td><td>2LMKL INF</td><td>kMLP-1</td></tr><tr><td>Breast</td><td>683/10</td><td>3.2± 1.0</td><td>3.5± 0.8</td><td>3.8± 0.7</td><td>3.0± 1.0</td><td>3.5± 0.7</td><td>2.9± 1.0</td><td>3.0± 1.0</td><td>3.1± 0.7</td><td>2.4± 0.7</td></tr><tr><td>Diabetes</td><td>768/8</td><td colspan="9">23.3± 1.8 24.2± 2.5 27.4± 2.5 33.6± 2.5 24.0± 3.0 24.2± 2.5 23.4± 1.6 23.4± 1.9 23.2± 1.9</td></tr><tr><td>Australian</td><td>690/14</td><td colspan="9">15.4± 1.4 15.0± 1.5 15.5± 1.6 20.0± 2.3 14.6± 1.2 14.7± 0.9 14.5± 1.6 14.3± 1.6 13.8± 1.7</td></tr><tr><td>Iono</td><td>351/33</td><td>7.2± 2.0</td><td>8.3± 1.9</td><td>7.4± 1.4</td><td>7.3± 1.8</td><td>6.3± 1.0</td><td>8.3± 2.7</td><td>7.7± 1.5</td><td>5.6± 0.9</td><td>5.0± 1.4</td></tr><tr><td>Ringnorm</td><td>400/20</td><td>1.5± 0.7</td><td>1.9± 0.8</td><td>3.3± 1.0</td><td>2.5± 1.0</td><td>1.5± 0.7</td><td>2.3± 1.0</td><td>2.1± 0.8</td><td>1.5± 0.8</td><td>1.5± 0.6</td></tr><tr><td>Heart</td><td>270/13</td><td colspan="9">17.9± 3.0 17.0± 2.9 23.3± 3.8 23.0± 3.6 16.7± 2.1 17.6± 2.5 16.9± 2.5 16.4± 2.1 15.5± 2.7</td></tr><tr><td>Thyroid</td><td>140/5</td><td>6.1± 2.9</td><td>7.1± 2.9</td><td>6.9± 2.2</td><td>5.4± 2.1</td><td>5.2± 2.0</td><td>7.4± 3.0</td><td>6.6± 3.1</td><td>5.2± 2.2</td><td>3.8± 2.1</td></tr><tr><td>Liver</td><td>345/6</td><td colspan="9">29.5± 4.1 37.7± 4.5 30.6± 2.9 36.4± 2.6 40.0± 2.9 29.9± 3.6 34.0± 3.4 37.3± 3.1 28.9± 2.9</td></tr><tr><td>German</td><td>1000/24</td><td colspan="9">24.8± 1.9 28.6± 2.8 25.7± 1.4 29.6± 1.6 30.0± 1.5 24.3± 2.3 25.2± 1.8 25.8± 2.0 24.0± 1.8</td></tr><tr><td>Waveform</td><td>400/21</td><td colspan="7">11.0± 1.8 11.8± 1.6 11.1± 2.0 11.8± 1.8 10.3± 2.3 10.0± 1.6 11.3± 1.9</td><td>9.6± 1.6</td><td>10.3± 1.9</td></tr><tr><td>Banana</td><td>400/2</td><td>10.3± 1.5</td><td>9.8± 2.0</td><td colspan="5">12.5± 2.6 16.6± 2.7 9.8± 1.8 19.5± 5.3 13.2± 2.1</td><td>9.8± 1.6</td><td>11.5± 1.9</td></tr><tr><td>Rank</td><td>-</td><td>4.2</td><td>6.3</td><td>7.0</td><td>6.9</td><td>4.3</td><td>5.4</td><td>5.0</td><td>2.8</td><td>1.6</td></tr></table>

Peer review (iclr2019):
 ****Reply to authors* rebuttal**** Dear Authors, Thank you very much for all the effort you have put into the rebuttal. Based on the improved theoretical and experimental results, I have decided to increase my score from 5 to 6. Best wishes, Rev 1 ****Original review**** This paper explores integration of kernel machines with neural networks based on replacing the non-linear function represented by each neuron with a function living in some pre-defined RKHS. From the theoretical standpoint, this work is a clear improvement upon the work of Zhang et al. (2017). Authors further propose a layer-wise training algorithm based on optimisation of a particular similarity measure between embeddings based on their class assignments at each layer, which eliminates necessity of gradient-based training. However, the experimental performance of the proposed algorithm is somewhat lacking in comparison, perhaps because the authors focus on kernelised equivalents of MLPs instead of CNNs as Zhang et al. My rating of the paper is mainly due to the lack of experimental evidence for usefulness of the layer-wise training, and absence of experimental comparison with several baselines (see details below). It is also unclear whether the structure of KNs is significantly better than that of NNs in terms of interpretability. Apart from the comments below, I would like to ask the authors to discuss relation to the following related papers: 1) Kulkarni & Karande, 2017: "Layer-wise training of deep networks using kernel similarity" https://arxiv.org/pdf/1703.07115.pdf 2) Scardapanea et al., 2017: "Kafnets: kernel-based non-parametric activation functions for neural networks" https://arxiv.org/pdf/1707.04035.pdf Detailed comments: Theory - (Sec 4.1) Backpropagation (BP) is being criticised: BP is only a particular implementation of gradient calculation. It seems to me that your criticisms are thus more related to use of iterative gradient-based optimisation algorithms, rather than to obtaining gradients through BP?! Regarding the criticism that BP forces intermediate layers to correct for "mistakes" made by layers higher up: it seems your layer-wise algorithm attempts to learn the best possible representation in first layer, and then progresses to the next layer where it tries to correct for the potential error of the first layer and so on. In other words it seems that the errors of layers are propagated from first to last, instead of last to first as in BP, but are still being propagated in a sense. I do not immediately see why propagation forward should be preferable. Can you please further explain this point? - It is proven in the appendix (Lemma B.3) that under certain conditions stacking additional layers never leads to degradation of training loss. Can you please clarify whether additional layers can be helpful even in the case where previous layers already succeeded in learning the optimal representation? - (Sec 4.1) Layer-wise vs. network-wise optimality: I find the claim that BP-based learner is not aware of the network-wise optimality confusing. BP explicitly optimises for network-wise optimality and the relative contribution to the network-wise error of each weight is propagated accordingly. I suppose my confusion stems from lack of a clear description of what defines a learner "aware" or "blind" to network-wise optimality. In general, I am not convinced layer-wise optimality is a useful criterion when what we want to achieve is network-wise optimality. As you show in the appendix, if layer-wise optimality is achieved then it implies network-wise optimality; however, layer-wise optimality is only a sufficient condition and likely not a necessary one (except for the simplified scenario studied in B.3). It is thus not clear to me why layer-wise training would always be preferable to network-wise training (e.g. using BP) especially because its greedy nature might intuitively prevent learning of hierarchical representations which are commonly claimed to be key to the success of neural networks. Can you please clarify? - (Sec 4.2) I think it would be beneficial to state in the introduction that the "risk" is with respect to the hinge loss which is common in the SVM/kernel literature but much less in the deep learning literature and thus could surprise a few people when they reach this point. Futher questions: - From Lemma 4.3, it seems that the derived representation is only optimal with respect to the **upper bound** on the empirical risk (which for .tau >= 2 will be an upper bound on the population risk). I got slightly confused at this point as my interpretation of the previous text was that the representation is optimal with respect to the population risk itself. Does the upper bound have the same set of optima? Please clarify. - (p.5) There are two assumptions that I find somewhat restrictive. Just before Lemma 4.3 you assume that the number of points in each class must be the same. Can you comment on whether you expect the same representation to be optimal for classification problems with significantly imbalanced number of samples per class? The second assumption is after Lemma 4.4 where you state that the stationary kernel k^{l-1} should attain its infinum for all x, y s.t. || x - y || greater than some threshold. This does not hold for many of the popular kernels like RBF, Matern, or inverse multiquadric. Do you think this assumption can be relaxed? - (p.5) Choice of the dissimilarity measure for G: Can you provide more intuition about why you selected L^1 distance and whether you would expect different results with L^2 or other common metrics? - (Sec 4.3) Can you please provide more detaild about the relation of the proposed objective (.hat(R)(F) + .tau max_j ||f_j||_H) to Lemmas 4.3 and 4.5 where the optimal representation was derived for functions that optimise an upper bound in terms of Gaussian complexity (e.g. is the representation that minimises risk w.r.t. the Gaussian bound also optimal with respect to functions that optimise this objective)? Experiments - I would appreciate addition of some standard baselines, like MLP combined with dropout or batch normalisation, and optimised with RMSProp (or similar). These would greatly help with assessing competitiveness with current SOTA results. - It would be nice to see the relative contribution of the two main components of the paper. Specifically, an experiment which would evaluate empirical performance of KNs optimised by some form of gradient descent vs. by your layer-wise training rule would be very insightful. Other - (p.2, 1st par in Sec 2) [minor] You state "a kernel machine is a universal function approximator". I suppose that might be true for a certain class of kernels but not in general?! Please clarify. - (p.2, 3rd par in Sec 2) [minor] Are you using a particular version of the representer theorem in the representation of f_j^{(i)} as linear combination of feature maps? Please clarify. - (p.2, end of 1st par in Sec 3) L^{(i)} is defined as sup over X_i. It is not clear to me that this constant is necessarily finite and I suspect it will not be in general (it will for the RBF kernel (and most stationary kernels) used in experiments though). Finiteness of L^{(i)} is necessary for the bound in Eq. (2) to be non-vacuous. Please clarify. - (p.3, after 1st display in Sec 4.2.1) [minor] Missing dot after "that we wish to minimise". Next sentence states "**the** optimal F" (emphasis mine) -- I am sorry if I overlooked it, but I did not notice a proof that a solution exists and is unique, and am not familiar enough with the literature to immediately see the answer. Perhaps a footnote clarifying the statement would help. - (p.4, 1st par in Sec 4) You say "A generalisation to regression is reserved for future work". I did not expect that based on the first few pages. On high-level, it seems that generalisation to regression need not be trivial as, for example, the optimal representation derived in Lemma 4.3 and Lemma 4.5 explicitly relies on the classification nature of the problem. Can you comment on expected difficulty of extension to regression? Possibly state in the introduction that only classification is considered in this paper. - (p.7, 1st par in Sec 6) [related] "However they did not extend the idea to any **arbitrary** NN" (emphasis mine). Can you please be more specific here? - (p.5-6) [minor] Last sentence in Lemmas 4.3 and 4.5 is slightly confusing. Can you rephrase please? - (p.6) [minor] You say "the learned decision boundary would generalise better to unseen data". Can you please clarify the last sentence (e.g. being more precise about the meaning of the word "simple" in the same sentence) and provide reference for why this is necessarily the case?