summaries: This paper considers a simplified model of the polynomial neural network, and shows that the eigenvalues of the induced tangent kernel have a slower decay rate, compared with the classical two-layer network. This also means the polynomial neural networks can learn the high-frequency components faster. The empirical evidence supports this tendency. The paper studies ReLU-activated two-layer neural networks. The main contribution is that unlike traditional deep neural networks (DNNs), this type of polynomial neural network (PNN) also learns high-frequency information quickly. This claim partially verifies why PNNs are suitable for tasks like face recognition, where high-frequency components are crucial for performance.  The paper is dedicated to the study of spectral bias in polynomial neural networks. For feed-forward neural networks with ReLU activation function it is well-known that learning high-frequency terms is made slower (Rahaman et al.).   Authors claim that polynomial networks are able to learn both low frequency and high-frequency terms almost equally fast. Theory, based on spectral analysis of an integral operator associated with the NTK kernel, is provided.  Theoretical approach was taken from arXiv: 1905.12173, but for NTK kernel of the network with an additional multiplicative interaction layer.  Then, experiments that support claims are provided.  It is claimed that multiplicative interactions in the architecture of a NN may be the reason for an ability to learn high frequencies.   The paper gives an analysis of polynomial networks in the NTK regime. They show a theoretical speed-up in learning higher frequencies when using polynomial networks. They also have some experiments that verify these properties. 
abstractText: Polynomial neural networks (PNNs) have been recently shown to be particularly effective at image generation and face recognition, where high-frequency information is critical. Previous studies have revealed that neural networks demonstrate a spectral bias towards low-frequency functions, which yields faster learning of low-frequency components during training. Inspired by such studies, we conduct a spectral analysis of the Neural Tangent Kernel (NTK) of PNNs. We find that the Π-Net family, i.e., a recently proposed parametrization of PNNs, speeds up the learning of the higher frequencies. We verify the theoretical bias through extensive experiments. We expect our analysis to provide novel insights into designing architectures and learning frameworks by incorporating multiplicative interactions via polynomials.
references.authors: Zeyuan Allen-Zhu Yuanzhi Li Zhao Song
references.title: A convergence theory for deep learning via over-parameterization
references.venue: CoRR, abs/1811.03962,
references.authors: Sanjeev Arora Simon Du Wei Hu Zhiyuan Li Ruosong Wang
references.title: Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks
references.venue: In International Conference on Machine Learning,
references.authors: Devansh Arpit Stanisław Jastrzebski Nicolas Ballas David Krueger Emmanuel Bengio Maxinder S. Kanwal Tegan Maharaj Asja Fischer Aaron Courville Yoshua Bengio Simon Lacoste-Julien
references.title: A closer look at memorization in deep networks
references.venue: URL http://arxiv.org/abs/1706.05394
references.authors: Francesca Babiloni Ioannis Marras Filippos Kokkinos Jiankang Deng Grigorios Chrysos Stefanos Zafeiriou
references.title: Poly-nl: Linear complexity non-local layers with polynomials
references.venue: In International Conference on Computer Vision (ICCV),
references.authors: Francis Bach
references.title: Breaking the Curse of Dimensionality with Convex Neural Networks. arXiv:1412.8690 [cs, math, stat], October 2016
references.venue: URL http://arxiv.org/abs/1412
references.authors: Francis R. Bach
references.title: Sharp analysis of low-rank kernel matrix approximations
references.venue: CoRR, abs/1208.2015,
references.authors: Dzmitry Bahdanau Kyunghyun Cho Yoshua Bengio
references.title: Neural Machine Translation by Jointly Learning to Align and Translate
references.venue: URL http: //arxiv.org/abs/1409.0473
references.authors: Ronen Basri David W. Jacobs Yoni Kasten Shira Kritchman
references.title: The convergence rate of neural networks for learned functions of different frequencies
references.venue: CoRR, abs/1906.00425,
references.authors: Ronen Basri Meirav Galun Amnon Geifman David W. Jacobs Yoni Kasten Shira Kritchman
references.title: Frequency bias in neural networks for input of non-uniform density
references.venue: CoRR, abs/2003.04560,
references.authors: Alberto Bietti Julien Mairal
references.title: On the Inductive Bias of Neural Tangent Kernels. arXiv:1905.12173 [cs, stat], October 2019
references.venue: URL http://arxiv.org/abs/1905
references.authors: Tom B. Brown Benjamin Mann Nick Ryder Melanie Subbiah Jared Kaplan Prafulla Dhariwal Arvind Neelakantan Pranav Shyam Girish Sastry Amanda Askell Sandhini Agarwal Ariel Herbert-Voss Gretchen Krueger Tom Henighan Rewon Child Aditya Ramesh Daniel M. Ziegler Jeffrey Wu Clemens Winter Christopher Hesse Mark Chen Eric Sigler Mateusz Litwin Scott Gray Benjamin Chess Jack Clark Christopher Berner Sam McCandlish Alec Radford Ilya Sutskever Dario Amodei
references.title: Language Models are Few-Shot Learners
references.venue: [cs],
references.authors: Yuan Cao Zhiying Fang Yue Wu Ding-Xuan Zhou Quanquan Gu
references.title: Towards Understanding the Spectral Bias of Deep Learning. arXiv:1912.01198 [cs, stat
references.venue: URL http://arxiv.org/abs/1912.01198
references.authors: L. Carlitz
references.title: The product of two ultraspherical polynomials
references.venue: Glasgow Mathematical Journal,
references.authors: Lenaic Chizat Edouard Oyallon Francis Bach
references.title: On Lazy Training in Differentiable Programming
references.venue: [cs, math],
references.authors: Youngmin Cho Lawrence Saul
references.title: Kernel Methods for Deep Learning
references.venue: In Advances in Neural Information Processing Systems,
references.authors: Grigorios G. Chrysos Stylianos Moschoglou Giorgos Bouritsas Yannis Panagakis Jiankang Deng Stefanos Zafeiriou
references.title: P-nets: Deep Polynomial Neural Networks
references.venue: URL https://openaccess.thecvf.com/content_CVPR_2020/ html/Chrysos_P-nets_Deep_Polynomial_Neural_Networks_CVPR_2020_ paper.html
references.authors: Simon Du Jason Lee
references.title: On the power of over-parametrization in neural networks with quadratic activation
references.venue: In International Conference on Machine Learning (ICML),
references.authors: Simon S. Du Xiyu Zhai Barnabas Poczos Aarti Singh
references.title: Gradient Descent Provably Optimizes Over-parameterized Neural Networks. arXiv:1810.02054 [cs, math, stat], February 2019
references.venue: URL http://arxiv.org/abs/1810.02054
references.authors: Christopher Frye Costas J Efthimiou
references.title: Spherical harmonics in p dimensions
references.venue: arXiv preprint arXiv:1205.3548,
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,
references.authors: Yanping Huang Youlong Cheng Ankur Bapna Orhan Firat Mia Xu Chen Dehao Chen HyoukJoong Lee Jiquan Ngiam Quoc V. Le Yonghui Wu Zhifeng Chen
references.title: GPipe: Efficient Training of Giant Neural Networks using Pipeline Parallelism
references.venue: [cs],
references.authors: A.G. Ivakhnenko
references.title: Polynomial Theory of Complex Systems
references.venue: IEEE Transactions on Systems, Man, and Cybernetics,
references.authors: Arthur Jacot Franck Gabriel Clément Hongler
references.title: Neural Tangent Kernel: Convergence and Generalization in Neural Networks
references.venue: URL http://arxiv.org/abs/1806.07572
references.authors: Siddhant M. Jayakumar Wojciech M. Czarnecki Jacob Menick Jonathan Schwarz Jack Rae Simon Osindero Yee Whye Teh Tim Harley Razvan Pascanu
references.title: Multiplicative interactions and where to find them
references.venue: In International Conference on Learning Representations (ICLR),
references.authors: Ziwei Ji Matus Telgarsky
references.title: The implicit bias of gradient descent on nonseparable data
references.venue: In Proceedings of the Thirty-Second Conference on Learning Theory,
references.authors: Ziwei Ji Justin D. Li Matus Telgarsky
references.title: Early-stopped neural networks are consistent. 2021
references.authors: Tero Karras Samuli Laine Timo Aila
references.title: A style-based generator architecture for generative adversarial networks
references.venue: In Conference on Computer Vision and Pattern Recognition
references.authors: Tero Karras Samuli Laine Timo Aila
references.title: A Style-Based Generator Architecture for Generative Adversarial Networks. arXiv:1812.04948 [cs, stat], March 2019
references.venue: URL http: //arxiv.org/abs/1812.04948
references.authors: Joe Kileel Matthew Trager Joan Bruna
references.title: On the expressive power of deep polynomial neural networks. Advances in neural information processing systems (NeurIPS)
references.authors: Yann LeCun Yoshua Bengio
references.title: Convolutional networks for images, speech, and time series. In The handbook of brain theory and neural networks, pp. 255–258
references.authors: Jaehoon Lee Lechao Xiao Samuel Schoenholz Yasaman Bahri Roman Novak Jascha Sohl- Dickstein Jeffrey Pennington
references.title: Wide neural networks of any depth evolve as linear models under gradient descent
references.venue: Advances in neural information processing systems,
references.authors: Victor Lempitsky Andrea Vedaldi Dmitry Ulyanov
references.title: Deep Image Prior
references.venue: IEEE/CVF Conference on Computer Vision and Pattern Recognition,
references.authors: Chien-Kuo Li
references.title: A Sigma-Pi-Sigma Neural Network (SPSNN)
references.venue: Neural Processing Letters,
references.authors: Stefano Sarao Mannelli Eric Vanden-Eijnden Lenka Zdeborová
references.title: Optimization and generalization of shallow neural networks with quadratic activation functions
references.venue: arXiv preprint arXiv:2006.15459,
references.authors: Atsushi Nitanda Taiji Suzuki
references.title: Optimal rates for averaged stochastic gradient descent under neural tangent kernel regime
references.venue: In International Conference on Learning Representations (ICLR),
references.authors: Nasim Rahaman Aristide Baratin Devansh Arpit Felix Draxler Min Lin Fred A. Hamprecht Yoshua Bengio Aaron Courville
references.title: On the Spectral Bias of Neural Networks
references.venue: URL http://arxiv.org/abs/1806.08734
references.authors: Nicolas Le Roux Yoshua Bengio
references.title: Continuous neural networks
references.venue: Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics,
references.authors: Meyer Scetbon Zaid Harchaoui
references.title: A spectral analysis of dot-product kernels, 2021
references.authors: Bernhard Schölkopf Alexander J Smola Francis Bach
references.title: Learning with kernels: support vector machines, regularization, optimization, and beyond
references.venue: MIT press,
references.authors: Y. Shin J. Ghosh
references.title: The pi-sigma network: an efficient higher-order neural network for pattern classification and function approximation
references.venue: In IJCNN-91-Seattle International Joint Conference on Neural Networks, volume i, pp. 13–18
references.authors: Alex Smola Zoltán Óvári Robert C Williamson
references.title: Regularization with dot-product kernels
references.venue: Advances in Neural Information Processing Systems,
references.authors: Daniel Soudry Elad Hoffer Mor Shpigel Nacson Suriya Gunasekar Nathan Srebro
references.title: The Implicit Bias of Gradient Descent on Separable Data. arXiv:1710.10345 [cs, stat], December 2018
references.venue: URL http://arxiv.org/abs/1710.10345
references.authors: Matthew Tancik Pratul Srinivasan Ben Mildenhall Sara Fridovich-Keil Nithin Raghavan Utkarsh Singhal Ravi Ramamoorthi Jonathan Barron Ren Ng
references.title: Fourier features let networks learn high frequency functions in low dimensional domains
references.venue: Advances in Neural Information Processing Systems,
references.authors: Guillermo Valle-Perez Chico Q. Camargo Ard A. Louis
references.title: Deep learning generalizes because the parameter-function map is biased towards simple functions
references.venue: In International Conference on Learning Representations,
references.authors: Yan Wang Lingxi Xie Chenxi Liu Siyuan Qiao Ya Zhang Wenjun Zhang Qi Tian Alan Yuille
references.title: Sort: Second-order response transform for visual recognition
references.venue: In International Conference on Computer Vision (ICCV),
references.authors: Zhi-Qin John Xu Yaoyu Zhang Tao Luo Yanyang Xiao Zheng Ma
references.title: Frequency principle: Fourier analysis sheds light on deep neural networks
references.venue: CoRR, abs/1901.06523,
references.authors: Greg Yang Hadi Salman
references.title: A fine-grained spectral perspective on neural networks
references.venue: arXiv preprint arXiv:1907.10599,
references.authors: Chiyuan Zhang Samy Bengio Moritz Hardt Benjamin Recht Oriol Vinyals
references.title: Understanding deep learning requires rethinking generalization. arXiv:1611.03530 [cs], February 2017
references.venue: URL http://arxiv.org/abs/1611.03530. arXiv:
sections.heading: 1 INTRODUCTION
sections.text: Deep neural networks (DNNs) have demonstrated remarkable success in different domains [22, 11]. DNNs can approximate complex functions or even datasets with randomized labels arbitrarily well [50, 3], which makes their ability to avoid over-fitting on real data surprising, since it seems to disagree with prior notions of model complexity. This has sparked the interest in investigating the notion of “implicit bias” in neural network training, which makes them favor low complexity solutions [43, 21, 26].
  The spectral analysis of deep networks offers one perspective on this implicit bias. Deep neural networks demonstrate a learning bias towards low frequency functions - i.e. functions that vary globally without local fluctuations are learned faster when training neural networks via gradient descent [37, 48]. The phenomenon, termed as the spectral bias of neural networks [37], has been explored from the perspective of the Neural Tangent Kernel (NTK) [24]. The eigenvalues of the obtained kernel influence important characteristics such as the approximation properties and rate of learning [10, 39, 36]. For standard two-layer ReLU networks within the NTK regime, the analysis supports the idea of a spectral bias by showing faster error convergence for information in lower frequencies [12].
  Recently, another class of models, Polynomial Neural Networks (PNNs), that express high order polynomial expansions, have demonstrated state-of-the-art performance in the challenging tasks of image generation [28] and face recognition [17]. Both tasks rely on fine-grained details, which correspond to high-frequency information. More generally, multiplicative interactions have demonstrated strong empirical performance on image-based applications [44, 47, 4]. Jayakumar et al. [25] highlight how multiplicative interactions, which construct second degree polynomials, can enlarge the hypothesis space and lead to faster learning for certain classes of functions.
  To understand this success of PNNs, we conduct a spectral analysis, inspired by the analysis of DNNs. We focus on one instance of polynomial networks called Π-Nets [16], where the output is a piece-wise polynomial function of the input obtained via multiplicative layers. The parameters of a Π-Net can be represented as high-order tensors, while polynomial expansions can offer increased representation power [17]. Our main contributions can be summarized as follows:
  1. We analyze two-layer polynomial networks in the NTK regime. By studying the spectral properties of the corresponding kernel, we prove a theoretical speed-up in learning higher frequencies over standard neural networks and validate the hypothesis in the approximate NTK regime, on the task of learning spherical harmonics.
  2. Beyond the NTK regime, we demonstrate this enhanced bias of Π-Nets towards higher frequencies in several experimental settings, beginning from synthetic learning tasks and then proceeding to state-of-art networks and inverse problems with 2D images.
  Aside from improving the understanding of polynomials in neural networks, our proposed analysis also sheds new light on the effect of multiplicative interactions in neural networks, prevalent in certain domains of machine learning including vision and natural language processing [29, 7].
sections.heading: 2 RELATED WORK
sections.text: Spectral Bias: Motivated by the empirical observations in [3, 37, 48, 46] that deep networks first learn “simple patterns”, several papers [49, 12, 8, 2] have conducted theoretical analyses to explain this bias towards lower frequencies. The work of [2], aiming to understand why random labels take longer to learn than natural labels, showed that alignment of the labels with the eigenvectors of the NTK Gram matrix determines the learning speed. Extending this result, Cao et al. [12] provide an explanation for the spectral bias by analyzing the decay rate of eigenvalues of the NTK when the input data is uniformly distributed on the sphere. Under the same assumption, Basri et al. [8] study training dynamics in the NTK setting for 2-layer ReLU networks with a fixed outer layer and an explicitely included linear bias term in the ReLU. [9] further extends this work to account for non-uniform data distributions. All these findings show that DNNs learn lower frequency functions faster which prompted Tancik et al. [45] to propose methods to mitigate this bias. Our paper aims to establish similar results for PNNs, to explain their good performance at learning higher frequencies.
  Polynomial neural networks (PNNs): The early papers that explore polynomials in the context of neural networks are mainly divided into two categories: 1) self-organizing networks with hard coded features [23], 2) Pi-Sigma networks [41, 34]. In both cases, the constructions did not scale well to higher dimensional inputs, and were not used for high-dimensional signals, such as images. More recent papers have used the Hadamard product to capture correlations between different branches of an architecture [4, 44, 47, 17]. Our goal is to analyze the properties of these polynomial neural networks that have shown sucess in practice.
  Polynomial activation functions (PAFs): It is important to note the distinction between PNNs and Polynomial activation functions. PAFs expand (element-wise) each feature to an rth degree, i.e., they assume a (deep) neural network where the element-wise activation functions are rth degree polynomials. This is substantially different from capturing higher-order correlations across input (or feature) elements like PNNs especially in the presence of non-linear activations. Theoretical work on over-parametrization [18], expressive power [30] and generalization of shallow nets [35] have emerged for PAFs. The aforementioned papers, however, do not conduct a spectral analysis and do not exhibit the benefits of PNNs for learning high-frequency information.
