Year: 2022

Abstract:
While the class of Polynomial Nets demonstrates comparable performance to neural networks (NN), it currently has neither theoretical generalization characterization nor robustness guarantees. To this end, we derive new complexity bounds for the set of Coupled CP-Decomposition (CCP) and Nested Coupled CP-decomposition (NCP) models of Polynomial Nets in terms of the `∞-operator-norm and the `2operator norm. In addition, we derive bounds on the Lipschitz constant for both models to establish a theoretical certificate for their robustness. The theoretical results enable us to propose a principled regularization scheme that we also evaluate experimentally in six datasets and show that it improves the accuracy as well as the robustness of the models to adversarial perturbations. We showcase how this regularization can be combined with adversarial training, resulting in further improvements.

heading: 1 INTRODUCTION
text: Recently, high-degree Polynomial Nets (PNs) have been demonstrating state-of-the-art performance in a range of challenging tasks like image generation (Karras et al., 2019; Chrysos and Panagakis, 2020), image classification (Wang et al., 2018), reinforcement learning (Jayakumar et al., 2020), non-euclidean representation learning (Chrysos et al., 2020) and sequence models (Su et al., 2020). In particular, in public benchmarks like the Face verification on MegaFace task1 (KemelmacherShlizerman et al., 2016), Polynomial Nets are currently the top performing model.
A major advantage of Polynomial Nets over traditional Neural Networks2 is that they are compatible with efficient Leveled Fully Homomorphic Encryption (LFHE) protocols (Brakerski et al., 2014). Such protocols allow efficient computation on encrypted data, but they only support addition or multiplication operations i.e., polynomials. This has prompted an effort to adapt neural networks by replacing typical activation functions with polynomial approximations (Gilad-Bachrach et al., 2016; Hesamifard et al., 2018). Polynomial Nets do not need any adaptation to work with LFHE.
Without doubt, these arguments motivate further investigation about the inner-workings of Polynomial Nets. Surprisingly, little is known about the theoretical properties of such high-degree polynomial expansions, despite their success. Previous work on PNs (Chrysos et al., 2020; Chrysos and Panagakis, 2020) have focused on developing the foundational structure of the models as well as their training, but do not provide an analysis of their generalization ability or robustness to adversarial perturbations.
In contrast, such type of results are readily available for traditional feed-forward Deep Neural Networks, in the form of high-probability generalization error bounds (Neyshabur et al., 2015; Bartlett et al., 2017; Neyshabur et al., 2017; Golowich et al., 2018) or upper bounds on their Lipschitz constant (Scaman and Virmaux, 2018; Fazlyab et al., 2019; Latorre et al., 2020). Despite their similarity in the compositional structure, theoretical results for Deep Neural Networks2 do not apply to Polynomial Nets, as they are essentialy two non-overlapping classes of functions.
Why are such results important? First, they provide key theoretical quantities like the sample complexity of a hypothesis class: how many samples are required to succeed at learning in the PAC-framework. Second, they provide certified performance guarantees to adversarial perturbations
1https://paperswithcode.com/sota/face-verification-on-megaface 2with non-polynomial activation functions.
(Szegedy et al., 2014; Goodfellow et al., 2015) via a worst-case analysis c.f. Scaman and Virmaux (2018). Most importantly, the bounds themselves provide a principled way to regularize the hypothesis class and improve their accuracy or robustness.
For example, Generalization and Lipschitz constant bounds of Deep Neural Networks that depend on the operator-norm of their weight matrices (Bartlett et al., 2017; Neyshabur et al., 2017) have layed out the path for regularization schemes like spectral regularization (Yoshida and Miyato, 2017; Miyato et al., 2018), Lipschitz-margin training (Tsuzuku et al., 2018) and Parseval Networks (Cisse et al., 2017), to name a few.
Indeed, such schemes have been observed in practice to improve the performance of Deep Neural Networks. Unfortunately, similar regularization schemes for Polynomial Nets do not exist due to the lack of analogous bounds. Hence, it is possible that PNs are not yet being used to their fullest potential. We believe that theoretical advances in their understanding might lead to more resilient and accurate models. In this work, we aim to fill the gap in the theoretical understanding of PNs. We summarize our main contributions as follows:
Rademacher Complexity Bounds. We derive bounds on the Rademacher Complexity of the Coupled CP-decomposition model (CCP) and Nested Coupled CP-decomposition model (NCP) of PNs, under the assumption of a unit `∞-norm bound on the input (Theorems 1 and 3), a natural assumption in image-based applications. Analogous bounds for the `2-norm are also provided (Appendices E.3 and F.3). Such bounds lead to the first known generalization error bounds for this class of models.
Lipschitz constant Bounds. To complement our understanding of the CCP and NCP models, we derive upper bounds on their `∞-Lipschitz constants (Theorems 2 and 4), which are directly related to their robustness against `∞-bounded adversarial perturbations, and provide formal guarantees. Analogous results hold for any `p-norm (Appendices E.4 and F.4).
Regularization schemes. We identify key quantities that simultaneously control both Rademacher Complexity and Lipschitz constant bounds that we previously derived, i.e., the operator norms of the weight matrices in the Polynomial Nets. Hence, we propose to regularize the CCP and NCP models by constraining such operator norms. In doing so, our theoretical results indicate that both the generalization and the robustness to adversarial perturbations should improve. We propose a Projected Stochastic Gradient Descent scheme (Algorithm 1), enjoying the same per-iteration complexity as vanilla back-propagation in the `∞-norm case, and a variant that augments the base algorithm with adversarial traning (Algorithm 2).
Experiments. We conduct experimentation in five widely-used datasets on image recognition and on dataset in audio recognition. The experimentation illustrates how the aforementioned regularization schemes impact the accuracy (and the robust accuracy) of both CCP and NCP models, outperforming alternative schemes such as Jacobian regularization and the L2 weight decay. Indeed, for a grid of regularization parameters we observe that there exists a sweet-spot for the regularization parameter which not only increases the test-accuracy of the model, but also its resilience to adversarial perturbations. Larger values of the regularization parameter also allow a trade-off between accuracy and robustness. The observation is consistent across all datasets and all adversarial attacks demonstrating the efficacy of the proposed regularization scheme.

heading: 2 RADEMACHER COMPLEXITY AND LIPSCHITZ CONSTANT BOUNDS FOR POLYNOMIAL NETS
text: Notation. The symbol ◦ denotes the Hadamard (element-wise) product, the symbol • is the facesplitting product, while the symbol ? denotes a convolutional operator. Matrices are denoted by uppercase letters e.g., V . Due to the space constraints, a detailed notation is deferred to Appendix C.
Assumption on the input distribution. Unless explicitly mentioned otherwise, we assume an `∞-norm unit bound on the input data i.e., ‖z‖∞ ≤ 1 for any input z. This is the most common assumption in image-domain applications in contemporary deep learning, i.e., each pixel takes values in [−1, 1] interval. Nevertheless, analogous results for `2-norm unit bound assumptions are presented in Appendices E.3, E.4, F.3 and F.4.
We now introduce the basic concepts that will be developed throughout the paper i.e., the Rademacher Complexity of a class of functions (Bartlett and Mendelson, 2002) and the Lipschitz constant.
Definition 1 (Empirical Rademacher Complexity). Let Z = {z1, . . . ,zn} ⊆ Rd and let {σj : j = 1, . . . , n} be independent Rademacher random variables i.e., taking values uniformly in {−1,+1}. Let F be a class of real-valued functions over Rd. The Empirical Rademacher complexity of F with respect to Z is defined as follows:
RZ(F) := Eσ [
sup f∈F
1
n n∑ j=1 σjf(zj) ] .
Definition 2 (Lipschitz constant). Given two normed spaces (X , ‖ · ‖X ) and (Y , ‖ · ‖Y ), a function f : X → Y is called Lipschitz continuous with Lipschitz constant K ≥ 0 if for all x1, x2 in X :
‖f(x1)− f(x2)‖Y ≤ K‖x1 − x2‖X .

heading: 2.1 COUPLED CP-DECOMPOSITION OF POLYNOMIAL NETS (CCP MODEL)
text: The Coupled CP-Decomposition (CCP) model of PNs (Chrysos et al., 2020) leverages a coupled CP Tensor decomposition (Kolda and Bader, 2009) to vastly reduce the parameters required to describe a high-degree polynomial, and allows its computation in a compositional fashion, much similar to a feed-forward pass through a traditional neural network. The CCP model was used in Chrysos and Panagakis (2020) to construct a generative model. CCP can be succintly defined as follows:
f(z) = C ◦ki=1 Uiz , (CCP)
where z ∈ Rd is the input data, f(z) ∈ Ro is the output of the model and Un ∈ Rm×d,C ∈ Ro×m are the learnable parameters, where m ∈ N is the hidden rank. In Fig. 1 we provide a schematic of the architecture, while in Appendix D.1 we include further details on the original CCP formulation (and how to obtain our equivalent re-parametrization) for the interested reader.
In Theorem 1 we derive an upper bound on the complexity of CCP models with bounded `∞-operatornorms of the face-splitting product of the weight matrices. Its proof can be found in Appendix E.1. For a given CCP model, we derive an upper bound on its `∞-Lipschitz constant in Theorem 2 and its proof is given in Appendix E.4.1.
Theorem 1. Let Z = {z1, . . . ,zn} ⊆ Rd and suppose that ‖zj‖∞ ≤ 1 for all j = 1, . . . , n. Let FkCCP := { f(z) = 〈 c, ◦ki=1Uiz 〉 : ‖c‖1 ≤ µ, ∥∥•ki=1Ui∥∥∞ ≤ λ} . The Empirical Rademacher Complexity of CCPk (k-degree CCP polynomials) with respect to Z is bounded as:
RZ(FkCCP) ≤ 2µλ √ 2k log(d)
n .
Proof sketch of Theorem 1 We now describe the core steps of the proof. For the interested reader, the complete and detailed proof steps are presented in Appendix E.1. First, Hölder’s inequality is used to bound the Rademacher complexity as:
RZ(FkCCP) = E sup f∈FkCCP
1
n
〈 c,
n∑ j=1 [σj◦ki=1(Uizj)]
〉 ≤ E sup
f∈FkCCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σj◦ki=1(Uizj)] ∥∥∥∥∥∥ ∞ .
(1)
This shows why the factor ‖c‖1 ≤ µ appears in the final bound. Then, using the mixed product property (Slyusar, 1999) and its extension to repeated Hadamard products (Lemma 7 in Appendix C.3), we can rewrite the summation in the right-hand-side of (1) as follows:
n∑ j=1 σj◦ki=1(Uizj) = n∑ j=1 σj •ki=1 (Ui) ∗ki=1 (zj) = •ki=1(Ui) n∑ j=1 σj ∗ki=1 (zj) .
This step can be seen as a linearization of the polynomial by lifting the problem to a higher dimensional space. We use this fact and the definition of the operator norm to further bound the term inside the `∞-norm in the right-hand-side of (1). Such term is bounded as the product of the `∞-operator norm of •ki=1(Ui), and the `∞-norm of an expression involving the Rademacher variables σj and the vectors ∗ki=1(zj). Finally, an application of Massart’s Lemma (Shalev-Shwartz and Ben-David (2014), Lemma 26.8) leads to the final result. Theorem 2. The Lipschitz constant (with respect to the `∞-norm) of the function defined in Eq. (CCP), restricted to the set {z ∈ Rd : ‖z‖∞ ≤ 1} is bounded as:
Lip∞(f) ≤ k‖C‖∞ k∏ i=1 ‖Ui‖∞ .

