{'fraction_non_alphanumeric': 0.042120332549678875, 'fraction_numerical': 0.014420993183519616, 'mean_word_length': 4.016801739474205, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 62, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Neural Spectrahedra and Semidefinite Lifts: Global Convex Optimization of Polynomial Activation Neural Networks in Fully Polynomial-Time

The training of two-layer neural networks with nonlinear activation functions is an important non-convex optimization problem with numerous applications and promising performance in layerwise deep learning. In this paper, we develop exact convex optimization formulations for two-layer neural networks with second degree polynomial activations based on semidefinite programming. Remarkably, we show that semidefinite lifting is always exact and therefore computational complexity for global optimization is polynomial in the input dimension and sample size for all input data. The developed convex formulations are proven to achieve the same global optimal solution set as their non-convex counterparts. More specifically, the globally optimal two-layer neural network with polynomial activations can be found by solving a semidefinite program (SDP) and decomposing the solution using a procedure we call Neural Decomposition. Moreover, the choice of regularizers plays a crucial role in the computational tractability of neural network training. We show that the standard weight decay regularization formulation is NP-hard, whereas other simple convex penalties render the problem tractable in polynomial time via convex programming. We extend the results beyond the fully connected architecture to different neural network architectures including networks with vector outputs and convolutional architectures with pooling. We provide extensive numerical simulations showing that the standard backpropagation approach often fails to achieve the global optimum of the training loss. The proposed approach is significantly faster to obtain better test accuracy compared to the standard backpropagation procedure.


Introduction
We study neural networks from the optimization perspective by deriving equivalent convex optimization formulations with identical global optimal solution sets. The derived convex problems have important theoretical and practical implications concerning the computational complexity of optimal training of neural network models. Moreover, the convex optimization perspective provides a more concise parameterization of neural network models that enables further analysis of their interesting properties.
In non-convex optimization, the choice of optimization method and its internal hyperparameters, such as initialization, mini-batching and step sizes, have a considerable effect on the quality of the learned model. This is in sharp contrast to convex optimization problems, where locally optimal solutions are globally optimal and optimizer parameters have no influence on the solution and therefore the model. Moreover, the solutions of convex optimization problems can be obtained in a very robust, efficient and reproducible manner thanks to the elegant and extensively studied structure of convex programs. Therefore, our convex optimization based globally optimal training procedure enables the study of the neural network model and the optimization procedure in a decoupled way. For instance, step sizes employed in the optimization can be considered hyperparameters of non-convex models, which affect the model quality and may require extensive tuning. For a classification task, in our convex optimization formulation, step sizes as well as the choice of the optimizers are no longer hyperparameters to obtain better classification accuracy. Any convex optimization solver can be applied to the convex problem to obtain a globally optimal model. Various types of activation functions were proposed in the literature as nonlinearities in neural network layers. Among the most widely adopted ones is the ReLU (rectified linear unit) activation given by σ(u) = max(0, u). A recently proposed alternative is the swish activation σ(u) = u(1+e −u ) −1 , which performs comparably well [39]. Another important class is the polynomial activation where the activation function is a scalar polynomial of a fixed degree. We focus on second degree polynomial activation functions, i.e., σ(u) = au 2 + bu + c. Although polynomial coefficients a, b, c can be regarded as hyperparameters, it is often sufficient to choose the coefficients in order to approximate a target nonlinear activation function such as the ReLU or swish activation. ReLU and swish activations are plotted in Figure 1 along with their second degree polynomial approximations.
Our derivation of the convex program for polynomial activations leverages convex duality and the S-procedure, and can be stated as a simple semidefinite program (SDP). We refer the reader to [38] for a survey of the S-procedure and applications in SDPs. In addition, another commonly used activation function in the literature, quadratic activation, is a special case of polynomial activations (b = c = 0) and we devote a separate section to this case (Section 5). The corresponding convex program is an SDP and takes a simpler form.
Main aspects of our work that differ from others in the literature that study the optimization landscape of two-layer neural networks (e.g. see section 1.2) are the following: Our results (1) provide global optimal solutions in fully polynomial time (polynomial in all problem parameters), (2) uncover an important role of the regularizer in computational tractability, (3) hold for arbitrary loss function and other network architectures such as vector output, convolutional and pooling, (4) are independent of the choice of the numerical optimizer and its parameters.
We summarize the types of neural network architectures considered in this work and the corresponding convex problems that we have derived to train them to global optimality in Table 1. The fourth column of Table 1 shows the upper bounds for critical width m * , i.e., the optimal number of neurons that one needs for global optimization of any problems with number of neurons m ≥ m * . The fifth column, named "construction algorithm", refers to the method for obtaining the optimal neural network weights from the solution of the associated convex program. The last column contains the references to the theorems for each result.

Overview of Our Contributions
• We show that the standard optimization formulation for training neural networks f θ (x) = m j=1 σ(x T u j )α j with trainable parameters θ = (u 1 , . . . , u m , α 1 , . . . , α m ) and degree two polynomial activations σ(u) = au 2 + bu + c, training data X = [x 1 , . . . , x n ] T ∈ R n×d , y ∈ R n , and is computationally intractable via a reduction to the NP-hard subset sum problem.
• Surprisingly, for quadratic activation networks, σ(u) = u 2 , we show that modifying the quadratic weight decay regularization to cubic regularization enables global optimization in fully polynomial time via convex semidefinite programming. The computational complexity is polynomial in all problem parameters (n, d, m).
• Furthermore, for any degree two polynomial activation σ, the non-convex neural network training problem min θ s.t. uj 2 =1 ,∀j∈ [m] (f θ (X), y) + β α 1 ( 3) can be equivalently stated as a convex semidefinite problem and globally solved in fully polynomial time. In fact, the cubic regularization strategy in (2) is a special case of this convex program. The result holds universally for all input data without any conditions and also holds when β → 0.
• In deriving the convex formulations, we identify a concise re-parameterization of the neural network parameters that enables exact convexification by removing the redundancy in the classical overparameterized formulation. This is similar in spirit to the semidefinite lifting procedure in relaxations of combinatorial optimization problems. In contrast to these relaxations, we show that our lifting is always exact as soon as the network width exceeds a critical threshold which can be efficiently determined.
• We develop a matrix decomposition procedure called Neural Decomposition to extract the optimal network parameters from the solution of convex optimization, which is guaranteed to produce an optimal neural network. Neural Decomposition transforms the convex re-parameterization to the overparameterized, i.e., redundant, formulation in a similar spirit to (a non-orthogonal version of) Eigenvalue Decomposition.
• In addition to the fully connected neural network architecture, we derive the equivalent convex programs for various other architectures such as convolutional, pooling and vector output architectures.
• We provide extensive numerical simulations showing that the standard backpropagation approach with or without regularization fails to achieve the global optimum of the training loss. Moreover, the test accuracy of the proposed convex optimization is considerably higher in standard datasets as well as random planted models. Our convex optimization solver is significantly faster in total computation time to achieve similar or better test accuracy.   Quadratic activation scalar output with cubic regularization, Quad (scalar, quad reg): Quadratic activation scalar output with quadratic regularization. K is the number of patches and f is the filter size for the convolutional architecture. C is the output dimension for the vector output case. P is the pool size for average pooling. σ(u) is defined as u 2 for quadratic activation, and au 2 + bu + c for polynomial activation.

Prior Work
A considerable fraction of recent works on the analysis of optimization landscape of neural networks focuses on explaining why gradient descent performs well. The works [12,43] consider the optimization landscape of a restricted class of neural networks with quadratic activation and quadratic regularization where the second layer weights are fixed. They show that when the neural network is overparameterized, i.e., m ≥ d, the non-convex loss function has benign properties: all local minima are global and all saddle points have a direction of negative curvature. However, in this paper we show that training both the first and second layer weights with quadratic regularization in fact makes global optimization NP-hard. In contrast, we provide a different formulation to obtain the global optimal solution via convex optimization in the more general case when the second layer weights are also optimized, the activation function is any arbitrary degree two polynomial, and global optimum is achieved for all values of m. The work in [31] similarly studies two-layer neural networks with quadratic activation function and squared loss and states results on both optimization and generalization properties. The authors in [19] focus on quadratic activation networks from the perspectives of optimization landscape and generalization performance, where the setting is based on a planted model with a full rank weight matrix. In [26,29] it was shown that sufficiently wide ReLU networks have a benign landscape when each layer is sufficiently wide, satisfying m ≥ n + 1. Another recent work analyzing the training of neural networks with quadratic-like activations for deeper architectures is [2]. Authors in [2] consider polynomial activation functions and investigate layerwise training and compare with end-to-end training of layers. It is demonstrated in [2] that the degree two polynomial activation function performs comparably to ReLU activation in deep networks. More specifically, it is reported in [2] that for deep neural networks, ReLU activation achieves a classification accuracy of 0.96 and a degree two polynomial activation yields an accuracy of 0.95 on the Cifar-10 dataset. Similarly for the Cifar-100 dataset, they obtain an accuracy of 0.81 for ReLU activation and 0.76 for the degree two activation. These numerical results are obtained for the activation σ(u) = u+0.1u 2 , which the authors prefer over the standard quadratic activation σ(u) = u 2 to make the neural network training stable. Moreover, the performance of layerwise learning with such activation functions is considerably high, although there is a gap between end-to-end trained models. These results verify that degree two polynomial activations are promising and worth studying from both theoretical and practical perspectives.
In a recent series of papers, the authors derived convex formulations for training ReLU neural networks to global optimality [37,15,16,14,40,41]. Our work takes a similar convex duality approach in deriving the convex equivalents of non-convex neural network training problems. In particular, the previous work in this area deals with ReLU activations while in this work we focus on polynomial activations. Hence, the mathematical techniques involved in deriving the convex programs and the resulting convex programs are substantially different. The convex program derived for ReLU activation in [37] has polynomial time trainability for fixed rank data matrices, whereas the convex programs developed in this work are all polynomial-time trainable with respect to all problem dimensions. More specifically, their convex program is given by where the neural network weights are constructed from v i ∈ R d and w i ∈ R d , i = 1, . . . , P . The matrices D i are diagonal matrices whose diagonal entries consist of 1 x T 1 u≥0 , 1 x T 2 u≥0 , . . . , 1 x T n u≥0 for all possible u ∈ R d . The number of distinct D i matrices, denoted by P is the number of hyperplane arrangements corresponding to the data matrix X. It is known that P is bounded by 2r e(n−1) r r where r = rank(X) (see [37] for the details). In particular, convolutional neural networks have a fixed value of r, for instance m filters of size 3 × 3 yield r = 9. This is an exponential improvement over previously known methods that train optimal ReLU networks which are exponential in the number of neurons m and/or the number of samples n [3,20,5].
The work in [6] presents formulations for convex factorization machines with nuclear norm regularization, which is known to obtain low rank solutions. Vector output extension for factorization machines and polynomial networks, which are different from polynomial activation networks, is developed in [7]. Polynomial networks are equivalent to quadratic activation networks with an addition of a linear neuron. In [7], the authors consider learning an infinitely wide quadratic activation layer by a greedy algorithm. However, this algorithm does not provide optimal finite width networks even in the quadratic activation case. Furthermore, [30] presents a greedy algorithm for training polynomial networks. The algorithm provided in [30] is based on gradually adding neurons to the neural network to reduce the loss. More recently, [42] considers applying lifting for quadratic activation neural networks and presents non-convex algorithms for low rank matrix estimation for two-layer neural network training.

Notation
Throughout the text, σ : R → R denotes the activation function of the hidden layer. We refer to the function σ(u) = u 2 as quadratic activation and σ(u) = au 2 + bu + c where a, b, c ∈ R as polynomial activation. We use X ∈ R n×d to denote the data matrix, where its rows x i ∈ R d correspond to data samples and columns are the features. In the text, whenever we have a function mapping from R to R with a vector argument (e.g., σ(v) or v 2 where v is a vector), this means the elementwise application of that function to all the components of the vector v. We denote a column vector of ones by1 and its dimension can be understood from the context. vec(·) denotes the vectorized version of its argument. In writing optimization problems, we use min and max to refer to "minimize" and "maximize". We use the notations [m] and 1, . . . , m interchangeably.
We use (ŷ, y) for convex loss functions throughout the text for both scalar and vector outputs. * (v) = sup z (v T z − (z, y)) denotes the Fenchel conjugate of the function (·, y). Furthermore, we assume * * = which holds when is a convex and closed function [8].
We use Z 0 for positive semidefinite matrices (PSD). S refers to the set of symmetric matrices. tr refers to matrix trace. ⊗ is used for outer product. The operator conv stands for the convex hull of a set.

Preliminaries on Semidefinite Lifting
We defer the discussion of semidefinite lifting for two-layer neural networks with polynomial activations to Section 2. We now briefly discuss a class of problems where SDP relaxations lead to exact solutions of the original problem and also instances where they fail to be exact. Let us consider the following quadratic objective problem with a single quadratic constraint: where Q 1 , Q 2 are indefinite, i.e., not assumed to be positive semidefinite. Due to the indefinite quadratics, this is a non-convex optimization problem. By introducing a matrix variable U = uu T , one can equivalently state this problem as This problem can be relaxed by replacing the equality by the matrix inequality U uu T . Re-writing the expression U uu T as a linear matrix inequality via the Schur complement formula yields the following SDP min U,u Remarkably, it can be shown that the original non-convex problem in (5) can be solved exactly by solving the convex SDP in (7) via duality, under the mild assumption that the original problem is strictly feasible (see [8]). This shows that the SDP relaxation is exact in this problem, returning a globally optimal solution when one exists. We note that there are alternative numerical procedures to compute the global optimum of quadratic programs with one quadratic constraint [8].
We also note that the lifting approach U = uu T and the subsequent relaxation U uu T for quadratic programs with more than two quadratic constraints is not tight in general [33,9]. A notable case with multiple constraints is the NP-hard Max-Cut problem and its SDP relaxation [21] max The SDP relaxation of Max-Cut is not tight since its feasible set contains the cut polytope conv uu T : u i ∈ {−1, +1} ∀i and other non-integral extreme points [27]. Nevertheless, an approximation ratio of 0.878 can be obtained via the Goemans-Williamson randomized rounding procedure [21]. It is conjectured that this is the best approximation ratio for Max-Cut [24], whereas it can be formally proven to be NPhard to approximate within a factor of 16 17 [22,44]. Hence, in general we cannot expect to obtain exact solutions to problems of combinatorial nature, such as Max-Cut and variants using SDP relaxations.
It is instructive to note that a naive application of the SDP lifting strategy is not immediately tractable for two-layer neural networks. For simplicity, consider a scalar output polynomial activation network f (x) = m j=1 σ(x T u j )α j where σ(u) = u 2 + u, and {u j , α j } m j=1 are trainable parameters. The corresponding training problem for a given loss function (·, y) and its SDP relaxation are as follows The above problem is non-convex due to the bilinear terms {U j α j } m j=1 . Moreover, a variable changê U j = U j α j does not respect semidefinite constraints U j u j u T j when α j ∈ R. Another limitation is the prohibitively high number of variables in the lifted space, which is d 2 m + dm + m as opposed to dm + m in the original problem. Therefore, a different convex analytic formulation is required to address all these concerns.
Although SDP relaxations are extensively studied for various non-convex problems (see e.g. [45] for a survey of applications), instances with exact SDP relaxations are exceptionally rare. As will be discussed in the sequel, our main result for two-layer neural networks is another instance of an SDP relaxation leading to exact formulations where the semidefinite lifting and relaxation is tight.
In convex geometry, a spectrahedron is a convex body that can be represented as a linear matrix inequality which are the feasible sets of semidefinite programs. An example is the elliptope defined as the feasible set of the Max-Cut relaxation given by U 0, U ii = 1 ∀i, which is a subset of n × n symmetric positive-definite matrices. Due to the existence of efficient projection operators and barrier functions of linear matrix inequalities, optimizing convex objectives over spectrahedra can be efficiently implemented, which renders SDPs tractable. We will show that polynomial activation neural networks can be represented via a class of simple linear matrix inequalities, dubbed neural spectrahedra (see Figure 2 for an example), and enables global optimization in fully polynomial time and elucidates their parameterization in convex analytic terms.

Paper Organization
Section 2 gives an overview of the theory developed in this work. Section 3 describes the convex optimization formulation via duality and S-procedure for polynomial activation neural networks. Section 4 establishes via the neural decomposition method that the convex problem developed in Section 3 can be used to train two-layer polynomial activation networks to global optimality. Quadratic activation neural networks and the hardness result are studied in Section 5 and 6. Vector output and convolutional neural network architectures are studied in Section 7 and 8, respectively and convolutional networks with average pooling is in Section 9. We discuss the implementation details for solving the convex programs and give experimental results in Section 10.

Lifted Representations of Networks with Polynomial Activations
Consider the network f (x) = m j=1 σ(x T u j )α j where the activation function σ is the degree two polynomial σ(u) = au 2 + bu + c. First, we note that the neural network output can be written as (14)).
where φ : R d → R d 2 +d+1 and ψ : R m(d+1) → R d 2 +d+1 are formally defined in the sequel. The above identity shows that the nonlinear neural network output is linear over the lifted features In turn, the nonlinear model f (x) is completely characterized by the lifted parameters which we define as the following matrix-vector-scalar triplet Optimizing over the lifted parameter space initially appears as hard as the original non-convex neural network training problem. This is due to the cubic and quadratic terms involving the weights of the hidden and output layer in the lifted parameters. Nevertheless, one of our main results shows that the lifted parameters can be exactly described using linear matrix inequalities. We begin by characterizing the lifted parameter space as a non-convex cone.
Definition 1 (Neural Cone of degree two). We define the non-convex cone C m 2 ⊆ R d 2 +d+1 as See Figure 2 (left) for a depiction of C 1 2 ⊆ R 3 corresponding to the case m = 1, d = 1.
Surprisingly, we will show that the original non-convex neural network problem is solved exactly to global optimality when the optimization is performed over a convex set which we define as the Neural Spectrahedron, given by the convex hull of the cone C 2 . In other words, every element of the convex hull can be associated with a neural network of the form f (x) = m j=1 σ(x T u j )α j through a special matrix decomposition procedure which we introduce in Section 4. Moreover, a Neural Spectrahedron can be described by a simple linear matrix inequality. Consequently, these two results enable global optimization of neural networks with polynomial activations of degree two in fully polynomial time with respect to all problem parameters: dimension d, number of samples n and number of neurons m. To the best of our knowledge, this is the first instance of a method that globally optimizes a standard neural network architecture with computational complexity polynomial in all problem dimensions. We refer the reader to the recent work [37] for a convex optimization formulation of networks with ReLU activation, where the worst case computational complexity is O(( n r ) r ) with r = rank(X). It is equally important that our results characterize neural networks as constrained linear learning methods φ(x), ψ in the lifted feature space φ(x), where the constraints on the lifted parameters ψ are precisely described by a Neural Spectrahedron via linear matrix inequalities. These constraints can be easily tackled with convex semidefinite programming or closed-form projections onto these sets in iterative first-order algorithms. We also investigate interesting regularization properties of this convex set, and draw similarities to 1 norm and nuclear norm. In contrast, Reproducing Kernel Hilbert Space methods and Neural Tangent Kernel approximations [23,10] are linear learning methods over lifted feature maps where the corresponding parameter constraints are ellipsoids. These approximations fall short of explaining the extraordinary power of finite width neural networks employed in practical applications.
We extend the definition of the Neural Cone to degree k activations as follows.
Definition 2 (Neural Cone of degree k). We define the non-convex cone C m k ⊆ R k i=0 d i as follows where we use the notation u ⊗k := u ⊗ · · · ⊗ u k times .
It is easy to see that two-layer neural networks with degree k polynomial activations can be represented linearly using the lifted parameter space C k and corresponding lifted features. Taking the closure of the union {C} ∞ k=0 , any analytic activation function can be represented in this fashion. In this paper we limit the analysis to the degree 2 case.
Next, we describe a compact set that we call neural spectrahedron which describes the lifted parameter space of networks with a constraint on the 1 norm of output layer weights.
We will show that a neural spectrahedron can be equivalently described as a linear matrix inequality via defining S m Z, Z ∈ S (d+1)×(d+1) , Z 1 , Z 1 ∈ S d×d , Z 2 , Z 2 ∈ R d×1 and Z 4 , Z 4 ∈ R + , and m * = m * (t) is a critical number of neurons that satisfies m * (0) = 0 and m * (t) ≤ 2(d + 1) ∀t, which will be explicitly defined in the sequel. Therefore, an efficient description of the set M(t) in terms of linear matrix inequalities enables efficient convex optimization methods in polynomial time. Moreover, it should be noted that in non-convex optimization, the choice of the optimization algorithm and its internal hyperparameters, such as initialization, mini-batching and step sizes have a substantial contribution to the quality of the learned neural network model. This is in stark contrast to convex optimization problems, where optimizer hyperparameters have no effect, and solutions can be obtained in a very robust, efficient and reproducible manner.