sections.heading: 3 ANALYSIS OF POLYNOMIAL NEURAL NETWORKS IN THE NTK REGIME
sections.text: In this section, we conduct a careful analysis of the kernel approximation of polynomial neural networks (PNNs) to gain insight on the effect of the multiplicative interactions. PNNs include multiplicative interactions and express high-degree polynomial expansions. The recent parametrization of the Π-Net family [16], which we summarize below, is used as a representative PNN. We review the tangent kernel approximation of neural networks and then we derive the tangent kernel of two-layer Π-Nets to study the spectral bias of Π-Nets.
  Notation: We denote by x¨, ¨y the standard inner-product on Rd. For two vectors x,y P Rd, x ˚ y denotes the element-wise or Hadamard product. We use ˆm to denote the mode-m vector product1.
  1The reader may refer to the appendix for more details on the mode-m product.
  We define the asymptotic notations Ωp¨q and Ω̃p¨q as follows: Let an and bn be two sequences. We write an “ Ωpbnq if lim infnÑ8 |an{bn| ą 0. We use Ω̃p¨q to hide the logarithmic factors in Ωp¨q.
sections.heading: 3.1 Π-NET FORMULATION
sections.text: A polynomial expansion of the input vector z P Rδ can be used to express the output x P Ro as an N th degree polynomial expansion as follows:
  x “ N ÿ
  n“1
  ˆ W rns n`1 ź
  j“2 ˆjz
  ˙
  ` β, (1)
  where β P Ro and W rns P Roˆ śn m“1 ˆmd (N
  n“1 are the learnable parameters, and ˆm is the modem vector product. Since the number of tensor parameters W rns grow exponentially with the degree of the polynomial, a coupled tensor decomposition with factor sharing can be used. The idea in Π-Nets is to propose such decompositions that can capture higher order correlations, and can be implemented in standard deep learning frameworks. One such decomposition uses the recursive formulation xn “ ´ ATrnsz ¯ ˚ ´ STrnsxn´1 ` B T rnsbrns ¯
  for n “ 1, . . . , N and expresses the output x as x “ CxN ` β. The term in the rightmost parenthesis is exactly the recursive form of a standard neural network (without activations). Therefore, Π-Nets augment standard neural network with multiplicative interactions via the Hadamard product. While Π-Nets can approximate the target function without activations, they achieve state-of-art performance with activation functions, wherein the output is a piece-wise polynomial. We include more details in the Appendix.
sections.heading: 3.2 THE NEURAL TANGENT KERNEL
sections.text: Consider the following two-layer ReLU neural network (without bias parameters) with width m that assumes the following form (defined as in [10]): fW pxq “ b 2 mW2σpW1xq, where W1 P Rmˆpd`1q, W2 P R1ˆm and we assume inputs txuni“1 follow some distribution τ on the unit sphere Sd P Rd`1; σp¨q denotes the element-wise ReLU operator. As the width of the network m goes to infinity, if the weights at initialization W p0q are independent and each follow the standard normal distribution, the inner product of the network gradient at initialization gives rise to a limiting kernel, namely the Neural Tangent Kernel (NTK) [24] κ defined as :
  κpx,x1q “ lim mÑ8 x∇W fW p0qpxq,∇W fW p0qpx1qy. (2)
  This kernel has been used to characterize the behavior of sufficiently wide networks fW during training. For instance, if the network is trained to minimize the `2 loss, then its training dynamics closely track those of kernel regression under κ. This holds in a particular training regime, referred to as “lazy training” [14], where the parameters hardly vary after initialization and the network can be approximated by its first-order Taylor expansion at initialization as:
  fW pxq « fW p0qpxq ` x∇W fW p0qpxq,W ´W p0qy. In practice however, the conditions of lazy training are often violated (notably, the requirement that the weights do not move) within the first few steps of gradient descent. Nevertheless, the NTK remains a useful theoretical tool for analyzing the neural network behavior as some of its predictions have been shown to hold in practice [12, 32].
  The NTK for the two-layer ReLU network fW pxq takes the following form [10, 14, 19] κpx,x1q “ 2xx,x1yκ1px,x1q ` 2κ2px,x1q, (3)
  where the kernels κ1 and κ2 are defined as follows: κ1px,x1q “ Ew„Np0,Iqrσ1pxw,xyqσ1pxw,x1yqs, κ2px,x1q “ Ew„Np0,Iqrσpxw,xyqσpxw,x1yqs,
  (4)
  where σp¨q denotes the ReLU operator and σ1p¨q denotes the indicator function 1r¨ ě 0s. The kernels κ1, κ2 have been explicitly computed for the ReLU activation in [38, 15, 10] where they were shown to be positive semi-definite dot-product kernels. The function value depends only on the value of the dot-product of the arguments or more formally, κ1px,x1q “ g1pxx,x1yq and κ2px,x1q “ g2pxx,x1yq for some functions g1, g2 : R Ñ R. These kernels will serve as building blocks for the tangent kernel of the Π-Net architecture studied in the next section.
sections.heading: 3.3 Π-KERNEL
sections.text: Following [8, 10, 12] that consider shallow two-layer ReLU networks, we derive the tangent kernel for a two-layer Π-Net by supplementing the standard network with a multiplicative interaction layer:
  fW pxq “ c 2
  m W3rσpW2xq ˚ σpW1xqs, (5)
  where W1,W2 P Rmˆpd`1q, W3 P R1ˆm, ˚ denotes the Hadamard product and we again assume inputs txuni“1 that follow some distribution τ on the unit sphere Sd Ă Rd`1. Remark 1 Note that formulation yields a piece-wise quadratic polynomial and that it is important to consider ReLU activations (or other non-linearities) since otherwise, the network yields a quadratic polynomial of the input and is no longer a universal approximator, even with infinite width.
  Theorem 1 Let κ1, κ2 as defined in Eq. (4). The limiting kernel for the two-layer Π-Net, called Π-kernel and denoted by κπpx,x1q, takes the following form:
  κπpx,x1q “ 2p2xx,x1yκ1px,x1q ` κ2px,x1qqκ2px,x1q. (6)
  The proof follows the standard NTK calculations and is presented in the Appendix. Importantly, the addition of a single multiplicative interaction induces a product of kernels form.
  Remark 2 The required width that ensures that the two-layer Π-Net stays close to initialization (the corresponding NTK is close to the limiting κπ) is slightly higher than standard networks, by a Ωp ?
  logmq factor. We provide a rough sketch of the proof in the Appendix (end of C). Remark 3 The advantage of using the 2-layer Π-Net formulation is that it allows us to directly contrast against prior results on standard 2-layer feed-forward network, by gauging the effect of the extra multiplicative layer. However, unlike feed-forward networks, the theory does not easily extend to polynomials of higher degree or depth, since even with a fixed degree and depth, the NTK varies with the placement of the multiplicative connections. In this work, we choose to focus primarily on the effect of the multiplicative layer and leave the extension to general polynomials for future work.
  Since the Π-kernel κπ can be expressed as a sum of products of continuous positive definite dotproduct kernels κ1 and κ2, it inherits their regularity properties and is itself a Mercer kernel, as a consequence of the Schur Product theorem.
sections.heading: 3.4 SPECTRAL ANALYSIS
sections.text: To properly characterize the approximation properties and the behavior during training of this newly derived kernel, we study its Mercer decomposition. Indeed, this approach was used in [10] to show that the 2-layer NTK (Eq. (3)) has better approximation properties than a fixed first layer network, and in [12] to study the spectral bias of feed forward ReLU networks.
  The Mercer decomposition of a kernel κ is derived from the eigenvalues and eigenfunction of the integral operator associated to the kernel (see [40] Theorem 2.10 and references therein). Taking X to be some compact set in Rd, recall that for any continuous kernel function κ : X ˆ X Ñ R and Borel measure τ on X , we can define an integral operator Lκ that “convolves” any square integrable function f P L2τ pX q with κ:
  Lκpfqpxq “ ż
  X
  κpx,yqfpyqdτpyq. (7)
  For a Mercer kernel, this linear operator admits countable, real, non-negative eigenvalues tµ1, µ2, . . . u and the associated eigenfunctions form an orthonormal basis of L2τ pX q. If the data is uniform on the sphere Sd, and κ is a dot-product kernel, then these eigenfunctions of Lκ are the spherical harmonics ([42], Lemma 4). Consequently, by applying Mercer’s Theorem, we can obtain the following decomposition:
  κpx,x1q “ 8 ÿ
  k“0 µk
  Npd,kq ÿ
  j“1 Yk,jpxqYk,jpx1q, (8)
  where Yk,j for j “ 1, 2...Npd, kq represent the spherical harmonics of degree k in d ` 1 variables (whose explicit formula is given in Appendix D.1 in [5]) , and Npd, kq :“ 2k`d´1k ` k`d´2 d´1 ˘ .
  From this, we can characterize the RKHS H associated to κ as follows:
  H “
  $
  &
  %
  f “ ÿ
  kě0,µk‰0
  Npd,kq ÿ
  j“1 ak,jYk,jp¨q subject to ‖f‖2H “
  ÿ
  kě0,µk‰0
  Npd,kq ÿ
  j“1
  a2k,j µk ă 8
  ,
  .
  -
  .
  The benefit of determining the decay rate of the eigenvalues tµ1, µ2, . . . u can be easily deduced from the previous set as the eigenvalues determine the size of H. Indeed, the slower the decay of tµ1, µ2, . . . u, the more sequences pak,jqk,j will verify ř kě0,µk‰0 řNpd,kq j“1 a2k,j µk
  ă 8. The results of [12] rely on this decay rate to show that gradient descent on wide networks picks up low-frequency information first2 (Theorem 4.2, [12]). This decay rate is therefore of central importance. Next, we briefly refer to the prior results on characterizing this decay rate for the kernel obtained for the twolayer feed forward ReLU network (Eq. (3)) and we derive the equivalent results for the Π-kernel.
  Proposition 1 (Theorem 4.3 in Cao et al. [12], Proposition 5 in Bietti & Mairal [10]) For the tangent kernel corresponding to two-layer feed forward ReLU network, the eigenvalues pµkqk satisfy the following:
  $
  &
  % µ0, µ1 “ Ωp1q, µk “ 0,when k is odd, µk “ Ωpk´d´1q when k " d.
  (9)
  The decay is exponentially fast in the input dimension d. For the Π-kernel, we show the following improvement. Theorem 2 Let tµπ,1, µπ,2, . . . u denote the eigenvalues of the linear operator Lκπ associated to the kernel κπ . For k " d pk ‰ 2 mod 4q, it holds that µπ,k “ Ωpk´d{2´2q. The proof of the theorem is provided in the Appendix. The main idea is to plug in the Mercer decomposition of each kernel, leading to a product of polynomials form. The expansion of this product then allows for isolating the dominant terms contributing to the eigenvalues for each frequency. The key take-away is the much slower rate of decay the kth eigenvalue (an order almost Ωpkpd{2qq) when compared to the standard two-layer NTK. An immediate consequence is a “larger” RKHS, which lends the Π-kernel superior approximation properties in the higher frequencies. Further, when combined with the findings of [12], it yields a speed-up in learning higher frequency harmonics.
  Remark 4 On the point of a “larger” RKHS, we note the two perspectives at play. First is the idea that a slower decay implies more functions contained in the RKHS, leading to better approximation properties. However, from the classical statistical learning point of view, while a larger RKHS makes the optimization problem easier, it may lead to a sub-optimal prediction performance [6], eventually relating to the traditional trade-off in machine learning. We note through our experiments that it is the perspective relating to better approximation and optimization for the Π-kernel that dominates.
  Remark 5 The k " d setting may not reflect the case for image-based applications. For instance, generation tasks wherein the output is an image of resolution dˆd, the higher frequencies of interest roughly correspond to the order of Ωpdq. Consequently, the exponential improvement in the decay rate for the k " d setting, and by extension in the approximation properties of the RKHS as well as the speed of learning higher frequency information, may not be as dramatic elsewhere. Nevertheless, our analysis provides sufficient intuition to expect some degree of improvement in realistic settings.
sections.heading: 4 NUMERICAL EVIDENCE
sections.text: The analysis in Section 3 reveals that polynomial networks in the NTK regime learn higher frequency information faster. In practice however, neural networks deviate from the near-initialization NTK conditions within just a few iterations of gradient descent. Therefore, to verify the analysis on the spectral bias of Π-Nets, we conduct a series of experiments that increasingly deviate from the NTK regime, including image-based datasets to further verify our theoretical analysis. We initially consider synthetic data (Section 4.1), then move onto realistic settings. For all experiments, we use Π-Nets based on the product of polynomials formulation (Appendix), in the same vein as [16].
  2We offer a brief characterization of the result in the Appendix.
sections.heading: 4.1 EXPERIMENTS ON SYNTHETIC DATA
sections.text: Our first experiment relates to learning spherical harmonics in the approximately infinite width NTK regime [12]. More precisely, the task comprises of learning linear combinations of spherical harmonics with data uniformly distributed on the unit sphere. We consider large-width two-layer networks to simulate the infinite-width settings and we compare the performance of Π-Nets against standard neural networks. The optimization method is full batch vanilla gradient descent, to avoid stochastic effects. The observations align with our theoretical findings, we defer the results to Appendix due to space constraints. Next, we assess whether the spectral bias manifests beyond the NTK regime. We consider the task of learning sinusoidal signals with smaller width networks and larger learning rates, so as to expressly violate the NTK assumptions.
  Learning Sinusoids The goal is to learn sinusoidal signals [37]. Given frequencies ki P K, with amplitudes Ai P α and phases φi P Φ, the target map f˚ : r0, 1s Ñ R is defined as f˚pxq “ Ai sinp2πkix`φiq. We approximate this map with two methods: i) a fully connected neural network and ii) a Π-Net. In the first setting, we compare a 256-unit wide, six-layer deep neural network against a 256-unit wide, six-layer deep Π-Net (which has five multiplicative layers). In each case, the network fθ (parametrized by θ) regresses over f˚ (K = p5, 10, ..., 45, 50q, φi „ Up0, 2πq and Ai “ 1 @ i), using N “ 200 evenly spaced input samples over r0, 1s and with a fixed learning rate (same for both networks); the spectrum of the network f̃θpkq is monitored during training by tracking the magnitude of learned components for each frequency, averaged over 5 runs (Fig. 1).
  Remark 6 The experimental setup for training the networks replicates the setup of Rahaman et al. [37], without any special initialization or hyper-parameter tuning for either network and the Π-Net is implemented exactly as the product of polynomials formulation specified in the Appendix.
  We observe that Π-Nets do speed up training of higher frequencies. To better substantiate our claim, we consider two additional variants, the results of which are included in the Appendix. In the first, we compare a deeper nine-layer feedforward network with a six-layer Π-Net with only three multiplicative layers, such that the feedforward network includes more parameters (Fig. 10). Since skip connections (additive) are known to speed-up training, we next compare the Π-Net with a feedforward network of same depth, but with additive skip connections instead of multiplicative (Fig. 11). In both cases, we verify that the speed-up in learning higher frequencies due to multiplicative layers is much more significant.
  Discussion: In the task of learning sinusoids, multiplicative interactions can speed up the learning of higher frequencies and are more effective at doing so than simply increasing the depth. In the Appendix, we conduct an experiment to evaluate the robustness of the networks in retaining high frequency information (Fig. 12) and we see that Π-Nets are more robust to perturbations in the higher frequencies. Remarkably, the Π-Net with more interactions is noticeably more robust than the lower degree polynomial network (Fig. 13 in the appendix). We hypothesize that Π-Nets enhance the
  representational space for higher frequencies, relating to our observations on the RKHS in the NTK regime, which makes them less susceptible than standard neural networks. This could also explain why more multiplicative interactions lead to more robustness.
sections.heading: 4.2 LEARNING IMAGES
sections.text: To further validate our claim in natural images, we adopt the convolutional layers often used in deep convolutional neural networks (DCNNs) [31]. DCNNs are ubiquitous in vision, partly due to the DCNN structure imposing a suitable prior for tasks with natural images, termed as the Deep Image Prior (DIP) [33]. We precisely assess how this prior changes with multiplicative interactions in Π-Nets in a denoising setup adapted from the DIP framework.
  Experiment - Denoising: We consider an image x P R3ˆHˆW and obtain its noisy version x0 “ x ` δ, perturbed with Gaussian noise. For the input, we sample a random tensor z P RNˆHˆW (N = 32 in our setup). We consider a neural network fθ (making the parametrization by θ explicit), and optimize the reconstruction loss, minθ ||fθpzq ´ x0||2, with respect to the noisy image. We can expect the DCNN structure to first learn the (“natural”) features corresponding to the true image x and pick up the noise only in the latter stages, which can then be avoided using early-stopping [33].
  Remark 7 Our goal is not to quantify the denoising performance but rather, to validate the experimental observation on sinusoids and assess whether Π-Nets can speed up learning of high frequencies in real-world applications. If the speed up is confirmed, Π-Nets can be early-stopped, i.e., require less iterations for achieving the target result.
  As in [37], we experiment with the U-net architecture, adapting it suitably for Π-Nets with an image of resolution 480 x 600. The implementation details are included in the Appendix. We train both networks for 2500 iterations, with the same learning rate and identical input tensor z. We monitor the Peak Signal-to-Noise Ratio (PSNR) curves while training, in Fig. 2. The peak-PSNR with the true image for both standard and Π-Nets are roughly the same, indicating similar denoising performance. However, Π-net achieves this peak PSNR in approximately half the number of iterations. Beyond this point, the Π-net PSNR with the true image decreases as it starts to pick up the high-frequency noise. We confirm this in the visual snapshot of the output of the two networks (Fig. 16 in the appendix). We conclude that Π-Nets do indeed speed up learning in images, by showing a reduced impedance towards high-frequency information. Next, we check how this bias affects learning high frequency information pertaining to natural images in the absence of noise.
  Experiment - Power Spectrum Analysis: We consider the identical setup as the denoising experiment but without the noise perturbation, setting δ “ 0. Therefore, the task for the network fθ learn θ˚ such that fθ˚pzq “ x. We allow 600 iterations of gradient descent and we monitor the radial power spectral density (p.s.d.) of the network output image against the ground truth p.s.d. It allows us to explicitly track the learning progress for each the frequency magnitude.
  We use the same U-net structure and suitably modify it for Π-Nets as before. Since the multiplicative interactions introduce more parameters, we reduce the depth/scale (by scale, we mean number of down/up-sampling operations in the U-net) of Π-Net to ensure roughly the same parameters („1.3 million). The learning speeds in different frequencies can be observed in Fig. 3. In the Appendix, we repeat the same experiment for equal depth networks.
  Discussion: From the preceding two experiments, we establish how Π-Nets pick up high frequency information faster than (larger) DCNNs. In the Appendix, we note the effect of multiplicative interactions on the deep image prior of the network. We also verify that this altered spectral bias remains relevant for standard learning tasks with natural images. Consequently, for tasks where capturing the finer details of the image is important (such as the denoising instance above), this potentially leads to a reduced number of iterations for Π-Nets as compared to standard networks. We expect the insights from our work to be useful for guiding network design with multiplicative interactions in large-scale experiments.
sections.heading: 4.3 EFFECT OF LABEL NOISE IN CLASSIFICATION
sections.text: Finally, we conduct an experiment in a classification setting, wherein the goal is to quantify the effect of this spectral bias in presence of label noise. The setup is adapted from Rahaman et al. [37], a fully connected 6-layer deep 256-unit wide network is trained on a binary classification on MNIST images (by only considering classes “3” and “8”). The labels are perturbed by label noise of different frequencies and the network is trained on the noisy labels with mean squared loss. Rahaman et al. [37] noted for feedforward networks that for a fixed amplitude, low frequency label noise degrades generalization performance (difference between training and validation losses) to a larger extent than high frequency noise. While the low frequency noise is learned instantly, the high-frequency noise is only fit later in the training. As a result, the network learns only true labels at first, and this corresponds to a “dip” in the true validation loss in the early stages. This “dip” becomes larger
  with increasing frequency of noise, indicating network impedance towards higher frequencies in the early iterations. Thus, it is only during the latter stages that the true validation loss degrades.
  We repeat the experiment with Π-Nets, by supplementing the feedforward network with exactly one multiplicative layer, and contrast the two networks. Since Π-Nets pick up high frequency variations faster, we consequently expect a smaller “dip”. We train for 5000 iterations with identical learning rates and observe the validation loss curves for different label noise frequencies. In Fig. 4, we zoom in on the first 1000 iterations for better visualization (complete curves included in the Appendix).
  Discussion: While the performance for the two networks is identical in absence of label noise (freq=0), the validation “dip” in Π-Net for higher frequencies (0.3, 0.5, 1) is visibly smaller and is negligible for lower frequencies (0.1, 0.2), indicating that Π-Nets pick up the high frequency label noise in the decision boundaries much faster. Additionally, we verify that increasing the number of multiplicative layers reduces the “dip” even further. The plots are deferred to the Appendix.
sections.heading: 5 DISCUSSION
sections.text: In this work, we focus on the spectral of polynomial neural networks. Our theoretical results in the NTK regime utilize a two-layer polynomial network and demonstrate a speed-up in the learning of higher frequencies over standard neural networks. We experimentally verify these properties, even outside of the NTK training regime. Additionally, the results offer intuition behind the success of networks that use multiplicative interactions, such as StyleGAN [28], whose connections to polynomials have been noted previously [17]. As a future direction, we aim to design controlled settings to study how polynomial and multiplicative interactions affect performance in state-of-the-art conditional generative models for image generation or deblurring, where high-frequency information is critical for photo-realistic outcomes. It should be noted that the spectral bias of standard neural networks towards low frequency (complexity) functions is believed to help in generalization. Therefore, another important direction is to explore whether this enhanced spectral bias of polynomials towards higher complexity functions translates to differences in their generalization properties. Finally, our results in the classification task yield a further research direction towards analyzing the smoothness of decision boundaries of Π-Nets with multiplicative interactions, especially focusing on areas wherein this effect becomes relevant such as adversarial susceptibility or knowledge distillation.
sections.heading: ACKNOWLEDGMENTS
sections.text: We are thankful to the anonymous reviewers for their constructive feedback. Research was sponsored by the Army Research Office and was accomplished under Grant Number W911NF-19-10404. This project has received funding from the European Research Council (ERC) under the
  European Union’s Horizon 2020 research and innovation programme (grant agreement n° 725594 - time-data). JM was supported by the ERC grant number 714381 (SOLARIS project) and by ANR 3IA MIAI@Grenoble Alpes, (ANR19-P3IA-0003).
sections.heading: A A PRIMER ON POLYNOMIAL NEURAL NETWORKS (Π-NETS)
sections.text: Polynomial neural networks (Π-Nets) [16] were recently introduced with the output being a high degree polynomial expansion of the input. The parameters of Π-Nets can be represented as higher order tensors. These networks are able to serve as function approximators even without the use of non-linear activations. When used in conjunction with activations, the models demonstrate better expressivity and achieve state-of-art results on a variety of learning tasks. We next describe the construction of Π-Nets.
  A.1 METHOD
  Mode-m vector product: Consider an M th order tensor X , with each of its element addressed by M indices, i.e., pX qi1,i2,...,iM “ xi1,i2,...,iM . An M th-order real-valued tensor X is defined over the tensor space RI1 ˆ I2 ˆ ... ˆ IM , where Im P Z for m “ 1, 2, ....M . The mode-m unfolding of a tensor X P RI1ˆI2ˆ...ˆIM maps X to a matrix Xpmq P RImˆÎm with Îm “ śM i“1,i‰m Ii. such that the tensor element X i1,i2,...,iM is mapped to the matrix element Xim,j where j “ 1 ` řM k“1,k‰mpik ´ 1qJk, where Jk “ śM n“1,n‰m In. The mode-m vector product of X with a vector u P RIm , denoted by X ˆm u P RI1ˆI2ˆ.....ˆIm´1ˆIm`1....IM , results in a tensor of order M ´ 1:
  pX ˆm uqi1,...,im´1,im`1,...,iM “ Im ÿ
  im“1 xi1,i2,...,iMuim .
  Furthermore, we denote X ˆ1 u1 ˆ2 u2.....ˆM uM “ X śm i“1ˆmum.
  Π-Net learns a function G : Rd Ñ Ro, such that each element of the output xj can be expressed as a polynomial of all the input elements zi, with i P r1, ds as follows:
  xj “ Gpzqj “ βj `wr1sj T z ` zTW r2sj z`
  W r3sj ˆ1 z ˆ2 z ˆ3 z ` ¨ ¨ ¨ `W rNs j
  N ź n“1 ˆnz,
  (10)
  where βj P R and W rnsj P R śn m“1 ˆmd (N
  n“1 are parameters for approximating the output xj . The correlations (of the input elements zi) up to N th order emerge in Eq. (10). More compactly:
  x “ Gpzq “ N ÿ
  n“1
  ˆ W rns n`1 ź
  j“2 ˆjz
  ˙
  ` β, (11)
  where β P Ro and W rns P Roˆ śn m“1 ˆmd (N
  n“1 are the learnable parameters. This form allows us to approximate any smooth function as per an extension of the Weierstrass Theorem. To prevent an exponential number of parameters, the authors propose using coupled tensor decompositions.
  A.2 TENSOR DECOMPOSITION FOR SINGLE POLYNOMIAL
  An appropriate tensor decomposition on the parameters in Eq. (11) allows for implementation with a neural network. Here, we briefly describe one such decomposition:
  Model: NCP: Next, we consider a joint hierarchical decomposition on the polynomial parameters. A Nested coupled CP decomposition (NCP) results in the following recursive relationship for N th order approximation:
  xn “ ´ ATrnsz ¯ ˚ ´ STrnsxn´1 `B T rnsbrns ¯ , (12)
  for n “ 2, . . . , N with x1 “ ´ ATrnsz ¯ ˚ ´ BTrnsbrns ¯
  and x “ CxN ` β. The parameters C P Roˆk,Arns P Rdˆk,Srns P Rkˆk,Brns P Rωˆk, brns P Rω for n “ 1, . . . , N , are learnable.
  brns P Rω (N
  n“1act as a scaling factor for each parameter tensor, whose role is illustrated in case of the third order approximation in Eq. (13):
  Gpzq “ β`W r1sˆ2 br1sˆ3 z`W r2sˆ2 br2sˆ3 zˆ4 z`W r3sˆ2 br3sˆ3 zˆ4 zˆ5 z. (13)
  Further, the joint factorization of the parameter tensors (in matrix form) for the third order NCP polynomial may be summarized as follows:
  • First order parameters : W r1sp1q “ CpAr3s dBr3sq T .
  • Second order parametes: W r2sp1q “ C " Ar3s d „ ´ Ar2s dBr2s ¯ Sr3s
  *T
  .
  • Third order parameters: W r3sp1q “ C " Ar3s d „ˆ Ar2s d !´ Ar1s dBr1s ¯ Sr2s ) ˙ Sr3s
  *T
  A.3 PRODUCT OF POLYNOMIALS
  The second scheme of implementation approximates the target function using a product of polynomials form, wherein the output of first polynomial is fed to the next and so on. The concept is visually depicted in 6; for instance, if each polynomial is degree two, then stacking N such polynomials results in an overall order of 2N .
  Remark: In practice, each matrix operation in the recursive formulation of Π-Nets represents an affine transform on a vector. Therefore, the implementation in standard libraries may be as simple as using a convolutional or a fully-connected layer. The other difference arises due to the multiplicative layers, which may be seen as special ’skip’ connections which are combined with the network output via the element-wise vector product (Hadamard product), as opposed to addition. Therefore, it is easy to see why these networks scale as well as standard deep networks.
sections.heading: B EIGENVALUES DETERMINE THE SPEED OF LEARNING
sections.text: As noted previously, the integral operator with respect to kernel function is defined as follows:
  Lκpfqpxq “ ż
  X
  κpx,yqfpyqdτpyq. (14)
  The key idea in [12] is to establish guarantees on the speed of convergence of gradient descent along different directions of Lκ, as defined by its eigenfunctions. Consider tλiuiě1 with λ1 ě λ2 ě .. be the strictly positive eigenvalues of Lκ with tφiuiě1 the respective eigenfunctions and define vi “ n´1pφpx1q, φpx2q.......φpxnqq. Since the eigenvalues may not be distinct, define rk as the sum of multiplicities of the first k distinct eigenvalues of Lκ. Define Vrk “ pv1..., vrkq. By definition, viiďrk are rescaled restrictions of orthonormal functions on the training examples. Therefore, they form a set of almost orthonormal bases in the vector space Rn.
  Finally, denote y “ py1, y2.....ynqT and ŷptq “ pfW ptqpx1q, .....fW ptqpxnqqT as the ground truth and predictions at time t, respectively, where fW ptq denotes the 2-layer ReLU network. Then the following result holds:
  Theorem 4.2 [12]: Suppose |φjpxq| ď M for j P rrks and x P Sd`1. For any , δ ě 0 and integer k, if n ě Ωp ´2.maxtpλrk ´ λrk`1q´2,M4r2kuq,m ě
  ˜ΩppolypT, λ´1rk , ´1qq, then with probability at least 1´ δ, gradient descent with step size η “ Õpm´1q satisfies:
  n´1{2 ¨ ||V Trkpy ´ ŷ ptqq||2 ď 2 ¨ p1´ λrkqT .n´1{2 ¨ ||V Trky||2 ` , (15)
  which is to say that the convergence rate of ||V Trkpy ´ ŷ ptqq||2, or alternatively the projection of the residual error on the space spanned by the first rk eigenvalues is controlled by the rthk eigenvalue λrk . With some additional work, this result can be extended to the 2-layer Π-Netby accounting for the extra width factor of Ωp ? logmq, as noted previously. The key takeaway in that with a wide enough network and large enough sample size, gradient descent first learns the target function along the eigen-directions with larger eigenvalues. Since the decay of eigenvalues is slower for the 2-layer Π-Net, it leads to a speed-up in learning higher frequencies.
sections.heading: C PROOF OF THEOREM 1: DERIVING THE Π-KERNEL
sections.text: In this section, we prove our first main result for the Π-kernel corresponding to theorem 1. Consider a two-layer Π-Net fW parametrized as follows :
  fW pxq “ c 2
  m W3rσpW2xq ˚ σpW1xqs, (16)
  where the weightsW1 P Rmˆd,W2 P Rmˆd,W3 P R1ˆm are initialized with independent identically distributed N p0, 1q coordinates.
  Recall that the neural tangent kernel κπ corresponds to the limit of the following inner product
  κπpx,x1q “ lim mÑ8 x∇W fW p0qpxq,∇W fW p0qpx1qy. (17)
  We can compute the inner-product Eq. (17) by computing the derivatives with respect to W1, W2, W3 of fW separately and sum up the inner products to obtain κπ , since the gradient can be split into three blocks.
  We denote by α̃ :“ W1x, and by β̃ :“ W2x the pre-activation vectors; and by α, β the postactivation vectors where the element-wise activation σ is applied to α̃ and β̃ respectively.
  First consider the derivative w.r.tW1 denoted BW1 :
  BW1fW pxq “ c 2
  m
  m ÿ i“1 Wi3σpβ̃ipxqqσ1pα̃ipxqqBW1 α̃ipxq. (18)
  Since α̃i “ eTi W1x “ xW1, eixT y, where ei P Rm is the i-th canonical basis vector of Rm, we have that
  BW1 α̃ipxq “ eixT .
  The contribution to the NTK (Eq. (17)) corresponds to xBW1fWpxq, BW1fW px1qy. Note that since the network is symmetric in tW1,W2u, we need only look at W1 to obtain the contribution for W2 as well. We find that
  xBW1fWpxq, BW1fW px1qy “ 2
  m
  m ÿ
  i,j“1 Wi3W j 3rσpβ̃ipxqqσpβ̃jpx1qqsrσ1pα̃ipxqqσ1pα̃jpx1qqsxeixT , ejx1T y
  “ 2 m
  m ÿ
  i,j“1 Wi3W j 3rσpβ̃ipxqqσpβ̃jpx1qqsrσ1pα̃ipxqqσ1pα̃jpx1qqsxTx1 δij
  “ 2 m
  m ÿ i“1 Wi3W i 3rσpβ̃ipxqqσpβ̃ipx1qqsrσ1pα̃ipxqqσ1pα̃ipx1qqsrxTx1s.
  (19)
  As limmÑ8, the above quantity converges to its expectation by the law of large numbers because the terms are independent and identically distributed with finite expectation. Consequently,
  lim mÑ8
  xBW1fWpxq, BW1fW px1qy “ 2xx,x1yκ1px,x1qκ2px,x1q (20)
  since $
  &
  % EW1rσ1pα̃ipxqqσ1pα̃ipx1qqs “ κ1px,x1q EW2rσpβ̃ipxqqσpβ̃ipx1qqs “ κ2px,x1q EW3rWi3Wi3s “ 1.
  By symmetry, we obtain with a similar term for W2 and their total contribution to Eq. (17) adds up to twice the term obtained forW1. Characterizing the kernel w.r.tW3 is more straightforward:
  BW3fW pxq “ c 2
  m σpα̃pxqq ˚ σpβ̃pxqq. (21)
  Therefore, its contribution to the NTK is
  xBW3fW pxq, BW3fW px1qy “ 2
  m
  m ÿ
  i
  σpα̃ipxqqσpβ̃ipxqqσpα̃ipx1qqσpβ̃ipx1qq. (22)
  Applying the law of large numbers again, as limmÑ8, this quantity tends to:
  EW1,W2„Np0,IqtrσpxW1,xyqqσpxW1,x 1yqsrσpxW2,xyqqσpxW2,x1yqsu
  “ EW1rσpxW1,xyqqσpxW1,x1yqsEW2rσpxW2,xyqqσpxW2,x1yqs “ 2κ2px,x1qκ2px,x1q.
  (23)
  The theorem follows by summing up 2ˆ Eq. (20) and Eq. (23). Width Requirements As noted previously, two-layer Π-Netneeds an the extra factor Ωp ?
  logmq in terms of width, to stay close to initialization The proof is largely derived from Lemma A.6 in [27], based on ideas first noted in [1] and therefore, we just provide a rough sketch here. For simplicity, consider the set of weights (W2,W3) fixed at initialization, i.e. at W p0q 2 ,W p0q 3 and note the first-order Taylor expansion for the Π-Net w.r.tW1 as follows:
  fW pxq « f0W pxq “ fW p0qpxq ` x∇W p0q1 fW p0qpxq,W1 ´W p0q 1 y
  Here f0W denotes the first order Taylor expansion of f at W0. We can then expand the RHS as: c
  2
  m
  ÿ
  j
  W p0q 3,j pσpx TW p0q 2,j qσpx TW p0q 1,j q ` σpx TW p0q 2,j qσ 1pxTWp0q1,j qx T pW1,j ´Wp0q1,j qq.
  Consequently, we can bound the approximation error as:
  |fW pxq´f0W pxq| ď c 2
  m
  ÿ
  j
  |σpxTWp0q2,j q||W p0q 3,j |pσpx TW1,jq´pσpxTWp0q1,j q`σ 1pxTWp0q1,j qx T pW1,j´Wp0q1,j qq|.
  Notice that aside from the |σpxTWp0q2,j q|, the error term is exactly the same as that for a standard network. Also note that P pxTW2,j ě τq “ P pσpxTW2,jq ą τq ď e´τ
  2{2, for τ ą 0 by the sub-gaussian tail bound. To ensure |σpxTWp0q2,j q| ď τ @ j with probability 1 ´ δ, we use the union bound to obtain the condition on τ as m.e´τ 2{2 ă δ. Consequently, we have that τ „ Ωp ?
  logmq, which essentially becomes an additional constant factor in the error term, over the usual error terms from the standard feed-forward network.
sections.heading: D PROOF OF THEOREM 2: CHARACTERIZING THE Π-KERNEL EIGENDECAY
sections.text: Here, we prove the Π-kernel eigenvalue decay rate stated in theorem 2. We first recall some connections between spherical harmonics and Gegenbauer polynomials.
  Definition 1 For a given α P R, Gegenbauer (or ultraspherical) polynomials denoted Cαk : r´1, 1s Ñ R are a family of orthogonal polynomials with respect to the weight function x ÞÑ p1´ x2qα´ 12 , i.e,
  ż 1
  ´1 Cαk pxqCα` pxqp1´ x2qα´ 1 2 “ 0,
  for k ‰ `.
  Remark: Gegenbauer polynomials are a generalization of Legendre polynomials which can be recovered by taking α “ 12 . The following addition formula expresses Gegenbauer polynomials in terms of spherical harmonics.
  Lemma 1 (Theorem 4.11 in [20]) (Addition formula) Let tYk,juNpd,kqj“1 denote spherical harmonics of degree k in d` 1 variables. It holds that, for any x,x1 P Sd,
  Npd,kq ÿ
  j“1 Yk,jpxqYk,jpx1q “ Npd, kqC
  p d´12 q k pxx,x 1yq.
  By making use of this Lemma, we can rewrite the Mercer decomposition of any admissible dot product kernel K when the data is uniform on the unit sphere. Indeed, by simplifying the Mercer’s decomposition in terms of spherical harmonics we have :
  Kpx,x1q “ 8 ÿ
  k“0 µk
  Npd,kq ÿ
  j“1 Yk,jpxqYk,jpx1q
  “ 8 ÿ
  k“0 µkNpd, kqC
  p d´12 q k pxx,x 1yq,
  (24)
  where pµkq8k“0 are its corresponding eigenvalues, and the second equality follows from Lemma 1. Now, since κ1 and κ2 are dot-product Mercer kernels, we can, akin to Eq. (24), obtain their respective decompositions in terms of Gegenbauer polynomials. There exist pµ1,kq8k“0 and pµ2,kq8k“0, both sequences of eigenvalues, such that
  xx,x1yκ1px,x1q “ 8 ÿ
  k“0 µ1,kNpd, kqC
  p d´12 q k pxx,x 1yq
  κ2px,x1q “ 8 ÿ
  k“0 µ2,kNpd, kqC
  p d´12 q k pxx,x 1yq. (25)
  The decay rate of the eigenvalues pµ1,kq8k“0 and pµ2,kq8k“0 is known [5, 10, 12] and is stated in the following Lemma.
  Lemma 2 (Appendix D.2 of [5], Lemma 17 of [10]) For k " d, k even, we have that µ1,k „ Apdqk´d´1 and µ2,k „ Bpdqk´d´2´1{2, where Apdq and Bpdq are constants only depending on the dimension d.
  Our goal is to establish the decay rate of the eigenvalues pµπ,kq8k“0 of the kernel
  κπpx,x1q “ 2p2xx,x1yκ1px,x1q ` κ2px,x1qqκ2px,x1q. (26)
  By plugging in the representations in Eq. (25) into the Π-kernel expression (Eq. (26)), we find that
  κπpx,x1q “ 2p 8 ÿ
  k“0 p2µ1,k ` µ2,kqNpd, kqC
  p d´12 q k pxx,x 1yqqp 8 ÿ
  k“0 µ2,kNpd, kqC
  p d´12 q k pxx,x 1yqq.
  (27) Moreover, since we can also write
  κπpx,x1q “ 8 ÿ
  k“0 µπ,kNpd, kqC
  p d´12 q k pxx,x 1yq,
  an appropriate factorisation of the product Eq. (27) will allow us to deduce the order of µπ,k by equating the coefficients appearing in front of the Gegenbauer polynomials.
  Developing the product in Eq. (27) leads to the appearance of products of Gegenbauer polynomials. The product of two Gegenbauer polynomials turns out to be a linear combination of other Gegenbauer polynomials.
  Lemma 3 (Equation 8 in [13]) Let α P R, for any m,n P N, there exists positive coefficients pλpm,nqs qminpm,nqs“0 such that
  Cpαqm pxqCpαqn pxq “ minpm,nq ÿ
  s“0 λpm,nqs C pαq m`n´2spxq. (28)
  By using this expansion, coefficients for each Gegenbauer polynomial may be identified.
  Take k P N to be an even number. We know from plugging Lemma 3 into Eq. (27) that the coefficient µπ,2kNpd, 2kq in front of the polynomial C p d´12 q 2k can be lower bounded as follows. µπ,2kNpd, 2kq ě Npd, kq2p2µ1,k ` µ2,kqµ2,kλpk,kq0 . (29) We obtain this lower bound by considering the contribution of a single term in Eq. (28) where m “ n “ k and s “ 0 to the coefficient in front of Cp d´1 2 q
  2k . This contribution is a lower bound because all coefficients involved in the product are non-negative.
  Consequently, in order to establish a decay rate for pµπ,kq it suffices to study the decay rate of λpk,kq0 . This rate is given by the following Lemma.
  Lemma 4 Let α “ d´12 be an integer, the coefficient λ pk,kq 0 defined in Lemma 3 admits the following expression λpk,kq0 “ ppα`k´1q!q2 pα´1q!ppkq!q2 p2kq! pα`2k´1q! . Moreover, for k " d, it holds that
  λ pk,kq 0 „ kpd{2q.
  See Appendix D.1 for a proof.
  We now dispose of all the necessary results to prove Theorem 2. Starting from Eq. (29), we find that for k " d, we have that
  µπ,2k ě Npd, kq2
  Npd, 2kq p2µ1,k ` µ2,kqµ2,kλ pk,kq 0
  „ k dkd
  p2kqd p2µ1,k ` µ2,kqµ2,kλ pk,kq 0 (by Stirling)
  “ k d
  2d Ωpk´2d´2qkpd{2q (by Lemma 2 and Lemma 4)
  “ Ωpp2kq´d{2´2q.
  (30)
  This concludes the proof that, for k divisible by 4, we have µπ,k “ Ωpk´d{2´2q. For k “ 1 mod 4, we conduct the same reasoning taking the coefficient λp2k`1,2kq0 . For k “ 3 mod 4, the coefficient to consider is λp2k`1,2k`2q0 . The equivalence derivation Lemma 4 proceeds in exactly the same manner for these coefficients.
  D.1 PROOF OF LEMMA 4
  Let us denote pαqk :“ αpα` 1qpα` 2q....pα` k ´ 1q and pαq0 :“ 1. From Equation 8 in [13], we have that
  λ pk,kq 0 “ 2k ` α 2k ` α. pαq0pαqkpαqk 0!pkq!pkq! . p2αq2k pαq2k . p2kq! p2αq2k
  “ pαqkpαqkpkq!pkq! p2kq! pαq2k
  “ ppα` k ´ 1q!q 2 ppα´ 1q!pkq!q2 p2kq!pα´ 1q! pα` 2k ´ 1q!
  “ ppα` k ´ 1q!q 2 pα´ 1q!ppkq!q2 p2kq! pα` 2k ´ 1q! .
  (31)
  We can apply Stirling’s approximation stating that n! „ ?
  2πnpne q n to find that:
  λ pk,kq 0 „ pα` k ´ 1qpα` k ´ 1q2pα`k´1q? α´ 1pα´ 1qα´1kpkq2k .
  ? 2kp2kq2k?
  α` 2k ´ 1pα` 2k ´ 1qα`2k´1 (32)
  Considering the case when k " α or equivalently, k " d since α “ pd ´ 1q{2, we obtain the following simplification:
  λ pk,kq 0 „
  kp2v`2k´1q k2k`1 . p2kq2k`0.5 p2kqv`2k´0.5
  „ pk 2 qα´1 „ kpd{2q.
  (33)
sections.heading: E LEARNING SPHERICAL HARMONICS WITH Π-KERNEL
sections.text: Following the setup in [12], we perform a similar experiment on learning combinations of spherical harmonics in the NTK regime. We define and initialize the Π-Net exactly as specified, with a width of 32768 neurons and using vanilla gradient descent, to approximate the kernel learning in the infinite width limit. We take n “ 1000 samples, txuni“1 from the uniform distribution on the unit sphere S10. We define our target function with integral k P K as follows:
  f˚pxq “ 1 NpKq ÿ
  kPK AkPkpxx, ζkyq, (34)
  where the Pkptq is the Gegenbauer polynomial with degree k, ζk are fixed vectors that are independently generated from uniform distribution on unit sphere in R10, andNpKq is a suitable normalizing constant to keep the order of magnitude of the target function approximately the same for different choices of K. As noted previously, f˚ may be seen as a linear combination of spherical harmonics and we compare the error residuals during the learning process in standard vs Π-Nets in the NTK regime, for varying K. We use a moving average of range 20 on these curves to smoothen out the heavy oscillations in the latter stages.
  Figure 8: The plots represent a comparison of log-scale convergence curve of error projection lengths for standard vs Π-kernel for different order harmonics with K “ t1, 2, 4u, indicating a clear improvement in the rate of convergence of error for higher harmonics
  For the first setting, we look at K “ t1, 2, 4u, with corresponding weight ratio as A1 : A2 : A4 “ 1 : 1 : 1. For each frequency, we look at the rate of convergence for the two kernels. We observe (Fig. 8) that the rate of error convergence for the Π-kernel is much faster than the standard two-layer NTK and especially for higher harmonics (clearest increase for the highest k “ 4), which is exactly what is expected from theoretical results.
  Next, we consider K “ t1, 3, 4, 5, 8, 12u i.e. we consider higher frequency harmonics and also introduce odd harmonics, with the respective Ak assigned the equal weights relative to each other. Recall that the eigenvalues corresponding to odd harmonics vanish for the standard kernel, meaning that we expect a much slower convergence rate for them. As before, we look at the rate of convergence for individual harmonics for the two kernels. The convergence curves are presented in Fig. 9. We verify the speed-up in higher frequencies, and an especially notable gap for odd harmonics other than k “ 1. Discussion: The empirical results strongly support our hypothesis that the Π-kernel can speed up learning higher harmonics faster than the standard two-layer kernel, even for settings outside of k " d.
sections.heading: F SYNTHETIC EXPERIMENTS WITH SINUSOIDS
sections.text: F.1 LEARNING SINUSOIDS
  F.2 ROBUSTNESS
  Robustness to Perturbations: Motivated by the observations of Rahaman et al. [37], we evaluate the susceptibility of higher frequencies to random perturbations for standard networks and Π-Nets. The setup relies on the learning task in the previous experiment. More precisely, following convergence of the network to a low-error approximation of target f˚ (denoted by fθ˚), random isotropic perturbations are introduced to the network parameters: θ “ θ˚ ` δθ̂. We monitor the effect of increasing the perturbation magnitude δ on the frequencies of interest. In the first setting (Fig. 12), we compare the effects of perturbations on the six-layer standard network to the six-layer Π-Net (with five multiplicative layers). In the second setting (Fig. 13), we directly compare the two variants of 6-layer deep Π-Nets i.e. one with five layers and the other with three layers. We observe that that Π-Nets are more robust to perturbations, especially in terms of retaining high frequency information, as compared to standard feedforward networks. We also note that the presence of more multiplicative interactions makes the network more robust.
sections.heading: G EXPERIMENTS WITH IMAGES
sections.text: G.1 IMPLEMENTATION DETAILS
  U-net type hourglass architectures provide the ideal inductive bias for a host of image restoration tasks in the DIP framework and for our experiments, we use the same architectural design as [33] for the standard networks. Since we require the output and input to have the same spatial dimensions, the number of downsampling blocks is equal to the upsampling blocks. We refer to this number as the ‘scale’ of the U-net. We provide a schematic illustration in Fig. 14.
  Note that the each down/upsampling block within itself contains convolutions, normalization and pooling but we abstract those details from the schematic for clarity. For the Π-network architecture, we consider a product of two polynomials model which modifies the standard network by introducing multiplicative connections. An illustrative schematic is presented in Fig. 15.
  Remark: In addition to the architecture, it is important to note another important detail about using Π-Nets. In practice, we recommend the use of two separate learning rates while training Π-Nets, wherein the learning rate for multiplicative connection parameters specifically is lower than the learning rate corresponding the parameters in the feed-forward part of the network. It leads to more stability while training.
  G.2 ADDITIONAL EXPERIMENTS
sections.heading: H ADDITIONAL EXPERIMENTS IN CLASSIFICATION
sections.text: Our results allow us to make a further overarching conclusion that Π-Nets, in addition to picking up higher frequencies w.r.t inputs faster (as demonstrated in the DIP settings), in classification settings can also pick up high frequency variations in the decision boundaries faster, thus making our general claims about the spectral bias of Π-Nets and multiplicative interactions stronger.