heading: 2.2 NESTED COUPLED CP-DECOMPOSITION (NCP MODEL)
text: The Nested Coupled CP-Decomposition (NCP) model leverages a joint hierarchical decomposition, which provided strong results in both generative and discriminative tasks in Chrysos et al. (2020). A slight re-parametrization of the model (Appendix D.2) can be expressed with the following recursive relation:
x1 = (A1z) ◦ (s1), xn = (Anz) ◦ (Snxn−1), f(z) = Cxk , (NCP) where z ∈ Rd is the input vector and C ∈ Ro×m,An ∈ Rm×d, Sn ∈ Rm×m and s1 ∈ Rm are the learnable parameters. In Fig. 1 we provide a schematic of the architecture.
In Theorem 3 we derive an upper bound on the complexity of NCP models with bounded `∞-operatornorm of a matrix function of its parameters. Its proof can be found in Appendix F.1. For a given NCP model, we derive an upper bound on its `∞-Lipschitz constant in Theorem 4 and its proof is given in Appendix F.4.1. Theorem 3. Let Z = {z1, . . . ,zn} ⊆ Rd and suppose that ‖zj‖∞ ≤ 1 for all j = 1, . . . , n. Define the matrix Φ(A1,S1, . . . ,An,Sn) := (Ak •Sk) ∏k−1 i=1 I⊗Ai •Si. Consider the class of functions:
FkNCP := {f(z) as in (NCP) : ‖C‖∞ ≤ µ, ‖Φ(A1,S1, . . . ,Ak,Sk)‖∞ ≤ λ} , where C ∈ R1×m (single output), thus, we will write it as c, and the corresponding bound also becomes ‖c‖1 ≤ µ. The Empirical Rademacher Complexity of NCPk (k-degree NCP polynomials) with respect to Z is bounded as:
RZ(FkNCP) ≤ 2µλ √ 2k log(d)
n .
Theorem 4. The Lipschitz constant (with respect to the `∞-norm) of the function defined in Eq. (NCP), restricted to the set {z ∈ Rd : ‖z‖∞ ≤ 1} is bounded as:
Lip∞(f) ≤ k‖C‖∞ k∏ i=1 (‖Ai‖∞‖Si‖∞) .

heading: 3 ALGORITHMS
text: By constraining the quantities in the upper bounds on the Rademacher complexity (Theorems 1 and 3), we can regularize the empirical loss minimization objective (Mohri et al., 2018, Theorem 3.3). Such method would prevent overfitting and can lead to an improved accuracy. However, one issue with the quantities involved in Theorems 1 and 3, namely∥∥•ki=1Ui∥∥∞ , ∥∥∥∥∥(Ak • Sk) k−1∏ i=1 I ⊗Ai • Si ∥∥∥∥∥ ∞ ,
is that projecting onto their level sets correspond to a difficult non-convex problem. Nevertheless, we can control an upper bound that depends on the `∞-operator norm of each weight matrix:
Lemma 1. It holds that ∥∥•ki=1Ui∥∥∞ ≤∏ki=1 ‖Ui‖∞.
Lemma 2. It holds that ∥∥∥(Ak • Sk)∏k−1i=1 I ⊗Ai • Si∥∥∥∞ ≤∏ki=1 ‖Ai‖∞‖Si‖∞.
The proofs of Lemmas 1 and 2 can be found in Appendix E.2 and Appendix F.2. These results mean that by constraining the operator norms of each weight matrix, we can control the overall complexity of the CCP and NCP models.
Projecting a matrix onto an `∞-operator norm ball is a simple task that can be achieved by projecting each row of the matrix onto an `1-norm ball, for example, using the well-known algorithm from Duchi et al. (2008). The final optimization objective for training a regularized CCP is the following:
min C,U1,...,Uk
1
n n∑ i=1 L(C,U1, . . . ,Uk;xi, yi) subject to ‖C‖∞ ≤ µ, ‖Ui‖∞ ≤ λ , (2)
where (xi, yi)ni=1 is the training dataset, L is the loss function (e.g., cross-entropy) and µ, λ are the regularization parameters. We notice that the constraints on the learnable parameters Ui and C have the effect of simultaneously controlling the Rademacher Complexity and the Lipschitz constant of the CCP model. For the NCP model, an analogous objective function is used.
To solve the optimization problem in Eq. (2) we will use a Projected Stochastic Gradient Descent method Algorithm 1. We also propose a variant that combines Adversarial Training with the projection step (Algorithm 2) with the goal of increasing robustness to adversarial examples.
Algorithm 1: Projected SGD Input: dataset Z, learning rate {γt > 0}T−1t=0 , iterations T , hyper-parameters R, f , Loss L. Output: model with parameters θ. Initialize θ. for t = 0 to T − 1 do
Sample (x, y) from Z θ = θ − γtOθL(θ;x, y). if t mod f = 0 then θ = ∏ {θ:‖θ‖∞≤R}(θ)
Algorithm 2: Projected SGD + Adversarial Training
Input: dataset Z, learning rate {γt > 0}T−1t=0 , iterations T and n, hyper-parameters R, f , and α, Loss L Output: model with parameters θ. Initialize θ. for t = 0 to T − 1 do
Sample (x, y) from Z for i = 0 to n− 1 do xadv = ∏ {x′ :||x′−x||∞≤ } {x+ α∇xL(θ;x, y)} θ = θ − γtOθL(θ;xadv, y) if t mod f = 0 then θ = ∏ {θ:‖θ‖∞≤R}(θ)
In Algorithms 1 and 2 the parameter f is set in practice to a positive value, so that the projection (denoted by Π) is made only every few iterations. The variable θ represents the weight matrices of the model, and the projection in the last line should be understood as applied independently for every weight matrix. The regularization parameter R corresponds to the variables µ, λ in Eq. (2).
Convolutional layers Frequently, convolutions are employed in the literature, especially in the image-domain. It is important to understand how our previous results extend to this case, and how the proposed algorithms work in that case. Below, we show that the `∞-operator norm of the convolutional layer (as a linear operator) is related to the `∞-operator norm of the kernel after a reshaping operation. For simplicity, we consider only convolutions with zero padding.
We study the cases of 1D, 2D and 3D convolutions. For clarity, we mention below the result for the 3D convolution, since this is relevant to our experimental validation, and we defer the other two cases to Appendix G.
Theorem 5. LetA ∈ Rn×m×r be an input image and letK ∈ Rh×h×r×o be a convolutional kernel with o output channels. For simplicity assume that k ≤ min(n,m) is odd. Denote byB = K ?A the output of the convolutional layer. Let U ∈ Rnmo×nmr be the matrix such that vec(B) = Uvec(A) i.e., U is the matrix representation of the convolution. Let M(K) ∈ Ro×hhr be the matricization of

heading: K, where each row contains the parameters of a single output channel of the convolution. It holds
text: that: ‖U‖∞ = ‖M(K)‖∞.
Thus, we can control the `∞-operator-norm of a convolutional layer during training by controlling that of the reshaping of the kernel, which is done with the same code as for fully connected layers. It can be seen that when the padding is non-zero, the result still holds.

heading: 4 NUMERICAL EVIDENCE
text: The generalization properties and the robustness of PNs are numerically verified in this section. We evaluate the robustness to three widely-used adversarial attacks in sec. 4.2. We assess whether the compared regularization schemes can also help in the case of adversarial training in sec. 4.3. Experiments with additional datasets, models (NCP models), adversarial attacks (APGDT, PGDT) and layer-wise bound (instead of a single bound for all matrices) are conducted in Appendix H due to the restricted space. The results exhibit a consistent behavior across different adversarial attacks, different datasets and different models. Whenever the results differ, we explicitly mention the differences in the main body below.

heading: 4.1 EXPERIMENTAL SETUP
text: The accuracy is reported as as the evaluation metric for every experiment, where a higher accuracy translates to better performance.
Datasets and Benchmark Models: We conduct experiments on the popular datasets of FashionMNIST (Xiao et al., 2017), E-MNIST (Cohen et al., 2017) and CIFAR-10 (Krizhevsky et al., 2014). The first two datasets include grayscale images of resolution 28 × 28, while CIFAR-10 includes 60, 000 RGB images of resolution 32 × 32. Each image is annotated with one out of the ten categories. We use two popular regularization methods from the literature for comparison, i.e., Jacobian regularization (Hoffman et al., 2019) and L2 regularization (weight decay).
Models: We report results using the following three models: 1) a 4th-degree CCP model named "PN-4", 2) a 10th-degree CCP model referenced as "PN-10" and 3) a 4th-degree Convolutional CCP model called "PN-Conv". In the PN-Conv, we have replaced all the Ui matrices with convolutional kernels. None of the variants contains any activation functions.
Hyper-parameters: Unless mentioned otherwise, all models are trained for 100 epochs with a batch size of 64. The initial value of the learning rate is 0.001. After the first 25 epochs, the learning rate is multiplied by a factor of 0.2 every 50 epochs. The SGD is used to optimize all the models, while the cross-entropy loss is used. In the experiments that include projection or adversarial training, the first 50 epochs are pre-training, i.e., training only with the cross-entropy loss. The projection is performed every ten iterations.
Adversarial Attack Settings: We utilize two widely used attacks: a) Fast Gradient Sign Method (FGSM) and b) Projected Gradient Descent (PGD). In FGSM the hyper-parameter represents the step size of the adversarial attack. In PGD there is a triple of parameters ( total, niters, iter), which represent the maximum step size of the total adversarial attack, the number of steps to perform for a single attack, and the step size of each adversarial attack step respectively. We consider the following hyper-parameters for the attacks: a) FGSM with = 0.1, b) PGD with parameters (0.1, 20, 0.01), c) PGD with parameters (0.3, 20, 0.03).

heading: 4.2 EVALUATION OF THE ROBUSTNESS OF PNS
text: In the next experiment, we assess the robustness of PNs under adversarial noise. That is, the method is trained on the train set of the respective dataset and the evaluation is performed on the test set perturbed by additive adversarial noise. That is, each image is individually perturbed based on the respective adversarial attack. The proposed method implements Algorithm 1.
The quantitative results in both Fashion-MNIST and E-MNIST using PN-4, PN-10 and PN-Conv under the three attacks are reported in Table 1. The column ‘No-proj’ exhibits the plain SGD training (i.e., without regularization), while the remaining columns include the proposed regularization, Jacobian and theL2 regularization respectively. The results without regularization exhibit a substantial decrease in accuracy for stronger adversarial attacks. The proposed regularization outperforms all methods consistently across different adversarial attacks. Interestingly, the stronger the adversarial
attack, the bigger the difference of the proposed regularization scheme with the alternatives of Jacobian and L2 regularizations.
Next, we learn the networks with varying projection bounds. The results on Fashion-MNIST and E-MNIST are visualized in Fig. 2, where the x-axis is plotted in log-scale. As a reference point, we include the clean accuracy curves, i.e., when there is no adversarial noise. Projection bounds larger than 2 (in the log-axis) leave the accuracy unchanged. As the bounds decrease, the results gradually improve. This can be attributed to the constraints the projection bounds impose into the Ui matrices.
Similar observations can be made when evaluating the clean accuracy (i.e., no adversarial noise in the test set). However, in the case of adversarial attacks a tighter bound performs better, i.e., the best accuracy is exhibited in the region of 0 in the log-axis. The projection bounds can have a substantial improvement on the performance, especially in the case of stronger adversarial attacks, i.e., PGD. Notice that all in the aforementioned cases, the intermediate values of the projection bounds yield an increased performance in terms of the test-accuracy and the adversarial perturbations.
Beyond the aforementioned datasets, we also validate the proposed method on CIFAR-10 dataset. The results in Fig. 3 and Table 2 exhibit similar patterns as the aforementioned experiments. Although the improvement is smaller than the case of Fashion-MNIST and E-MNIST, we can still obtain about 10% accuracy improvement under three different adversarial attacks.

heading: 4.3 ADVERSARIAL TRAINING (AT) ON PNS
text: Adversarial training has been used as a strong defence against adversarial attacks. In this experiment we evaluate whether different regularization methods can work in conjunction with adversarial training that is widely used as a defence method. Since multi-step adversarial attacks are computationally intensive, we utilize the FGSM attack during training, while we evaluate the trained model in all three
adversarial attacks. For this experiment we select PN-10 as the base model. The proposed model implements Algorithm 2.
The accuracy is reported in Table 3 with Fashion-MNIST on the top and E-MNIST on the bottom. In the FGSM attack, the difference of the compared methods is smaller, which is expected since similar attack is used for the training. However, for stronger attacks the difference becomes pronounced with the proposed regularization method outperforming both the Jacobian and the L2 regularization methods.
The limitations of the proposed work are threefold. Firstly, Theorem 1 relies on the `∞-operator norm of the face-splitting product of the weight matrices, which in practice we relax in Lemma 1 for performing the projection. In the future, we aim to study if it is feasible to compute the non-convex projection onto the set of PNs with bounded `∞-norm of the face-splitting product of the weight matrices. This would allow us to let go off the relaxation argument and directly optimize the original tighter Rademacher Complexity bound (Theorem 1).
Secondly, the regularization effect of the projection differs across datasets and adversarial attacks, a topic that is worth investigating in the future.
Thirdly, our bounds do not take into account the algorithm used, which corresponds to a variant of the Stochastic Projected Gradient Descent, and hence any improved generalization properties due to possible uniform stability (Bousquet and Elisseeff, 2002) of the algorithm or implicit regularization properties (Neyshabur, 2017), do not play a role in our analysis.

heading: 5 CONCLUSION
text: In this work, we explore the generalization properties of the Coupled CP-decomposition (CCP) and nested coupled CP-decomposition (NCP) models that belong in the class of Polynomial Nets (PNs). We derive bounds for the Rademacher complexity and the Lipschitz constant of the CCP and the NCP models. We utilize the computed bounds as a regularization during training. The regularization terms have also a substantial effect on the robustness of the model, i.e., when adversarial noise is added to the test set. A future direction of research is to obtain generalization bounds for this class of functions using stability notions. Along with the recent empirical results on PNs, our derived bounds can further explain the benefits and drawbacks of using PNs.

heading: ACKNOWLEDGEMENTS
text: We are thankful to Igor Krawczuk and Andreas Loukas for their comments on the paper. We are also thankful to the reviewers for providing constructive feedback. Research was sponsored by the Army Research Office and was accomplished under Grant Number W911NF-19-1-0404. This work is funded (in part) through a PhD fellowship of the Swiss Data Science Center, a joint venture between EPFL and ETH Zurich. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement number 725594 - time-data).

heading: A APPENDIX INTRODUCTION
text: The Appendix is organized as follows:
• The related work is summarized in Appendix B.
• In Appendix C the notation and the core Lemmas from the literature are described.
• Further details on the Polynomial Nets are provided in Appendix D.
• The proofs on the CCP and the NCP models are added in Appendix E and Appendix F respectively.
• The extensions of the theorems for convolutional layers and their proofs are detailed in Appendix G.
• Additional experiments are included in Appendix H.

heading: B RELATED WORK
text: Rademacher Complexity: Known bounds for the class of polynomials are a consequence of more general result for kernel methods (Mohri et al., 2018, Theorem 6.12). Support Vector Machines (SVMs) with a polynomial kernel of degree k effectively correspond to a general polynomial with the same degree. In contrast, our bound is tailored to the parametric definition of the CCP and the NCP models, which are a subset of the class of all polynomials. Hence, they are tighter than the general kernel complexity bounds.
Bounds for the class of neural networks were stablished in (Bartlett et al., 2017; Neyshabur et al., 2017), but they require a long and technical proof, and in particular it assumes an `2-bound on the input, which is incompatible with image-based applications. This bound also depend on the product of spectral norms of each layer. In contrast, our bounds are more similar in spirit to the path-norm-based complexity bounds (Neyshabur et al., 2015), as they depend on interactions between neurons at different layers. This interaction precisely corresponds to the face-splitting product between weight matrices that appears in Theorem 1.
Lipschitz constant: A variety of methods have been proposed for estimating the Lipschitz constant of neural networks. For example, Scaman and Virmaux (2018) (SVD), Fazlyab et al. (2019) (Semidefinite programming) and Latorre et al. (2020) (Polynomial Optimization) are expensive optimization methods to compute tighter bounds on such constant. These methods are unusable in our case as they would require a non-trivial adaptation to work with Polynomial Nets. In contrast we find an upper bound that applies to such family of models, and it can be controlled with efficient `∞-operator-norm projections. However, our bounds might not be the tightest. Developing tighter methods to bound and control the Lipschitz constant for Polynomial Nets is a promising direction of future research.

heading: C BACKGROUND
text: Below, we develop a detailed notation in Appendix C.1, we include related definitions in Appendix C.2 and Lemmas required for our proofs in Appendix C.3. The goal of this section is to cover many of the required information for following the proofs and the notation we follow in this work. Readers familiar with the different matrix/vector products, e.g., Khatri-Rao or face-splitting product, and with basic inequalities, e.g., Hölder’s inequality, can skip to the next section.

heading: C.1 NOTATION
text: Different matrix products and their associated symbols are referenced in Table 4, while matrix operations on a matrix A are defined on Table 5. Every matrix product, e.g., Hadamard product, can be used in two ways: a)A ◦B, which translates to Hadamard product of matricesA andB, b) ◦Ni=1Ai abbreviates the Hadamard productsA1 ◦A2 ◦ . . .AN︸ ︷︷ ︸
N products
.
The symbol xji refers to j th element of vector xi.

heading: C.2 DEFINITIONS
text: For thoroughness, we include the definitions of the core products that we use in this work. Specifically, the definitions of the Hadamard product (Definition 3), Kronecker product (Definition 4), the KhatriRao product (Definition 5), column-wise Khatri-Rao product (Definition 6) and the face-splitting product (Definition 7) are included. Definition 3 (Hadamard product). For two matrices A and B of the same dimension m × n, the Hadamard productA ◦B is a matrix of the same dimension as the operands, with elements given by
(a ◦ b)i,j = ai,jbi,j .
Definition 4 (Kronecker product). If A is an m × n matrix and B is a p × q matrix, then the Kronecker productA ⊗B is the pm × qn block matrix, given as follows:
A⊗B = a1,1B · · · a1,nB... . . . ... am,1B · · · am,nB  . Example: the Kronecker product of the matricesA ∈ R2×2 andB ∈ R2×2 is computed below:
[ a1,1 a1,2 a2,1 a2,2 ] ︸ ︷︷ ︸
A
⊗ [ b1,1 b1,2 b2,1 b2,2 ] ︸ ︷︷ ︸
B
= a1,1b1,1 a1,1b1,2 a1,2b1,1 a1,2b1,2a1,1b2,1 a1,1b2,2 a1,2b2,1 a1,2b2,2a2,1b1,1 a2,1b1,2 a2,2b1,1 a2,2b1,2 a2,1b2,1 a2,1b2,2 a2,2b2,1 a2,2b2,2  ︸ ︷︷ ︸
A⊗B
.
Definition 5 (Khatri–Rao product). The Khatri–Rao product is defined as:
A ∗B = (Ai,j ⊗Bi,j)i,j ,
in which the (i, j)-th block is the mipi × njqj sized Kronecker product of the corresponding blocks ofA andB, assuming the number of row and column partitions of both matrices is equal. The size of the product is then ( ∑ imipi)× ( ∑ i njqj).
Example: ifA andB both are 2 × 2 partitioned matrices e.g.:
A = [ A1,1 A1,2 A2,1 A2,2 ] =  a1,1 a1,2 a1,3a2,1 a2,2 a2,3 a3,1 a3,2 a3,3  ,
B = [ B1,1 B1,2 B2,1 B2,2 ] =  b1,1 b1,2 b1,3b2,1 b2,2 b2,3 b3,1 b3,2 b3,3  , then we obtain the following:
A ∗B = [ A1,1 ⊗B1,1 A1,2 ⊗B1,2 A2,1 ⊗B2,1 A2,2 ⊗B2,2 ] =  a1,1b1,1 a1,2b1,1 a1,3b1,2 a1,3b1,3a2,1b1,1 a2,2b1,1 a2,3b1,2 a2,3b1,3a3,1b2,1 a3,2b2,1 a3,3b2,2 a3,3b2,3 a3,1b3,1 a3,2b3,1 a3,3b3,2 a3,3b3,3  . Definition 6 (Column-wise Khatri–Rao product). A column-wise Kronecker product of two matrices may also be called the Khatri–Rao product. This product assumes the partitions of the matrices are their columns. In this case m1 = m, p1 = p, n = q and for each j: nj = pj = 1. The resulting product is a mp × n matrix of which each column is the Kronecker product of the corresponding columns ofA andB.
Example: the Column-wise Khatri–Rao product of the matrices A ∈ R2×3 and B ∈ R3×3 is computed below:
[ a1,1 a1,2 a1,3 a2,1 a2,2 a2,3 ] ︸ ︷︷ ︸
A
∗ [ b1,1 b1,2 b1,3 b2,1 b2,2 b2,3 b3,1 b3,2 b3,3 ] ︸ ︷︷ ︸
B
=  a1,1b1,1 a1,2b1,2 a1,3b1,3 a1,1b2,1 a1,2b2,2 a1,3b2,3 a1,1b3,1 a1,2b3,2 a1,3b3,3 a2,1b1,1 a2,2b1,2 a2,3b1,3 a2,1b2,1 a2,2b2,2 a2,3b2,3 a2,1b3,1 a2,2b3,2 a2,3b3,3  ︸ ︷︷ ︸
A∗B
.
From here on, all ∗ refer to Column-wise Khatri–Rao product. Definition 7 (Face-splitting product). The alternative concept of the matrix product, which uses row-wise splitting of matrices with a given quantity of rows. This matrix operation was named the face-splitting product of matrices or the transposed Khatri–Rao product. This type of operation is based on row-by-row Kronecker products of two matrices.
Example: the Face-splitting product of the matricesA ∈ R3×2 andB ∈ R3×2 is computed below:
[ a1,1 a1,2 a2,1 a3,2 a3,1 a3,2 ] ︸ ︷︷ ︸
A
• [ b1,1 b1,2 b2,1 b2,2 b3,1 b3,2 ] ︸ ︷︷ ︸
B
= [ a1,1b1,1 a1,2b1,1 a1,1b1,2 a1,2b1,2 a2,1b2,1 a2,2b2,1 a2,1b2,2 a2,2b2,2 a3,1b3,1 a3,2b3,1 a3,1b3,2 a3,2b3,2 ] ︸ ︷︷ ︸
A•B
.

heading: C.3 WELL-KNOWN LEMMAS
text: In this section, we provide the details on the Lemmas required for our proofs along with their proofs or the corresponding citations where the Lemmas can be found as well.
Lemma 3. (Federer, 2014) Let g, h be two composable Lipschitz functions. Then g ◦ h is also Lipschitz with Lip(g ◦ h) ≤ Lip(g)Lip(h). Here and only here ◦ represents function composition. Lemma 4. (Federer, 2014) Let f : X ⊆ Rn → Rm be differentiable and Lipschitz continuous. Let Jf (x) denote its total derivative (Jacobian) at x. Then Lipp(f) = supx∈X ‖Jf (x)‖p where ‖Jf (x)‖p is the induced operator norm on Jf (x).
Lemma 5 (Hölder’s inequality). (Cvetkovski, 2012) Let (S,Σ, µ) be a measure space and let p, q ∈ [1,∞] with 1p + 1 q = 1. Then, for all measurable real-valued functions f and g on S, it holds that: ‖fg‖1 ≤ ‖f‖p ‖g‖q .
Lemma 6 (Mixed Product Property 1). (Slyusar, 1999) The following holds:
(A1B1) ◦ (A2B2) = (A1 •A2)(B1 ∗B2) . Lemma 7 (Mixed Product Property 2). The following holds:
◦Ni=1(AiBi) = •Ni=1(Ai) ∗Ni=1 (Bi) .
Proof. We prove this lemma by induction on N .
Base case (N = 1): A1B1 = A1B1.
Inductive step: Assume that the induction hypothesis holds for a particular k, i.e., the case N = k holds. That can be expressed as:
◦ki=1(AiBi) = •ki=1(Ai) ∗ki=1 (Bi) . (3)
Then we will prove that it holds for N = k + 1:
◦k+1i=1 (AiBi) = [◦ki=1(AiBi)] ◦ (Ak+1Bk+1) = [•ki=1(Ai) ∗ki=1 (Bi)] ◦ (Ak+1Bk+1) use inductive hypothesis (Eq. (3)) = [•ki=1(Ai) •Ak+1][∗ki=1(Bi) ∗Bk+1] Lemma 6 [Mixed product property 1] = •k+1i=1 (Ai) ∗ k+1 i=1 (Bi) .
That is, the case N = k + 1 also holds true, establishing the inductive step.
Lemma 8 (Massart Lemma. Lemma 26.8 in Shalev-Shwartz and Ben-David (2014)). Let A ={a1, · · · ,aN} be a finite set of vectors in Rm. Define ā = 1N ∑N i=1 ai. Then:
R(A) ≤ max a∈A ‖a− ā‖
√ 2 logN
m .
Definition 8 (Consistency of a matrix norm). A matrix norm is called consistent on Cn,n, if ‖AB‖ ≤ ‖A‖ ‖B‖ .
holds forA,B ∈ Cn,n. Lemma 9 (Consistency of the operator norm). (Lyche, 2020) The operator norm is consistent if the vector norm ‖·‖α is defined for all m ∈ N and ‖·‖β = ‖·‖α
Proof.
‖AB‖ = max Bx6=0 ‖ABx‖α ‖x‖α = max Bx6=0 ‖ABx‖α ‖Bx‖α ‖Bx‖α ‖x‖α
≤ max y 6=0 ‖Ay‖α ‖y‖α max x 6=0 ‖Bx‖α ‖x‖α = ‖A‖ ‖B‖ .
Lemma 10. (Rao, 1970) (AC) ∗ (BD) = (A⊗B)(C ∗D) .

heading: D DETAILS ON POLYNOMIAL NETWORKS
text: In this section, we provide further details on the two most prominent parametrizations proposed in Chrysos et al. (2020). This re-parametrization creates equivalent models, but enables us to absorb the bias terms into the input terms. Firstly, we provide the re-parametrization of the CCP model in Appendix D.1 and then we create the re-parametrization of the NCP model in Appendix D.2.

heading: D.1 RE-PARAMETRIZATION OF CCP MODEL
text: The Coupled CP-Decomposition (CCP) model of PNs (Chrysos et al., 2020) leverages a coupled CP Tensor decomposition. A k-degree CCP model f(ζ) can be succinctly described by the following recursive relations:
y1 = V1ζ, yn = (Vnζ) ◦ yn−1 + yn−1, f(ζ) = Qyk + β , (4) where ζ ∈ Rδ is the input data with δ ∈ N, f(ζ) ∈ Ro is the output of the model and Vn ∈ Rµ×δ, Q ∈ Ro×µ and β ∈ Ro are the learnable parameters, where µ ∈ N is the hidden rank. In order to simplify the bias terms in the model, we will introduce a minor re-parametrization in Lemma 11 that we will use to present our results in the subsequent sections. Lemma 11. Let z = [ζ>, 1]> ∈ Rd, xn = [y>n , 1]> ∈ Rm, C = [Q,β] ∈ Ro×m, d = δ + 1,m = µ+ 1. Define:
U1 = [ V1 0 0> 1 ] ∈ Rm×d, Ui = [ Vi 1 0> 1 ] ∈ Rm×d (i > 1) .
where the boldface numbers 0 and 1 denote all-zeros and all-ones column vectors of appropriate size, respectively. The CCP model in Eq. (4) can be rewritten as f(z) = C ◦ki=1 Uiz, which is the one used in Eq. (CCP).
As a reminder before providing the proof, the core symbols for this proof are summarized in Table 6.
Table 6: Core symbols in the proof of Lemma 11.
Symbol Dimensions Definition ◦ - Hadamard (element-wise) product. ζ Rδ Input of the polynomial expansion.
f(ζ) Ro Output of the polynomial expansion. k N Degree of polynomial expansion. m N Hidden rank of the expansion. Vn Rµ×δ Learnable matrices of the expansion. Q Ro×µ Learnable matrix of the expansion. β Ro Bias of the expansion. z Rd Re-parametrization of the input. C Ro×m C = (Q,β).
Proof. By definition, we have:
x1 = [y > 1 , 1] > = [ y1 1 ] = [ V1ζ 1 ] = [ V1 0 0> 1 ] [ ζ 1 ] = [ V1 0 0> 1 ] [ζ>, 1]> = U1z .
xn = [y > n , 1] > = [ yn 1 ] = [ (Vnζ) ◦ yn−1 + yn−1 1 ] = [ Vnζ + 1 1 ] ◦ [ yn−1 1 ] = [ Vn 1 0> 1 ] [ ζ 1 ] ◦ [ yn−1 1 ] = [ Vn 1 0> 1 ] [ζ>, 1]> ◦ [y>n−1, 1]> = Unz ◦ xn−1 .
Hence, it holds that:
f(z) = Qyk + β = (Q,β) [ yk 1 ] = Cxk
= C Ukz ◦ xk−1 = C Ukz ◦ (Uk−1z) ◦ xk−2 = · · · = C Ukz ◦ (Uk−1z) ◦ · · · ◦ (U2z) ◦ x1 = C Ukz ◦ (Uk−1z) ◦ · · · ◦ (U2z) ◦ (U1z) = C ◦ki=1 (Uiz) .

heading: D.2 REPARAMETRIZATION OF THE NCP MODEL
text: The nested coupled CP decomposition (NCP) model of PNs (Chrysos et al., 2020) leverages a joint hierarchical decomposition. A k-degree NCP model f(ζ) is expressed with the following recursive relations:
y1 = (V1ζ) ◦ (b1), yn = (Vnζ) ◦ (Unyn−1 + bn), f(ζ) = Qyk + β . (5)
where ζ ∈ Rδ is the input data with δ ∈ N, f(ζ) ∈ Ro is the output of the model and Vn ∈ Rµ×δ, bn ∈ Rµ, Un ∈ Rµ×µ, Q ∈ Ro×µ and β ∈ Ro are the learnable parameters, where µ ∈ N is the hidden rank. In order to simplify the bias terms in the model, we will introduce a minor re-parametrization in Lemma 12 that we will use to present our results in the subsequent sections.
Lemma 12. Let z = [ζ>, 1]> ∈ Rd, xn = [y>n , 1]> ∈ Rm, C = [Q,β] ∈ Ro×m, d = δ + 1,m = µ+ 1. Let:
s1 = [b > 1 , 1] > ∈ Rm, Si = [ Ui bi 0> 1 ] ∈ Rm×m(i > 1), Ai = [ Vi 0 0> 1 ] ∈ Rm×d .
where the boldface numbers 0 and 1 denote all-zeros and all-ones column vectors of appropriate size, respectively. The NCP model in Eq. (5) can be rewritten as
x1 = (A1z) ◦ (s1), xn = (Anz) ◦ (Snxn−1), f(z) = Cxk . (6)
In the aforementioned Eq. (6), we have written Sn even for n = 1, when s1 is technically a vector, but this is done for convenience only and does not change the end result.

heading: E RESULT OF THE CCP MODEL

heading: E.1 PROOF OF THEOREM 1: RADEMACHER COMPLEXITY BOUND OF CCP UNDER `∞ NORM
text: To facilitate the proof below, we include the related symbols in Table 7. Below, to avoid cluttering the notation, we consider that the expectation is over σ and omit the brackets as well.

heading: Proof.
text: RZ(FkCCP) = E sup f∈FkCCP
1
n n∑ j=1 σjf(zj)
= E sup f∈FkCCP
1
n n∑ j=1 ( σj 〈 c, ◦ki=1(Uizj) 〉) = E sup
f∈FkCCP
1
n
〈 c,
n∑ j=1 [σj◦ki=1(Uizj)]
〉
≤ E sup f∈FkCCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σj◦ki=1(Uizj)] ∥∥∥∥∥∥ ∞
Lemma 5 [Hölder’s inequality]
= E sup f∈FkCCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σj •ki=1 (Ui) ∗ki=1 (zj)] ∥∥∥∥∥∥ ∞
Lemma 7 [Mixed product property]
= E sup f∈FkCCP
1 n ‖c‖1 ∥∥∥∥∥∥•ki=1(Ui) n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ ∞
≤ E sup f∈FkCCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ ∞
∥∥•ki=1(Ui)∥∥∞ ≤ µλ
n E ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ ∞ .
(7)
Next, we compute the bound of E ∥∥∥∑nj=1[σj ∗ki=1 (zj)]∥∥∥∞.
Let Zj = ∗ki=1(zj) ∈ Rd k . For each l ∈ [dk], let vl = (Zl1, . . . ,Zln) ∈ Rn. Note that ‖vl‖2 ≤√ n maxj ‖Zj‖∞. Let V = {v1, . . . ,vdk}. Then, it is true that:
E ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ ∞ = E ∥∥∥∥∥∥ n∑ j=1 σjZj ∥∥∥∥∥∥ ∞ = E d k max l=1 ∣∣∣∣∣∣ n∑ j=1 σj(vl)j ∣∣∣∣∣∣ = nR(V ) . (8)
Using Lemma 8 [Massart Lemma] we have that:
R(V ) ≤ 2 max j ‖Zj‖∞
√ 2 log (dk)/n . (9)
Then, it holds that:
RZ(FkCCP) = E sup f∈FkCCP
1
n n∑ j=1 σjf(zj)
≤ µλ n E ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ ∞
Eq. (7)
= µλ
n nR(V ) Eq. (8)
≤ 2µλmax j ‖Zj‖∞
√ 2 log(dk)/n Eq. (9)
= 2µλmax j ∥∥∗ki=1(zj)∥∥∞√2 log (dk)/n ≤ 2µλ(max
j ‖zj‖∞)
k √
2 log (dk)/n ≤ 2µλ √ 2k log (d)/n .
(10)