A geometric description of the Neural Spectrahedron for the special case of nonnegative output layer weights
Here we describe a simpler case with the restriction α j ≥ 0 ∀j ∈ [m] in the Neural Cone C m 2 and we will suppose that m ≥ d + 1. In this special case, let us define the one-sided positive Neural Spectrahedron as We observe that Next, we note that as soon as the network width 1 satisfies m ≥ d + 1, we have where 0 is the zero matrix, since m j=1 is a positive semidefinite matrix, and hence can be factorized 2 as a convex combination of at most d+1 rank-one matrices of the Note that the zero matrix is included to account for the inequality m j=1 α j ≤ 1 in (16). This important observation enables us to represent the convex hull of the non-convex Neural Cone (an example is shown in Figure 2), via the simple convex body + M(t) given in (17).
Most importantly, the positive Neural Spectrahedron set + M(t) provides a representation of the non-convex Neural Cone C m 2 via its extreme points. Furthermore, + M(t) has a simple description as a linear matrix inequality provided in the following lemma (the proof can be found in the appendix).
Therefore the positive Neural Spectrahedron can be represented as the intersection of the positive semidefinite cone and linear inequalities. Moreover, every element of + M(t) can be factorized as m j=1 which can be identified as an element of the non-convex Neural Cone C m 2 and a neural network in the lifted parameter space as shown in (10).
The assumption m ≥ d + 1 is not required and only used here to illustrate this simpler special case. In the more general case of arbitrary output layer weights α j ∈ R, ∀j ∈ [m], we have the more general linear matrix inequality representation in (14), which is in terms of two positive semidefinite cones and three linear inequalities. In general, such a restriction on the number of neurons m in terms of the dimension d is not necessary. In the next sections, we only require m ≥ m * , where m * can be determined via a convex program. Furthermore, the regularization parameter directly controls the number of neurons m * . We illustrate the effect of the regularization parameter on m * in the numerical experiments section, and show that m * can be made arbitrarily small.

