summaries: In this paper, the authors study some simple convolutional kernels. They show that two or three layers convolutional kernels with polynomial kernels on the higher level and Gaussian pooling provide similar performance as the much more complicated state-of-the-art convolutional kernels (e.g., Myrtle kernel). Motivated by these good performance, they proceed to characterize the RKHS of these kernels, and describe how extra layers and pooling allows to capture interaction between patches with more or less spatial dependency. They then use the RKHS norm and some standard bound on the generalization error to show how choosing an architecture adapted to the target function can improve the statistical efficiency. The paper studies the RKHS and generalization properties of convolutional kernel networks, the NNGP corresponding to CNNs and the class of kernels that achieve state of the art performance on image classification. Among the theoretical contributions are analysis of the regularization induced by pooling, and interactions between patches captured with iterated convolutions. The paper also includes experimental results on CIFAR10 matching the prior state of the art for a kernel method while using a shallower architecture, as well as some ablations on the size of the convolutional kernels, the size of the pooling filters, and the number of layers.  This paper analyzed convolutional kernels. Experimentally, the authors showed that convolutional kernel formed by shallow convolutional neural networks (2 - 3 layers) performs as well as the state of art deep convolutional kernels (e.g., Shankar et. al.). The authors also provide an exact description of the RKHS functions and their norm, of the one layer convolutional kernel and two layer convolutional kernel with low degree polynomial activation function on the top layer. Finally, assuming that the target function has a specific form, then theoretical results can give the generalization upper bound of kernel ridge regression in terms of the sample complexity; this shows that using a proper architecture can same the sample complexity by a polynomial factor of the input size.   The paper aims to study properties of deep convolutional models via the surrogate of simple hierarchical kernels with convolution and pooling layers. The authors characterise the underlying RKHS and their norms, and provide generalization bounds. The results imply that convolution operations such as pooling and patches, if present in the data, lead to improve guarantees. An empirical study both justifies studying these simple kernels in the first place, and illustrates the obtained theoretical results. 
abstractText: The empirical success of deep convolutional networks on tasks involving highdimensional data such as images or audio suggests that they can efficiently approximate certain functions that are well-suited for such tasks. In this paper, we study this through the lens of kernel methods, by considering simple hierarchical kernels with two or three convolution and pooling layers, inspired by convolutional kernel networks. These achieve good empirical performance on standard vision datasets, while providing a precise description of their functional space that yields new insights on their inductive bias. We show that the RKHS consists of additive models of interaction terms between patches, and that its norm encourages spatial similarities between these terms through pooling layers. We then provide generalization bounds which illustrate how pooling and patches yield improved sample complexity guarantees when the target function presents such regularities.
references.authors: Zeyuan Allen-Zhu Yuanzhi Li
references.title: Backward feature correction: How deep learning performs deep learning
references.venue: arXiv preprint arXiv:2001.04413,
references.authors: Sanjeev Arora Simon S Du Wei Hu Zhiyuan Li Russ R Salakhutdinov Ruosong Wang
references.title: On exact computation with an infinitely wide neural net
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Francis Bach
references.title: Breaking the curse of dimensionality with convex neural networks
references.venue: Journal of Machine Learning Research (JMLR),
references.authors: Francis Bach
references.title: On the equivalence between kernel quadrature rules and random feature expansions
references.venue: Journal of Machine Learning Research (JMLR),
references.authors: Francis Bach
references.title: Learning Theory from First Principles (draft)
references.venue: URL https://www.di.ens. fr/~fbach/ltfp_book.pdf
references.authors: Gregory Beylkin Martin J Mohlenkamp
references.title: Numerical operator calculus in higher dimensions
references.venue: Proceedings of the National Academy of Sciences,
references.authors: Alberto Bietti Francis Bach
references.title: Deep equals shallow for ReLU networks in kernel regimes
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Alberto Bietti Julien Mairal
references.title: Group invariance, stability to deformations, and complexity of deep convolutional representations
references.venue: Journal of Machine Learning Research (JMLR),
references.authors: Alberto Bietti Julien Mairal
references.title: On the inductive bias of neural tangent kernels
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Joan Bruna Stéphane Mallat
references.title: Invariant scattering convolution networks
references.venue: IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
references.authors: Andrea Caponnetto Ernesto De Vito
references.title: Optimal rates for the regularized least-squares algorithm
references.venue: Foundations of Computational Mathematics,
references.authors: Lin Chen Sheng Xu
references.title: Deep neural tangent kernel and laplace kernel have the same rkhs
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Minshuo Chen Yu Bai Jason D Lee Tuo Zhao Huan Wang Caiming Xiong Richard Socher
references.title: Towards understanding hierarchical learning: Benefits of neural representations
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Lenaic Chizat Francis Bach
references.title: Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss
references.venue: In Conference on Learning Theory,
references.authors: Lenaic Chizat Edouard Oyallon Francis Bach
references.title: On lazy training in differentiable programming
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Youngmin Cho Lawrence K Saul
references.title: Kernel methods for deep learning
references.venue: In Advances in Neural Information Processing Systems (NIPS),
references.authors: Carlo Ciliberto Francis Bach Alessandro Rudi
references.title: Localized structured prediction
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Nadav Cohen Amnon Shashua
references.title: Convolutional rectifier networks as generalized tensor decompositions
references.venue: In Proceedings of the International Conference on Machine Learning (ICML),
references.authors: Nadav Cohen Amnon Shashua
references.title: Inductive bias of deep convolutional networks through pooling geometry
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Amit Daniely Roy Frostig Yoram Singer
references.title: Toward deeper understanding of neural networks: The power of initialization and a dual view on expressivity
references.venue: In Advances in Neural Information Processing Systems (NIPS),
references.authors: Simon S Du Yining Wang Xiyu Zhai Sivaraman Balakrishnan Ruslan Salakhutdinov Aarti Singh
references.title: How many samples are needed to estimate a convolutional neural network
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Costas Efthimiou Christopher Frye
references.title: Spherical harmonics in p dimensions
references.venue: World Scientific,
references.authors: Alessandro Favero Francesco Cagnetta Matthieu Wyart
references.title: Locality defeats the curse of dimensionality in convolutional teacher-student scenarios
references.venue: arXiv preprint arXiv:2106.08619,
references.authors: Adrià Garriga-Alonso Laurence Aitchison Carl Edward Rasmussen
references.title: Deep convolutional networks as shallow gaussian processes
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Amnon Geifman Abhay Yadav Yoni Kasten Meirav Galun David Jacobs Ronen Basri
references.title: On the similarity between the laplace and neural tangent kernels
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Behrooz Ghorbani Song Mei Theodor Misiakiewicz Andrea Montanari
references.title: When do neural networks outperform kernel methods
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Suriya Gunasekar Jason D Lee Daniel Soudry Nati Srebro
references.title: Implicit bias of gradient descent on linear convolutional networks
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Wolfgang Hackbusch Stefan Kühn
references.title: A new scheme for the tensor representation
references.venue: Journal of Fourier analysis and applications,
references.authors: Kaiming He Xiangyu Zhang Shaoqing Ren Jian Sun
references.title: Deep residual learning for image recognition
references.venue: In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition
references.authors: Reinhard Heckel Mahdi Soltanolkotabi
references.title: Denoising and regularization via exploiting the structural bias of convolutional generators
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Arthur Jacot Franck Gabriel Clément Hongler
references.title: Neural tangent kernel: Convergence and generalization in neural networks
references.venue: In Advances in Neural Information Processing Systems (NIPS),
references.authors: Hervé Jégou Florent Perronnin Matthijs Douze Jorge Sánchez Patrick Pérez Cordelia Schmid
references.title: Aggregating local image descriptors into compact codes
references.venue: IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),
references.authors: Jaehoon Lee Samuel Schoenholz Jeffrey Pennington Ben Adlam Lechao Xiao Roman Novak Jascha Sohl-Dickstein
references.title: Finite versus infinite neural networks: an empirical study
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Zhiyuan Li Ruosong Wang Dingli Yu Simon S Du Wei Hu Ruslan Salakhutdinov Sanjeev Arora
references.title: Enhanced convolutional neural tangent kernels
references.authors: Zhiyuan Li Yi Zhang Sanjeev Arora
references.title: Why are convolutional nets more sample-efficient than fully-connected nets
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Yi Lin
references.title: Tensor product space anova models
references.venue: Annals of Statistics,
references.authors: David G Lowe
references.title: Object recognition from local scale-invariant features
references.venue: In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR),
references.authors: Julien Mairal
references.title: End-to-End Kernel Learning with Supervised Convolutional Kernel Networks
references.venue: In Advances in Neural Information Processing Systems (NIPS),
references.authors: Julien Mairal Piotr Koniusz Zaid Harchaoui Cordelia Schmid
references.title: Convolutional kernel networks
references.venue: In Advances in Neural Information Processing Systems (NIPS),
references.authors: Eran Malach Shai Shalev-Shwartz
references.title: Computational separation between convolutional and fullyconnected networks
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Stéphane Mallat
references.title: Group invariant scattering
references.venue: Communications on Pure and Applied Mathematics,
references.authors: Song Mei Theodor Misiakiewicz Andrea Montanari
references.title: Learning with invariances in random features and kernel models
references.venue: In Conference on Learning Theory (COLT),
references.authors: Hrushikesh N Mhaskar Tomaso Poggio
references.title: Deep vs. shallow networks: An approximation theory perspective
references.venue: Analysis and Applications,
references.authors: Ha Quang Minh Partha Niyogi Yuan Yao
references.title: Mercer’s theorem, feature maps, and smoothing
references.venue: In Conference on Learning Theory (COLT),
references.authors: Roman Novak Lechao Xiao Yasaman Bahri Jaehoon Lee Greg Yang Jiri Hron Daniel A Abolafia Jeffrey Pennington Jascha Sohl-Dickstein
references.title: Bayesian deep convolutional networks with many channels are gaussian processes
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Tomaso Poggio Hrushikesh Mhaskar Lorenzo Rosasco Brando Miranda Qianli Liao
references.title: Why and when can deep-but not shallow-networks avoid the curse of dimensionality: a review
references.venue: International Journal of Automation and Computing,
references.authors: Saburou Saitoh
references.title: Integral transforms, reproducing kernels and their applications, volume 369
references.venue: CRC Press,
references.authors: Jorge Sánchez Florent Perronnin Thomas Mensink Jakob Verbeek
references.title: Image classification with the fisher vector: Theory and practice
references.venue: International Journal of Computer Vision (IJCV),
references.authors: Meyer Scetbon Zaid Harchaoui
references.title: Harmonic decompositions of convolutional networks
references.venue: In Proceedings of the International Conference on Machine Learning (ICML),
references.authors: Johannes Schmidt-Hieber
references.title: Nonparametric regression using deep neural networks with relu activation function
references.venue: Annals of Statistics,
references.authors: Bernhard Schölkopf Alexander J Smola
references.title: Learning with kernels: support vector machines, regularization, optimization, and beyond
references.authors: Vaishaal Shankar Alex Fang Wenshuo Guo Sara Fridovich-Keil Jonathan Ragan-Kelley Ludwig Schmidt Benjamin Recht
references.title: Neural kernels without tangents
references.venue: In Proceedings of the International Conference on Machine Learning (ICML),
references.authors: Winfried Sickel Tino Ullrich
references.title: Tensor products of sobolev–besov spaces and applications to approximation from the hyperbolic cross
references.venue: Journal of Approximation Theory,
references.authors: Alex J Smola Zoltan L Ovari Robert C Williamson
references.title: Regularization with dot-product kernels
references.venue: In Advances in Neural Information Processing Systems (NIPS),
references.authors: Louis Thiry Michael Arbel Eugene Belilovsky Edouard Oyallon
references.title: The unreasonable effectiveness of patches in deep convolutional kernels methods
references.venue: In Proceedings of the International Conference on Learning Representations (ICLR),
references.authors: Ulrike von Luxburg Olivier Bousquet
references.title: Distance-based classification with lipschitz functions
references.venue: Journal of Machine Learning Research (JMLR),
references.authors: Grace Wahba
references.title: Spline models for observational data, volume 59
references.authors: Martin J Wainwright
references.title: High-dimensional statistics: A non-asymptotic viewpoint, volume 48
references.authors: Colin Wei Jason Lee Qiang Liu Tengyu Ma
references.title: Regularization matters: Generalization and optimization of neural nets vs their induced kernel
references.venue: In Advances in Neural Information Processing Systems (NeurIPS),
references.authors: Thomas Wiatowski Helmut Bölcskei
references.title: A mathematical theory of deep convolutional neural networks for feature extraction
references.venue: IEEE Transactions on Information Theory,
references.authors: Greg Yang
references.title: Scaling limits of wide neural networks with weight sharing: Gaussian process behavior, gradient independence, and neural tangent kernel derivation
references.authors: Matthew D Zeiler Rob Fergus
references.title: Visualizing and understanding convolutional networks
references.venue: In Proceedings of the European Conference on Computer Vision (ECCV),
references.authors: Y. Zhang P. Liang M.J. Wainwright
references.title: Convexified convolutional neural networks
references.venue: In International Conference on Machine Learning (ICML),
references.authors: e.g Thiry
references.title: One simplified example that encodes such a prior is if the target function f∗ only depends on the input image through a localized part of the input such as a patch xu = (x[u + v])v∈S ∈ Rp|S|, where S is a small box centered around
references.authors: He
references.title: 2016), as the feature selection capabilities of SGD may alleviate the dependence on dimension, e.g., by finding Gabor-like filters. Gaussian vs average pooling. Table 7 shows the differences in performance between two or three layer architectures considered in Table 2, when Gaussian pooling filters
references.authors: Shankar
references.title: patches from each image, and reconstructing the image from whitened patches via averaging. In contrast, other approaches use global whitening of the entire image Lee et al
references.venue: Finite networks and comparison to Shankar et al
references.authors: Chen
references.title: 87.9%). This may suggest that for convolutional networks, deeper networks may have additional advantages when using optimization algorithms, in terms of adapting to possibly relevant structure of the problem, such as hierarchical representations (see, e.g., Allen-Zhu
references.authors: Minh
references.title: For instance, for the exponential kernel (or Gaussian on the sphere) κ(〈x, y
sections.heading: 1 INTRODUCTION
sections.text: Deep convolutional models have been at the heart of the recent successes of deep learning in problems where the data consists of high-dimensional signals, such as image classification or speech recognition. Convolution and pooling operations have notably contributed to the practical success of these models, yet our theoretical understanding of how they enable efficient learning is still limited.
  One key difficulty for understanding such models is the curse of dimensionality: due to the highdimensionality of the input data, it is hopeless to learn arbitrary functions from samples. For instance, classical non-parametric regression techniques for learning generic target functions typically require either low dimension or very high degrees of smoothness in order to obtain good generalization (e.g., Wainwright, 2019), which makes them impractical for dealing with high-dimensional signals. Thus, further assumptions on the target function are needed to make the problem tractable, and we seek assumptions that make convolutions a useful modeling tool. Various works have studied approximation benefits of depth with models that resemble deep convolutional architectures (Cohen & Shashua, 2017; Mhaskar & Poggio, 2016; Schmidt-Hieber et al., 2020). Nevertheless, while such function classes may provide improved statistical efficiency in theory, it is unclear if there exist efficient algorithms to learn such models, and hence, whether they might correspond to what convolutional networks learn in practice. To overcome this issue, we consider instead function classes based on kernel methods (Schölkopf & Smola, 2001; Wahba, 1990), which are known to be learnable with efficient (polynomial-time) algorithms, such as kernel ridge regression or gradient descent.
  We consider “deep” structured kernels known as convolutional kernels, which yield good empirical performance on standard computer vision benchmarks (Li et al., 2019; Mairal, 2016; Mairal et al., 2014; Shankar et al., 2020), and are related to over-parameterized convolutional networks (CNNs) in so-called “kernel regimes” (Arora et al., 2019; Bietti & Mairal, 2019b; Daniely et al., 2016; GarrigaAlonso et al., 2019; Jacot et al., 2018; Novak et al., 2019; Yang, 2019). Such regimes may be seen as providing a first-order description of what common deep models trained with gradient methods may learn. Studying the corresponding function spaces (reproducing kernel Hilbert spaces, or RKHS) may then provide insight into the benefits of various architectural choices. For fully-connected architectures, such kernels are rotation-invariant, and the corresponding RKHSs are well understood in terms of regularity properties on the sphere (Bach, 2017a; Smola et al., 2001), but do not show any major differences between deep and shallow kernels (Bietti & Bach, 2021; Chen & Xu, 2021; Geifman et al., 2020). In contrast, in this work we show that even in the kernel setting, multiple layers
  of convolution and pooling operations can be crucial for efficient learning of functions with specific structures that are well-suited for natural signals. Our work paves the way for further studies of the inductive bias of optimization algorithms on deep convolutional networks beyond kernel regimes, for instance by incorporating adaptivity to low-dimensional structure (Bach, 2017a; Chizat & Bach, 2020; Wei et al., 2019) or hierarchical learning (Allen-Zhu & Li, 2020).
  We make the following contributions: • We revisit convolutional kernel networks (Mairal, 2016), finding that simple two or three layers
  models with Gaussian pooling and polynomial kernels of degree 2-4 at higher layers provide competitive performance with state-of-the-art convolutional kernels such as Myrtle kernels (Shankar et al., 2020) on Cifar10.
  • For such kernels, we provide an exact description of the RKHS functions and their norm, illustrating representation benefits of multiple convolutional and pooling layers for capturing additive and interaction models on patches with certain spatial regularities among interaction terms.
  • We provide generalization bounds that illustrate the benefits of architectural choices such as pooling and patches for learning additive interaction models with spatial invariance in the interaction terms, namely, improvements in sample complexity by polynomial factors in the size of the input signal.
  Related work. Convolutional kernel networks were introduced by Mairal et al. (2014); Mairal (2016). Empirically, they used kernel approximations to improve computational efficiency, while we evaluate the exact kernels in order to assess their best performance, as in (Arora et al., 2019; Li et al., 2019; Shankar et al., 2020). Bietti & Mairal (2019a;b) show invariance and stability properties of its RKHS functions, and provide upper bounds on the RKHS norm for some specific functions (see also Zhang et al., 2017); in contrast, we provide exact characterizations of the RKHS norm, and study generalization benefits of certain architectures. Scetbon & Harchaoui (2020) study statistical properties of simple convolutional kernels without pooling, while we focus on the role of architecture choices with an emphasis on pooling. Cohen & Shashua (2016; 2017); Mhaskar & Poggio (2016); Poggio et al. (2017) study expressivity and approximation with models that resemble CNNs, showing benefits thanks to hierarchy or local interactions, but such models are not known to be learnable with tractable algorithms, while we focus on (tractable) kernels. Regularization properties of convolutional models were also considered in (Gunasekar et al., 2018; Heckel & Soltanolkotabi, 2020), but in different regimes or architectures than ours. Li et al. (2021); Malach & Shalev-Shwartz (2021) study benefits of convolutional networks with efficient algorithms, but do not study the gains of pooling. Du et al. (2018) study sample complexity of learning CNNs, focusing on parametric rather than non-parametric models. Mei et al. (2021) study statistical benefits of global pooling for learning invariant functions, but only consider one layer with full-size patches. Concurrently to our work, Favero et al. (2021) study benefits of local patches, but focus on one-layer models.
sections.heading: 2 DEEP CONVOLUTIONAL KERNELS
sections.text: In this section, we recall the construction of multi-layer convolutional kernels on discrete signals, following most closely the convolutional kernel network (CKN) architectures studied by Mairal (2016); Bietti & Mairal (2019a). These architectures rely crucially on pooling layers, typically with Gaussian filters, which make them empirically effective even with just two convolutional layers. These kernels define function spaces that will be the main focus of our theoretical study of approximation and generalization in the next sections. In particular, when learning a target function of the form f∗(x) = ∑ i fi(x), we will show that they are able to efficiently exploit two useful properties of f
  ∗: (locality) each fi may depend on only one or a few small localized patches of the signal; (invariance) many different terms fi may involve the same function applied to different input patches. We provide further background and motivation in Appendix A.
  For simplicity, we will focus on discrete 1D input signals, though one may easily extend our results to 2D or higher-dimensional signals. We will assume periodic signals in order to avoid difficulties with border effects, or alternatively, a cyclic domain Ω = Z/|Ω|Z. A convolutional kernel of depth L may then be defined for input signals x, x′ ∈ L2(Ω,Rp) by KL(x, x′) = 〈Ψ(x),Ψ(x′)〉, through the explicit feature map Ψ(x) = ALMLPL · · ·A1M1P1x. (1) Here, P`,M` and A` are linear or non-linear operators corresponding to patch extraction, kernel mapping and pooling, respectively, and are described below. They operate on feature maps in L2(Ω`)
  (with Ω0 =Ω) with values in different Hilbert spaces, starting fromH0 := Rp, and are defined below. An illustration of this construction is given in Figure 1.
  Patch extraction. Given a patch shape S` ⊂ Ω`−1, such as S` = [−1, 0, 1] for one-dimensional patches of size 3, the operator P` is defined for x ∈ L2(Ω`−1,H`−1) by1
  P`x[u] = (x[u+ v])v∈S` ∈ H |S`| `−1.
  Kernel mapping. The operators M` perform a nonlinear embedding of patches into a new Hilbert space using dot-product kernels. We consider homogeneous dot-product kernels given for z, z′ ∈ H|S`|`−1 by
  k`(z, z ′) = ‖z‖‖z′‖κ` ( 〈z, z′〉 ‖z‖‖z′‖ ) = 〈ϕ`(z), ϕ`(z′)〉H` , (2)
  where ϕ` : H|S`|`−1 → H` is a feature map for the kernel. The kernel functions take the form κ`(u) =∑
  j≥0 bju j with bj ≥ 0. This includes the exponential kernel κ(u) = eα(u−1) (i.e., Gaussian kernel on the sphere) and the arc-cosine kernel arising from random ReLU features (Cho & Saul, 2009), for which our construction is equivalent to that of the conjugate or NNGP kernel for an infinite-width random ReLU network with the same architecture. The operator M` is then defined pointwise by
  M`x[u] = ϕ`(x[u]). (3) At the first layer on image patches, these kernels lead to functional spaces consisting of homogeneous functions with varying degrees of smoothness on the sphere, depending on the properties of the kernel (Bach, 2017a; Smola et al., 2001). At higher layers, our theoretical analysis will also consider simple polynomial kernels such as k`(z, z′) = (〈z, z′〉)r, in which case the feature map may be explicitly written in terms of tensor products. For instance, r = 2 gives ϕ`(z) = z ⊗ z and H` = (H|S`|`−1)⊗2 = (H`−1 ⊗ H`−1)|S`|×|S`|. See Appendix A.2 for more background on dot-product kernels, their tensor products, and their regularization properties.
  Pooling. Finally, the pooling operators A` perform local averaging through convolution with a filter h`[u], which we may consider to be symmetric (h`[−u]=: h̄`[u]=h`[u]). In practice, the pooling operation is often followed by downsampling by a factor s`, in which case the new signal A`x is defined on a new domain Ω` with |Ω`| = |Ω`−1|/s`, and we may write for x ∈ L2(Ω`−1) and u ∈ Ω`,
  A`x[u] = ∑
  v∈Ω`−1 h`[s`u− v]x[v]. (4)
  Our experiments consider Gaussian pooling filters with a size and bandwidth proportional to the downsampling factor s`, following Mairal (2016), namely, size 2s` + 1 and bandwidth √ 2s`. In Section 3, we will often assume no downsampling for simplicity, in which case we may see the filter bandwidth as increasing with the layers.
  Links with other convolutional kernels. We note that our construction closely resembles kernels derived from infinitely wide convolutional networks, known as conjugate or NNGP kernels (GarrigaAlonso et al., 2019; Novak et al., 2019), and is also related to convolutional neural tangent kernels (Arora et al., 2019; Bietti & Mairal, 2019b; Yang, 2019). The Myrtle family of kernels (Shankar et al., 2020) also resembles our models, but they use small average pooling filters instead of Gaussian filters, which leads to deeper architectures due to smaller receptive fields.
sections.heading: 3 APPROXIMATION WITH (DEEP) CONVOLUTIONAL KERNELS
sections.text: In this section, we present our main results on the approximation properties of convolutional kernels, by characterizing functions in the RKHS as well as their norms. We begin with the one-layer case, which does not capture interactions between patches but highlights the role of pooling, before moving multiple layers, where interaction terms play an important role. Proofs are given in Appendix E.
  1L2(Ω,H) denotes the space of H-valued signals x such that ‖x‖2L2(Ω,H) := ∑ u∈Ω ‖x[u]‖ 2 H <∞.
sections.heading: 3.1 THE ONE-LAYER CASE
sections.text: We begin by considering the case of a single convolutional layer, which can already help us illustrate the role of patches and pooling. Here, the kernel is given by
  K1(x, x ′) = 〈AΦ(x), AΦ(x′)〉L2(Ω,H),
  with Φ(x)[u] = ϕ(xu), where we use the shorthand xu = Px[u] for the patch at position u. We now characterize the RKHS of K1, showing that it consists of additive models of functions inH defined on patches, with spatial regularities among the terms, induced by the pooling operator A. (Notation: A∗ and A† denote the adjoint and pseudo-inverse of an operator A, respectively.) Proposition 1 (RKHS for 1-layer CKN.). The RKHS of K1 consists of functions f(x) = 〈G,Φ(x)〉L2(Ω,H) = ∑ u∈ΩG[u](xu), with G ∈ Range(A∗), and with RKHS norm
  ‖f‖2HK1 = infG∈L2(Ω,H)‖A †∗G‖2L2(Ω,H) s.t. f(x) = ∑ u∈Ω G[u](xu) (5)
  Note that if A∗ is not invertible (for instance in the presence of downsampling), the constraint G ∈ Range(A∗) is active and A†∗ is its pseudo-inverse. In the extreme case of global average pooling, we have A = (1, . . . , 1)⊗ Id : L2(Ω,H)→ H, so that G ∈ Range(A∗) is equivalent to G[u] = g for all u, for some fixed g ∈ H. In this case, the penalty in (5) is simply the squared RKHS norm ‖g‖2H. In order to understand the norm (5) for general pooling, recall that A is a convolution operator with filter h1, hence its inverse (which we now assume exists for simplicity) may be easily computed in the Fourier basis. In particular, for a patch z ∈ Rp|S1|, defining the scalar signal gz[u] = G[u](z), we may write the following using the reproducing property and linearity:
  A†∗G[u](z) = (A−1)>gz[u] = F−1 diag(F h̄1)−1Fgz[u], where h̄1[u] := h1[−u] arises from transposition, F is the discrete Fourier transform, and both F and A are viewed here as |Ω| × |Ω| matrices. From this expression, we see that by penalizing the RKHS norm of f , we are implicitly penalizing the high frequencies of the signals gz[u] for any z, and this regularization is stronger when the pooling filter h1 has a fast spectral decay. For instance, as the spatial bandwidth of h1 increases (approaching a global pooling operation), Fh1 decreases more rapidly, which encourages gz[u] to be more smooth as a function of u, and thus prevents f from relying too much on the location of patches. If instead h1 is very localized in space (e.g., a Dirac filter, which corresponds to no pooling), gz[u] may vary much more rapidly as a function of u, which then allows f to discriminate differently depending on the spatial location. This provides a different perspective on the invariance properties induced by pooling. If we denote G̃[u] = A†∗G[u], the penalty writes
  ‖G̃‖2L2(Ω,H) = ∑ u∈Ω ‖G̃[u]‖2H.
  Here, the RKHS norm ‖ · ‖H also controls smoothness, but this time for functions G̃[u](·) defined on input patches. For homogeneous dot-product kernels of the form (2), the norm takes the form ‖g‖H = ‖T− 12 g‖L2(Sd−1), where the regularization operator T− 1 2 is the self-adjoint inverse square root of the integral operator for the patch kernel restricted to L2(Sd−1). For instance, when the eigenvalues of T decay polynomially, as for arc-cosine kernels, T− 1 2 behaves like a power of the spherical Laplacian (see Bach, 2017a, and Appendix A.2). Then we may write
  ‖G̃‖2L2(Ω,H) = ‖((A−1)> ⊗ T− 1 2 )G‖2L2(Ω)⊗L2(Sd−1),
  which highlights that the norm applies two regularization operators (A−1)> and T− 1 2 independently on the spatial variable and the patch variable of (u, z) 7→ G[u](z), viewed here as an element of L2(Ω)⊗ L2(Sd−1).
sections.heading: 3.2 THE MULTI-LAYER CASE
sections.text: We now study the case of convolutional kernels with more than one convolutional layer. While the patch kernels used at higher layers are typically similar to the ones from the first layer, we show empirically on Cifar10 that they may be replaced by simple polynomial kernels with little loss in accuracy. We then proceed by studying the RKHS of such simplified models, highlighting the role of depth for capturing interactions between different patches via kernel tensor products.
  An empirical study. Table 1 shows the performance of a given 2-layer convolutional kernel architecture, with different choices of patch kernels κ1 and κ2. The reference model uses exponential kernels in both layers, following the construction in Mairal (2016). We find that replacing the second layer kernel by a simple polynomial kernel of degree 3, κ2(u) = u3, leads to roughly the same test accuracy. By changing κ2 to κ2(u) = u2, the test accuracy is only about 1% lower, while doing the same for the first layer decreases it by about 3%. The shallow kernel with a single non-linear convolutional layer (shown in the last line of Table 1) performs significantly worse. This suggests that the approximation properties described in Section 3.1 may not be sufficient for this task, while even a simple polynomial kernel of order 2 at the second layer may substantially improve things by capturing interactions, in a way that we describe below.
  Two-layers with a quadratic kernel. Motivated by the above experiments, we now study the RKHS of a two-layer kernel K2(x, x′) = 〈Ψ(x),Ψ(x′)〉L2(Ω2,H2) with Ψ as in (1) with L = 2, where the second-layer uses a quadratic kernel2 on patches k2(z, z′) = (〈z, z′〉)2. An explicit feature map for k2 is given by ϕ2(z) = z ⊗ z. Denoting by H the RKHS of k1, the patches z lie in H|S2|, thus we may view ϕ2 as a feature map into a Hilbert space H2 = (H ⊗H)|S2|×|S2| (by isomorphism to H|S2| ⊗ H|S2|). The following result characterizes the RKHS of such a 2-layer convolutional kernel, showing that it consists of additive models of interaction functions inH⊗H on pairs of patches, with different spatial regularities on the interaction terms induced by the two pooling operations. (Notation: we use the notations diag(M)[u] = M [u, u] for M ∈ L2(Ω2), diag(x)[u, v] = 1{u=v}x[u] for x ∈ L2(Ω), and Lc is the translation operator Lcx[u] = x[u− c].) Proposition 2 (RKHS of 2-layer CKN with quadratic k2). Let Φ(x) = (ϕ1(xu)⊗ ϕ1(xv))u,v∈Ω ∈ L2(Ω2,H⊗H). The RKHS of K2 when k2(z, z′) = (〈z, z′〉)2 consists of functions of the form
  f(x) = ∑ p,q∈S2 〈Gpq,Φ(x)〉 = ∑ p,q∈S2 ∑ u,v∈Ω Gpq[u, v](xu, xv),
  where Gpq ∈ L2(Ω2,H ⊗ H) obeys the constraints Gpq ∈ Range(Epq) and diag((LpA1 ⊗ LqA1)
  †∗Gpq) ∈ Range(A∗2). Here, Epq : L2(Ω1)→ L2(Ω2) is a linear operator given by Epqx = (LpA1 ⊗ LqA1)∗ diag(x). (6)
  The squared RKHS norm ‖f‖2HK2 is then equal to the minimum over such decompositions of the quantity ∑
  p,q∈S2 ‖A†∗2 diag((LpA1 ⊗ LqA1)†∗Gpq)‖2L2(Ω2,H⊗H). (7)
  As discussed in the one-layer case, the inverses should be replaced by pseudo-inverses if needed, e.g., when using downsampling. In particular, if A∗2 is singular, the second constraint plays a similar role to the one-layer case. In order to understand the first constraint, we show in Figure 2 the outputs of Epqx for Dirac delta signals x[v] = δu[v]. We can see that if the pooling filter h1 has a small support of size m, then Gpq[u− p, v − q] must be zero when |u− v| > m, which highlights that the functions in Gpq may only capture interactions between pairs of patches where the (signed) distance between the first and the second is close to p− q. The penalty then involves operators LpA1 ⊗ LqA1, which may be seen as separable 2D convolutions on the “images” Gpq[u, v]. Then, if z, z′ ∈ Rp|S1| are two fixed patches, defining gz,z′ [u, v] =
  2For simplicity we study the quadratic kernel instead of the homogeneous version used in the experiments, noting that it still performs well (78.0% instead of 79.4% on 10k examples).
  Gpq[u, v](z, z ′), we have, assuming A1 is invertible and symmetric,
  (LpA1 ⊗ LqA1)†∗Gpq[u, v](z, z′) = (A1 ⊗A1)−1gz,z′ [u− p, v − q] = F−12 diag(F2(h1 ⊗ h1))−1F2gz,z′ [u− p, v − q],
  where F2 = F ⊗ F is the 2D discrete Fourier transform. Thus, this penalizes the variations of gz,z′ in both dimensions, encouraging the interaction functions to not rely too strongly on the specific positions of the two patches. This regularization is stronger when the spatial bandwidth of h1 is large, since this leads to a more localized filter in the frequency domain, with stronger penalties on high frequencies. In addition to this 2D smoothness, the penalty in Proposition 2 also encourages smoothness along the p− q diagonal of this resulting 2D image using the pooling operator A2. This has a similar behavior to the one-layer case, where the penalty prevents the functions from relying too much on the absolute position of the patches. Since A2 typically has a larger bandwidth than A1, interaction functionsGpq[u, u+r] are allowed to vary with r more rapidly than with u. The regularity of the resulting “smoothed” interaction terms as a function of the input patches is controlled by the RKHS norm of the tensor product kernel k1 ⊗ k1 as described in Appendix A.2.
  Extensions. When using a polynomial kernel k2(z, z′) = (〈z, z′〉)α with α > 2, we obtain a similar picture as above, with higher-order interaction terms. For example, if α = 3, the RKHS contains functions with interaction terms of the form Gpqr[u, v, w](xu, xv, xw), with a penalty∑
  p,q,r∈S2 ‖A†∗2 diag((A1p ⊗A1q ⊗A1r)†∗Gpqr)‖2L2(Ω2,H⊗3),
  where A1c = LcA1. Similarly to the quadratic case, the first-layer pooling operator encourages smoothness with respect to relative positions between patches, while the second-layer pooling penalizes dependence on the global location. One may extend this further to higher orders to capture more complex interactions, and our experiments suggest that a two-layer kernel of this form with a degree-4 polynomial at the second layer may achieve state-of-the-art accuracy for kernel methods on Cifar10 (see Table 2). We note that such fixed-order choices for κ2 lead to convolutional kernels that lower-bound richer kernels with, e.g., an exponential kernel at the second layer, in the Loewner order on positive-definite kernels. This imples in particular that the RKHS of these “richer” kernels also contains the functions described above. For more than two layers with polynomial kernels, one similarly obtains higher-order interactions, but with different regularization properties (see Appendix D).
sections.heading: 4 GENERALIZATION PROPERTIES
sections.text: In this section, we study generalization properties of the convolutional kernels studied in Section 3, and show improved sample complexity guarantees for architectures with pooling and small patches when the problem exhibits certain invariance properties.
  Learning setting. We consider a non-parametric regression setting with data distribution ρ over (x, y), where the goal is to minimize R(f) = E(x,y)∼ρ[(y − f(x))2]. We denote by f∗ = Eρ[y|x] = arg minf R(f) the regression function, and assume f∗ ∈ H for some RKHS H with kernel K. Without any further assumptions on the kernel, we have the following generalization bound
  on the excess risk for the kernel ridge regression (KRR) estimator, denoted f̂n (see Proposition 7 in Appendix E): E[R(f̂n)−R(f∗)] ≤ C‖f∗‖H √ τ2ρ Ex∼ρX [K(x, x)]
  n , (8)
  where ρX is the marginal distribution of ρ on inputs x, C is an absolute constant, and τ2ρ is an upper bound on the conditional noise variance Var[y|x]. We note that this 1/√n rate is optimal if no further assumptions are made on the kernel (Caponnetto & De Vito, 2007). The quantity Ex∼ρX [K(x, x)] corresponds to the trace of the covariance operator, and thus provides a global control of eigenvalues through their sum, which will already highlight the gains that pooling can achieve. Faster rates can be achieved, e.g., when assuming certain eigenvalue decays on the covariance operator, or when further restricting f∗. We discuss in Appendix F.2 how similar gains to those described in this section can extend to fast rate settings under specific scenarios.
  One-layer CKN with invariance. As discussed in Section 3, the RKHS of 1-layer CKNs consists of sums of functions that are localized on patches, each belonging to the RKHS H of the patch kernel k1. The next result illustrates the benefits of pooling when f∗ is translation invariant.
  Proposition 3 (Generalization for 1-layer CKN.). Assume f∗(x) = ∑ u∈Ω g(xu) with g ∈ H of minimal norm, and assume Ex∼ρX [k1(xu, xv)] ≤ σ2u−v for some (σ2r)r∈Ω. For a 1-layer CKN K1 with any pooling filter h ≥ 0 with ‖h‖1 = 1, we have ‖f∗‖HK1 = √ |Ω|‖g‖H, KRR satisfies
  ER(f̂n)−R(f∗) . |Ω|‖g‖H√
  n √ τ2ρ ∑ r∈Ω 〈h, Lrh〉σ2r . (9)
  The quantities σ2r can be interpreted as auto-correlations between patches at distance r from each other. Note that if h is a Dirac filter, then 〈h, Lrh〉 = 1{r = 0} (recall Lrh[u] = h[u − r]), thus only σ20 plays a role in the bound, while if h is an average pooling filter, we have 〈h, Lrh〉 = 1/|Ω|, so that σ20 is replaced by the average σ̄ 2 := ∑ r σ 2 r/|Ω|. Natural signals commonly display a decay in their auto-correlation functions, suggesting that a similar decay may be present in σ2r as a function of r. In this case, σ̄2 may be much smaller than σ20 , which in turn yields an improved sample complexity guarantee for learning such an f∗ with global pooling, by a factor up to |Ω| in the extreme case where σ2r vanishes for r ≥ 1. In Appendix F.1, we provide simple models where this can be quantified. For more general filters, such as local averaging or Gaussian filters, and assuming σ2r ≈ 0 for r 6= 0, the bound interpolates between no pooling and global pooling through the quantity ‖h‖22. While this yields a worse bound than global pooling on invariant functions, such filters enable learning functions that are not fully invariant, but exhibit some smoothness along the translation group, more efficiently than with no pooling. It should also be noted that the requirement that g belongs to an RKHSH is much weaker when the patches are small, as this typically implies that g admits more than p|S|/2 derivatives, a condition which becomes much stronger as the patch size grows. In the fast rate setting that we study in Appendix F.2, this also leads to better rates that only depend on the dimension of the patch instead of the full dimension.
  Two layers. When using two layers with polynomial kernels at the second layer, we saw in Section 3 that the RKHS of CKNs consists of additive models of interaction terms of the order of the polynomial kernel used. The next proposition illustrates how pooling filters and patch sizes at the second layer may affect generalization on a simple target function consisting of order-2 interactions.
  Proposition 4 (Generalization for 2-layer CKN.). Consider a 2-layer CKN K2 with quadratic k2, as in Proposition 2, and pooling filters h1, h2 with ‖h1‖1 = ‖h2‖1 = 1. Assume that ρX satisfies Ex∼ρX [k1(xu, xu′)k1(xv, xv′)] ≤ if u 6= u′ or v 6= v′, and ≤ 1 otherwise. We have
  Ex∼ρX [K2(x, x)] ≤ |S2|2|Ω| (∑
  v
  〈h2, Lvh2〉〈h1, Lvh1〉2 + ) . (10)
  As an example, consider f∗(x) = ∑ u,v g(xu, xv) for g ∈ H ⊗H of minimal norm. The following table illustrates the obtained generalization bounds R(f̂n)−R(f∗) for KRR with various two-layer architectures (δ: Dirac filter; 1: global average pooling):
  h1 h2 |S2| ‖f∗‖K2 Ex∼ρX [K2(x, x)] Bound ( = 0, τ2ρ = 1) δ δ |Ω| |Ω|‖g‖ |Ω|3 + |Ω|3 ‖g‖|Ω|2.5/√n δ 1 |Ω| |Ω|‖g‖ |Ω|2 + |Ω|3 ‖g‖|Ω|2/√n 1 1 |Ω| √ |Ω|‖g‖ |Ω|+ |Ω|3 ‖g‖|Ω|/√n
  1 δ or 1 1 √ |Ω|‖g‖ |Ω|−1 + |Ω| ‖g‖/√n
  The above result shows that the two-layer model allows for a much wider range of behaviors than the one-layer case, between approximation (through the norm ‖f∗‖K2) and estimation (through Ex∼ρX [K2(x, x)]), depending on the choice of architecture. Choosing the right architecture may lead to large improvements in sample complexity when the target functions has a specific structure, for instance here by a factor up to |Ω|2.5. In Appendix F.1, we discuss simple possible models where we may have a small . Note that choosing filters that are less localized than Dirac impulses, but more than global average pooling, will again lead to different “variance” terms (10), while providing more flexibility in terms of approximation compared to global pooling. This result may be easily extended to higher-order polynomials at the second layer, by increasing the exponents on |S2| and 〈h1, Lvh1〉 to the degree of the polynomial. Other than the gains in sample complexity due to pooling, the bound also presents large gains compared to a “fully-connected” architecture, as in the one-layer case, since it only grows with the norm of a local interaction function inH⊗H that depends on two patches, which may then be small even when this function has low smoothness.
sections.heading: 5 NUMERICAL EXPERIMENTS
sections.text: In this section, we provide additional experiments illustrating numerical properties of the convolutional kernels considered in this paper. We focus here on the Cifar10 dataset, and on CKN architectures based on the exponential kernel. Additional results are given in Appendix B.
  Experimental setup on Cifar10. We consider classification on Cifar10 dataset, which consists of 50k training images and 10k test images with 10 different output categories. We pre-process the images using a whitening/ZCA step at the patch level, which is commonly used for such kernels on images (Mairal, 2016; Shankar et al., 2020; Thiry et al., 2021). This may help reduce the effective dimensionality of patches, and better align the dominant eigen directions to the target function, a property which may help kernel methods (Ghorbani et al., 2020). Our convolutional kernel evaluation code is written in C++ and leverages the Eigen library for hardware-accelerated numerical computations. The computation of kernel matrices is distributed on up to 1000 cores on a cluster consisting of Intel Xeon processors. Computing the full Cifar10 kernel matrix typically takes around 10 hours when running on all 1000 cores. Our results use kernel ridge regression in a one-versus-all approach, where each class uses labels 0.9 for the correct label and −0.1 for the other labels. We report the test accuracy for a fixed regularization parameter λ = 10−8 (we note that the performance typically remains the same for smaller values of λ). The exponential kernel always refers to κ(u) = e 1 σ2 (u−1) with σ = 0.6.
  Varying the kernel architecture. Table 2 shows test accuracies for different architectures compared to Table 1, including 3-layer models and 2-layer models with larger patches. In both cases, the full models with exponential kernels outperform the 2-layer architecture of Table 1, and provide comparable accuracy to the Myrtle10 kernel of Shankar et al. (2020), with an arguably simpler architecture. We also see that using degree-3 or 4 polynomial kernels at the second second layer of the two-layer model essentially provides the same performance to the exponential kernel, and that degree-2 at the second and third layer of the 3-layer model only results in a 0.3% accuracy drop. The two-layer model with degree-2 at the second layer loses about 1% accuracy, suggesting that certain Cifar10 images may require capturing interactions between at least 3 different patches in the image for good classification, though even with only second-order interactions, these models significantly outperform single-layer models. While these results are encouraging, computing such kernels is prohibitively costly, and we found that applying the Nyström approach of Mairal (2016) to these kernels with more layers or larger patches requires larger models than for the architecture of Table 1 for a similar accuracy. Figure 3(left) shows learning curves for different architectures, with slightly better convergence rates for more expressive models involving higher-order kernels or more layers; this suggests that their approximation properties may be better suited for these datasets.
  Role of pooling. Figure 3 shows the spectral decays of the empirical kernel matrix on 1000 Cifar images, which may help assess the “effective dimensionality” of the data, and are related to generalization properties (Caponnetto & De Vito, 2007). While multi-layer architectures with pooling seem to provide comparable decays for various depths, removing pooling leads to significantly slower decays, and hence much larger RKHSs. In particular, the “strided pooling” architecture (i.e., with Dirac pooling filters and downsampling) shown in Figure 3(right), which resembles the kernel considered in (Scetbon & Harchaoui, 2020), obtains less than 40% accuracy on 10k examples. This suggests that the regularization properties induced by pooling, studied in Section 3, are crucial for efficient learning on these problems, as shown in Section 4. Appendix B provides more empirics on different pooling configurations.
sections.heading: 6 DISCUSSION AND CONCLUDING REMARKS
sections.text: In this paper, we studied approximation and generalization properties of convolutional kernels, showing how multi-layer models with convolutional architectures may effectively break the curse of dimensionality on problems where the input consists of high-dimensional natural signals, by modeling localized functions on patches and interactions thereof. We also show how pooling induces additional smoothness constraints on how interaction terms may or may not vary with global and relative spatial locations. An important question for future work is how optimization of deep convolutional networks may further improve approximation properties compared to what is captured by the kernel regime presented here, for instance by selecting well-chosen convolution filters at the first layer, or interaction patterns in subsequent layers, perhaps in a hierarchical manner.
sections.heading: A FURTHER BACKGROUND
sections.text: This section provides further background on the problem of approximation of functions defined on signals, as well as on the kernels considered in the paper. We begin by introducing and motivating the problem of learning functions defined on signals such as images, which captures tasks such as image classification where deep convolutional networks are predominant. We then recall properties of dot-product kernels and kernel tensor products, which are key to our study of approximation.
sections.heading: A.1 NATURAL SIGNALS AND CURSE OF DIMENSIONALITY
sections.text: We consider learning problems consisting of labeled examples (x, y) ∼ ρ from a data distribution ρ, where x is a discrete signal x[u] with u ∈ Ω denoting the position (e.g., pixel location in an image) in a domain Ω, x[u] ∈ Rp (e.g., p = 3 for RGB pixels), and y ∈ R is a target label. In a non-parametric setup, statistical learning may be framed as trying to approximate the regression function
  f∗(x) = Eρ[y|x] using samples from the data distribution ρ. If f∗ is only assumed to be Lipschitz, learning requires a number of samples that scales exponentially in the dimension (see, e.g., von Luxburg & Bousquet (2004); Wainwright (2019)), a phenomenon known as the curse of dimensionality. In the case of natural signals, the dimension d = p|Ω| scales with the size of the domain |Ω| (e.g., the number of pixels), which is typically very large and thus makes this intractable. One common way to alleviate this is to assume that f∗ is smooth, however the order of smoothness typically needs to be of the order of the dimension in order for the problem to become tractable, which is a very strong assumption here when d is very large. This highlights the need for more structured assumptions on f∗ which may help overcome the curse of dimensionality.
  Insufficiency of invariance and stability. Two geometric properties that have been successful for studying the benefits of convolutional architectures are (near-)translation invariance and stability to deformations. Various works have shown that certain convolutional models f yield good invariance and stability Mallat (2012); Bruna & Mallat (2013); Bietti & Mairal (2019a), in the sense that when x̃ is a translation or a small deformation of x, then |f(x̃)− f(x)| is small. Nevertheless, one can show that for band-limited signals (such as discrete signals), ‖x̃ − x‖2 can be controlled in a similar way (though with worse constants, see (Wiatowski & Bölcskei, 2018, Proposition 5)), so that Lipschitz functions on such signals obey such stability properties. Thus, deformation stability is not a much stronger assumption than Lipschitzness, and is insufficient by itself to escape the curse of dimensionality.
  Spatial localization. One successful strategy for learning image recognition models which predates deep learning is to rely on simple aggregations of local features. These may be extracted using hand-crafted procedures (Lowe, 1999; Sánchez et al., 2013; Jégou et al., 2011), or using learned feature extractors, either through learned filters in the early layers of a CNN (Zeiler & Fergus, 2014), or other procedures (e.g., Thiry et al. (2021)). One simplified example that encodes such a prior is if the target function f∗ only depends on the input image through a localized part of the input such as a patch xu = (x[u + v])v∈S ∈ Rp|S|, where S is a small box centered around 0, that is, f∗(x) = g∗(xu). Then, if g∗ is assumed to be Lipschitz, we would like a sample complexity that only scales exponentially in the dimension of a patch p|S|, which is much smaller than dimension of the entire image p|Ω|. This is indeed the case if we use a kernel defined on such patches, such as
  K(x, x′) = ∑ u k(xu, x ′ u),
  where k is a “simple” kernel such as a dot-product kernel, as discussed in Appendix C. In contrast, if K is a dot-product kernel on the entire image, corresponding to an infinite-width limit of a fullyconnected network, then approximation is more difficult and is generally cursed by the full dimension (see Appendix C). While some models of wide fully-connected networks provide some adaptivity to low-dimensional structures such as the variables in a patch (Bach, 2017a), no tractable algorithms are currently known to achieve such behavior provably, and it is reasonable to instead encode such prior information in a convolutional architecture.
  Modeling interactions. Modeling interactions between elements of a system at different scales, possibly hierarchically, is important in physics and complex systems, in order to efficiently handle systems with large numbers of variables (Beylkin & Mohlenkamp, 2002; Hackbusch & Kühn, 2009). As an example, one may consider target functions f∗(x) that consist of interaction functions of the form g(xp, xq), where p, q denote locations of the corresponding patches, and higher-order interactions may also be considered. In the context of image recognition, while functions of a single patch may capture local texture information such as edges or color, such an interaction function may also respond to specific spatial configurations of relevant patches, which could perhaps help identify properties related to the “shape” of an object, for instance. If such functions g are too general, then the curse of dimensionality may kick in again when one considers more than a handful of patches. Certain idealized models of approximation may model such interactions more efficiently through hierarchical compositions (e.g., Poggio et al. (2017)) or tensor decompositions (Cohen & Shashua, 2016; 2017), though no tractable algorithms are known to find such models. In this work, we tackle this in a tractable way using multi-layer convolutional kernels. We show that they can model interactions through kernel tensor products, which define functional spaces that are typically much smaller and more structured than for a generic kernel on the full vector (xp, xq).
sections.heading: A.2 DOT-PRODUCT KERNELS AND THEIR TENSOR PRODUCTS
sections.text: In this section, we review some properties of dot-product kernels, their induced RKHS and regularization properties. We then recall the notion of tensor product of kernels, which allows us to describe the RKHS of products of kernels in terms of that of individual kernels.
  Dot-product kernels. The rotation-invariance of dot-product kernels provides a natural description of their RKHS in terms of harmonic decompositions of functions on the sphere using spherical harmonics (Smola et al., 2001; Bach, 2017a). This leads to natural connections with regularity properties of functions defined on the sphere. For instance, if the kernel integral operator on L2(Sd−1) has a polynomially decaying spectral decay, as is the case for kernels arising from the ReLU activation (Bach, 2017a; Bietti & Mairal, 2019b), then the RKHS contains functions g ∈ L2(Sd−1) with an RKHS norm equivalent to ‖∆β/2Sd−1g‖L2(Sd−1), (11) for some β that depends on the decay exponent and must be larger than (d − 1)/2, with ∆Sd−1 the Laplace-Beltrami operator on the sphere. This resembles a Sobolev norm of order β, and the RKHS contains functions with bounded derivatives up to order β. When d is small (e.g., at the first layer with small images patches), the space contains functions that need not be too regular, and may thus be quite discriminative, while for large d (e.g., for a fully-connected network), the functions must be highly smooth in order to be in the RKHS, and large norms are necessary to approach
  non-smooth functions. For kernels with decays faster than polynomial, such as the Gaussian kernel, the RKHS contains smooth functions, but may still provide good approximation to non-smooth functions, particularly with small d and when using small bandwidth parameters. The homogeneous case (2) leads to functions f(x) = ‖x‖g( x‖x‖ ) with g defined on the sphere, with a norm given by the same penalty (11) on the function g (Bietti & Mairal, 2019b).
  Kernel tensor products. For more than one layer, the convolutional kernels we study in Section 3 can be expressed in terms of products of kernels on patches, of the form
  K((x1, . . . , xm), (x ′ 1, . . . , x ′ m)) = m∏ j=1 k(xj , x ′ j), (12)
  where x1, . . . , xm, x′1, . . . , x ′ m ∈ Rd are patches which may come from different signal locations. If ϕ : Rd → H is the feature map into the RKHSH of k, then ψ(x1, . . . , xm) = ϕ(x1)⊗ · · · ⊗ ϕ(xm)
  is a feature map for K, and the corresponding RKHS, denoted H⊗m = H⊗ · · · ⊗ H, contains all functions
  f(x1, . . . , xm) = n∑ i=1 gi,1(x1) . . . gi,m(xm),
  for some n, with gi,m ∈ H for i ∈ [n] and j ∈ [m] (see,e.g., (Wainwright, 2019, Section 12.4.2) for a precise construction). The resulting RKHS is often much smaller than for a more generic kernel on Rd×m; for instance, if H is a Sobolev space of order β in dimension d, then H⊗m is much smaller than the Sobolev space of order β in d×m dimensions, and corresponds to stronger, mixed regularity conditions (see, e.g., Bach, 2017b; Sickel & Ullrich, 2009, for the d = 1 case). This can yield improved generalization properties if the target function has such a structure (Lin, 2000). Kernels of the form (12) and sums of such kernels have been useful tools for avoiding the curse of dimensionality by encoding interactions between variables that are relevant to the problem at hand (Wahba, 1990, Chapter 10). In what follows, we show how patch extraction and pooling operations shape the properties of such interactions between patches in convolutional kernels through additional spatial regularities.
sections.heading: B ADDITIONAL EXPERIMENTS
sections.text: In this section, we provide additional experiments to those presented in Section 5, using different patch kernels, patch sizes, pooling filters, preprocessings, and datasets.
  Three-layer architectures with different patch kernels. Table 3 provides more results on 3-layer architectures compared to Table 2, including different changes in the degrees of polynomial kernels at the second and third layer. In particular we see that the architecture with degree-2 kernels at both layers, which captures interactions of order 4, also outperforms the simpler ones using degree-4 kernels at either layer, suggesting that a deeper architecture may better model relevant interactions terms on this problem.
  Varying the second layer patch size. Table 4 shows the variations in test performance when changing the size of the second patches at the second layer. We see that intermediate sizes between 3x3 and 9x9 work best, but that performance degrades when using patches that are too large or too small. For very large patches, this may be due to the large variance in (10), or perhaps instability (Bietti & Mairal, 2019a). For |S2| =1x1, note that while pooling after the first layer allows even 1x1 patches to capture interactions across different input image patches, these may be limited to short range interactions when the pooling filter is localized (see Proposition 2), which may limit the expressivity of the model.
  Arc-cosine kernel. In Table 5, we consider 2-layer convolutional kernels with a similar architecture to those considered in Table 1, but where we use arc-cosine kernels arising from ReLU activations instead of the exponential kernel used in Section 5, given by
  κ(u) = 1
  π (u · (π − arccos(u)) +
  √ 1− u2).
  The obtained convolutional kernel then corresponds to the conjugate kernel or NNGP kernel arising from an infinite-width convolutional network with the ReLU activation (Daniely et al., 2016; GarrigaAlonso et al., 2019; Novak et al., 2019). We may also consider the neural tangent kernel (NTK) for the same architecture, which additionally involves arc-cosine kernels of degree 0, which correspond to random feature kernels for step activations u 7→ 1{u ≥ 0}. We find that the NTK performs slightly better than the conjugate kernel, but both kernels achieve lower accuracy compared to the Exponential kernel shown in Table 1. Nevertheless, we observe a similar pattern regarding the use of polynomial kernels at the second layer, namely, the drop in accuracy is much smaller when using a quadratic kernel compared to a linear kernel, suggesting that non-linear kernels on top of the first layer, and the interactions they may capture, are crucial on this dataset for good accuracy.
  One-layer architectures and larger initial patches. Table 6 shows the accuracy for one-layer convolutional kernels with 6x6 patches3 and various pooling sizes, with a highest accuracy of 75.8% for a pooling size of 8. While this improves on the accuracy obtained with 3x3 patches (slightly
  3Note that in this case the ZCA/whitening step is applied on these larger 6x6 patches.
  above 74% for the architectures in Tables 1 and 3 with a single non-linear kernel at the first layer), these accuracies remain much lower than those achieved by two-layer architectures with even quadratic kernels at the second layer. While using larger patches may allow capturing patterns that are less localized compared to small 3x3 patches, the neighborhoods that they model need to remain small in order to avoid the curse of dimensionality when using dot-product kernels, as discussed in Section 2. Instead, the multi-layer architecture may model information at larger scales with a much milder dependence on the size of the neighborhood, thanks to the structure imposed by tensor product kernels (see Section A.2) and the additional regularities induced by pooling.
  We also found that larger patches at the first layer may hurt performance in multi-layer models: when considering the architecture of Table 1 with exponential kernels, using 5x5 patches instead of 3x3 at the first layer yields an accuracy of 79.6% instead of 80.5% on Cifar10 when training on the same 10k images. This again reflects the benefits of using small patches at the first layer for allowing better approximation on small neighborhoods, while modeling larger scales using interaction models according to the structure of the architecture. We note nevertheless that for standard deep networks, larger patches are often used at the first layer (e.g., He et al., 2016), as the feature selection capabilities of SGD may alleviate the dependence on dimension, e.g., by finding Gabor-like filters.
  Gaussian vs average pooling. Table 7 shows the differences in performance between two or three layer architectures considered in Table 2, when Gaussian pooling filters are replaced by average pooling filters. For both architectures considered, average pooling leads to a significant performance drop. This suggests that one may need deeper architectures in order for such average pooling filters to work well, as in (Shankar et al., 2020), either with multiple 3x3 convolutional layers before applying pooling, or by applying multiple average pooling layers in a row as in certain Myrtle kernels. Note that iterating multiple average pooling layers in a row is equivalent to using a larger and more smooth pooling filter (with one more order of smoothness at each layer), which may then be more comparable to our Gaussian pooling filters.
  SVHN dataset. We now consider the SVHN dataset, which consists of 32x32 images of digits from Google Street View images, 73 257 for training and 26 032 for testing. Due to the larger dataset size, we only consider the kernel approximation approach of Mairal (2016) based on the Nyström method, which projects the patch kernel feature maps at each layer to finite-dimensional subspaces generated by a set of anchor points (playing the role of convolutional filters), themselves computed via a K-means clustering of patches.4 We train one-versus-all classifiers on the resulting finite-dimensional representations using regularized ERM with the squared hinge loss, and simply report the best test accuracy over a logarithmic grid of choices for the regularization parameter, ignoring model selection issues in order to assess approximation properties. We use the same ZCA preprocessing as on Cifar10 and the same architecture as in Table 1, with a relatively small number of filters (256 at the first layer,
  4We use the PyTorch implementation available at https://github.com/claying/ CKN-Pytorch-image.
  4096 at the second layer, leading to representations of dimension 65 536), noting that the accuracy can further improve when increasing this number. Our observations are similar to those for the Cifar10 dataset: using a degree-3 polynomial kernel at the second layer reaches very similar accuracy to the exponential kernel; using a degree-2 polynomial leads to a slight drop, but a smaller drop than when making this same change at the first layer; using a linear kernel at the second layer leads to a much larger drop. This again highlights the importance of using non-linear kernels on top of the first layer in order to capture interactions at larger scales than the scale of a single patch.
  Local versus global whitening. Recall that our pre-processing is based on a patch-level whitening or ZCA on each image, following Mairal (2016). In practice, this is achieved by whitening extracted patches from each image, and reconstructing the image from whitened patches via averaging. In contrast, other approaches use global whitening of the entire image Lee et al. (2020); Shankar et al. (2020). For the 2-layer model shown in Table 2 with 5x5 patches at the second layer, we found global ZCA to provide significantly worse performance, with a drop from 88.3% to about 80%.
  Finite networks and comparison to Shankar et al. (2020). The work Shankar et al. (2020) introduces Myrtle kernels but also consider similar architectures for usual CNNs with finite-width, trained with stochastic gradient descent. Obtaining competitive architectures for the finite-width case is not the goal of our work, which focuses on good architectures for the kernel setup, yet it remains interesting to consider this question. In the case of Shankar et al. (2020), training the finite-width networks yields better accuracy compared to their “infinite-width” kernel counterparts, a commonly observed phenomenon which may be due to better “adaptivity” of optimization algorithms compared to kernel methods, which have a fixed representation and thus may not learn representations adapted to the data (see, e.g., Allen-Zhu & Li (2020); Bach (2017a); Chizat et al. (2019)). Nevertheless, we found that for the two-layer architecture considered in Table 1, which has many fewer layers compared to the Myrtle architectures of Shankar et al. (2020), using a finite-width ReLU network yields poorer performance compared to the kernel (around 83% at best, compared to 87.9%). This may suggest that for convolutional networks, deeper networks may have additional advantages when using optimization algorithms, in terms of adapting to possibly relevant structure of the problem, such as hierarchical representations (see, e.g., Allen-Zhu & Li (2020); Chen et al. (2020); Poggio et al. (2017) for theoretical justifications of the benefits of depth in non-kernel regimes).
sections.heading: C COMPLEXITY OF SPATIALLY LOCALIZED FUNCTIONS
sections.text: In this section, we briefly elaborate on our discussion in Section A.1 on how simple convolutional structure may improve complexity when target functions are spatially localized. We assume f∗(x) = g∗(xu) with xu = (x[u+ v])v∈S ∈ Rp|S| a patch of size |S|, where g∗ is a Lipschitz function. If we define the kernel Ku(x, x′) = k(xu, x′u), where k is a dot-product kernel arising from a one-hidden layer network with positively-homogeneous activation such as the ReLU, and further assume patches to be bounded and g∗ to be bounded, then the uniform approximation error bound of Bach (2017a, Proposition 6) together with a simple O(1/ √ n) Rademacher complexity bound on estimation error shows that we may achieve a generalization bound with a rate that only depends on the patch dimension p|S| rather than p|Ω| in this setup (i.e., a sample complexity that is exponential in p|S|, which is much smaller than p|Ω|). If we consider the kernel K(x, x′) = ∑ u∈Ω k(xu, x ′ u), the RKHS contains all functions in the RKHS of Ku for all u ∈ Ω, with the same norm (this may be seen as an application of Theorem 6 with a feature map given by concatenating the kernel maps of each Ku), so that we may achieve the same approximation error as above, and thus a similar generalization bound that is not cursed by dimension. This kernel also allows us to obtain similar generalization guarantees when f∗ consists of linear combinations of such spatially localized functions on different patches within the image.
  In contrast, when using a similar dot-product kernel on the full signal, corresponding to using a fully-connected network in a kernel regime, one may construct functions f∗(x) = g∗(xu) with g∗ Lipschitz where an RKHS norm that is exponentially large in the (full) dimension p|Ω| is needed for a small approximation error (see (Bach, 2017a, Appendix D.5)).
  Related to this, Malach & Shalev-Shwartz (2021) show a separation in the different setting of learning certain parity functions on the hypercube using gradient methods; their upper bound for convolutional
  networks is based on a similar kernel regime as above. We note that kernels that exploit such a localized structure have also been considered in the context of structured prediction for improved statistical guarantees (Ciliberto et al., 2019).
sections.heading: D EXTENSIONS TO MORE LAYERS
sections.text: In this section, we study the RKHS for convolutional kernels with more than 2 convolutional layers, by considering the simple example of a 3-layer convolutional kernel K3 defined by the feature map
  Ψ(x) = A3M3P3A2M2P2A1M1P1x,
  with quadratic kernels at the second and third layer, i.e., k2(z, z′) = (〈z, z′〉)2 and k3(z, z′) = (〈z, z′〉)2. By isomorphism, we may consider the sequence of Hilbert spaces H` to be H1 = H, H2 = (H ⊗ H)|S2|×|S2|, and H3 = (H⊗4)(|S3|×|S2|×|S2|) 2
  . For some domain Ω, we define the operators diag2 : L
  2(Ω4)→ L2(Ω2) and its adjoint diag2 : L2(Ω2)→ L2(Ω4) by diag2(M)[u, v] = M [u, u, v, v] for M ∈ L2(Ω4)
  diag2(M)[u1, u2, u3, u4] = 1{u1 = u2}1{u3 = u4}M [u1, u3] for M ∈ L2(Ω2). We may then describe the RKHS as follows. Proposition 5 (RKHS of 3-layer CKN with quadratic k2/3). The RKHS of K3 when k2 and k3 are quadratic kernels (〈·, ·〉)2 consists of functions of the form
  f(x) = ∑
  α∈(S3×S2×S2)2 ∑ u1,u2,u3,uv∈Ω Gα[u1, u2, u3, u4](xu1 , xu2 , xu3 , xu4), (13)
  where Gα ∈ L2(Ω4,H⊗4) obeys the constraint Gα ∈ Range(Eα), (14)
  where the linear operator Eα : L2(Ω3) → L2(Ω4) for α = (p, q, r, p′, q′, r′) (with p, p′ ∈ S3 and q, r, q′, r′ ∈ S2) is defined by
  Eαx = A ∗ 1,α diag2(A ∗ 2,α diag(A ∗ 3x)).