heading: E.2 PROOF OF LEMMA 1
text: Proof. ∥∥•ki=1(Ui)∥∥∞ = mmaxj=1 ∥∥[•ki=1(Ui)]j∥∥1 =
m max j=1 ∥∥∥⊗ki=1[U ji ]∥∥∥ 1 Definition of Face-splitting product
= m
max j=1 [ k∏ i=1 ∥∥∥U ji ∥∥∥ 1 ] Multiplicativity of absolute value
≤ k∏ i=1 [ m max j=1 ∥∥∥U ji ∥∥∥ 1 ]
= k∏ i=1 ‖Ui‖∞ .

heading: E.3 RADEMACHER COMPLEXITY BOUND UNDER `2 NORM
text: Theorem 6. Let Z = {z1, . . . ,zn} ⊆ Rd and suppose that ‖zj‖∞ ≤ 1 for all j = 1, . . . , n. Let FkCCP := { f(z) = 〈 c, ◦ki=1Uiz 〉 : ‖c‖2 ≤ µ, ∥∥•ki=1Ui∥∥2 ≤ λ} .
The Empirical Rademacher Complexity of CCPk (k-degree CCP polynomials) with respect to Z is bounded as:
RZ(FkCCP) ≤ µλ√ n .
To facilitate the proof below, we include the related symbols in Table 9. Below, to avoid cluttering the notation, we consider that the expectation is over σ and omit the brackets as well.
Table 9: Core symbols for proof of Theorem 6.
Symbol Dimensions Definition ◦ - Hadamard (element-wise) product. • - Face-splitting product. ∗ - Column-wise Khatri–Rao product. z Rd Input of the polynomial expansion.
f(z) R Output of the polynomial expansion. k N Degree of polynomial expansion. m N Hidden rank of the expansion. Ui Rm×d Learnable matrices. c R1×m Learnable matrix. µ R ‖c‖2 ≤ µ. λ R ∥∥•ki=1(Ui)∥∥2 ≤ λ.

heading: Proof.
text: RZ(FkCCP) = E sup f∈FkCCP
1
n n∑ j=1 σjf(zj)
= E sup f∈FkCCP
1
n n∑ j=1 ( σj 〈 c, ◦ki=1(Uizj) 〉) = E sup
f∈FkCCP
1
n
〈 c,
n∑ j=1 [σj◦ki=1(Uizj)]
〉
≤ E sup f∈FkCCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σj◦ki=1(Uizj)] ∥∥∥∥∥∥ 2
Lemma 5 [Hölder’s inequality]
= E sup f∈FkCCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σj •ki=1 (Ui) ∗ki=1 (zj)] ∥∥∥∥∥∥ 2
Lemma 7 [Mixed product property]
= E sup f∈FkCCP
1 n ‖c‖2 ∥∥∥∥∥∥•ki=1(Ui) n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2
≤ E sup f∈FkCCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2 ∥∥•ki=1(Ui)∥∥2 .
E ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2 = E √√√√√ ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2 2
≤ √√√√√E ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2
2
Jensen’s inequality
= √√√√E n∑ s,j [σsσj 〈 ∗ki=1(zs), ∗ki=1(zj) 〉 ]
= √√√√ n∑ j=1 [ ∥∥∗ki=1(zj)∥∥22]
= √√√√ n∑ j=1 ( k∏ i=1 ‖zj‖22) ≤ √ n .
(11)
So:
RZ(FkCCP) ≤ E sup f∈FkCCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2
∥∥•ki=1(Ui)∥∥2 ≤ supf∈FkCCP ‖c‖2 ∥∥•ki=1(Ui)∥∥2√ n
≤ µλ√ n .

heading: E.4 LIPSCHITZ CONSTANT BOUND OF THE CCP MODEL
text: We will first prove a more general result about the `p-Lipschitz constant of the CCP model. Theorem 7. The Lipschitz constant (with respect to the `p-norm) of the function defined in Eq. (CCP), restricted to the set {z ∈ Rd : ‖z‖p ≤ 1} is bounded as:
Lipp(f) ≤ k‖C‖p k∏ i=1 ‖Ui‖p .
Proof. Let g(x) = Cx and h(z) = ◦ki=1(Uiz). Then it holds that f(z) = g(h(z)). By Lemma 3, we have: Lipp(f) ≤ Lipp(g)Lipp(h). We will compute an upper bound of each function individually.
Let us first consider the function g(x) = Cx. By Lemma 4, because g is a linear map represented by a matrix C, its Jacobian is Jg(x) = C. So:
Lipp(g) = ‖C‖p := sup ‖x‖p=1 ‖Cx‖p .
where ‖C‖p is the operator norm on matrices induced by the vector p-norm.
Now, let us consider the function h(z) = ◦ki=1Uiz. Its Jacobian is given by:
dh dz = k∑ i=1 diag(◦j 6=iUjz)Ui .
Using Lemma 4 we have:
Lipp(h) ≤ sup z:‖z‖p≤1 ∥∥∥∥∥ k∑ i=1 [diag(◦j 6=i(Ujz))Ui] ∥∥∥∥∥ p
≤ sup z:‖z‖p≤1 k∑ i=1 ‖diag(◦j 6=i(Ujz))Ui‖p Triangle inequality
≤ sup z:‖z‖p≤1 k∑ i=1 ‖diag(◦j 6=i(Ujz))‖p ‖Ui‖p Lemma 9 [consistency]
≤ sup z:‖z‖p≤1 k∑ i=1 ‖◦j 6=i(Ujz)‖p ‖Ui‖p
≤ sup z:‖z‖p≤1 k∑ i=1 ∏ j 6=i (‖Ujz‖p) ‖Ui‖p
≤ sup z:‖z‖p≤1 k∑ i=1 ∏ j 6=i (‖Uj‖p ‖z‖p) ‖Ui‖p
≤ sup z:‖z‖p≤1 k∑ i=1 k∏ j=1 (‖Uj‖p)
= k k∏ j=1 ‖Uj‖p .
So:
Lipp(FL) ≤ Lipp(g)Lipp(h)
≤ k‖C‖p k∏ i=1 ‖Ui‖p .

heading: E.4.1 PROOF OF THEOREM 2
text: Proof. This is a particular case of Theorem 7 when p =∞.

heading: F RESULT OF THE NCP MODEL

heading: F.1 PROOF OF THEOREM 3: RADEMACHER COMPLEXITY OF NCP UNDER `∞ NORM

heading: Proof.
text: RZ(FkNCP) = E sup f∈FkNCP
1
n n∑ j=1 σjf(zj)
= E sup f∈FkNCP
1
n n∑ j=1 (σj 〈c,xk(zj)〉)
= E sup f∈FkNCP
1
n
〈 c,
n∑ j=1 [σjxk(zj)]
〉
≤ E sup f∈FkNCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σjxk(zj)] ∥∥∥∥∥∥ ∞
Lemma 5 [Hölder’s inequality]
= E sup f∈FkNCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σj((Akzj) ◦ (Skxk−1(zj)))] ∥∥∥∥∥∥ ∞
= E sup f∈FkNCP
1 n ‖c‖1 ∥∥∥∥∥∥ n∑ j=1 [σj((Ak • Sk)(zj ∗ xk−1(zj)))] ∥∥∥∥∥∥ ∞
Lemma 7 [Mixed product property]
= E sup f∈FkNCP
1 n ‖c‖1 ∥∥∥∥∥∥(Ak • Sk) n∑ j=1 [σj(zj ∗ xk−1(zj))] ∥∥∥∥∥∥ ∞ .
(12)
Now, because of the recursive definition of the Eq. (NCP), we obtain:
n∑ j=1 σj(zj ∗ xk−1(zj)) = n∑ j=1 σj(zj ∗ (Ak−1zj) ◦ (Sk−1xk−2(zj)))
= n∑ j=1 σj(zj ∗ ((Ak−1 • Sk−1)(zj ∗ xk−2(zj)))) Lemma 7
= n∑ j=1 σj(I ⊗ (Ak−1 • Sk−1))(zj ∗ (zj ∗ xk−2(zj))) Lemma 10
= I ⊗ (Ak−1 • Sk−1) n∑ j=1 [σj(zj ∗ (zj ∗ xk−2(zj))) .
(13)
recursively applying this argument we have:
n∑ j=1 σj(zj ∗ xk−1(zj)) = ( k−1∏ i=1 I ⊗Ai • Si ) n∑ j=1 σj ∗ki=1 (zj) . (14)
Combining the two previous equations (Eqs. (13) and (14)) inside Eq. (12) we finally obtain
RZ(FkNCP) ≤ E sup f∈FkNCP
1 n ‖c‖1 ∥∥∥∥∥∥(Ak • Sk) ( k−1∏ i=1 I ⊗Ai • Si ) n∑ j=1 σj ∗ki=1 (zj) ∥∥∥∥∥∥ ∞
≤ E sup f∈FkNCP
1 n ‖c‖1 ∥∥∥∥∥(Ak • Sk) ( k−1∏ i=1 I ⊗Ai • Si )∥∥∥∥∥ ∞ ∥∥∥∥∥∥ n∑ j=1 σj ∗ki=1 (zj) ∥∥∥∥∥∥ ∞
≤ µλ n E ∥∥∥∥∥∥ n∑ j=1 σj ∗ki=1 (zj) ∥∥∥∥∥∥ ∞
= µλ
n nR(V ) = µλR(V ) . Eq. (8) .
following the same arguments as in Eq. (10) it follows that:
RZ(FkNCP) ≤ 2µλ √ 2k log (d)
n .

heading: F.2 PROOF OF LEMMA 2
text: Proof.
∥∥∥∥∥(Ak • Sk) k−1∏ i=1 I ⊗Ai • Si ∥∥∥∥∥ ∞ ≤ ‖Ak • Sk‖∞ k−1∏ i=1 ‖I ⊗Ai • Si‖∞ Lemma 9 [consistent]
= k∏ i=1 ‖Ai • Si‖∞
≤ k∏ i=1 ‖Ai‖∞‖Si‖∞ Lemma 1 .

heading: F.3 RADEMACHER COMPLEXITY UNDER `2 NORM
text: Theorem 8. Let Z = {z1, . . . ,zn} ⊆ Rd and suppose that ‖zj‖∞ ≤ 1 for all j = 1, . . . , n. Define the matrix Φ(A1,S1, . . . ,An,Sn) := (Ak •Sk) ∏k−1 i=1 I⊗Ai •Si. Consider the class of functions:
FkNCP := {f(z) as in (NCP) : ‖C‖2 ≤ µ, ‖Φ(A1,S1, . . . ,Ak,Sk)‖2 ≤ λ} ,
where C ∈ R1×m (single output case). The Empirical Rademacher Complexity of NCPk (k-degree NCP polynomials) with respect to Z is bounded as:
RZ(FkNCP) ≤ µλ√ n .

heading: Proof.
text: RZ(FkNCP) = E sup f∈FkNCP
1
n n∑ j=1 σjf(zj)
= E sup f∈FkNCP
1
n n∑ j=1 (σj 〈c,xk(zj)〉)
= E sup f∈FkNCP
1
n
〈 c,
n∑ j=1 [σjxk(zj)]
〉
≤ E sup f∈FkNCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σjxk(zj)] ∥∥∥∥∥∥ 2
Lemma 5 [Hölder’s inequality]
= E sup f∈FkNCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σj((Akzj) ◦ (Skxk−1(zj)))] ∥∥∥∥∥∥ 2
= E sup f∈FkNCP
1 n ‖c‖2 ∥∥∥∥∥∥ n∑ j=1 [σj((Ak • Sk)(zj ∗ xk−1(zj)))] ∥∥∥∥∥∥ 2
Lemma 7 [Mixed product property]
= E sup f∈FkNCP
1 n ‖c‖2 ∥∥∥∥∥∥(Ak • Sk) n∑ j=1 [σj(zj ∗ xk−1(zj))] ∥∥∥∥∥∥ 2
= E sup f∈FkNCP
1 n ‖c‖2 ∥∥∥∥∥∥(Ak • Sk) ( k−1∏ i=1 I ⊗Ai • Si ) n∑ j=1 σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2
Eq. (14)
≤ E sup f∈FkNCP
1 n ‖c‖2 ∥∥∥∥∥(Ak • Sk) ( k−1∏ i=1 I ⊗Ai • Si )∥∥∥∥∥ 2 ∥∥∥∥∥∥ n∑ j=1 σj ∗ki=1 (zj)] ∥∥∥∥∥∥ 2
≤ µλ n E ∥∥∥∥∥∥ n∑ j=1 σj ∗ki=1 (zj) ∥∥∥∥∥∥ 2
≤ µλ√ n . Eq. (11) .

heading: F.4 LIPSCHITZ CONSTANT BOUND OF THE NCP MODEL
text: Theorem 9. Let FL be the class of functions defined as FL := { x1 = (A1z) ◦ (S1),xn = (Anz) ◦ (Snxn−1), f(z) = Cxk :
‖C‖p ≤ µ, ‖Ai‖p ≤ λi, ‖Si‖p ≤ ρi, ‖z‖p ≤ 1 } .
The Lipschitz Constant of FL (k-degree NCP polynomial) under `p norm restrictions is bounded as:
Lipp(FL) ≤ kµ k∏ i=1 (λiρi) .
Proof. Let g(x) = Cx, h(z) = (Anz) ◦ (Snxn−1(z)). Then it holds that f(z) = g(h(z)).
By Lemma 3, we have: Lip(f) ≤ Lip(g)Lip(h). This enables us to compute an upper bound of each function (i.e., g, h) individually.
Let us first consider the function g(x) = Cx. By Lemma 4, because g is a linear map represented by a matrix C, its Jacobian is Jg(x) = C. So:
Lipp(g) = ‖C‖p := sup ‖x‖p=1 ‖Cx‖p = { σmax(C) if p = 2 maxi ∑ j
∣∣C(i,j)∣∣ if p =∞ . where ‖C‖p is the operator norm on matrices induced by the vector p-norm, and σmax(C) is the largest singular value of C.
Now, let us consider the function xn(z) = h(z) = (Anz) ◦ (Snxn−1(z)). Its Jacobian is given by:
Jxn = diag(Anz)SnJxn−1 + diag(Snxn−1)An, Jx1 = diag(S1)A1 .
Lipp(h) = sup z:‖z‖p≤1 ‖Jxn‖p
= sup z:‖z‖p≤1 ‖ diag(Anz)SnJxn−1 + diag(Snxn−1)An‖p
≤ sup z:‖z‖p≤1 ‖ diag(Anz)SnJxn−1‖p + ‖diag(Snxn−1)An‖p Triangle inequality
≤ sup z:‖z‖p≤1 ‖diag(Anz)‖p‖Sn‖p‖Jxn−1‖p + ‖ diag(Snxn−1)‖p‖An‖p Lemma 9 [consistent]
≤ sup z:‖z‖p≤1 ‖Anz‖p‖Sn‖p‖Jxn−1‖p + ‖Snxn−1‖p‖An‖p
≤ sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p‖Jxn−1‖p + ‖Sn‖p‖xn−1‖p‖An‖p
= sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p‖Jxn−1‖p + ‖Sn‖p‖(An−1z) ◦ (Sn−1xn−2)‖p‖An‖p
≤ sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p‖Jxn−1‖p + ‖Sn‖p‖An−1z‖p‖Sn−1xn−2‖p‖An‖p
≤ sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p‖Jxn−1‖p + ‖Sn‖p‖An−1‖p‖z‖p‖Sn−1‖p‖xn−2‖p‖An‖p
= sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p(‖Jxn−1‖p + ‖An−1‖p‖Sn−1‖p‖xn−2‖p)
≤ sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p(‖Jxn−1‖p + n−1∏ i=1 (‖Si‖p‖Ai‖p)‖z‖n−2p ) .
Then we proof the result by induction.
Inductive hypothesis:
sup z:‖z‖p≤1 ‖Jxn‖p ≤ n n∏ i=1 (‖Si‖p‖Ai‖p) .
Case k = 1:
Lipp(h) = sup z:‖z‖p≤1 ‖Jx1‖p
= ‖ diag(S1)A1‖p ≤ ‖diag(S1)‖p‖A1‖p ≤ ‖S1‖p‖A1‖p .
Case k = n:
Lipp(h) = sup z:‖z‖p≤1 ‖Jxn‖p
≤ sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p(‖Jxn−1‖p + n−1∏ i=1 (‖Si‖p‖Ai‖p)‖z‖n−2p )
≤ sup z:‖z‖p≤1 ‖An‖p‖z‖p‖Sn‖p((n− 1) n−1∏ i=1 (‖Si‖p‖Ai‖p) + n−1∏ i=1 (‖Si‖p‖Ai‖p)‖z‖n−2p )
≤ n n∏ i=1 (‖Si‖p‖Ai‖p) .
So:
Lipp(FL) ≤ Lipp(g)Lipp(h)
≤ k‖C‖p k∏ i=1 (‖Si‖p‖Ai‖p) .

heading: F.4.1 PROOF OF THEOREM 4
text: Proof. This is particular case of Theorem 9 with p =∞.

heading: G RELATIONSHIP BETWEEN A CONVOLUTIONAL LAYER AND A FULLY CONNECTED LAYER
text: In this section we discuss various cases of input/output types depending on the dimensionality of the input tensor and the output tensor. We also provide the proof of Theorem 5. Theorem 10. LetA ∈ Rn. LetK ∈ Rh be a 1-D convolutional kernel. For simplicity, we assume h is odd and h ≤ n. Let B ∈ Rn, B = K ? A be the output of the convolution. Let U be the convolutional operator i.e., the linear operator (matrix) U ∈ Rn×n such thatB = K ?A = UA. It holds that ‖U‖∞ = ‖K‖1. Theorem 11. Let A ∈ Rn×m, and let K ∈ Rh×h be a 2-D convolutional kernel. For simplicity assume h is odd number and h ≤ min(n,m). Let B ∈ Rn×m, B = K ?A be the output of the convolution. Let U be the convolutional operator i.e., the linear operator (matrix) U ∈ Rnm×nm such that vec(B) = Uvec(A). It holds that ‖U‖∞ = ‖vec(K)‖1.

heading: G.1 PROOF OF THEOREM 10
text: Proof. From,B = K ?A = UA we can obtain the following:
 u1,1 u1,2 · · · u1,n u2,1 u2,2 · · · u2,n
... ...
. . . ...
un,1 un,2 · · · un,n  (A1, · · · ,An)> = (K1, · · · ,Kh)> ? (A1, · · · ,An)> . We observe that:
ui,j =
{ K
h+1 2 +j−i if |i− j| ≤ h−12 ; 0 if |i− j| > h−12 .
Then, it holds that:
‖U‖∞ = n
max i=1 n∑ j=1 |ui,j | ≤ n max i=1 h∑ j=1 ∣∣Kj∣∣ = ‖K‖1 .

heading: G.2 PROOF OF THEOREM 11
text: Proof. We partition U into n × n partition matrices of shape m ×m. Then the (i, j)th partition matrix U(i,j) describes the relationship between theBi and theAj . So, U(i,j) is also similar to the Toeplitz matrix in the previous result.
U(i,j) =
{ Mh+1
2 +j−i if |i− j| ≤ h−12 ; 0 if |i− j| > h−12 . Meanwhile, the matrix M satisfies:
mi(s,l) =
{ k(i,h+12 +l−s)
if |s− l| ≤ h−12 ; 0 if |s− l| > h−12 .
Then, we have the following:
‖U‖∞ = n×m max i=1 n×m∑ j=1 |ui,j | ≤ n max i=1 n∑ j=1 ‖Ui,j‖∞ ≤ n max i=1 h∑ j=1 ∥∥Kj∥∥ 1 = h∑ i=1 ∥∥Ki∥∥ 1 = ‖vec(K)‖1 .
In addition, by h ≤ n, h ≤ m we have:∥∥∥U h+12 +M(h−12 )∥∥∥ 1 = h∑ i=1 ∥∥Ki∥∥ 1 . (15)
Then, it holds that ‖U‖∞ = ∑h i=1 ∥∥Ki∥∥ 1 .

heading: G.3 PROOF OF THEOREM 5
text: Proof. We partition U into o× r partition matrices of shape nm× nm. Then the (i, j)th partition of the matrix U(i,j) describes the relationship between the ith channel of B and the jth channel of A. Then, the following holds: ∥∥U(i,j)∥∥∞ = ∑hi=1 ∥∥Kiij∥∥1, whereKij means the two-dimensional tensor obtained by the third dimension ofK takes j and the fourth dimension ofK takes i.
‖U‖∞ = max n×m×o i=1 n×m×r∑ j=1 ∣∣u(i,j)∣∣ =
o−1 max l=0 n×m max i=1 r−1∑ s=0 n×m∑ j=1 ∣∣u(i+nml,j+nms)∣∣ ≤ o−1max
l=0 r−1∑ s=0 ( n×m max i=1 n×m∑ j=1 ∣∣u(i+nml,j+nms)∣∣) =
o−1 max l=0 r−1∑ s=0 ∥∥U(l+1,s+1)∥∥∞ =
o max l=1 r∑ s=1 ∥∥U(l,s)∥∥∞ =
o max l=1 r∑ s=1 h∑ i=1 ∥∥Kils∥∥1 ≤ omax
l=1 ∥∥∥K̂l∥∥∥ 1
= ∥∥∥K̂∥∥∥
∞ .
Similar to Eq. (15): for every nm rows, we choose k+12 th
row. Then its 1-norm is equal to this nm rows of the K̂’s∞-norm. So the equation holds.

heading: H AUXILIARY NUMERICAL EVIDENCE
text: A number of additional experiments are conducted in this section. Unless explicitly mentioned otherwise, the experimental setup remains similar to the one in the main paper. The following experiments are conducted below:
1. The difference between the theoretical and the algorithmic bound and their evolution during training is studied in Appendix H.1.
2. An ablation study on the hidden size is conducted in Appendix H.2.
3. An ablation study is conducted on the effect of adversarial steps in Appendix H.3.
4. We evaluate the effect of the proposed projection into the testset performance in Appendix H.4.
5. We conduct experiments on four new datasets, i.e., MNIST, K-MNIST, E-MNIST-BY, NSYNTH in Appendix H.5. These experiments are conducted in addition to the datasets already presented in the main paper.
6. In Appendix H.6 experiments on three additional adversarial attacks, i.e., FGSM-0.01, APGDT and TPGD, are performed.
7. We conduct an experiment using the NCP model in Appendix H.7.
8. The layer-wise bound (instead of a single bound for all matrices) is explored in Appendix H.8.
9. The comparison with adversarial defense methods is conducted in Appendix H.9.

heading: H.1 THEORETICAL AND ALGORITHMIC BOUND
text: As mentioned in Section 3, projecting the quantity θ = ∥∥•ki=1Ui∥∥∞ onto their level set corresponds to a difficult non-convex problem. Given that we have an upper bound
θ = ∥∥•ki=1Ui∥∥∞ ≤ Πki=1‖Ui‖∞ =: γ .
we want to understand in practice how tight is this bound. In Fig. 4 we compute the ratio γθ for PN-4. In Fig. 5 the ratio is illustrated for randomly initialized matrices (i.e., untrained networks).

heading: H.2 ABLATION STUDY ON THE HIDDEN SIZE
text: Initially, we explore the effect of the hidden rank of PN-4 and PN-10 on Fashion-MNIST. Fig. 6 exhibits the accuracy on both the training and the test-set for both models. We observe that PN-10 has a better accuracy on the training set, however the accuracy on the test set is the same in the two models. We also note that increasing the hidden rank improves the accuracy on the training set, but not on the test set.

heading: H.3 ABLATION STUDY ON THE EFFECT OF ADVERSARIAL STEPS
text: Our next experiment scrutinizes the effect of the number of adversarial steps on the robust accuracy. We consider in all cases a projection bound of 1, which provides the best empirical results. We vary the number of adversarial steps and report the accuracy in Fig. 7. The results exhibit a similar performance both in terms of the dataset (i.e., Fashion-MNIST and K-MNIST) and in terms of the network (PN-4 and PN-Conv). Notice that when the adversarial attack has more than 10 steps the performance does not vary significantly from the performance at 10 steps, indicating that the projection bound is effective for stronger adversarial attacks.
‖•ki=1Ui‖∞
for
different log bound values for PN-4. In (b), (d) the exact values of the two bounds are computed over the course of the unregularized training. Notice that there is a gap between the two bounds, however importantly the two bounds are increasing at the same rate, while their ratio is close to 1.
2 3 4 5 6 7 8 Log of hidden rank
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
Ra tio
/
PN-10
(a)
2 4 6 8 10 12 14 Depth
1.0
1.1
1.2
1.3
1.4
1.5
1.6
Ra tio
/
hidden rank = 16
(b)
Figure 5: Visualization of the ratio
∏k
i=1 ‖Ui‖∞ ‖•ki=1Ui‖∞ in a randomly initialized network (i.e., using normal
distribution random matrices). Specifically, in (a) we visualize the ratio for different log hidden rank values for PN-10. In (b) we visualize the ratio for different depth values for hidden rank = 16. Neither of two plots contain any regularization.

heading: H.4 EVALUATION OF THE ACCURACY OF PNS
text: In this experiment, we evaluate the accuracy of PNs. We consider three networks, i.e., PN-4, PN-10 and PN-Conv, and train them under varying projection bounds using Algorithm 1. Each model is evaluated on the test set of (a) Fashion-MNIST and (b) E-MNIST.
The accuracy of each method is reported in Fig. 8, where the x-axis is plotted in log-scale (natural logarithm). The accuracy is better for bounds larger than 2 (in the log-axis) when compared to tighter bounds (i.e., values less than 0). Very tight bounds stifle the ability of the network to learn, which explains the decreased accuracy. Interestingly, PN-4 reaches similar accuracy to PN-10 and PN-Conv in Fashion-MNIST as the bound increases, while in E-MNIST it cannot reach the same performance as the bound increases. The best bounds for all three models are observed in the intermediate values, i.e., in the region of 1 in the log-axis for PN-4 and PN-10.
We scrutinize further the projection bounds by training the same models only with cross-entropy loss (i.e., no bound regularization). In Table 10, we include the accuracy of the three networks with and without projection. Note that projection consistently improves the accuracy, particularly in the case of larger networks, i.e., PN-10.

heading: H.5 EXPERIMENTAL RESULTS ON ADDITIONAL DATASETS
text: To validate even further we experiment with additional datasets. We describe the datasets below and then present the robust accuracy in each case. The experimental setup remains the same as in Section 4.2 in the main paper. As a reminder, we are evaluating the robustness of the different models under adversarial noise.
Dataset details: There are six datasets used in this work:
1. Fashion-MNIST (Xiao et al., 2017) includes grayscale images of clothing. The training set consists of 60, 000 examples, and the test set of 10, 000 examples. The resolution of each image is 28× 28, with each image belonging to one of the 10 classes.
2. E-MNIST (Cohen et al., 2017) includes handwritten character and digit images with a training set of 124, 800 examples, and a test set of 20, 800 examples. The resolution of each image is 28× 28. E-MNIST includes 26 classes. We also use the variant EMNIST-BY that includes 62 classes with 697, 932 examples for training and 116, 323 examples for testing.
3. K-MNIST (Clanuwat et al., 2018) depicts grayscale images of Hiragana characters with a training set of 60, 000examples, and a test set of 10, 000 examples. The resolution of each image is 28× 28. K-MNIST has 10 classes.
4. MNIST (Lecun et al., 1998) includes handwritten digits images. MNIST has a training set of 60, 000 examples, and a test set of 10, 000 examples. The resolution of each image is 28× 28.
5. CIFAR-10 (Krizhevsky et al., 2014) depicts images of natural scenes. CIFAR-10 has a training set of 50, 000 examples, and a test set of 10, 000 examples. The resolution of each RGB image is 32× 32.
6. NSYNTH (Engel et al., 2017) is an audio dataset containing 305, 979 musical notes, each with a unique pitch, timbre, and envelope.
We provide a visualization3 of indicative samples from MNIST, Fashion-MNIST, K-MNIST and E-MNIST in Fig. 9.
We originally train PN-4, PN-10 and PN-Conv without projection bounds. The results are reported in Table 11 (columns titled ‘No proj’) for MNIST and K-MNIST, Table 13 (columns titled ‘No proj’) for E-MNIST-BY and Table 14 (columns titled ‘No proj’) for NSYNTH. Next, we consider the performance under varying projection bounds; the accuracy in each case is depicted in Fig. 10 for K-MNIST, MNIST and E-MNIST-BY and Fig. 11 for NSYNTH. The figures (and the tables) depict the same patterns that emerged in the two main experiments, i.e., the performance can be vastly improved for intermediate values of the projection bound. Similarly, we validate the performance when using adversarial training. The results in Table 12 demonstrate the benefits of using projection bounds even in the case of adversarial training.

heading: H.6 EXPERIMENTAL RESULTS OF MORE TYPES OF ATTACKS
text: To further verify the results of the main paper, we conduct experiments with three additional adversarial attacks: a) FGSM with = 0.01, b) Projected Gradient Descent in Trades (TPGD) (Zhang et al., 2019), c) Targeted Auto-Projected Gradient Descent (APGDT) (Croce and Hein, 2020). In TPGD and APGDT, we use the default parameters for a one-step attack.
The quantitative results are reported in Table 15 for four datasets and the curves of Fashion-MNIST and E-MNIST are visualized in Fig. 12 and the curves of K-MNIST and MNIST are visualized in Fig. 13. The results in both cases remain similar to the attacks in the main paper, i.e., the proposed projection improves the performance consistently across attacks, types of networks and adversarial attacks.

heading: H.7 EXPERIMENTAL RESULTS IN NCP MODEL
text: To complement, the results of the CCP model, we conduct an experiment using the NCP model. That is, we use a 4th degree polynomial expansion, called NCP-4, for our experiment. We conduct an
3The samples were found in https://www.tensorflow.org/datasets/catalog.
experiment in the K-MNIST dataset and present the result with varying bound in Fig. 14. Notice that the patterns remain similar to the CCP model, i.e., intermediate values of the projection bound can increase the performance significantly.

heading: H.8 LAYER-WISE BOUND
text: To assess the flexibility of the proposed method, we assess the performance of the layer-wise bound. In the previous sections, we have considered using a single bound for all the matrices, i.e., ‖Ui‖∞ ≤ λ, because the projection for a single matrix has efficient projection algorithms. However, Lemma 1 enables each matrix Ui to have a different bound λi. We assess the performance of having different bounds for each matrix Ui.
We experiment on PN-4 that we set a different projection bound for each matrix Ui. Specifically, we use five different candidate values for each λi and then perform the grid search on the FashionMNIST FGSM-0.01 attack. The results on Fashion-MNIST in Table 16 exhibit how the layer-wise bounds outperform the previously used single bound4. The best performing values for PN-4 are λ1 = 1.5, λ2 = 2, λ3 = 1.5, λ4 = 2, µ = 0.8. The values of λi in the first few layers are larger, while the value in the output matrix C is tighter.
To scrutinize the results even further, we evaluate whether the best performing λi can improve the performance in different datasets and the FGSM-0.1 attack. In both cases, the best performing λi can improve the performance of the single bound.
4The single bound is mentioned as ‘Our method’ in the previous tables. In this experiment both ‘single bound’ and ‘layer-wise bound’ are proposed.

heading: H.9 ADVERSARIAL DEFENSE METHOD
text: One frequent method used against adversarial perturbations are the so called adversarial defense methods. We assess the performance of adversarial defense methods on the PNs when compared with the proposed method.
We experiment on PN-4 in Fashion-MNIST. We chose three different methods: gaussian denoising, median denoising and guided denoising (Liao et al., 2018). Gaussian denoising and median denoising are the methods of using gaussian filter and median filter for feature denoising (Xie et al., 2019). The results in Table 17 show that in both attacks our method performs favourably to the adversarial defense methods.

Review summaries:
This paper analyzed the empirical Rademacher complexity of polynomial nets, the Coupled CP-Decomposition (CCP), and Nested Coupled CP-decomposition(NCP) nets. The corresponding Lipschitz constants for both models are derived. The Projected SGD and the Projected SGD + Adversarial Training algorithms are used to optimize the CCP and the NCP with the Lipschitz constants. Experiments on classification and adversarial attack tasks are performed to evaluate the proposed method. Results show that the proposed projection regularization is better than the Jacobian regularization and the L2 regularization methods in adversarial attack tasks. 

The paper provides theoretical guarantees for polynomial neural networks, specifically bounds on Empirical Rademacher Complexity (which implies generalization guarantees) and on the Lipschitz constant of the networks (implies robustness). Polynomial nets have shown state-of-the-art results for several machine learning tasks. Specifically CCP and NCP models of Chrysos et al. (2020) are considered, and the bounds assume bounded l-infinity-norm which is meaningful for image inputs. Further an algorithm is proposed which adds regularization to polynomial networks using the terms in the theoretical upper bounds (corresponding to weight norms) which leads to gain in accuracy, especially in adversarial settings (where the regularization is combined with adversarial training).

The paper introduces Rademacher Complexity bounds and Lipschitz constant bounds (both $\ell_2$ and $\ell_\infty$) for two families of polynomial networks. While prior theoretical results on complexity and Lipschitz constant bounds exist for standard neural networks that consist of linear layers and activation functions, these results do not directly apply to Polynomial networks due to their vastly different structures from standard neural networks.   On the practical side, the authors found that imposing $L_\infty$ norm constraint on the weight matrices of the polynomial networks is effective against adversarial perturbation. In particular, the paper focuses on robust image and audio recognition tasks and found out that imposing $L_\infty$ norm constraint via projected gradient descent is more effective than other commonly deployed regularization techniques (in standard networks) such as Jacobian regularization and weight decay. 

The paper derives empirical Rademacher complexity bounds for two variants of PN, along with their Lipschitz constants. The derived bounds motivate practical regularization schemes that are implemented based on projected gradient descent. Empirical results demonstrate the usefulness of the proposed method.