Convex Duality for Polynomial Activation Networks
We consider the non-convex training of a two-layer fully connected neural network with polynomial activation and derive a convex dual optimization problem. The input-output relation for this architecture is where σ is the degree two polynomial σ(u) = au 2 +bu+c. This neural network has m neurons with the first layer weights u j ∈ R d and second layer weights α j ∈ R. We refer to this case where f : R d → R as the scalar output case. Section 7 extends the results to the vector output case. It is relatively easy to obtain a weak dual that provides a lower-bound via Lagrangian duality. However, in non-convex problems, a duality gap may exist since strong duality does not hold in general. Remarkably, we show that strong duality holds as soon as the network width exceeds a critical threshold which can be easily determined.
We will assume 1 norm regularization on the second layer weights as regularization and include constraints that the first layer weights are unit norm. We note that 1 norm regularization on the second layer weights results in a special dual problem and hence is crucial in the derivations. We show in Section 5 that this formulation is equivalent to cubic regularization when the activation is quadratic. For the standard 2 2 , i.e., weight decay regularization, we will in fact show that the problem is NPhard (see Section 6). The training of a network under this setting requires solving the non-convex optimization problem given by Theorem 3.1 states the main result for polynomial activation neural networks that the non-convex optimization problem in (20) can be solved globally optimally via a convex problem. Before we state Theorem 3.1, we briefly describe the numerical examples shown in Figure 3 and 4 which compare the solution of the non-convex problem via backpropagation and the solution of the corresponding convex program via a convex solver (see Section 10 for details on the solver). Figure 3 shows the training and test costs on a regression task with randomly generated data for the two-layer quadratic activation neural network. We observe that convex SDP takes a much shorter time to optimize and obtains a globally optimal solution while the SGD algorithm converges to local minima in some of the trials where the initialization is different. Furthermore, Figure 4 compares the classification accuracies for the two-layer vector output polynomial activation network on a multiclass classification problem with real data. The exact statement of the vector output extension of the main result is provided in Section 7. In Section 10, we present additional numerical results verifying all of the theoretical results on various datasets. Figure 5 compares the accuracy of the non-convex polynomial activation model when it is trained with different optimizers (SGD and Adam) for a range of step sizes. Figure 5 shows that the convex formulations outperform the non-convex solution via SGD and Adam. The extension of the main result to convolutional neural networks is discussed in Section 8 and 9.
provides a global optimal solution for the non-convex problem in (20) when the number of neurons satisfies m ≥ m * where Here Z * and Z * denote the solution of (21). The variables Z ∈ S (d+1)×(d+1) and Z ∈ S (d+1)×(d+1) are defined in (29). It follows that the optimal number of neurons is upper bounded by m * ≤ 2(d + 1).
The proof of Theorem 3.1 is established in this section and the next. In this section we show that the solution of the convex program (21) provides a lower bound for the solution of the non-convex problem (20). In the next section, we prove, via the method of neural decomposition, that the solution of the convex problem provides also an upper bound, which concludes the proof of Theorem 3.1.
In proving the lower bound, we leverage duality. Minimizing over first α j 's and then u j 's, we can restate the problem in (20) as The dual problem for the inner minimization problem is given by Next, let us call the optimal solution of the following problem d * By changing the order of the minimization and maximization operations, we obtain the following bound We note that the constraints |v T σ(Xu j )| ≤ β can equivalently be written as two quadratic (in u j ) inequalities for each j = 1, . . . , m, Next, we use the S-procedure given in Corollary 3.3 to reformulate the quadratic inequality constraints as linear matrix inequality constraints. Corollary 3.3 is based on Lemma 3.2 which characterizes the solvability of a quadratic system. The proof of Corollary 3.3 is given in the appendix.
Lemma 3.2 (Proposition 3.1 from [38]). Let f 1 and f 2 be quadratic functions where f 2 is strictly concave (or strictly convex) and assume that f 2 takes both positive and negative values. Then, the following two statements are equivalent: Corollary 3.3 allows us to write the maximization problem in (26) as the equivalent problem given by where we note the two additional variables ρ 1 , ρ 2 ∈ R are introduced. Next, we will find the dual of the problem in (28). Let us first define the following Lagrange multipliers where Z, Z ∈ S (d+1)×(d+1) are symmetric matrices, and the dimensions for each block matrix are . We note that because of the symmetry of Z and Z , we have Z T 2 = Z 3 and Z 2 T = Z 3 . The Lagrangian for the problem in (28) is Maximizing the Lagrangian with respect to v, ρ 1 , ρ 2 , we obtain the problem in (21), which concludes the lower bound part of the proof. In the next section, we introduce a method for decomposing the solution of this convex program (i.e. Z * and Z * ) into feasible neural network weights to prove the upper bound.