sections.heading: The operators A1,α and A2,α denote:
sections.text: A1,α = LqA1 ⊗ LrA1 ⊗ Lq′A1 ⊗ Lr′A1 A2,α = LpA2 ⊗ Lp′A2.
  The squared RKHS norm ‖f‖2HK3 is then equal to the minimum over decompositions (13) of the quantity ∑
  α
  ‖A†∗3 diag(A†∗2,α diag2(A†∗1,αGα))‖2L2(Ω3). (15)
  The constraint (14) and penalty (15) resemble the corresponding constraint/penalty in the twolayer case for an order-4 polynomial kernel at the second layer, but provide more structure on the interactions, using a multi-scale structure that may model interactions between certain pairs of patches ((xu1 , xu2) and (xu3 , xu4) in (13)) more strongly than those between all four patches. In addition to localizing the interactions Gα around certain diagonals, the kernel also promotes spatial regularities: assuming that the spatial bandwidths of A` increase with `, the functions (u, v, w1, w2) 7→ Gα[u, u+ w1, u+ v, u+ v+w2] may vary quickly with w1 or w2 (distances between patches in each of the two pairs), but should vary more slowly with v (distance between the two pairs) and even more slowly with u (a global position).
sections.heading: E PROOFS
sections.text: We recall the following result about reproducing kernel Hilbert spaces, which characterizes the RKHS of kernels defined by explicit Hilbert space features maps (see, e.g., Saitoh, 1997, §2.1).
  Theorem 6 (RKHS from explicit feature map). Let H be some Hilbert space, ψ : X → H a feature map, and K(x, x′) = 〈ψ(x), ψ(x′)〉H a kernel on X . The RKHS H of K consists of functions f = 〈g, ψ(·)〉H , with norm
  ‖f‖H = inf{‖g′‖H : g′ ∈ H s.t. f = 〈g′, ψ(·)〉H} (16)
  We also state here the generalization bound for kernel ridge regression used in Section 4, adapted from Bach (2021, Proposition 7.1).
  Proposition 7 (Generalization bound for kernel ridge regression). Denote f∗(x) = Eρ[y|x]. Assume f∗ ∈ H, Varρ[y|x] ≤ σ2, K(x, x) ≤ 1 a.s., and define
  f̂λ := arg min f∈H
  1
  n n∑ i=1 (yi − f(xi))2 + λ‖f‖2H, (17)
  for i.i.d. data (xi, yi) ∼ ρ, i = 1, . . . , n. Let n ≥ max {‖f∗‖∞
  σ2 , ‖f∗‖2H σ2 EρX [K(x, x)] , σ2 EρX [K(x, x)] ‖f∗‖2H
  } .
  For λ = √ σ2 EρX [K(x, x)]/n‖f∗‖2H, we have
  E[R(f̂λ)−R(f∗)] ≤ C‖f∗‖H √ σ2 EρX [K(x, x)]
  n , (18)
  where C is an absolute constant.
  Proof. Under the conditions of the theorem, we may apply (Bach, 2021, Proposition 7.1), which states that for λ ≤ 1 and n ≥ 5λ (1 + log(1/λ)), we have
  E[R(f̂λ)−R(f∗)] ≤ 16 σ2
  n Tr((Σ + λI)−1Σ) + 16λ〈f∗, (Σ + λI)−1Σf∗〉H +
  24 n2 ‖f∗‖2∞,
  where Σ = EρX [K(x, ·)⊗K(x, ·)] is the covariance operator. We conclude by using the inequalities
  Tr((Σ + λI)−1Σ) ≤ Tr(Σ) λ = EρX [K(x, x)] λ
  〈f∗, (Σ + λI)−1Σf∗〉H ≤ ‖f∗‖2H, and optimizing for λ.