Neural Decomposition
We have shown that a lower bound on the optimal value of the non-convex problem in (20) is obtained via the solution of the convex program in (21) that we have derived using Lagrangian duality. Now we show that this lower bound is in fact identical to the optimal value of the non-convex problem, thus proving strong duality. Our approach is based on proving an upper bound by constructing neural network weights from the solution of the convex problem such that the convex objective achieves the same objective as the non-convex objective. Suppose that (Z * , Z * ) is a solution to (21). Let us denote the rank of Z * by r and the rank of Z * by r . We will discuss the decomposition for Z * and then complete the picture by considering the same decomposition for Z * . We begin by noting that Z * satisfies the constraints of (21), i.e., Suppose that we have a decomposition of Z * as a sum of rank-1 matrices such that Z * = r j=1 p j p T j where p j ∈ R d+1 and tr(p j p T j G) = p T j Gp j = 0 for j = 1, . . . , r. We show how this can always be done in subsection 4.1 by introducing a new matrix decomposition method, dubbed the neural decomposition procedure. Letting We may assume p j = 0, ∀j in the decomposition (otherwise we can simply remove zero components), implying c j 2 2 > 0, ∀j. Furthermore, this expression for p j 's allows us to establish that As a result, we have the following decompositions: where we have introduced the normalized weights u j = cj cj 2 , j = 1, . . . , r. If d j ≤ 0 for some j, we redefine the corresponding p j as p j ← −p j , which does not modify the decomposition j p j p T j and the equality p T j Gp j = 0. Hence, without loss of generality, we can assume that d j ≥ 0 for all j = 1, . . . , r, which leads to Similarly for Z * , we will form the following decompositions: Considering the decompositions for both Z * and Z * , finally we obtain a neural network with first layer weights as {u 1 , . . . , u r , u 1 , . . . , u r }, and second layer weights as We note that this corresponds to a neural network with r + r neurons. If both Z * and Z * are full rank, then we will have 2(d + 1) neurons, which is the maximum.
To see why we can use the decompositions of Z * and Z * to construct neural network weights, we plug-in the expressions (36) and (37) in the objective of the convex program in (21): We note that this expression exactly matches the optimal value of the non-convex objective in (20) for a neural network with r + r neurons. Also, the unit norm constraints on the first layer weights are satisfied (hence feasible) since u j 's and u j 's are normalized. This establishes that the neural network weights obtained from the solution of the convex program provide an upper bound for the minimum value of the original non-convex problem. Consequently, we have shown that the optimal solution of the convex problem (21) provides a global optimal solution to the non-convex problem (20) and this concludes the proof of Theorem 3.1.

Neural Decomposition Procedure
Here we describe the procedure for computing the decomposition Z * = r j=1 p j p T j 0 such that p T j Gp j = 0, j = 1, . . . , r. This algorithm is inspired by the constructive proof of the S-procedure given in Lemma 2.4 of [38] with modifications to account for the equalities p T j Gp j = 0.
Neural Decomposition for Symmetric Matrices: 0. Compute a rank-1 decomposition Z * = r j=1 p j p T j . This can be done with the eigenvalue decomposition Z * = r j=1 q j q T j λ j . Since Z * 0, we have λ j > 0, for j = 1, . . . , r. Then we can obtain the desired rank-1 decomposition Z * = r j=1 p j p T j by defining p j = λ j q j , j = 1, . . . , r.
We know such j exists since tr(Z * G) = r j=1 p T j Gp j = 0 (this is true since it is one of the constraints of the convex program), and p T 1 Gp 1 = 0. Hence, for at least one j ∈ {2, . . . , r}, p T j Gp j must have the opposite sign as p T 1 Gp 1 .

Return y
We know that such α exists since the quadratic equation has real solutions since the discriminant 4(p T 1 p j ) 2 − 4(p T 1 Gp 1 )(p T j Gp j ) is positive due to step 1 where we picked j such that (p T 1 Gp 1 )(p T j Gp j ) < 0. To find α, we simply solve the quadratic equation for α.
3. Update r ← r − 1, and then the vectors p 1 , . . . , p r as follows: Remove p 1 and p j and insert u = pj −αp1 √ 1+α 2 . Consequently, we will be dealing with the updated matrix Z * ← Z * − yy T in the next iteration, which is of rank r − 1: Note that Step 0 is carried out only once and then steps 1 through 3 are repeated r − 1 times. At the end of r − 1 iterations, we are left with the rank-1 matrix p 1 p T 1 which satisfies p T 1 Gp 1 = 0 since initial Z * satisfies tr(Z * G) = 0 and the following r − 1 updates are of the form yy T which satisfies y T Gy = 0. If we denote the returned y vectors as y i for the iteration i and y r is the last one we are left with, then y i 's satisfy the desired decomposition that Z * = r i=1 y i y T i and y T i Gy i = 0, i = 1, . . . , r. Figure 6 is an illustration of the neural decomposition procedure for a toy example with d = 2 where the eigenvectors of Z * and the vectors p j are plotted together. Due to the constraints  Figure 6: Illustration of the neural decomposition procedure for d = 2 (i.e. Z * ∈ R 3×3 ). The dashed red arrows correspond to the eigenvectors of Z * (q 1 , q 2 , q 3 ) and the solid blue arrows show the decomposed vectors p 1 and p 2 . In this example, the rank of Z * is 2 where q 1 and q 2 are its two principal eigenvectors. The eigenvalue corresponding to the eigenvector q 1 is zero. The light blue colored surface shows the Lorentz cones z = x 2 + y 2 and z = − x 2 + y 2 . We observe that the decomposed vectors p 1 and p 2 lie on the boundary of Lorentz cones.
p T j Gp j = 0, j = 1, 2, the vectors p j have to lie on the boundary of Lorentz cones 3 z = x 2 + y 2 and z = − x 2 + y 2 . Decomposing the solution of the convex problem Z * and Z * onto these cones, i.e., neural decomposition, enables the construction of neural network weights from Z * and Z * .

Quadratic Activation Networks
In this section, we derive the corresponding convex program when the activation function is quadratic, i.e., σ(u) = u 2 . The resulting convex problem takes a simpler form than the polynomial activation case. We start by noting that the bound in (26) holds for any activation function. The inequalities |v T σ(Xu j )| ≤ β however lead to different constraints than the polynomial activation case. Note that |v T (Xu j ) 2 | ≤ β is equivalent to the inequalities The constraint max uj : uj 2=1 |v T (Xu j ) 2 | ≤ β can be expressed as largest eigenvalue inequalities where λ max denotes the maximum eigenvalue. Next, representing the largest eigenvalue constraints as linear matrix inequality constraints, we arrive at the following maximization problem Writing the Lagrangian for (43) as with Z 1 , Z 2 ∈ S d×d and maximizing with respect to v, we obtain the following convex problem min Z1,Z2 0 Replacing Z = Z 1 −Z 2 , where Z 1 0, Z 2 0, we recall that any matrix Z can be uniquely decomposed in this form thanks to the Moreau decomposition onto the cone of positive definite matrices and its polar dual, which is the set of negative semidefinite matrices. In particular, suppose that the eigenvalue decomposition of Z is Z = j λ j z j z T j . Then, Z 1 and Z 2 are uniquely determined by Z 1 = j:λj >0 λ j z j z T j and Z 2 = − j:λj <0 λ j z j z T j . Note that tr(Z 1 + Z 2 ) = j:λj >0 λ j + j:λj <0 (−λ j ) = j |λ j | = Z * is the sum of the absolute values of the eigenvalues of Z, which is equivalent to the nuclear norm for symmetric matrices. Consequently, this leads to the following simplified problem with nuclear norm regularization: The optimal neural network weights are determined from the solution of the convex problem via eigenvalue decomposition of Z * and the rescaling given in (51). The optimal number of neurons is upper bounded by m * ≤ d since rank(Z * ) ≤ d.

Strong Duality for Quadratic Activation
We have shown that a lower bound on the non-convex problem for quadratic activation is given by the nuclear norm regularized convex objective. Now we show that this lower bound is in fact identical to the non-convex problem. Suppose that Z * is a solution to (45). Let us decompose Z * via eigenvalue decomposition as Z * = j λ j z j z T j . We can generate an upper bound on the non-convex problem by constructing neural network parameters as α j = λ j , and u j = z j with objective value j (Xz j ) 2 λ j , y +β j |λ j |. Noting that this value exactly matches the optimal value of the convex objective in (45), we conclude that the optimal solution of (45) provides a global optimal solution to the non-convex problem.

Equivalent Non-convex Problem: Quadratic Activation with Cubic Regularization
We now show that the non-convex problem with unit norm first layer weights and the 1 norm regularized second layer weights is in fact equivalent to the non-convex problem with cubic regularization on all the weights. Let us consider the unconstrained problem with cubic regularization: where c = 2 Noting the regularization term is convex in t j for t j > 0 and optimizing it with respect to t j , we obtain t j = 2 2/9 |αj | uj 2

2/3
. Plugging the expression for t j in yields Now we define the scaled second layer weights α j = α j u j 2 2 . Noting that (Xu j ) 2 α j = (X uj uj 2 ) 2 α j and defining u j = u j / u j 2 , we obtain the equivalent problem with the 1 norm of the second layer weights as the regularization term

Rescaling
We note that the weights α j and u j that the eigenvalue decomposition of the solution of (45) gives are scaled versions of the weights of the problem with cubic regularization in (47). The solution to the problem in (47) can be constructed by rescaling the weights as This concludes the proof of Theorem 5.1.

Comparison with Polynomial Activation Networks
In this subsection, we list the important differences between the results for quadratic activation and polynomial activation neural networks. The convex program for the quadratic activation network does not have the equality constraints that appear in the convex program for the polynomial activation.
In addition, for the quadratic activation, the upper bound on the critical width m * is d while it is 2(d + 1) for the polynomial activation case. We note that in the case of quadratic activation, the optimal neural network weights are determined from eigenvalue decomposition of Z * . This results in the first layer weights to be orthonormal because they can be chosen as the eigenvectors of the real and symmetric matrix Z * . In contrast, we do not have this property for polynomial activations as the associated optimal weights are determined via neural decomposition. In this case, the resulting hidden neurons are not necessarily orthogonal, which shows that the Neural Decomposition is a type of non-orthogonal matrix decomposition. This can also be seen in Figure 6.

Constructing Multiple Globally Optimal Solutions in the Neural Network Parameter Space
Once we find an optimal Z * using the SDP in (45), we can transform it to the neural network parameter space with at most d neurons using the eigenvalue decomposition of Z * as Z * = d j=1 u j u T j α j . However, we can also generate a neural network with an arbitrary number of neurons, which is also optimal. We now describe this construction below for an arbitrary number of neurons m ≥ 2d. Let us pick an arbitrary m/2 × d matrix H with orthonormal columns, i.e., where h 1 , . . . , h m/2 are the rows of H and we assume m/2 ≥ d. One can generate such matrices using randomized Haar ensemble, or partial Hadamard matrices. Then, we can represent Z * using Since Z * is a symmetric matrix, m/2 j=1 Z * h j h T j is also symmetric, and we can write Finally, for each term in the above summation, we employ the symmetrization identity valid for any x, y ∈ R d . We arrive at the representation where u j = Z * h j + h j , α j = 1/4 for j = 1, . . . , m/2 and u j = Z * h j − h j , α j = −1/4 for j = m/2 + 1, . . . , m.
Since the matrix H is arbitrary, one can map an optimal Z * matrix from the convex semidefinite program to infinitely many optimal solutions in the neural network parameterization space.

Standard Weight Decay Formulation is NP-Hard
In Section 5, we have studied two-layer neural networks with quadratic activation and cubic regularization and derive a convex program whose solution globally optimizes the non-convex problem. In this section, we show that if, instead of cubic regularization, we have quadratic regularization (i.e. weight decay), the resulting optimization problem is an NP-hard problem.
The remainder of this section breaks down the proof of Theorem 6.1. At the core of the proof is the polynomial-time reduction of the problem to the NP-hard problem of phase retrieval.

Reduction to an Equivalent Problem
The optimization problem for training a two-layer fully connected neural network with quadratic activation and quadratic regularization can be stated as where the scaling factor c is the same as before (i.e. c = 2 1 3 + 2 − 2 3 ≈ 1.88988). Rescaling u j ← u j t 1/2 j and α j ← α j /t j for t j > 0, j = 1, . . . , m, we obtain the following equivalent optimization problem Note that the regularization term is convex in t j for t j > 0. Optimizing the regularization term with respect to t j leads to t j = 2 1/3 |αj | Defining scaled weights α j = α j u j 2 2 and u j = u j / u j 2 , we obtain the equivalent problem This shows that solving the standard weight decay formulation is equivalent to solving a 2/3-norm penalized problem with unit norm first layer weights.

Hardness Result
We design a data matrix such that the solution coincides with solving the phase retrieval problem which is NP-hard (see [18]). We consider the equality constrained version of (58), i.e., β → 0, which is given by (59)