sections.heading: E.1 PROOF OF PROPOSITION 1 (RKHS OF ONE-LAYER CONVOLUTIONAL KERNEL)
sections.text: Proof. From Theorem 6, the RKHS contains functions of the form
  f(x) = 〈F,AΦ(x)〉L2(Ω1,H), with RKHS norm equal to the minimum of ‖F‖L2(Ω1,H) over such decompositions. We may alternatively write f(x) = 〈G,Φ(x)〉L2(Ω,H) with G = A∗F . The mapping from F to G is one-to-one if G ∈ Range(A∗). Then, we obtain that equivalently, the RKHS contains functions of this form, with G ∈ Range(A∗), and with RKHS norm equal to the minimum of ‖A†∗G‖L2(Ω1,H) over such decompositions.
  E.2 PROOF OF PROPOSITION 2 (RKHS OF 2-LAYER CKN WITH QUADRATIC k2)
  Proof. From Theorem 6, the RKHS contains functions of the form
  f(x) = 〈F,A2M2P2A1Φ1(x)〉L2(Ω2,H2), with RKHS norm equal to the minimum of ‖F‖L2(Ω2,H2) over such decompositions. Here, Φ1(x) ∈ L2(Ω,H) is given by Φ1(x)[u] = ϕ1(xu), so that Φ(x) in the statement is given by Φ(x) =
  Φ1(x)⊗ Φ1(x). We also have thatH2 = (H⊗H)|S2|×|S2|, so that we may write F = (Fpq)p,q∈S2 with Fpq ∈ L2(Ω2,H⊗H). For p, q ∈ S2, denoting by Lc the translation operator Lcx[u] = x[u− c], we have
  (M2P2A1Φ1(x)[u])pq = LpA1Φ1(x)[u]⊗ LqA1Φ1(x)[u] = diag(LpA1Φ1(x)⊗ LqA1Φ1(x))[u] = diag((LpA1 ⊗ LqA1)Φ(x))[u].
  Then, we have
  〈Fpq, (A2M2P2A1Φ1(x))pq〉L2(Ω2,H⊗H) = 〈Fpq, A2 diag((LpA1 ⊗ LqA1)Φ(x))〉L2(Ω2,H⊗H) = 〈A∗2Fpq,diag((LpA1 ⊗ LqA1)Φ(x))〉L2(Ω1,H⊗H) = 〈diag(A∗2Fpq), (LpA1 ⊗ LqA1)Φ(x)〉L2(Ω21,H⊗H) = 〈(LpA1 ⊗ LqA1)∗ diag(A∗2Fpq),Φ(x)〉L2(Ω2,H⊗H).
  We may then write this as 〈Gpq,Φ(x)〉L2(Ω2,H⊗H) with
  Gpq = (LpA1 ⊗ LqA1)∗ diag(A∗2Fpq),
  and the mapping between Fpq ∈ L2(Ω2,H ⊗ H) and Gpq ∈ L2(Ω2,H ⊗ H) is one-to-one if Gpq ∈ Range((LpA1 ⊗ LqA1)∗), and diag((LpA1 ⊗ LqA1)†∗Gpq) ∈ Range(A∗2). We may then equivalently write the RKHS norm as the minimum over Gpq satisfying such constraints for all p, q ∈ S2, of the quantity∑
  p,q∈S2 ‖Fpq‖2 = ∑ p,q∈S2 ‖A†∗2 diag((LpA1 ⊗ LqA1)†∗Gpq)‖2L2(Ω2,H⊗H).
  E.3 PROOF OF PROPOSITION 5 (RKHS OF 3-LAYER CKN WITH QUADRATIC k2/3)
  Proof. Let Φ(x) = (ϕ1(xu))u ∈ L2(Ω,H), so that we may write∑ u1,u2,u3,uv∈Ω G[u1, u2, u3, u4](xu1 , xu2 , xu3 , xu4) = 〈G,Φ(x)⊗4〉L2(Ω4,H⊗4),
  for some G ∈ L2(Ω4,H⊗4). From Theorem 6, the RKHS contains functions of the form
  f(x) = 〈F,A3M3P3A2Φ2(x)〉L2(Ω3,H3), (19)
  with RKHS norm equal to the minimum of ‖F‖L2(Ω3,H3) over such decompositions. Here, Φ2(x) ∈ L2(Ω1,H2) = L2(Ω1, (H⊗H)|S2|×|S2|) is given as in the proof of Proposition 2, by
  Φ2,qr(x)[u] = diag((LqA1 ⊗ LrA1)(Φ(x)⊗ Φ(x)))[u], for q, r ∈ S2. A patch P3A2Φ2(x)[u] is then given by
  P3A2Φ2(x)[u] = (LpA2Φ2,qr(x)[u])p∈S3,q,r∈S2 ∈ (H⊗H)|S3|×|S2|×|S2|.
  Applying the quadratic feature map given by ϕ3(z) = z ⊗ z ∈ (H⊗4)(|S3|×|S2|×|S2|) 2 for z ∈ (H⊗H)|S3|×|S2|×|S2|, we obtain for α = (p, q, r, p′, q′, r′) ∈ (S3 × S2 × S2)2,
  (M3P3A2Φ2(x)[u])α = LpA2Φ2,qr(x)[u]⊗ Lp′A2Φ2,q′r′(x)[u] = diag(A2,α(Φ2,qr(x)⊗ Φ2,q′r′(x)))[u],
  where A2,α = LpA2 ⊗ Lp′A2.
  Now, one can check that we have the following relation:
  Φ2,qr(x)⊗ Φ2,q′r′(x) = diag((LqA1 ⊗ LrA1)(Φ(x)⊗ Φ(x)))⊗ diag((Lq′A1 ⊗ Lr′A1)(Φ(x)⊗ Φ(x))) = diag2(A1,αΦ(x) ⊗4),
  with A1,α = LqA1 ⊗ LrA1 ⊗ Lq′A1 ⊗ Lr′A1.
  Since H3 = (H⊗4)(|S3|×|S2|×|S2|) 2 , we may write F = (Fα)α∈(S3×S2×S2)2 , with each Fα ∈ L2(Ω3,H⊗4). We then have
  〈Fα,(A3M3P3A2Φ2(x))α〉L2(Ω3) = 〈Fα, A3 diag(A2,α diag2(A1,αΦ(x)⊗4))〉L2(Ω3) = 〈diag(A∗3Fα), A2,α diag2(A1,αΦ(x)⊗4)〉L2(Ω22) = 〈diag2(A∗2,α diag(A∗3Fα)), A1,αΦ(x)⊗4〉L2(Ω41) = 〈A∗1,α diag2(A∗2,α diag(A∗3Fα)),Φ(x)⊗4〉L2(Ω4).
  We may write this as 〈Gα,Φ(x)⊗4〉L2(Ω4,H⊗4), with Gα = A ∗ 1,α diag2(A ∗ 2,α diag(A ∗ 3Fα)).
  The mapping from Fα to Gα is bijective if Gα is constrained to the lie in the range of the operator Eα. If Gα satisfies this constraint, we may write
  Fα = A †∗ 3 diag(A †∗ 2,α diag2(A †∗ 1,αGα)).
  Then, the resulting penalty on Gα is as desired.
sections.heading: E.4 PROOF OF PROPOSITION 3 (GENERALIZATION FOR ONE-LAYER CKN)
sections.text: Proof. Note that we have Ex[K1(x, x)] = ∑ u ∑ v,r h[u− v]h[u− v − r]Ex[k1(xv, xv−r)]
  ≤ ∑ v,r 〈h, Lrh〉σ2r
  = |Ω| ∑ r 〈h, Lrh〉σ2r .
  It remains to verify that ‖f∗‖HK1 = √ |Ω|‖g‖H. Note that if we denote G = (g)u∈Ω ∈ L2(Ω,H),
  then we haveA∗G = G, since ∑ v h[v−u]G[v] = ( ∑ v h[v−u])g = g. This implies thatA∗†G = G, regardless of which pooling filter is used. Then we have, by (5) that ‖f∗‖2 ≤ |Ω|‖g‖2. Further, since g is of minimal norm, no other G ∈ L2(Ω,H) may lead to a smaller norm, so that we can conclude ‖f∗‖2 = |Ω|‖g‖2.
sections.heading: E.5 PROOF OF PROPOSITION 4 (GENERALIZATION FOR TWO-LAYER CKN WITH
sections.text: QUADRATIC k2)
  Proof. We begin by studying the “variance” quantity Ex[K2(x, x)]. By expanding the construction of the kernel K2, we may write
  K2(x, x) = ∑ p,q∈S2 ∑ u,v,v′
  h2[u−v]h2[u−v′]×∑ w1,w2 w′1,w ′ 2 h1[v−w1]h1[v−w2]h1[v′−w′1]h1[v′−w′2]k1(xw1−p, xw′1−p)k1(xw2−q, xw′2−q).
  Upper bounding the quantity E[k1(xw1−p, xw′1−p)k1(xw2−q, xw′2−q)] by 1 when w1 = w ′ 1 and w2 = w′2, and by otherwise, the sum of the coefficients in front of 1 can be bounded as follows:∑ p,q∈S2 ∑ u,v,v′,w1,w2 h2[u−v]h2[u−v′]h1[v−w1]h1[v−w2]h1[v′−w1]h1[v′−w2]
  = |S2|2 ∑ u,v,v′ ∑ v h2[u−v]h2[u−v]〈Lvh1, Lv′h1〉2
  = |S2|2 ∑ v,v′ 〈Lvh2, Lv′h2〉〈Lvh1, Lv′h1〉2
  = |Ω||S2|2 ∑ v 〈h2, Lvh2〉〈h1, Lvh1〉2.
  while the sum of coefficients bounded by is upper bounded by |Ω||S2|2. Overall, this yields Ex[K2(x, x)] ≤ |S2|2|Ω| (∑
  v
  〈h2, Lvh2〉〈h1, Lvh1〉2 + ) . (20)
  The bounds obtained for the example function f∗(x) = ∑ u,v g(xu, xv) rely on plugging in the values of the Dirac filter h[u] = δ(u = 0) or average pooling filters h[u] = 1/|Ω| in the expression of Ex[K2(x, x)], and on computing the norm ‖f∗‖K2 for different architectures using Eq. (7) in Proposition 2. Computing Ex[K2(x, x)] is immediate using the expression above. Computing ‖f∗‖K2 is more involved, and requires finding appropriate decompositions of the form f∗(x) = ∑ p,q ∑ u,v Gpq[u, v](xu, xv) in order to leverage Proposition 2:
  • When using a Dirac filter at the first layer, we need |S2| = |Ω| in order to capture logrange interaction terms, and we represent f∗ as a sum of Gp that are non-zero and equal to g only on the p-th diagonal, i.e., Gp[u, v] = g when v = u + p, and zero otherwise. We then verify ∑ p ∑ u,v Gp[u, v](xu, xv) = ∑ p ∑ u g(xu, xu+p) = f
  ∗(x). Then, the expression (7) on this decomposition yields |S2||Ω|‖g‖2H⊗H = |Ω|2‖g‖2H⊗H, for any choice of pooling h2 such that ‖h2‖1 = 1 (using similar arguments to the proof of Proposition 3). This is then equal to the squared norm of f∗ due to the minimality of ‖g‖H⊗H.
  • When using average pooling at the first layer and |S2| = |Ω|, we may use a single termG[u, v] with all entries equal to g, i.e., a decomposition f∗(x) = ∑ u,v G[u, v](xu, xv).
  Using (7), we obtain |Ω|‖g‖2 for the squared norm. The same decomposition can be used when |S2| = 1, leading to the same norm.