Addition of a Simplex Constraint
Let the first d rows of the data matrix X be e T 1 , . . . , e T d and let the first d entries of y be 1/d. Then, the constraint m j=1 (Xu j ) 2 = y implies m j=1 u 2 jk α j = 1/d for k = 1, . . . , d .
Summing the above for all k = 1, . . . , d, and noting that d k=1 u 2 jk = 1 lead to the constraint m j=1 α j = 1.

Reduction to the NP-Hard Phase Retrieval and Subset Sum Problem
We let X = [I;X] and y = [ 1 d1 ;ỹ] to obtain the simplex constraint m j=1 α j = 1 as shown in the previous subsection. In this case, the optimization problem reduces to Suppose that there exists a feasible solution {α * j , u * j } m j=1 , which satisfies α * 0 = 1, where α * 1 = 1 and u * 1 T u * 1 = 1 with only one nonzero neuron. Then, it follows from Lemma 6.2 that this solution is strictly optimal. Consequently, the problem in (61) is equivalent to Lemma 6.2 ( p minimization recovers 1-sparse solutions when 0 < p < 1).
Consider the optimization problem min α1,...,αm where C is a convex set and p ∈ (0, 1). Suppose that there exists a feasible solution α * ∈ C and i α * i = 1 such that α * 0 = 1. Then, α * is strictly optimal with objective value 1. More precisely, any solution with cardinality strictly greater than 1 has objective value strictly larger than 1.

NP-hardness Proof
Subset sum problem given in Definition 4 is a decision problem known to be NP-complete (e.g. [18]). The decision version of the problem in (62) can be stated as follows: Does there exist a feasible u 1 ?
We show that this decision problem is NP-hard via a polynomial-time reduction to the subset sum problem.
Definition 4 (Subset sum problem). Given a set of integers A, does there exist a subset A S whose elements sum to z? Lemma 6.3 establishes the reduction of the decision version of (62) to the subset sum problem. The proof is provided in the appendix and follows the same approach used in the proof for the NPhardness of phase retrieval in [18], with the main difference being the additional constraints u 2 1k = 1/d, k = 1, . . . , d in (62). Finally, Lemma 6.3 concludes the proof of Theorem 6.1. only if the answer for the subset sum problem with A = {a 1 , . . . , a d } is 'yes'.
Remark 6.1. It follows from Theorem 6.1 that the two-layer neural network training problem with polynomial activation and unit norm first layer weights and j |α j | p as the regularization term with p < 1 is also NP-hard for β → 0 since it reduces to the quadratic activation case for the polynomial coefficients a = 1, b = 0, c = 0.

Vector Output Networks
The derivations until this point have been for neural network architectures with scalar outputs, i.e., y i ∈ R. In this section, we turn to the vector output case y i ∈ R C where C is the output dimension, and derive a convex problem that has the same optimal value as the non-convex neural network optimization problem. We exploit the same techniques described in the scalar output case except for the part for constructing the vector second layer weights from the solution of the convex program.
In the scalar output case, the convex problem is over the symmetric matrices Z, Z and in the vector output case, the optimization is over C such matrix pairs Z k , Z k , k = 1, . . . , C.
We begin our treatment of the vector output case by considering the neural network defined by where α j ∈ R C , j = 1, . . . , m are the vector second layer weights. Note that in the scalar output case, the second layer weights α j were scalars. Taking the regularization to be the 1 norm of the second layer weights, the neural network training requires solving the following non-convex optimization problem where Y ∈ R n×C is the output matrix. Equivalently, The dual problem for the inner minimization problem is given by where v ∈ R n×C is the dual variable and v k ∈ R n is the k'th column of v. Theorem 7.1 gives the main result of this section.
The optimal neural network weights are determined from the solution of the convex problem via the neural decomposition procedure for each Z * k and Z k * and the construction given in (75). The optimal number of neurons is upper bounded by m * ≤ 2(d + 1)C.
Proof of Theorem 7.1. Applying the S-procedure for the constraints in the dual problem (67), we obtain the following maximization problem Next, let us introduce the following Lagrange multipliers Then, the Lagrangian is Finally maximizing the Lagrangian leads to the following convex SDP: We construct the neural network weights from the optimal solution of the convex program as follows. We follow the neural decomposition procedure from Section 4 for extracting neurons from each of the matrices Z * k and Z k * , k = 1, . . . , C. The decompositions for Z * k will be of the form Then, the weights due to Z * k , k = 1, . . . , C are determined as follows: First layer weights: {u 1,1 , u 1,2 , . . . , u 1,r1 }, . . . , {u C,1 , u C,2 , . . . , u C,r C } Second layer weights: where e k denotes the k'th C-dimensional unit vector, and r k is the rank of the matrix Z * k . In short, the matrix Z * k with rank r k leads to the first layer weights {u k,1 , u k,2 , . . . , u k,r k } and the second layer weights {d 2 k,1 e T k , d 2 k,2 e T k , . . . , d 2 k,r k e T k }. The weights due to Z k * , k = 1, . . . , C are determined the same way. Then, we reach the following neural network construction: Finally, the total number of neurons that the convex problem finds is C k=1 (r k +r k ). The maximum number of neurons occurs if all Z * k and Z k * are full rank, and this corresponds to a maximum total of 2(d + 1)C neurons.
We plug the decomposition expressions given in (73) in the convex program in (72) to conclude that the optimal value of the convex program is an upper bound for the non-convex optimization problem (65). The k'th entry of the estimate for the i'th training sample iŝ It follows that the output vector for the i'th sample iŝ We note that this output is of the same form as the non-convex case (66). We also need to check that the regularization term is equivalent to the sum of 1 norms of the second layer weights: which is of the form β m j=1 α j 1 . Hence, the neural network weights that we obtain via the neural decomposition procedure lead to an upper bound for the original non-convex optimization problem. This concludes the proof that the optimal solution of the convex problem (72) provides a global optimal solution to the non-convex problem (65).

Convolutional Neural Networks
In this section, we consider two-layer convolutional networks with a convolutional first layer and a fully connected second layer. We will denote the filter size by f . Let us denote the patches of a data sample x by x 1 , . . . , x K where the patches have the same dimension as the filters, i.e., x k ∈ R f . The stride and padding do not affect the below derivations as they can be readily handled when forming the patches. The output of this network is expressed as: where u j ∈ R f denotes the j'th filter. We will take the regularization to be the 1 norm of the second layer weights α j = α j1 . . . α jK T ∈ R K , j = 1, . . . , m: where we use X k ∈ R n×f to denote the matrix with the k'th patch of all the data samples. The dual for the inner minimization problem is given by We state the main result of this section in Theorem 8.1.
Theorem 8.1 (Globally optimal convex program for polynomial activation convolutional neural networks). The solution of the convex problem in (85) provides a global optimal solution for the nonconvex convolutional neural network problem in (79) when the number of filters is at least (rank(Z * k ) + rank(Z k * )) and equivalently, the number of neurons satisfies m ≥ m * where The optimal neural network weights are determined from the solution of the convex problem via the neural decomposition procedure for each Z * k and Z k * . The optimal number of filters is upper bounded by 2(f + 1)K and the optimal number of neurons is upper bounded by m * ≤ 2(f + 1)K 2 .
Proof of Theorem 8.1. We apply the S-procedure to replace the constraints of (81) with equivalent LMI constraints and this yields where x i,k ∈ R f denotes the k'th patch of the i'th data sample. The Lagrangian is as follows where Z k , Z k are (f + 1) × (f + 1) dimensional symmetric matrices. Maximizing the Lagrangian with respect to v, ρ k,1 , ρ k,2 , k = 1, . . . , K yields the convex SDP tr(Z k,1 ) = Z k,4 , tr(Z k,1 ) = Z k,4 , k = 1, . . . , K We now show that the convex program in (85) provides an upper bound for the non-convex problem via the same strategy that we have used for the vector output case in Section 7. We construct the neural network weights from each of the matrices Z * k and Z k * , k = 1, . . . , K via neural decomposition: and the weights due to each Z * k are First layer filters: u k,1 , u k,2 , . . . , u k,r k Second layer weights: To clarify, for each filter u k,j , we have K (scalar) weights in the second layer because we apply the same filter to K different patches and the resulting K numbers (after being input to the activation function) each are multiplied by a different second layer weight. The second layer weights associated with the filter u k,j will be these K numbers: {0, . . . , 0, d 2 k,j , 0 . . . , 0}, where the only nonzero entry is the j'th one. Consequently, each Z * k matrix produces rank(Z * k ) filters and K rank(Z * k ) neurons. Including the weights due to Z k * as well, we will have K k=1 (r k + r k ) filters and K K k=1 (r k + r k ) neurons in total. The optimal number of filters is upper bounded by 2(f + 1)K and the optimal number of neurons is upper bounded by 2(f + 1)K 2 .
We omit the details of plugging the weights into the convex objective to show that it becomes equivalent to the non-convex objective. The details are similar to the vector output case.

Average Pooling
In this section we will consider convolutional neural networks with average pooling. We will denote the pool size by P . Let us consider a two-layer neural network where the first layer is a convolutional layer with filter size f . The convolutional layer is followed by the polynomial activation, average pooling, and a fully connected layer. We will denote the number of patches per sample by K. The output of this architecture can be expressed as We note that the number of parameters in the second layer (i.e. α jk 's) is equal to m K P . The optimization problem for this architecture can be written as σ(X (k−1)P +l u j ) α jk , where α j = α j1 . . . α j,K/P T , j = 1, . . . , m. The dual of the inner minimization problem is given Theorem 9.1 states our result for CNN with average pooling.
Theorem 9.1 (Globally optimal convex program for polynomial activation convolutional neural networks with average pooling). The solution of the convex problem in (95) provides a global optimal solution for the non-convex problem for the convolutional neural network with average pooling in (89) when the number of neurons satisfies m ≥ m * where The optimal neural network weights are determined from the solution of the convex problem via the neural decomposition procedure for each Z * k and Z k * . The optimal number of neurons is upper bounded by m * ≤ 2(f + 1) K 2 P 2 . Proof of Theorem 9.1. We rewrite the constraints of the dual problem (90) as follows: S-procedure allows us to write this problem equivalently as The Lagrangian is as follows where Z k , Z k are (f + 1) × (f + 1) dimensional symmetric matrices. Maximizing the Lagrangian with respect to v, ρ k,1 , ρ k,2 , k = 1, . . . , K/P yields the following convex SDP: tr(Z k,1 ) = Z k,4 , tr(Z k,1 ) = Z k,4 , k = 1, . . . , K/P Z k 0, Z k 0, k = 1, . . . , K/P .
We omit the details of constructing the neural network weights from the solution of the convex SDP Z * k , Z k * , k = 1, . . . , K/P which follows in a similar fashion as the proof of Theorem 8.1.
We note that when we pick the pool size as P = 1, this is the same as not having average pooling, and the corresponding convex program is the same as (85), derived in Section 8. The other extreme for the pool size is when P = K and this corresponds to what is known as global average pooling in which case the convex SDP simplifies to We note that the problem (96) has only two variables Z and Z . This should be contrasted with the convolutional architecture with no pooling (85) which has 2K variables.

Numerical Results
In this section, we present numerical results that verify the presented theory of the convex formulations along with experiments comparing the test set performance of the derived formulations. All experiments have been run on a MacBook Pro with 16GB RAM. Solvers: We have used CVXPY [11,1] for solving the convex SDPs. In particular, we have used the open source solver SCS (splitting conic solver) [34,35] in CVXPY, which is a scalable first order solver for convex cone problems.
Furthermore, we have solved the non-convex problems via backpropagation for which we have used PyTorch [36]. We have used the SGD algorithm for the non-convex models. For all the experiments involving SGD in this section, we show only the results corresponding to the best learning rate that we select via an offline hyperparameter search. The momentum parameter is 0.9. In the plots, the non-convex models are either labeled as 'Backpropagation (GD)' or 'Backpropagation (SGD)'. The first one, short for gradient descent, means that the batch size is equal to the number of samples n, and the second one, short for stochastic gradient descent, means that the batch size is not n and the exact batch size is explicitly stated in the figure captions.
Polynomial approximation of activation functions: To obtain the degree-2 polynomial approximation of a given activation function σ(u) such as the ReLU activation, one way is to select the polynomial coefficients a, b, c that minimize the 2 norm objective T a b c where t i 's are linearly spaced in [L, U ]. The lower and upper limits L and U specify the range in which we would like to approximate the given activation function. For instance, when L = −5, U = 5, N = 1000 and σ(u) is the ReLU activation, the optimal polynomial coefficients are a = 0.09, b = 0.5, c = 0.47. When we change the approximation range to a slightly narrower one with L = −4, U = 4, the coefficients then become a = 0.12, b = 0.5, c = 0.38. Note that the training data can be normalized appropriately to confine the range of the input to the neurons and control the approximation error.

Results for Verifying the Theoretical Formulations
The first set of numerical results in Figure 7 is for verifying that the derived convex problems have the same optimal value as their non-convex counterparts. The plots in Figure 7 show the non-convex cost against time when 1) the non-convex problem is solved in PyTorch and 2) the corresponding convex problem (see Table 1) is solved using CVXPY. The number of neurons for the non-convex models in all of the plots in Figure 7 is set to the optimal number of neurons m * found by the convex problem. Figure 7 demonstrates that solving the convex SDP takes less time than solving the associated non-convex problem using backpropagation for all of the neural network architectures. Figure 7 also shows that the training of the non-convex models via the backpropagation algorithm does not always yield the global optimal but instead may converge to local minima. In addition, we note that the plots do not reflect the time it takes to tune the learning rate for the non-convex models, which was performed offline.

Experiments on UCI datasets
We now show how the derived convex programs perform in the context of classification datasets. The datasets used in this subsection are from the UCI machine learning repository [13]. The plots in Figure 8 show the training and test set costs and classification accuracies for binary classification datasets and the plots in Figure 9 are for multiclass classification datasets. The convex program used for solving the binary classification problem is the scalar output polynomial activation problem given in (21) and for the multiclass problem it is the vector output version given in (72).
We note that the training cost plots of Figure 8 and 9 are consistent with the theoretical results. The accuracy plots show that the convex programs achieve the same final accuracy of the non-convex models or higher accuracies in shorter amounts of time. Table 2 shows the classification accuracies of various fully connected neural network architectures on binary classification UCI datasets. For each dataset, the training and validation partitions are as pre-processed in [17]. The training and validation partitions are used to select the best hyperparameters. The hyperparameter search for the non-convex models includes searching for the best regularization coefficient β and learning rate. Gradient descent has been used to optimize the nonconvex models and the number of epochs is 1000. After determining the best hyperparameters, we compute the 4-fold cross validation accuracy and report it in this table. The partitions for the 4-fold cross validation are also the same as those pre-processed by [17]. Furthermore, for the results shown in Table 2, the number of neurons for all the non-convex models is set to 2(d + 1), which is the maximum number of neurons that the polynomial activation convex SDP could output (see Theorem 3.1). Table  2 shows that the convex SDP achieves better or similar accuracy values compared to the non-convex models on most of the datasets. (a) Quad act (10,20,9)

Comparison with ReLU Networks
We compare the classification accuracies for polynomial activation and ReLU activation in Figure 10 on three different binary classification UCI datasets. The regularization coefficient has been picked separately for polynomial activation and ReLU activation networks to maximize the accuracy. Figure  10 demonstrates that the convex SDP shows competitive accuracy performance and faster run times compared to ReLU activation networks. Figure 11 shows the binary classification accuracy performance of the CNN architecture with global average pooling on MNIST [28], Fashion MNIST [46], and Cifar-10 [25] datasets. Figure 11 compares the non-convex tractable problem, the corresponding convex formulation, and the non-convex weight decay formulation. By the weight decay formulation, we mean quadratic regularization on both the first layer filters and the second layer weights. We observe that the accuracy of the convex SDP is slightly better or the same as SGD while the run time for the convex SDP solution is consistently shorter than the time it takes for SGD to converge. Figure 12 shows how the accuracy changes as a function of the regularization coefficient β for the convex problem for two-layer polynomial activation networks. Figure 12 highlights that the choice of the regularization coefficient is critical in the accuracy performance. In plot a, we see that the value of β that maximizes the test set accuracy is β = 10 for which the optimal number of neurons m * is near 20. We note that for the dataset in plot a, the optimal number of neurons is upper bounded by m * ≤ 2(d + 1) = 32. Similarly for plot b, the best choice for the regularization coefficient is β = 1 and the optimal number of neurons for β = 1 is near 40. Furthermore, we observe that a higher value for β tends to translate to a lower optimal number of neurons m * (plotted on the right vertical axis). Even though the convex optimization problem in (21) has a fixed number of variables (in this case, 2(d + 1) 2 ) for a given dataset, a low number of neurons is still preferable for many reasons such as inference speed. We observe that the number of neurons can be controlled via the regularization coefficient β.