sections.heading: F GENERALIZATION GAINS UNDER SPECIFIC DATA MODELS
sections.text: In this section, we consider simple models of architectures and data distribution where we may quantify more precisely the improvements in sample complexity guarantees thanks to pooling.
  We consider architectures with non-overlapping patches, and a data distribution where the patches are independent and uniformly distributed on the sphere Sd−1 sphere in d := p|S1| dimensions. Further, we consider a dot-product kernel on patches of the form k1(z, z′) = κ(〈z, z′〉) for z, z′ ∈ Sd−1, with the common normalization κ(1) = 1.
  In the construction of Section 2, using non-overlapping patches corresponds to taking a downsampling factor s1 = |S1| at the first layer, and a corresponding pooling filter such that h1[u] = 0 for u 6= 0 mod |S1|. Alternatively, we may more simply denote patches by xu with u ∈ Ω by considering a modified signal with more channels (xu = x[u] of dimension pe instead of p, where e is the patch size) so that extracting patches of size 1 actually corresponds to a patch of size e of the underlying signal. We then have that the patches xu in a signal x are i.i.d., uniformly distributed on the sphere Sd−1. We denote the uniform measure on Sd−1 by dτ .
  We note that when patches are in high dimension, overlapping patches may become near-orthogonal, which could allow extensions of our arguments below to the case with overlap, yet this may require different tools similar to Mei et al. (2021). We leave these questions to future work.
  F.1 QUANTIFYING THE TRACE OF THE COVARIANCE OPERATOR E[K(x, x)]
  In this section, we focus on the “variance” term E[K(x, x)] which is used in the generalization results of Section 4.
  One layer. In the one-layer case, we clearly have σ20 := E[κ(〈xu, xu〉)] = κ(1) = 1. For u 6= v, since patches xu and xv are independent and i.i.d., σ2u−v is a constant independent of u, v, which may be computed by integration on the sphere as:
  σ2u−v = E[κ(〈xu, xv〉)] = Exu∼τ [Exv∼τ [κ(〈xu, xv〉)|xu]] = ωd−2 ωd−1 ∫ 1 −1 κ(t)(1− t2) d−32 dt, (21)
  where we used a standard change of variable t = 〈xu, xv〉 when integrating xv over Sd−1 (see, e.g., Efthimiou & Frye (2014)), with ωp−1 = 2πp/2/Γ(p/2) the surface measure of the sphere Sp−1. Note that the integral in (21) corresponds to the constant component in the Legendre decomposition of κ, which is known for common kernels as consider in various works studying spectral properties of dot-product kernels Bach (2017a); Minh et al. (2006). For instance, for the exponential kernel (or Gaussian on the sphere) κ(〈x, y〉) = e− ‖x−y‖2 2σ2 = e 1 σ2
  (〈x,y〉−1), which is used in most of our experiments with σ = 0.6, we have (Minh et al., 2006, Theorem 2):
  ωd−2 ωd−1 ∫ 1 −1 κ(t)(1− t2) d−32 dt = e−1/σ2(2σ2)(d−2)/2Id/2−1(1/σ2)Γ(d/2),
  where I denotes the modified Bessel function of the first kind. For arc-cosine kernels, it may be obtained by leveraging the random feature expansion of the kernel Bach (2017a). More generally, we also note that when the patch dimension d is large, we have
  ωd−2 ωd−1 ∫ 1 −1 κ(t)(1− t2) d−32 dt→ κ(0), as d→∞,
  since ωd−2ωd−1 (1 − t 2) d−3 2 is a probability density that converges weakly to a Dirac mass at 0. In particular, for the exponential kernel with σ = 0.6, we have κ(0) = e−1/σ 2 ≈ 0.06. When learning a translation-invariant function, the bound in Prop. 3 then shows that global average pooling yields an improvement w.r.t. no pooling of order |Ω|/(1+0.06|Ω|). Note that removing the constant component of κ, i.e., using the kernel κ(u)−κ(0)κ(1)−κ(0) , may further improve this bound, leading to a denominator very close to 1 when d is large, and hence an improvement in sample complexity of order |Ω|. We also remark that the dependence on κ(0) may be removed by using a finer generalization analysis beyond uniform convergence that leverages spectral properties of the kernel (see Section F.2).
  Two layers with quadratic k2. For the two-layer case, we may obtain expressions of E[k1(xu, xu′)k1(xv, xv′)] as above. Denote := Ez,z′∼τ [k1(z, z′)] where z, z′ are independent, which is given in (21). We may have the following cases:
  • If u = u′ and v = v′, we have, trivially, E[k1(xu, xu′)k1(xv, xv′)] = 1. • If u = u′ and v 6= v′, we have E[k1(xu, xu′)k1(xv, xv′)] = E[k1(xv, xv′)] = , since xv and xv′
  are independent. The same holds if u 6= u′ and v = v′. • If |{u, v, u′, v′}| = 4, we have
  E[k1(xu, xu′)k1(xv, xv′)] = E[k1(xu, xu′)]E[k1(xv, xv′)] = 2.
  • If u = v and |{u, u′, v′}| = 3, then we have E[k1(xu, xu′)k1(xv, xv′)] = Exu [Exu′ [k1(xu, xu′)|xu]Ex′v [k1(xu, xv′)|xu]] = 2,
  by using Exu′ [k1(xu, xu′)|xu] = Ez,z′∼τ [k1(z, z′)], which holds by rotational invariance.
  • If u = v′ and u′ = v, we have E[k1(xu, xu′)k1(xv, xv′)] = Ez,z′∼τ [k1(z, z′)2] =: ̃. This takes the same form as (21), but with a different kernel function κ2 instead of κ. Note that in the case of the Exponential kernel, κ2 is also an exponential kernel with different bandwidth.
  Overall, when and ̃ are small compared to 1, we can see that the quantity E[k1(xu, xu′)k1(xv, xv′)] is small compared to 1 unless u = u′ and v = v′, thus satisfying the assumptions in Prop. 4. As described above, we may obtain expressions of and ̃ in various cases, and in particular these vanish in high dimension when using a kernel with κ(0) = 0.