Other Losses
We have so far evaluated the performance of the derived convex programs for squared loss, i.e. (ŷ, y) = ŷ − y 2 2 . We reiterate that the derived convex programs are general in the sense the formulations hold for any convex loss function . To verify this numerically, we now present results for additional loss functions such as Huber loss and 1 norm loss in Figure 13. More concretely, Huber loss is defined as (ŷ, y) = n i=1 Huber(ŷ i − y i ) where Huber(x) = 2|x| − 1 for |x| > 1 and Huber(x) = x 2 for |x| ≤ 1. The 1 norm loss is (ŷ, y) = ŷ − y 1 . We observe that in the case of 1 norm loss, backpropagation takes longer to converge.

The Effect of Polynomial Coefficients
The plots in Figure 14 show the classification accuracy against the polynomial coefficients a, b, c for the polynomial activation convex problem. In each plot, we vary one of the coefficients and fix the other two coefficients as 1. We observe that the coefficient of the quadratic term a plays the most important role in the accuracy performance. The accuracy is not affected by the choice of the coefficient c.

Discussion
In this paper, we have studied the optimization of two-layer neural networks with degree two polynomial activations. We have shown that regularization plays an important role in the tractability of the problems associated with neural network training. We have developed convex programs for the cases where the regularization leads to tractable formulations. Convex formulations are useful since they have many well-known advantages over non-convex optimization such as having to optimize fewer hyperparameters and no risk of getting stuck at local minima. The methods presented in this work optimize the neural network parameters in a higher dimensional space in which the problem becomes convex. For fully connected neural networks with quadratic activation, the standard non-convex problem requires optimizing m neurons (i.e. a d-dimensional first layer weight and a 1-dimensional second layer weight per neuron). The convex program for this neural network finds the optimal network parameters in the lifted space S d×d . For polynomial activations, convex optimization takes place for Z and Z in S (d+1)×(d+1) . We note that the dimensions of the convex programs are polynomial with respect to all problem dimensions. In contrast, the convex program of [37] has 2dP variables where P grows exponentially with respect to the rank of the data matrix.
We have used the SCS solver with CVXPY for solving the convex problems in the numerical experiments. It is important to note that there is room for future work in terms of which solvers to use. Solvers specifically designed for the presented convex programs could enjoy faster run times.
The scope of this work is limited to two-layer neural networks. We note that it is a promising direction to consider the use of our convex programs for two-layer neural networks as building blocks in learning deep neural networks. Many recent works such as [2] and [4] investigate layerwise learning algorithms for deep neural networks. The training of individual layers in layerwise learning could be improved by the presented convex programs since the convex programs can be efficiently solved and eliminate much of the hyperparameter tuning involved in standard neural network training.

A.1 Constrained Least Squares Form for the Squared Loss
Let us consider the polynomial activation scalar output case. In the case of squared loss (ŷ, y) = ŷ − y 2 2 , the convex program takes the following form: Noting that ax T i (Z 1 − Z 1 )x i = vec(x i x T i ) T vec(Z 1 − Z 1 ), we can write the squared loss term as where we have defined X V ∈ R n×(d 2 +d+1) and z ∈ R (d 2 +d+1) . The squared loss term is equal to z T X T V X V z − 2y T X V z + y 2 2 . If we pre-compute X T V X V ∈ R (d 2 +d+1)×(d 2 +d+1) and X T V y ∈ R (d 2 +d+1) , then the objective no longer has dependence on the number of samples n. We note that the pre-computation of X T V X V and X T V y is useful when one is performing hyperparameter tuning for the regularization coefficient β.

B Proofs
Proof of Lemma 2.1. We will denote the set in (17) as S 1 and the set in (18) as S 2 to simplify the notation. We will prove S 1 = S 2 by showing S 1 ⊆ S 2 and S 2 ⊆ S 1 . We first show S 1 ⊆ S 2 . Let us take a point S ∈ S 1 . This implies that S is a matrix of the form t m j=1 u j 1 with j α j ≤ 1 and u j 2 = 1 for all j. We note that tr(t m j=1 u j u T j α j ) = t m j=1 tr(u j u T j α j ) = t m j=1 tr(u T j u j )α j = t m j=1 α j ≤ t. This shows that S satisfies the equality condition in the definition (18). Now, we show that S is a PSD matrix. Note that each of the rank-1 matrices u j 1 u j 1 T is a PSD matrix and since the coefficients α j 's and t are nonnegative, it follows that S is PSD. This proves that S ∈ S 2 . We next show S 2 ⊆ S 1 . Let us take a point S ∈ S 2 . This implies that S is PSD and tr(S 1 ) = S 4 = t 0 ≤ t. We show in Section 4 that it is possible to decompose S via the neural decomposition procedure to obtain the expressions given in (36). It follows that we can write S in the following form where the scaling factor t0 m j=1 d 2 j is to ensure that tr(S 1 ) = S 4 = t 0 ≤ t. It is obvious to see that S is in S 1 when t 0 = t by the definition of S 1 given in (17). When t 0 < t, we still have that S is in S 1 which can be seen by noting that S 1 is defined as the convex hull of rank-1 matrices and the zero matrix. We can scale all the rank-1 matrices in the convex combination with t0 t and change the weight of the zero matrix accordingly.
Proof of Lemma 6.3. Let us define the set A = {a 1 , a 2 , . . . , a d } where a i are integers. We need to show that the problem (62) finds a feasible solution u 1 if and only if there exists a subset A S of the set A that satisfies a∈A S a = z. We assume n = 2d + 1 and henceX is (d + 1) × d andỹ is (d + 1) dimensional. LetX D ∈ R d×d denote the matrix with the first d rows ofX, andx d+1 is the last sample inX. Let us defineỹ i as where w =X −T Dx d+1 ∈ R d . Direction 1: Suppose there exists u 1 ∈ R d such that (x T i u 1 ) 2 =ỹ i for every i = 1, . . . , d + 1 and u 2 1k = 1/d for every k = 1, . . . , d. Then there exists a subset A S with a∈A S a = z. Proof of direction 1: AssumingX D is invertible, it follows thatXX −1 D = I d w T where I d is the d × d identity matrix. Let us consider a feasible u 1 . Then, v =X D u 1 satisfies ((X −T Dx i ) T v) 2 =ỹ i for i = 1, . . . , d + 1. Consequently, we haveX −T Dx i = e i for i = 1, . . . , d, andX −T Dx d+1 = w. As a result, we obtain the following relation between v andỹ: Next, because of (102), we have |v i | = a i w i , i = 1, . . . , d, and |w T v| = 2z − d j=1 a j .
Let us define ε i such that v i = ε i a i /w i for i = 1, . . . , d. Note that ε i ∈ {−1, 1}. Then, This means we either have z = i:εi=1 a i or z = − i:εi=1 a i + d i=1 a i = i:εi=−1 a i . This shows that the sum of the elements of A S is equal to z when A S is either equal to {a i |ε i = 1} or {a i |ε i = −1}.
In proving direction 2, it is straightforward to show the existence of u 1 that satisfies the constraint (x T i u 1 ) 2 =ỹ i . To show that there is a u 1 that satisfies the constraint u 2 1k = 1 d , we pickX in a certain way that we discuss now: To prove direction 2, we will need to make sure |X −1 We pickX D to be any diagonal matrix with arbitrary −1's and +1's on the diagonal and pick x d+1 = √ d a 1 . . . a d T . Since w =X −T Dx d+1 , we will have |w i | = |x d+1,i | = √ d|a i | for i = 1, . . . , d. This choice forX D andx d+1 ensures that |X −1 D v| =1 1 √ d . Direction 2: Suppose there is a subset A S with a∈A S a = z. Then there exists a feasible u 1 ∈ R d .
Proof of direction 2: Define ε i such that for a i in A S , it is equal to 1, and otherwise it is equal to −1. Next,