sections.heading: F.2 FAST RATES
sections.text: In this section, we derive spectral compositions of 1-layer CKN architectures with non-overlapping patches under the product of spheres distribution described above. This allows us to derive fast rates that depend on the complexity of the target functions on patches, and shows similar improvement factors to those derived in Section F.1, without the κ(0) term, which in fact turns out to only be due to a single eigenspace, namely constant functions. We note that our derivation extends (Favero et al., 2021) to the case of generic pooling filters, and considers a different data distribution.
  Before studying the 1-layer case, we remark that while it may seem natural to extend such decompositions to the 2-layer case using tensor products of spherical harmonics, as done by Scetbon & Harchaoui (2020) in the case without pooling, it appears that pooling may make it more challenging to find an eigenbasis since subspaces consisting of tensor products of spherical harmonics with fixed total degree are no longer left stable by the kernel.5 We thus leave such a study to future work.
  We begin by considering the following Mercer decomposition of the patch kernel
  k1(z, z ′) = κ(〈z, z′〉) = ∞∑ k=0 µk N(d,k)∑ j=1 Yk,j(z)Yk,j(z ′),
  where Yk,j for k ≥ 0 and j = 1, . . . , N(d, k) are spherical harmonic polynomials of degree k forming an orthonormal basis of L2(dτ). The µk here are Legendre/Gegenbauer coefficients of the function κ (see, e.g., Bach, 2017a; Smola et al., 2001).
  Note that the 1-layer kernel may be written as K(x, y) = ∑ u,v∈Ω h ∗ h̄[u− v]κ(〈xu, yv〉),
  where ∗ denotes convolution and we denote h̄[u] = h[−u]. Denoting by ew[u], k = 0, . . . , |Ω| − 1 the orthonormal DFT basis vectors, normalized such that 〈ew, ēw〉 = 1, we have
  h ∗ h̄[u− v] = |Ω|−1∑ w=0 ĥ ∗ h̄[w]ew[u]ēw[v] = |Ω|−1∑ w=0 λwew[u]ēw[v],
  with λw := ĥ ∗ h̄[w] = |ĥ[w]|2. Then, we have
  K(x, y) = |Ω|−1∑ w=0 ∑ k≥0 λwµk N(d,k)∑ j=1 (∑ u ew[u]Yk,j(xu) )(∑ u ēw[u]Yk,j(yu) ) .
  Note that when k = 0, we have Y0,1(xu) = 1 for all u, hence(∑ u ew[u]Y0,1(xu) )(∑ u ēw[u]Y0,1(yu) ) = ( ∑ u ew[u])( ∑ v ēw[v]) = |Ω|1{w = 0},
  since e0 = |Ω|−1/2(1, . . . , 1) and ∑ u ew[u] = 0 for w > 0. Then, we may write
  K(x, y) = |Ω|λ0µ0φ0(x)φ̄0(y) + |Ω|−1∑ w=0 ∑ k≥1 λwµk N(d,k)∑ j=1 φw,k,j(x)φ̄w,k,j(y),
  5For instance, the term k1(xw, yu)k1(xw, yv), which may only appear in the presence of pooling, maps the polynomial Yk(xu)Yk(xv) of degree 2k to a polynomial µ2kYk(xw)
  2 which is not necessarily orthogonal to all spherical harmonics tensor products of degree smaller than 2k.
  with φ0(x) = 1, φw,k,j(x) = ∑ u ew[u]Yk,j(xu), and φ̄ denotes the complex conjugate of φ.
  It is then easy to check that the φ0 and φw,k,j form an orthonormal basis of L2(dτ⊗|Ω|). We thus have obtained a Mercer decomposition of the kernel K w.r.t. the data distribution, so that its eigenvalues are also the eigenvalues of the covariance operator (Caponnetto & De Vito, 2007), and control generalization performance, typically the degrees of freedom
  N (λ) = Tr((Σ + λI)−1Σ) = ∑ m≥0 ξm λ+ ξm ,
  where Σ is the covariance operator, and (ξm)m is the collection of its eigenvalues. In particular, if we have a decay ξm m−α (with α < 1), then we have N (λ) ≤ O(λ−1/α), which then leads to a fast rate of n−α/(α+1) on the excess risk when optimizing for λ in kernel ridge regression (Bach, 2021; Caponnetto & De Vito, 2007). In our case, we have
  N (λ) = |Ω|λ0µ0 λ+ |Ω|λ0µ0 + |Ω|−1∑ w=0 ∑ k≥1 N(d, k) λwµk λ+ λwµk
  We make a few remarks:
  • Given that the number of spatial frequencies w is fixed, the asymptotic decay rate of eigenvalues associated to the φw,k,j (that is, λwµk, each with multiplicity N(d, k)) is the same as that of the eigenvalues associated to φw0,k,j for some fixed w0, which in turn corresponds to the decay for the corresponding dot-product kernel on the sphere. For instance, if κ is the arc-cosine kernel, we have ξm m−α with α = d+2d−1 , and more generally α = 2sd−1 for a kernel resembling a Sobolev space with s bounded derivatives. This then leads to a rate n2s/(2s+(d−1)), which only depends on the dimension d of patches, rather than the full dimension d|Ω|. One may also add more general assumption on the smoothness of the localized components of f∗ (such as g in Proposition 3) in order to get rates that depend explicitly on the order of smoothness s of such components (as in Caponnetto & De Vito, 2007).
  • The eigenvalue associated to φ0 plays a minor role as by itself as it only contributes at most τ2ρ/n to the excess risk, which is negligible compared to the rest of the eigenvalues which lead to a slower n−α/(α+1) rate.
  • With no pooling (h is a Dirac delta), we have λw = |ĥ[w]|2 = 1/|Ω| for all w. We then have N (λ) ≤ 1 + |Ω|Nκ(|Ω|λ),
  where we defined Nκ(λ) := ∑ k≥1N(d, k)µk/(λ+ µk).
  • With global pooling (h = 1/|Ω| is constant), we have λ0 = 1/|Ω|, and λw = 0 for w > 0. This yields
  N (λ) ≤ 1 +Nκ(|Ω|λ). This then yields an improvement by a factor |Ω| in sample complexity guarantees compared to the scenario above with no pooling, namely the dominant term in the excess risk bound will be C(1/n)α/(α+1) compared to C(|Ω|/n)α/(α+1).
  • For more general pooling, one may exploit specific decays of λw to obtain finer bounds. We may also obtain the following bound by Jensen’s inequality
  N (λ) ≤ 1 + |Ω| ∑ k≥1 N(d, k) λ̄µk λ+ λ̄µk
  ≤ 1 + |Ω|Nκ ( |Ω|λ ‖h‖22 ) where we used that λ̄ = ( ∑ w λw)/|Ω| = ‖h‖22/|Ω|. When Nκ(λ) ≤ Cλ−1/α, we get a bound N (λ) ≤ 1 + C|Ω|1−1/α‖h‖2/α2 λ−1/α instead of a bound N (λ) ≤ 1 + C|Ω|1−1/αλ−1/α for the case of no pooling, that is, the improvement is again controlled by ‖h‖22, which goes from 1/|Ω| for global pooling, to 1 for no pooling.