Convex Formulations for Training Two-Layer ReLU Neural Networks

## Abstract

Abstract Solving non-convex, NP-hard optimization problems is crucial for training machine learning models, including neural networks. However, non-convexity often leads to black-box machine learning models with unclear inner workings. While convex formulations have been used for verifying neural network robustness, their application to training neural networks remains less explored. In response to this challenge, we reformulate the problem of training infinite-width two-layer ReLU networks as a convex completely positive program in a finite-dimensional (lifted) space. Despite the convexity, solving this problem remains NP-hard due to the complete positivity constraint. To overcome this challenge, we introduce a semidefinite relaxation that can be solved in polynomial time. We then experimentally evaluate the tightness of this relaxation, demonstrating its competitive performance in test accuracy across a range of classification tasks.

## 1 Introduction

The outstanding performance of deep neural networks has driven significant changes in the research directions of optimization for machine learning over the past decade, replacing traditional convex optimization techniques with non-convex methods.
Nonetheless, non-convexity introduces significant challenges to the analysis, resulting in the use of models that lack a comprehensive explanation of their inner workings. While convex formulations have been studied for enhancing reliability through applications like estimating the Lipschitz constant of neural networks (Fazlyab et al., ; Chen et al., ; Latorre et al., ; Pauli et al., ) or verifying their robustness (Raghunathan et al., ; Fazlyab et al., ; Zhang, ; Lan et al., ; Chiu & Zhang, ), their application to the training of neural networks remains less explored.

In this paper, we study convex optimization representations for training a two-layer neural network with rectifier linear unit (ReLU) activations and a sufficiently wide hidden layer.
While the relationships between classical matrix factorization problems and two-layer linear neural networks, as well as those between matrix factorization and semidefinite programming, are well-established, our understanding of the connections between neural networks and convex optimization becomes less clear when non-linear activation functions are involved.
Our paper advances these connections by establishing links between convex optimization and neural networks with ReLU activations.

In light of this brief introduction, we summarize our key contributions:

- $\blacksquare$
We present a novel copositive program for training a two-layer ReLU neural network, providing an exact reformulation of the training problem when the network is sufficiently wide. Specifically, for fixed input and output dimensions and a given number of datapoints, there exists a finite critical width beyond which the network’s expressivity saturates, making the training problem equivalent to the proposed copositive formulation. We further extend this equivalence to the more general infinite-width regime, where the neural network is represented as a Lebesgue integral over a probability measure on the weights.
- $\blacksquare$
Although copositive programs are convex, solving them is NP-hard (Bomze et al., ). To mitigate this obstacle, we propose a semidefinite programming relaxation of the original formulation. We then numerically evaluate its tightness on two synthetic examples. We also assess its performance on real-data classification tasks, where, combined with a rounding heuristic, our approach achieves competitive test accuracy compared to Neural Network Gaussian Process (NNGP) (Lee et al., ) and Neural Tangent Kernel (NTK) (Jacot et al., ) methods.

We make our implementation available under
[https://github.com/KarthikPrakhya/SDPNN-IW](https://github.com/KarthikPrakhya/SDPNN-IW).

## 2 Background

###### Definition 1 (Positive semidefinite cone) .

A symmetric matrix $\bm{W}\in\mathbb{R}^{n\times n}$ is said to be positive semidefinite if $\bm{u}^{\top}\bm{W}\bm{u}\geq 0$ for all $\bm{u}\in\mathbb{R}^{n}$.
The set of all positive semidefinite matrices forms a self-dual convex cone, known as the positive semidefinite cone, which can be defined by

$$ $\mathcal{PSD}_{n}:=\mathrm{conv}\{\bm{w}\bm{w}^{\top}:\bm{w}\in\mathbb{R}^{n}\},$ (1) $$

where $\mathrm{conv}$ represents the convex hull.
We omit the subscript when the size is clear from the context.

###### Definition 1 (Positive semidefinite cone) .

#### Semidefinite programming

(SDP) is a powerful framework in convex optimization focused on minimization of a convex objective over the positive semidefinite cone.
SDPs can be solved in polynomial time under mild technical assumptions using interior-point methods (Vandenberghe & Boyd, ).
We refer to (Majumdar et al., ; Yurtsever et al., ) for an overview of recent advances in more scalable SDP solvers.

SDPs have been applied in neural networks for various tasks (often as a convex relaxation), including the estimation of Lipschitz constant (Chen et al., ; Fazlyab et al., ; Latorre et al., ; Pauli et al., ), verification of neural networks (Raghunathan et al., ; Zhang, ; Chiu & Zhang, ; Lan et al., ; Fazlyab et al., ), and for stability guarantees (Pauli et al., , ; Revay et al., ; Yin et al., ).
A recent line of work, initiated by Pilanci & Ergen (), investigates the use of SDPs for training neural networks.
We provide a detailed discussion and comparison with these approaches in Section 6.

###### Definition 2 (Copositive cone) .

A matrix $\bm{W}\in\mathbb{R}^{n\times n}$ is said to be copositive if $\bm{u}^{\top}\bm{W}\bm{u}\geq 0$ for all $\bm{u}\in\mathbb{R}^{n}_{+}$.
The cone of copositive matrices is called the copositive cone ($\mathcal{COP}$).

###### Definition 3 (Completely positive cone) .

The dual of copositive cone is known as the completely positive cone, defined as

$$ $\mathcal{CP}_{\!n}:=\mathrm{conv}\{\bm{w}\bm{w}^{\top}:\bm{w}\in\mathbb{R}^{n} _{+}\}.$ (2) $$

It is easy to see that completely positive cone is a subset of positive semidefinite cone, and positive semidefinite cone is a subset of copositive cone, as shown in Figure 1.

Figure: Figure 1: Cones.
Refer to caption: x1.png

Copositive programming (CP) is a subfield of convex optimization concerned with minimization of a convex objective function over copositive or completely positive matrices.
Despite its convexity, solving a CP is NP-Hard (Bomze et al., ).
Numerous NP-Hard problems are formulated as a CP, including the binary quadratic problems (Burer, ), problems of finding stability and chromatic numbers of a graph (De Klerk & Pasechnik, ; Dukanovic & Rendl, ), 3-partitioning problem (Povh & Rendl, ), and the quadratic assignment problem (Povh & Rendl, ).
We refer to the excellent surveys (Dür, ; Dür & Rendl, ) and references therein for more examples.

Recent research has broadened the scope of CP into data science.
For example, CP formulations have been utilized in machine learning for verifying neural networks (Brown et al., ), and problems such as graph matching and permutation synchronization in computer vision (Yurtsever et al., ).
CP formulations have also been proposed for training two-layer ReLU networks, though under specific data-related assumptions (Sahiner et al., ).

###### Definition 4 (Non-negative cone) .

The non-negative cone $\mathcal{N}$ comprises of entry-wise non-negative matrices.

###### Definition 5 (Doubly non-negative cone) .

The doubly non-negative cone is the set of matrices that are both positive semidefinite and elementwise non-negative, defined by $\mathcal{DNN}:=\mathcal{PSD}\cap\mathcal{N}.$
It is an important subset of the positive semidefinite cone that contains the cone of completely positive matrices and is frequently used in relaxations of CP formulations.

###### Definition 6 (Rank) .

For a matrix $\bm{W}\in\mathcal{PSD}_{n}$, the rank is the minimum number of vectors $\bm{w}_{i}\in\mathbb{R}^{n}$ required to express $\bm{W}$ as a convex combination of $\bm{w}_{i}\bm{w}_{i}^{\top}$.

###### Definition 7 (CP-rank) .

We refer to the rank restricted to decompositions involving non-negative vectors as the CP-rank:
The CP-rank of a matrix $\bm{W}\in\mathcal{CP}_{\!n}$ is the minimum number of non-negative vectors $\bm{w}_{i}\in\mathbb{R}_{+}^{n}$ that are needed to express $\bm{W}$ as a convex combination of $\bm{w}_{i}\bm{w}_{i}^{\top}$.

It is well known that the rank of any $n\times n$ matrix $\bm{W}$ cannot exceed $n$.
The following result from (Barioli & Berman, , Theorem 3.1) characterizes the maximal CP-rank.

###### Lemma 1 .

Let $\bm{W}$ be a completely positive matrix of rank $r\geq 2$.
Then, the CP-rank of $\bm{W}$ cannot exceed ${r(r+1)}/{2}-1$.
As an immediate corollary, the CP-rank of any $n\times n$ completely positive matrix for $n\geq 2$ is bounded above by ${n(n+1)}/{2}-1$, which we denote as the maximal CP-rank.

###### Definition 2 (Copositive cone) .

###### Definition 3 (Completely positive cone) .

###### Definition 4 (Non-negative cone) .

###### Definition 5 (Doubly non-negative cone) .

###### Definition 6 (Rank) .

###### Definition 7 (CP-rank) .

###### Lemma 1 .

## 3 CP Formulation for Training Two-Layer ReLU Networks

We now present our CP formulation for training two-layer RELU networks. We start by the *wide* regime with finite width before moving onto the infinite case that carries theoretical importance.

###### Definition 8 (Two-Layer ReLU Network) .

A two layer neural network $\psi:\mathbb{R}^{d}\rightarrow\mathbb{R}^{c}$ with $m$ hidden neurons and ReLU activations is given by

$$ $\psi(\bm{x})=\sum_{j=1}^{m}(\bm{x}^{\top}\bm{u}_{j})_{+}\bm{v}_{j}^{\top},$ ( Nn ) $$

where the function $(\cdot)_{+}=\max(0,\cdot)$ denotes the ReLU activation, $\{\bm{u}_{j}\in\mathbb{R}^{d}\}_{j=1}^{m}$ are the weights for the first layer and $\{\bm{v}_{j}\in\mathbb{R}^{c}\}_{j=1}^{m}$ are the second layer weights.

Assuming that the data matrix $\bm{X}\in\mathbb{R}^{n\times d}$ and labels $\bm{Y}\in\mathbb{R}^{n\times c}$ are given, and using the mean squared error loss and Tikhonov regularization (also known as weight-decay), the problem for training this neural network can be formulated as follows:

$$ $\min_{\begin{subarray}{c}\bm{u}_{j}\in\mathbb{R}^{d}\\ \bm{v}_{j}\in\mathbb{R}^{c}\end{subarray}}~{}\Big{\|}\sum_{j=1}^{m}(\bm{X}\bm{ u}_{j})_{+}\bm{v}_{j}^{\top}-\bm{Y}\Big{\|}_{F}^{2}+\frac{\gamma}{2}\sum_{j=1} ^{m}(\|\bm{u}_{j}\|_{2}^{2}+\|\bm{v}_{j}\|_{2}^{2}),$ ( Nn -Train) $$

where $\gamma\geq 0$ is the regularization parameter.

In the following, our main result is an equivalence between the above training problem and the following convex optimization problem with $\mathcal{PSD}$ and $\mathcal{CP}$ constraints when $m$ is sufficiently large:

###### Theorem 1 (CP Formulation) .

Consider the two-layer ReLU neural network in (Nn).
For any fixed input dimension $d$, output dimension $c$, and number of data points $n$, there exists a finite critical width $R\leq\max\{p,2n^{2}+n-1\}$, where $p=c+d+2n$, such that, for $m\geq R$ hidden neurons, the training problem (Nn-Train) is equivalent to the convex optimization problem (Cp-Nn), meaning that they have identical optimal values:

$$ $\displaystyle\min_{\bm{\Lambda}\in\mathbb{R}^{p\times p}}\,~{}\left\|\bm{P}_{ \alpha}\bm{\Lambda}\bm{P}_{v}^{\top}-\bm{Y}\right\|_{F}^{2}+\frac{\gamma}{2} \left(\mathrm{trace}(\bm{P}_{u}\bm{\Lambda}\bm{P}_{u}^{\top})+\mathrm{trace}( \bm{P}_{v}\bm{\Lambda}\bm{P}_{v}^{\top})\right)$ ( Cp - Nn ) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathrm{trace}( \bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})+\mathrm{trace}(\bm{M}\bm{ \Lambda}\bm{M}^{\top})=0$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{\Lambda} \in\mathcal{PSD}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{P}_{ \alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{CP},$ $$

where $\bm{M}=-\bm{P}_{\alpha}+\bm{P}_{\beta}+\bm{X}\bm{P}_{u}$ and the selection matrices are defined as follows:

$$ $\begin{aligned} &\bm{P}_{u}=\begin{bmatrix}\bm{0}_{d\times n}&\bm{0}_{d\times n }&\bm{I}_{d\times d}&\bm{0}_{d\times c}\end{bmatrix}&&~{}\bm{P}_{\alpha}= \begin{bmatrix}\bm{I}_{n\times n}&\bm{0}_{n\times n}&\bm{0}_{n\times d}&\bm{0} _{n\times c}\end{bmatrix}\\ &\bm{P}_{v}=\begin{bmatrix}\bm{0}_{c\times n}&\bm{0}_{c\times n}&\bm{0}_{c \times d}&\bm{I}_{c\times c}\end{bmatrix}&&~{}\bm{P}_{\beta}=\begin{bmatrix} \bm{0}_{n\times n}&\bm{I}_{n\times n}&\bm{0}_{n\times d}&\bm{0}_{n\times c} \end{bmatrix}\end{aligned}~{}~{}~{}\text{and}~{}~{}~{}\bm{P}_{\alpha\beta}= \begin{bmatrix}\bm{P}_{\alpha}\\ \bm{P}_{\beta}\end{bmatrix}$ $$

It is noteworthy that this formulation is independent of $m$ and hence does not scale in size with the number of hidden neurons.

###### Proof sketch.

The first challenge in (Nn-Train) arises from the non-linear ReLU activation function, which we address by decomposing $\bm{X}\bm{u}$ into its positive and negative parts, inspired by (Brown et al., ).
Specifically, for non-negative vectors $\bm{\alpha}$ and $\bm{\beta}$ such that $\mathrm{trace}(\bm{\alpha}\bm{\beta}^{\top})=0$, we have $\bm{X}\bm{u}=\bm{\alpha}-\bm{\beta}$ if and only if $\bm{\alpha}=(\bm{X}\bm{u})_{+}$ and $\bm{\beta}=(\bm{X}\bm{u})_{+}-\bm{X}\bm{u}$, corresponding to the positive and negative components, respectively.

Introducing the variable $\bm{\lambda}_{j}=[\bm{\alpha}_{j}^{\top}~{}\bm{\beta}_{j}^{\top}~{}\bm{u}^{
\top}~{}\bm{v}^{\top}]^{\top}\in\mathbb{R}^{p}$, we can reformulate (Nn-Train) as:

$$ $\displaystyle\min_{\bm{\lambda}_{j}\in\mathbb{R}^{p}}~{}~{}~{}\Big{\|}\sum_{j= 1}^{m}\bm{P}_{\alpha}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\bm{P}_{v}^{\top}- \bm{Y}\Big{\|}_{F}^{2}+\frac{\gamma}{2}\sum_{j=1}^{m}\big{(}\|\bm{P}_{u}\bm{ \lambda}_{j}\|_{2}^{2}+\|\bm{P}_{v}\bm{\lambda}_{j}\|_{2}^{2}\big{)},$ (3) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}\bm{P}_{\alpha\beta}\bm{ \lambda}_{j}\geq\bm{0}_{2n\times 1},~{}~{}~{}\mathrm{trace}(\bm{P}_{\alpha}\bm {\lambda}_{j}\bm{\lambda}_{j}^{\top}\bm{P}_{\beta}^{\top})=0,~{}~{}~{}\bm{M} \bm{\lambda}_{j}=\bm{0}_{n\times 1}.$ $$

Our goal is to show that (Cp-Nn) is equivalent to the problem (3) when $m$ is sufficiently large.
To this end, we establish that any matrix $\bm{\Lambda}\in\mathcal{PSD}$ satisfying $\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{CP}$ can be factorized as $\smash{\bm{\Lambda}=\sum_{j=1}^{R}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}}$, where $\bm{P}_{\alpha\beta}\bm{\lambda}_{j}\geq\bm{0}$, and $R$ is a finite number bounded by $\smash{R\leq\max\{p,2n^{2}+n-1\}}$ (see Lemma 2 in the supplementary material).
Here, the first term ($p$) represents the maximal rank of $\bm{\Lambda}$, and the second term ($2n^{2}+n-1$) is the maximal CP-rank of $\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}$.

Thus, for any $m\geq R$, we can reformulate (Cp-Nn) by expressing $\bm{\Lambda}$ as $\sum_{i=1}^{m}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}$, with the constraint $\bm{P}_{\alpha\beta}\bm{\lambda}_{j}\geq\bm{0}_{2n\times 1}$, without changing the global solution:

$$ $\displaystyle\min_{\bm{\lambda}_{j}\in\mathbb{R}^{p}}~{}~{}~{}\Big{\|}\sum_{j= 1}^{m}\bm{P}_{\alpha}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\bm{P}_{v}^{\top}- \bm{Y}\Big{\|}_{F}^{2}+\frac{\gamma}{2}\sum_{j=1}^{m}\big{(}\|\bm{P}_{u}\bm{ \lambda}_{j}\|_{2}^{2}+\|\bm{P}_{v}\bm{\lambda}_{j}\|_{2}^{2}\big{)},$ (4) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}\bm{P}_{\alpha\beta}\bm{ \lambda}_{j}\geq\bm{0}_{2n\times 1},~{}~{}~{}\sum_{j=1}^{m}\big{(}\mathrm{ trace}(\bm{P}_{\alpha}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\bm{P}_{\beta}^{ \top})+\|\bm{M}\bm{\lambda}_{j}\|_{2}^{2}\big{)}=0,$ $$

where we used the fact that $\mathrm{trace}(\bm{a}\bm{b}^{\top})=\mathrm{trace}(\bm{b}^{\top}\!\bm{a})$, which further simplifies to $\|\bm{a}\|_{2}^{2}$ when $\bm{a}=\bm{b}$.
We have omitted the $\mathcal{PSD}$ and $\mathcal{CP}$ constraints since they are inherently satisfied by the factorization.

Problems (3) and (4) share the same objective function.
Next, we show that their feasible sets are also identical.
Note that since each term in the summation constraint of problem (4) is non-negative, the sum can be zero only if each individual term is zero.
We complete the proof by noting $\|\bm{M}\bm{\lambda}_{j}\|_{2}^{2}=0$ if and only if $\bm{M}\bm{\lambda}_{j}=\bm{0}_{n\times 1}$.
For more details on the proof, we refer the reader to Appendix A.
∎

###### Remark 1 .

Although we focused on the MSE loss for simplicity, our results extend to any convex loss function $\ell:\mathbb{R}^{d}\times\mathbb{R}^{c}\to\mathbb{R}$.

###### Definition 8 (Two-Layer ReLU Network) .

###### Theorem 1 (CP Formulation) .

###### Proof sketch.

###### Remark 1 .

### 3.1 The Infinite Width Regime

###### Definition 9 (Infinite-width Two Layer RELU Networks) .

An infinite-width fully connected two-layer ReLU network can be expressed as a Lebesgue integral over a signed measure $\nu(\bm{u},\bm{v})$ defined on the weights $\bm{u}$ and $\bm{v}$, subject to the condition that $(\bm{x}^{\top}\bm{u})_{+}\bm{v}$ is integrable w.r.t. $\nu$ for any $\bm{x}\in\mathbb{R}^{d}$:

$$ $\psi(\bm{x})=\int_{\mathbb{R}^{d}\times\mathbb{R}^{c}}(\bm{x}^{\top}\bm{u})_{+ }\bm{v}^{\top}d\nu(\bm{u},\bm{v}).$ ( Nn $\!\smallint$ ) $$

A similar characterization of infinite-width neural networks can be found in (Mhaskar, ; Bach, ).
For our analysis, we restrict $\nu(\bm{u},\bm{v})$ to be a probability measure induced over the product space, $\nu:\mathcal{B}_{d}\times\mathcal{B}_{c}\rightarrow[0,1]$, where $\mathcal{B}_{n}$ is the Borel $\sigma$-algebra over $\mathbb{R}^{n}$.

While the case of countably infinite number of hidden neurons, which is represented by an infinite sum in (Nn), is captured as a special case of (Nn$\!\smallint$) with a discrete probability measure, the (Nn$\!\smallint$) model also accommodates other notions of infinite-width networks, as the underlying probability measure can be continuous or a combination of discrete and continuous components.

The training of an infinite-width ReLU neural network in (Nn$\!\smallint$) with the mean squared error loss can be expressed as follows:

$$ $\min_{\nu:\mathcal{B}_{d}\times\mathcal{B}_{c}\rightarrow[0,1]}\left\|\int_{ \mathbb{R}^{d}\times\mathbb{R}^{c}}(\bm{X}\bm{u})_{+}\bm{v}d\nu(\bm{u},\bm{v}) -\bm{Y}\right\|_{F}^{2}+\frac{\gamma}{2}\int_{\mathbb{R}^{d}\times\mathbb{R}^{ c}}\big{(}\|\bm{u}\|_{2}^{2}+\|\bm{v}\|_{2}^{2}\big{)}d\nu(\bm{u},\bm{v}),$ ( Nn $\!\smallint$ -Train) $$

where optional Tikhonov regularization is controlled by the parameter $\gamma\geq 0$.
The following theorem presents our second main result, which shows that the training problem in (Nn$\!\smallint$-Train) is also equivalent to our CP formulation.

###### Theorem 2 .

Consider the two-layer infinite-width ReLU neural network defined in (Nn$\!\smallint$).
For any input dimension $d$ and output dimension $c$, the training problem in (Nn$\!\smallint$-Train) is equivalent to the convex optimization problem in (Cp-Nn), implying that they share the same optimal values.

We refer the reader to Appendix B for the proof.
This result implies that the infinite-width ReLU network training problem in (Nn$\!\smallint$-Train) can be solved to global optimality by training a finite-width ReLU network with a width exceeding the critical threshold defined in Theorem 1.
The resulting solution can then be represented as a discrete probability measure with finite support.

###### Definition 9 (Infinite-width Two Layer RELU Networks) .

###### Theorem 2 .

## 4 SDP Relaxation for Training Two-Layer ReLU Networks

Despite its convexity, solving (Cp-Nn) remains intractable due to the complete positivity constraint.
To tackle this challenge, we propose an SDP relaxation.

###### Proposition 1 (SDP-NN) .

The following SDP is a relaxation of the problem (Cp-Nn), obtained by replacing the completely positive cone with the doubly non-negative cone, resulting in the constraints ${\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{N}}$ and $\smash{\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{
PSD}}$.
The latter is omitted as it is inherently satisfied for all $\bm{\Lambda}\in\mathcal{PSD}$:

$$ $\displaystyle\min_{\bm{\Lambda}\in\mathbb{R}^{p\times p}}\,~{}\left\|\bm{P}_{ \alpha}\bm{\Lambda}\bm{P}_{v}^{\top}-\bm{Y}\right\|_{F}^{2}+\frac{\gamma}{2} \left(\mathrm{trace}(\bm{P}_{u}\bm{\Lambda}\bm{P}_{u}^{\top})+\mathrm{trace}( \bm{P}_{v}\bm{\Lambda}\bm{P}_{v}^{\top})\right)$ ( Sdp - Nn ) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}~{}\,\,\mathrm{trace}( \bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})+\mathrm{trace}(\bm{M}\bm{ \Lambda}\bm{M}^{\top})=0$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{\Lambda} \in\mathcal{PSD}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{P}_{ \alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{N}.$ $$

###### Remark 2 .

This is essentially the zeroth-order Sum-of-Squares (SoS) relaxation of (Cp-Nn).
While tighter relaxations could be constructed using the SoS hierarchy (Parrilo, ), the computational requirements often make these higher-order formulations impractical.

After solving (LABEL:eqn:sdp-nn-1), a rounding procedure is required to extract the weights of the trained neural network from the lifted solution $\bm{\Lambda}_{\star}$.

###### Proposition 2 (Rounding) .

Selecting the rounding width $R$ as the critical width of the network, we can formulate the rounding problem as follows:

$$ $\displaystyle\min_{\bm{\lambda}\in\mathbb{R}^{p\times R}}\,~{}\,~{}~{}~{}\phi( \bm{\lambda})~{}:=~{}\big{\|}\bm{\Lambda}_{\star}-\bm{\lambda}\bm{\lambda}^{ \top}\big{\|}_{F}^{2}$ (5) $\displaystyle~{}~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bm{M}\bm{ \lambda}=\bm{0}_{n\times R},\quad\bm{P}_{\alpha}\bm{\lambda}\odot\bm{P}_{\beta }\bm{\lambda}=\bm{0}_{n\times R},\quad\bm{P}_{\alpha\beta}\bm{\lambda}\geq\bm{ 0}_{2n\times R},$ $$

where $\odot$ denotes the Hadamard product.
The first and second layer weights can be recovered from the rounded solution as $\bm{P}_{u}\bm{\lambda}$ and $\bm{P}_{v}\bm{\lambda}$, respectively.

It is straightforward to verify that the feasible set of problem (5) coincides with that in Eq. 3, which is, in turn, equivalent to the feasible set of (Cp-Nn).
Consequently, this formulation aims to identify the closest point $\bm{\lambda}\bm{\lambda}^{\top}$ to the SDP solution $\bm{\Lambda}_{\star}$ within the feasible set of (Cp-Nn).
While finding a global solution is intractable due to the non-convexity of the objective function and constraints, we propose an effective heuristic based on the three-operator splitting (TOS) method (Davis & Yin, ).
We define the two sets comprising the constraints as:

$$ $\displaystyle\mathcal{D}_{1}$ $\displaystyle:=\{\bm{\lambda}\in\mathbb{R}^{p\times R}:\bm{P}_{\alpha}\bm{ \lambda}\odot\bm{P}_{\beta}\bm{\lambda}=\bm{0}_{n\times R},~{}\bm{P}_{\alpha \beta}\bm{\lambda}\geq\bm{0}_{2n\times R}\}$ $\displaystyle\mathcal{D}_{2}$ $\displaystyle:=\{\bm{\lambda}\in\mathbb{R}^{p\times R}:\bm{M}\bm{\lambda}=\bm{ 0}_{n\times R}\}$ $$

Then, we can apply TOS for the rounding, starting from an initial estimate $\bar{\bm{\lambda}}^{0}\in\mathbb{R}^{p\times R}$ and iteratively updating it by the following formula:

$$ $\displaystyle\bm{\lambda}^{k}$ $\displaystyle=\mathrm{proj}_{\mathcal{D}_{1}}\big{(}\bar{\bm{\lambda}}^{k}\big {)}$ (TOS) $\displaystyle\hat{\bm{\lambda}}^{k}$ $\displaystyle=\mathrm{proj}_{\mathcal{D}_{2}}\big{(}2\bm{\lambda}^{k}-\bar{\bm {\lambda}}^{k}-\eta\nabla\phi(\bm{\lambda}^{k})\big{)}$ $\displaystyle\bar{\bm{\lambda}}^{k+1}$ $\displaystyle=\bar{\bm{\lambda}}^{k}-\bm{\lambda}^{k}+\hat{\bm{\lambda}}^{k}$ $$

where $\eta>0$ is the step-size parameter.
We can compute the gradient by $\nabla\phi(\bm{\lambda})=4(\bm{\lambda}\bm{\lambda}^{\top}-\bm{\Lambda}_{\star
})\bm{\lambda}$.
$\mathrm{proj}_{\mathcal{D}_{2}}(\bm{\lambda})=(\bm{I}-\bm{M}^{\dagger}\bm{M})
\bm{\lambda}$ projects $\bm{\lambda}$ onto the nullspace of $\bm{M}$.
And $\mathrm{proj}_{\mathcal{D}_{1}}$ can be computed as

$$ $\mathrm{proj}_{\mathcal{D}_{1}}(\bm{\lambda})_{i,j}=\begin{cases}0&1\leq i\leq 2 n,~{}\bm{\lambda}_{i,j}<0\\ 0&1\leq i\leq n,~{}\bm{\lambda}_{i+n,j}>\bm{\lambda}_{i,j}\\ 0&n+1\leq i\leq 2n,~{}\bm{\lambda}_{i-n,j}>\bm{\lambda}_{i,j}\\ \bm{\lambda}_{i,j}&\mbox{otherwise}\end{cases}$ (6) $$

###### Remark 3 .

We employ the rounding step primarily to demonstrate that the solution obtained from the SDP relaxation is reasonable.
Therefore, our main objective is to evaluate the quality of the SDP relaxation, rather than to guarantee convergence in the rounding phase.
To the best of our knowledge, there are no known convergence guarantees for TOS in our setting, as existing results for TOS in non-convex optimization apply only to smooth non-convex terms (Bian & Zhang, ; Yurtsever et al., ) and do not extend to non-convex constraint sets.
Nonetheless, TOS has shown to be an effective heuristic in our experiments.

**Table 1: Summary of datasets, their sizes, and dimensions of the corresponding SDP relaxations.**
| Dataset | $\#$ Instances | $\#$ Inp. Feat | $\#$ Out. Feat | $\#$ Variables | $\#$ Constraints |
| --- | --- | --- | --- | --- | --- |
| Random | 25 | 2 | 5 | 3249 | 8999 |
| Spiral | 60 | 2 | 3 | 15625 | 45651 |
| Iris | 75 | 4 | 3 | 24649 | 71799 |
| Ionosphere | 175 | 34 | 2 | 148996 | 420493 |
| Pima Indians | 383 | 8 | 2 | 602176 | 1791109 |
| Bank Notes | 685 | 4 | 2 | 1893376 | 5663653 |

###### Proposition 1 (SDP-NN) .

###### Remark 2 .

###### Proposition 2 (Rounding) .

###### Remark 3 .

## 5 Numerical Experiments

We conduct a series of experiments over synthetic and real datasets to investigate the empirical tightness of our SDP relaxation.

#### Datasets

We use the following datasets in our experiments:

Random:
A synthetic dataset with a data matrix $\bm{X}\in\mathbb{R}^{25\times 2}$ and labels $\bm{Y}\in\mathbb{R}^{25\times 5}$, generated using a 2-layer neural network with 100 hidden neurons.
We create 100 random datasets by initializing the entries of $\bm{X}$ and the generator network’s weights using a standard normal distribution.

Spiral:
An artificially-generated 3-class classification dataset, described by (Sahiner et al., ).
It includes 60 samples (20 per class), with a data matrix $\bm{X}\in\mathbb{R}^{60\times 2}$ containing 2 input features and one-hot encoded labels $\bm{Y}\in\mathbb{R}^{60\times 3}$.

Iris (Fisher, ):
The dataset consists of 150 samples of iris flowers from three different classes, each sample is described by four features.
The training partition includes a data matrix $\bm{X}\in\mathbb{R}^{75\times 4}$ and one-hot encoded labels $\bm{Y}\in\mathbb{R}^{75\times 3}$.

Ionosphere (Sigillito et al., ):
A radar dataset with 351 instances and 34 input features for binary classification, with class imbalance.
The training partition includes a data matrix $\bm{X}\in\mathbb{R}^{175\times 34}$ and one-hot encoded labels $\bm{Y}\in\mathbb{R}^{175\times 2}$.

Pima Indians Diabetes (Smith et al., ):
A diabetes prediction dataset that consists of 768 patients, with 8 medical predictors as features, and a binary classification task, with class imbalance.
Out of 768 instances, we have 268 positive instances.
The training partition includes a data matrix $\bm{X}\in\mathbb{R}^{383\times 8}$ and one-hot encoded labels $\bm{Y}\in\mathbb{R}^{383\times 2}$.

Bank Notes Authentication (Lohweg, ):
A binary classification dataset with 1372 instances, where features extracted using wavelet transforms are used to determine whether a banknote is genuine or forged.
The training data includes $\bm{X}\in\mathbb{R}^{685\times 4}$ and one-hot encoded labels $\bm{Y}\in\mathbb{R}^{685\times 2}$.

Table 1 provides an overview of the dataset sizes and the dimensionality of the corresponding problem size of (LABEL:eqn:sdp-nn-1).
All real datasets were split 50-50 into train and test partitions.

#### Computational environment.

The experiments were conducted on a Intel Xeon Gold 6132 with 192 GB of RAM and 2x14 cores.
SDP formulations were solved by CVXPY (Diamond & Boyd, ), employing either the interior point method (IPM) solver by MOSEK (Andersen & Andersen, ), or the Splitting Cone Solver (SCS) (O’donoghue et al., ), depending on the problem size.

### 5.1 Empirical Approximation Ratio of the SDP Relaxation

We evaluate the tightness of our SDP relaxation by comparing its optimal objective value to the training loss obtained using Stochastic Gradient Descent (SGD) on (Nn-Train).
Due to the difficulty of approximating the global optimum with SGD on real datasets, we conduct this experiment using the synthetic Random and Spiral datasets.

The results are summarized in Table 2.
For each dataset, we solve the training problem using regularization parameters $\gamma=0.1$ and $\gamma=0.01$.
The table columns show the training loss achieved by SGD for various hidden layer sizes, ranging from $m=5$ to $m=300$, along with the training loss achieved by the SDP-NN formulation.
For the Random dataset, the results are averaged over 100 trials, with error bars indicating standard deviations.

The last column reports the approximation ratio (AR) of the SDP relaxation, computed as the ratio of the loss obtained by SDP-NN to the loss achieved by SGD with $300$ hidden neurons.
Empirically, we find that the approximation ratio is above $76.90\%$ in these experiments.
We also observe that the ratio improves for higher regularization parameters, which is expected, as stronger regularization tends to smooth the optimization landscape and reduces the impact of non-convexity.

###### Remark 4 .

Ratio reported in Table 2 represents only a lower bound on the actual approximation ratio; since, (i) SGD may converge to a local solution, and (ii) $300$ hidden neurons used in this experiment is smaller than the actual critical width, which scales quadratically with the number of data points, and quickly exceeds practical limits.

To partially address these challenges, we ran SGD five times with different random initializations for each configuration and selected the best outcome.
Even so, there is no guarantee of reaching a global optimum.
Similarly, although $m=300$ is below the critical width, we observe that the results tend to saturate, and increasing $m$ further does not lead to any significant improvements for SGD.
Despite these limitations, SGD training can provide an upper bound on the optimal value for infinite-width neural network training.
To complement this, we compute a lower bound on the optimal values of the SDP formulations using the dual objectives in CVXPY, utilizing the interior-point solver from MOSEK.
Together, these bounds allow us to establish a reliable approximation lower-bound for the relaxation.

The training loss curves of SGD for each experiment, along with implementation details such as learning rate and number of epochs, are provided in Appendix D.

**Table 2: SGD Loss and Approximation Ratio (AR) (%) of SDP-NN.**
| Dataset | $\gamma$ | Training Objective | AR |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  | SGD-5 | SGD-10 | SGD-100 | SGD-200 | SGD-300 | SDP-NN | SDP-NN |
| Random | 0.1 | 18.27 ±5.76 | 8.60 ±1.16 | 8.09 ±1.04 | 8.09 ±1.04 | 8.09 ±1.04 | 7.28 ±0.98 | 89.93 |
| 0.01 | 11.78 ±5.53 | 1.32 ±0.25 | 0.94 ±0.12 | 0.94 ±0.12 | 0.94 ±0.12 | 0.76 ±0.10 | 80.66 |  |
| Spiral | 0.1 | 16.64 | 16.59 | 16.59 | 16.59 | 16.59 | 16.24 | 97.84 |
| 0.01 | 15.56 | 15.21 | 15.16 | 15.16 | 15.16 | 11.66 | 76.90 |  |

###### Remark 4 .

### 5.2 Prediction Performance of the SDP Relaxation

We assess the prediction quality of our SDP relaxation combined with the TOS-based rounding procedure on real datasets using regularization parameters $\gamma=0.1$ and $\gamma=0.01$.
Due to computational constraints, we set the rank $R$ in the rounding procedure to $300$, which is lower than the actual critical width.
The rounding procedure is initialized using the square root of the SDP-NN solution, obtained via SVD decomposition.
We then adjust the resulting factor matrices by either truncating or expanding them to dimensions $\mathbb{R}^{p\times 300}$.
The TOS step size is set as $\eta={1}/{\|\bm{\Lambda}_{\star}\|_{2}}$, where $\bm{\Lambda}_{\star}$ denotes the solution to the SDP-NN, and the algorithm is run for $1000$ iterations.

We benchmark our test accuracy against the NNGP and NTK kernel methods and include results from SGD with $300$ hidden neurons as a baseline.
The results of this experiment are summarized in Table 3, where the weighted average F1 score and test accuracy are used as the performance metrics.
For classification, the threshold was varied from 0 0 to $1$ in increments of $0.01$, and the threshold that maximized the overall test accuracy was selected.
Overall, the predictions obtained from the SDP relaxation are competitive across most settings.

## 6 Discussion of Related Work

We now review the literature on (in)finite width neural networks as well as convex optimization techniques for training them.

#### Convex optimization for neural network training

Early discussions of convex neural networks can be found in (Bengio et al., ; Bach, ). These initial studies primarily focus on networks with infinite width, which leads to optimization problems in an infinite-dimensional space.

Sahiner et al. () developed convex equivalents for the non-convex neural network training problem with ReLU activations, based on a theory of convex semi-infinite duality.
Different than ours, their approach involves an explicit summation over all possible sign patterns in the ReLU layer.
Since the number of distinct sign patterns grows exponentially with the rank of the data matrix, their analysis is primarily focused on data matrices of fixed small rank or the spike-free structures.
Building on this duality theory, Ergen & Pilanci () derived convex optimization formulations for two- or three-layer convolutional neural networks (CNNs) with ReLU activations.
Bartan & Pilanci () explored semidefinite lifting for training neural networks with polynomial activations, which is further extended to quantized neural networks with polynomial activations in (Bartan & Pilanci, ), and to deep neural networks with polynomial and ReLU activations in (Bartan & Pilanci, ).
Ergen & Pilanci () demonstrated that training multiple three-layer ReLU regularized sub-networks can be equivalently cast as a convex optimization problem in a higher-dimensional space.
Ergen & Pilanci (, ) showed that optimal hidden-layer weights for two-layer and deep ReLU neural networks are the extreme points of a convex set.
They further demonstrated the application of these methods in transfer learning tasks with large language models.
Sahiner et al. () developed convex optimization problems for vision transformers, focusing on a single self-attention block with ReLU activation.
Finally, Wang et al. () introduced randomized algorithms from a geometric algebra perspective to efficiently address the enumeration of sign patterns.
Sahiner et al. () applied the Burer-Monteiro factorization (Burer & Monteiro, ) to efficiently solve these formulations at a larger scale.

Our mathematical techniques differ from prior work.
We do not require enumeration or sampling of sign patterns, nor do we rely on the aforementioned duality theory.
We formulate the training of (in)finite-width ReLU neural networks as a finite-dimensional convex optimization problem, assuming fixed input/output dimensions and a set number of data points.
Since our exact CP formulation is intractable, we propose SDP relaxations that can be solved in polynomial time.

**Table 3: Prediction quality metrics.**
| Method | Iris | Ionosphere | Pima Indians | Bank Notes |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\gamma=0.1$ | $\gamma=0.01$ | $\gamma=0.1$ | $\gamma=0.01$ | $\gamma=0.1$ | $\gamma=0.01$ | $\gamma=0.1$ | $\gamma=0.01$ |  |
| Weighted Average F1 Score |  |  |  |  |  |  |  |  |
| SGD | 0.96 | 0.987 | 0.915 | 0.898 | 0.626 | 0.594 | 0.993 | 0.992 |
| NNGP | 0.909 | 0.987 | 0.912 | 0.910 | 0.650 | 0.709 | 0.986 | 0.999 |
| NTK | 0.993 | 0.993 | 0.924 | 0.920 | 0.666 | 0.699 | 0.992 | 0.997 |
| SDP-NN | 0.987 | 0.987 | 0.924 | 0.927 | 0.679 | 0.703 | 0.930 | 0.893 |
| Overall Test Accuracy |  |  |  |  |  |  |  |  |
| SGD | 0.960 | 0.987 | 0.891 | 0.88 | 0.583 | 0.557 | 0.988 | 0.985 |
| NNGP | 0.827 | 0.973 | 0.891 | 0.897 | 0.534 | 0.602 | 0.972 | 0.999 |
| NTK | 0.987 | 0.987 | 0.903 | 0.903 | 0.586 | 0.620 | 0.984 | 0.994 |
| SDP-NN | 0.987 | 0.987 | 0.920 | 0.909 | 0.646 | 0.625 | 0.860 | 0.767 |

#### Infinite-Width Neural Networks (IWNNs)

Investigations of IWNNs address fundamental questions about the limits of their learning ability and the types of functions they can approximate.
This focus is also relevant as overparameterized networks, which are commonly used in practice, often demonstrate strong generalization to unseen data (Zhang et al., ).
Another motivation for studying IWNNs is their connection to other areas of machine learning, such as Gaussian processes (Neal, ).
There are two main approaches for training infinite-width fully connected neural networks in the literature: weakly-trained networks and fully-trained networks.

Weakly-trained networks are those where only the top classification layer is trained, while all other layers remain randomly initialized.
Lee et al. (); Matthews et al. () studied fully-connected weakly-trained IWNNs, while Novak et al. (); Garriga-Alonso et al. () focused on convolutional variants.
Yang () extended this analysis to other weakly-trained infinite-width architectures.
All these works adopt a Gaussian process (GP) perspective of weakly-trained IWNNs.

In contrast, fully-trained networks are those in which all the layers are trained, using variations of gradient descent.
Jacot et al. () introduced Neural Tangent Kernel (NTK), which captures the behavior of fully-connected deep neural networks in the infinite-width limit, trained using gradient descent.
The NTK is different from the GP kernels: it is defined by the gradient of a fully connected network’s output with respect to its weights, based on random initialization.
Arora et al. () gave a rigorous, non-asymptotic proof that the NTK captures the behavior of a fully-trained wide neural network under weaker conditions and introduced the concept of convolutional neural tangent kernel (CNTK).
They also developed an exact and efficient dynamic programming algorithm to compute CNTKs for ReLU activation.
Lee et al. () conducted an empirical study comparing the Neural Network Gaussian Process (NNGP) kernel (Lee et al., ) and the NTK in both finite-width and IWNNs. Unlike these, we do not adopt a GP-based view or kernel methods. Instead, we model fully-trained fully-connected IWNNs using a CP formulation and corresponding computationally tractable SDP relaxations.

#### Other related works

We address the non-convex problem of training two-layer ReLU neural networks by transforming it into a convex optimization problem in a higher-dimensional space.
This approach belongs to a class of techniques known as *lifting*, where the decision variable $\bm{w}$ is replaced with a quadratic term $\bm{W}=\bm{w}\bm{w}^{\top}$, enabling the convex reformulation or relaxation of an otherwise non-convex problem, see (Bomze et al., ; Burer, ; Bao et al., ; Anstreicher, ) and references therein.
In the context of neural networks, Brown et al. () applied lifting for the verification of neural network outputs.
A key distinction between our formulation and the prior work in (Burer, ; Brown et al., ) lies in the structure of the objective function.
In previous work, the objective function is quadratic in the original space and simplifies to a linear objective in the lifted space, which plays a crucial role in the analysis, as solutions to linear minimization problems always appear at extreme points of the feasible set.
In our problem, the objective function is convex in the lifted space (rather than linear), and we establish an exact correspondence between the feasible sets in the original and lifted spaces for both wide and infinite-width networks.

Our work should not be confused with the line of research introduced by Askari et al. () known as lifted neural networks, as the notion of lifting in our work differs from theirs.
Both the goals and approaches are distinct.
They reformulate the problem as a biconvex optimization, enabling layer-wise training via block coordinate descent.
They do not introduce the quadratic terms that are central to our approach, and their resulting problem formulation remains non-convex.

## 7 Conclusions and Future Directions

We derived a convex optimization formulation for the training problem of a two-layer ReLU neural network with a sufficiently wide hidden layer.
Despite its convexity, the problem is intractable using classical computational methods due to the completely positive cone constraints.
To address this, we proposed an SDP relaxation that can be solved in polynomial time using off-the-shelf SDP solvers, combined with a rounding heuristic.
Notably, the size of our formulation is independent of the network width.
We evaluated the tightness and predictive performance of this relaxation through numerical experiments.
The results are promising, as our approach yields competitive performance compared to established methods such as NNGP and NTK kernels on test datasets.

#### Limitations

An obvious limitation of the proposed framework is its lack of scalability, a common issue for many SDP formulations in machine learning applications.
Since the problem size scales quadratically with the number of data points, solving these formulations at the scale required for real-world neural network applications poses a significant computational challenge.
Nevertheless, our findings offer valuable theoretical insights into ReLU networks, for instance, allowing us to interpret the network width as the (CP-)rank of the corresponding completely positive program.

#### Future work

We proposed an exact CP formulation, which is computationally intractable for classical computers.
However, rapid progress in quantum computing poses a promising alternative (Birdal et al., ).
In particular, hybrid classical-quantum algorithms have demonstrated potential for solving specific classes of copositive programs (Yurtsever et al., ; Brown et al., ).
Developing and implementing algorithms—even as proof-of-concept experiments—to train ReLU networks on these platforms represents a compelling research direction.

Over-parameterized networks are often easier to optimize using traditional methods such as stochastic gradient descent (Livni et al., ).
Although it is possible to construct compact neural networks with fewer parameters that generalize comparably, these compressed models are significantly more challenging to optimize in the reduced-dimensional space (Arora et al., ).
A similar phenomenon is observed in the SDP literature: Waldspurger & Waters () demonstrated that even when an SDP has a unique rank-1 solution, factorizing the variable with a rank below a certain critical threshold may introduce spurious local minima.
Given these parallels, a natural question arises: is there a corresponding bound for completely positive programs, and could this provide insights into the required width of ReLU networks in practice?

## Ethics, Reproducibility & Acknowledgments

#### Ethics statement

This work adheres to the Imperial College London Code of Ethics, ensuring responsible use of publicly available datasets and maintaining full transparency in experimental design and reporting. The study complies with ethical guidelines on research integrity and considers potential societal impacts. Since the contributions are primarily theoretical, the focus is on advancing our understanding of the theory of artificial intelligence rather than targeting specific applications that may have direct societal impacts.

#### Reproducibility statement

All datasets used in this work are referenced and publicly accessible.
The experimental configurations and computational environment are outlined in detail within the main text and the appendix.
We make our implementation available under [https://github.com/KarthikPrakhya/SDPNN-IW](https://github.com/KarthikPrakhya/SDPNN-IW).
Furthermore, all random seeds are provided to facilitate precise reproducibility.

#### Acknowledgments

AY and KP were supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.
We thank the High Performance Computing Center North (HPC2N) at Umeå University for providing computational resources and valuable support during test and performance runs.
The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS), partially funded by the Swedish Research Council through grant agreement no. 2022-06725.
TB was supported by a UKRI Future Leaders Fellowship [grant number MR/Y018818/1].
We acknowledge the use of OpenAI’s ChatGPT for editorial assistance in preparing this manuscript.

## Appendix

## Appendix A Proof of Theorem 1

Our formulation in Theorem 1 involves a $(p\times p)$ matrix decision variable $\bm{W}\in\mathcal{PSD}$, which contains a $(q\times q)$ principal submatrix $\bm{P}\bm{W}\bm{P}^{\top}\in\mathcal{CP}$, where $\bm{P}=\left[\bm{I}_{q\times q}~{}~{}\bm{0}_{q\times(p-q)}\right]$ is the row-selection operator.
The following lemma plays a crucial role in our analysis, in showing that there exists a finite critical width for a two-layer ReLU network, beyond which increasing the number of hidden neurons does not improve the network’s expressive power.

###### Lemma 2 .

Let $\bm{W}\in\mathbb{R}^{p\times p}$ be a positive semidefinite matrix of rank $k$, with a completely positive principal submatrix $\widehat{\bm{W}}\in\mathbb{R}^{q\times q}$ of CP-rank $K$.
Without loss of generality, we assume $\widehat{\bm{W}}$ is the leading principal submatrix of $\bm{W}$ given by $\widehat{\bm{W}}=\bm{P}\bm{W}\bm{P}^{\top}$ with $\bm{P}=\left[\bm{I}_{q\times q}~{}~{}\bm{0}_{q\times(p-q)}\right]$.
Then, $\bm{W}$ can be decomposed as

$$ $\bm{W}=\sum_{j=1}^{R}\bm{w}_{j}\bm{w}_{j}^{\top},\quad\mathrm{such\,that}\quad \bm{w}_{j}\in\mathbb{R}^{p}~{}~{}\mathrm{and}~{}~{}\bm{P}\bm{w}_{j}\in\mathbb{ R}_{+}^{q}$ (7) $$

for some $R\leq\max\{k,K\}$.
As an immediate corollary, for any $(p\times p)$ PSD matrix with a $(q\times q)$ dimensional CP submatrix (without any restriction on their rank or CP-rank), we have $R\leq\max\{p,R_{q}\}$, where $R_{q}$ is the maximal CP-rank $R_{q}:=q(q+1)/2-1$.

We will use the following result (see Lemma 1 in (Xu, )) in the proof of Lemma 2.

###### Fact 1 .

If a PSD matrix $\bm{A}\in\mathbb{R}^{n\times n}$ admits two distinct factorizations, $\bm{A}=\bm{C}\bm{C}^{\top}$ and $\bm{A}=\bm{D}\bm{D}^{\top}$, where $\bm{C},\bm{D}\in\mathbb{R}^{n\times N}$, then we can find an orthogonal matrix $\bm{Q}\in\mathbb{R}^{N\times N}$ such that $\bm{D}=\bm{C}\bm{Q}$.

###### Proof of 1 .

Let $\bm{c}_{i}\in\mathbb{R}^{N}$ and $\bm{d}_{i}\in\mathbb{R}^{N}$ be the $i$th rows of $\bm{C}$ and $\bm{D}$, respectively.
Since $\{\bm{c}_{i}\}_{i=1}^{n}$ and $\{\bm{d}_{i}\}_{i=1}^{n}$ span the same subspace, there exists a linear transformation $\bm{Q}$ mapping $\mathrm{span}(\bm{c}_{1},\dots,\bm{c}_{n})$ to $\mathrm{span}(\bm{d}_{1},\dots,\bm{d}_{n})$ such that $\bm{d}_{i}=\bm{c}_{i}\bm{Q}$ for each $i=1,\dots,n$.
By construction, we have $\langle\bm{d}_{i},\bm{d}_{j}\rangle=\langle\bm{c}_{i}\bm{Q},\bm{c}_{j}\bm{Q}
\rangle=\langle\bm{c}_{i},\bm{c}_{j}\rangle$, so $\bm{Q}$ preserves inner products on this subspace.
Therefore, by extending $\bm{Q}$ to all of $\mathbb{R}^{N}$, we obtain an orthogonal matrix.
∎

###### Proof of Lemma 2 .

For the PSD matrix $\bm{W}$ and any non-negative integer $R\geq k$, we have a factorization of the form $\bm{W}=\bm{B}\bm{B}^{\top}$, where $\bm{B}\in\mathbb{R}^{p\times R}$.
This factorization is not necessarily unique, and the columns of $\bm{B}$ are not required to be linearly independent.
Similarly, for the CP matrix $\widehat{\bm{W}}$, we have a factorization $\widehat{\bm{W}}=\bm{D}\bm{D}^{\top}$, where $\bm{D}\in\mathbb{R}_{+}^{q\times R}$ for any $R\geq K$.
As before, this factorization is not necessarily unique, and the columns of $\bm{D}$ need not be linearly independent.

Consider these factorizations with $R=\max\{k,K\}$.
Define $\bm{C}$ as the first $q$ rows of $\bm{B}$, given by $\bm{C}=\bm{P}\bm{B}$.
We have $\widehat{\bm{W}}=\bm{C}\bm{C}^{\top}$.
Then, by 1, since any CP matrix is also PSD, there exists an orthogonal matrix $\bm{Q}\in\mathbb{R}^{R\times R}$ such that $\bm{D}=\bm{C}\bm{Q}$.
Finally, let us define $\bm{E}=\bm{B}\bm{Q}$.
It is easy to verify that $\bm{E}\bm{E}^{\top}=\bm{W}$, and $\bm{P}\bm{E}=\bm{D}$ is entrywise non-negative.
Hence, if we set $\bm{w}_{j}$ equal to the $j$th column of $\bm{E}$, we obtain the desired decomposition.
∎

The following lemma will be used to handle the non-convexity introduced by ReLU activation:

###### Lemma 3 .

For real-valued vector $\hat{\bm{\alpha}}\in\mathbb{R}^{n}$ and non-negative real-valued vectors $\bm{\alpha},\bm{\beta}\in\mathbb{R}_{+}^{n}$ with $\bm{\alpha}:=(\hat{\bm{\alpha}})_{+}$ and $\bm{\beta}:=\bm{\alpha}-\hat{\bm{\alpha}}$, the following conditions are equivalent:

$$ $\bm{\alpha}=(\hat{\bm{\alpha}})_{+}\quad\text{if and only if}\quad\mathrm{ trace}({\bm{\alpha}}\bm{\beta}^{\top})=0.$ (8) $$

###### Proof of Lemma 3 .

First, observe that $\mathrm{trace}(\bm{\alpha}\bm{\beta}^{\top})=\mathrm{trace}(\bm{\beta}^{\top}
\bm{\alpha})=\bm{\beta}^{\top}\bm{\alpha}=\sum_{i=1}^{n}\alpha_{i}\beta_{i}.$

Assume that the left-hand side in Eq. 8 holds.
Note that if $\bm{\alpha}$ represents the positive part of $\hat{\bm{\alpha}}$, then $\bm{\beta}$ is the negative part.
Consequently, for each index $i$, either $\alpha_{i}$ or $\beta_{i}$ is zero, hence $\mathrm{trace}(\bm{\alpha}\bm{\beta}^{\top})=0$.

Now, suppose that the right-hand side in Eq. 8 holds.
Given that the elements of $\bm{\alpha}$ and $\bm{\beta}$ are non-negative, we have $\mathrm{trace}(\bm{\alpha}\bm{\beta}^{\top})=0$ if and only if $\alpha_{i}\beta_{i}=\alpha_{i}(\alpha_{i}-\hat{\alpha}_{i})=0$ for all $i=1,\ldots,n$.
This implies that if $\alpha_{i}>0$, then $\hat{\alpha}_{i}=\alpha_{i}$.
Otherwise, if $\alpha_{i}=0$, we have $\hat{\alpha}_{i}=\alpha_{i}-\beta_{i}\leq 0$.
Thus, we conclude that $\bm{\alpha}=(\hat{\bm{\alpha}})_{+}$.
∎

###### Proof of Theorem 1 .

We start by rewriting Eq. Nn-Train in the constrained form as

$$ $\displaystyle\min_{\begin{subarray}{c}\bm{u}_{j}\in\mathbb{R}^{d},\bm{v}_{j} \in\mathbb{R}^{c}\\ \bm{\alpha}_{j}\in\mathbb{R}^{n},\hat{\bm{\alpha}}_{j}\in\mathbb{R}^{n}\end{ subarray}}~{}~{}~{}~{}\Big{\|}\sum_{j=1}^{m}\bm{\alpha}_{j}\bm{v}_{j}^{\top}- \bm{Y}\Big{\|}_{F}^{2}+\frac{\gamma}{2}\sum_{j=1}^{m}(\|\bm{u}_{j}\|_{2}^{2}+ \|\bm{v}_{j}\|_{2}^{2})$ (9) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}\hat{\bm{\alpha}}_{j}= \bm{X}\bm{u}_{j}\quad\text{and}\quad\bm{\alpha}_{j}=(\hat{\bm{\alpha}}_{j})_{+ },~{}~{}~{}~{}~{}\mathrm{for}~{}~{}~{}j=1,\ldots,m.$ $$

We introduce non-negative variables $\bm{\beta}_{j}:=\bm{\alpha}_{j}-\hat{\bm{\alpha}}_{j}$.
As both $\bm{\alpha}_{j}$ and $\bm{\beta}_{j}$ are non-negative valued, it can be shown that $\bm{\alpha}_{j}=(\hat{\bm{\alpha}}_{j})_{+}$ if and only if $\mathrm{trace}({\bm{\alpha}_{j}}\bm{\beta}_{j}^{\top})=0$, see Lemma 3.
As a result, the constraints in Eq. 9 can be equivalently rewritten in terms of $\bm{\alpha}_{j}$ and $\bm{\beta}_{j}$ as:

$$ $\bm{\alpha}_{j},\bm{\beta_{j}}\in\mathbb{R}_{+}^{n},\quad\mathrm{trace}({\bm{ \alpha}_{j}}\bm{\beta}_{j}^{\top})=0,\quad\text{and}\quad\bm{X}\bm{u}_{j}-\bm{ \alpha}_{j}+\bm{\beta}_{j}=\bm{0}_{n\times 1},~{}~{}~{}~{}~{}\mathrm{for}~{}~{ }~{}j=1,\ldots,m.$ (10) $$

The quadratic terms in the objective function in Eq. 9 and the trace constraint in Eq. 10 introduce challenges due to non-convexity.
We address this by lifting the problem into a higher-dimensional space, introducing new variables $\bm{\lambda}_{j}\in\mathbb{R}^{p}$ and $\bm{\Lambda}_{j}\in\mathbb{R}^{p\times p}$, where $p=2n+c+d$:

$$ $\displaystyle\bm{\lambda}_{j}\!:=\!\begin{bmatrix}\bm{\alpha}_{j}\[3.50006pt] \bm{\beta}_{j}\[3.50006pt] \bm{u}_{j}\[3.50006pt] \bm{v}_{j}\end{bmatrix}~{}~{}~{}\text{and}~{}~{}~{}\bm{\Lambda}_{j}\!:=\! \begin{bmatrix}\bm{\Lambda}_{\bm{\alpha}_{j}\bm{\alpha}_{j}^{\top}}&\bm{ \Lambda}_{\bm{\alpha}_{j}\bm{\beta}_{j}^{\top}}&\bm{\Lambda}_{\bm{\alpha}_{j} \bm{u}_{j}^{\top}}&\bm{\Lambda}_{\bm{\alpha}_{j}\bm{v}_{j}^{\top}}\[2.5pt] \bm{\Lambda}_{\bm{\beta}_{j}\bm{\alpha}_{j}^{\top}}&\bm{\Lambda}_{\bm{\beta}_{ j}\bm{\beta}_{j}^{\top}}&\bm{\Lambda}_{\bm{\beta}_{j}\bm{u}_{j}^{\top}}&\bm{ \Lambda}_{\bm{\beta}_{j}\bm{v}_{j}^{\top}}\[2.5pt] \bm{\Lambda}_{\bm{u}_{j}\bm{\alpha}_{j}^{\top}}&\bm{\Lambda}_{\bm{u}_{j}\bm{ \beta}_{j}^{\top}}&\bm{\Lambda}_{\bm{u}_{j}\bm{u}_{j}^{\top}}&\bm{\Lambda}_{ \bm{u}_{j}\bm{v}_{j}^{\top}}\[2.5pt] \bm{\Lambda}_{\bm{v}_{j}\bm{\alpha}_{j}^{\top}}&\bm{\Lambda}_{\bm{v}_{j}\bm{ \beta}_{j}^{\top}}&\bm{\Lambda}_{\bm{v}_{j}\bm{u}_{j}^{\top}}&\bm{\Lambda}_{ \bm{v}_{j}\bm{v}_{j}^{\top}}\end{bmatrix}~{}~{}~{}\text{such that}~{}~{}~{}\bm {\Lambda}_{j}=\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}.$ (11) $$

This allows us to express the objective as a convex function of $\bm{\Lambda}_{j}$:

$$ $\Big{\|}\sum_{j=1}^{m}\bm{P}_{\alpha}\bm{\Lambda}_{j}\bm{P}_{v}^{\top}-\bm{Y} \Big{\|}_{F}^{2}+\frac{\gamma}{2}\sum_{j=1}^{m}\langle\bm{P}_{u}^{\top}\bm{P}_ {u}+\bm{P}_{v}^{\top}\bm{P}_{v},\bm{\Lambda}_{j}\rangle.$ (12) $$

Similarly, we can reformulate the constraints in Eq. 10 as

$$ $\bm{P}_{\alpha\beta}\bm{\lambda}_{j}\geq 0,\quad\mathrm{trace}(\bm{P}_{\alpha} \bm{\Lambda}_{j}\bm{P}_{\beta}^{\top})=0,\quad\text{and}\quad\bm{M}\bm{\lambda }_{j}=\bm{0}_{n\times 1},~{}~{}~{}~{}~{}\mathrm{for}~{}~{}~{}j=1,\ldots,m.$ (13) $$

So far, we have shown that Eq. Nn-Train is equivalent to the following problem:

$$ $\displaystyle\min_{\bm{\lambda}_{j}\in\mathbb{R}^{p}}\,~{}\Big{\|}\sum_{j=1}^{ m}\bm{P}_{\alpha}\bm{\Lambda}_{j}\bm{P}_{v}^{\top}-\bm{Y}\Big{\|}_{F}^{2}+ \frac{\gamma}{2}\sum_{j=1}^{m}\langle\bm{P}_{u}^{\top}\bm{P}_{u}+\bm{P}_{v}^{ \top}\bm{P}_{v},\bm{\Lambda}_{j}\rangle$ ( Nn -Train-Lifted) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}\mathrm{for}~{}j=1,\ldots,m:$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{M}\bm {\lambda}_{j}=\bm{0}_{n\times 1}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\mathrm{ trace}(\bm{P}_{\alpha}\bm{\Lambda}_{j}\bm{P}_{\beta}^{\top})=0$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{ \Lambda}_{j}=\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\quad\text{and}\quad\bm{P} _{\alpha\beta}\bm{\lambda}_{j}\geq\bm{0}_{2n\times 1}.$ $$

Next, we will show that Eq. Cp-Nn is also equivalent to Eq. Nn-Train-Lifted when $m$ is sufficiently large.
Since $\bm{\Lambda}\in\mathcal{PSD}$, it follows that $\mathrm{trace}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=0$ and $\mathrm{trace}(\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})=0$ is equivalent to $\mathrm{trace}(\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})+\mathrm{trace
}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=0$.
By Lemma 2, we can find a positive integer $R$ such that any matrix $\bm{\Lambda}\in\mathcal{PSD}$ such that $\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{CP}$ can be factorized as

$$ $\bm{\Lambda}=\sum_{j=1}^{R}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}~{}~{}~{} \text{for some $\bm{\lambda}_{j}\in\mathbb{R}^{p}~{}~{}~{}$such that$~{}~{}~{} \bm{P}_{\alpha\beta}\bm{\lambda}_{j}\geq\bm{0}_{2n\times 1}$.}$ (14) $$

Given this decomposition, we can write

$$ $\mathrm{trace}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=\sum_{j=1}^{R}\mathrm{trace} \left(\bm{M}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\bm{M}^{\top}\right)\\ =\sum_{j=1}^{R}\|\bm{M}\bm{\lambda}_{j}\|^{2}.$ (15) $$

Therefore, $\mathrm{trace}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=0$ if and only if $\bm{M}\bm{\lambda}_{j}=\bm{0}_{n\times 1}$ for all $r=1,\ldots,R$.
Similarly,

$$ $\displaystyle\mathrm{trace}(\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})= \sum_{j=1}^{R}\mathrm{trace}\left(\bm{P}_{\alpha}\bm{\lambda}_{j}\bm{\lambda}_ {j}^{\top}\bm{P}_{\beta}^{\top}\right)=\sum_{j=1}^{R}\langle\bm{P}_{\alpha}\bm {\lambda}_{j},\bm{P}_{\beta}\bm{\lambda}_{j}\rangle.$ (16) $$

Note that $\langle\bm{P}_{\alpha}\bm{\lambda}_{j},\bm{P}_{\beta}\bm{\lambda}_{j}\rangle\geq
0$ for all $j$ since both $\bm{P}_{\alpha}\bm{\lambda}_{j}$ and $\bm{P}_{\beta}\bm{\lambda}_{j}$ are non-negative valued.
As a result, $\mathrm{trace}(\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})=0$ if and only if $\mathrm{trace}\left(\bm{P}_{\alpha}\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\bm{
P}_{\beta}^{\top}\right)=\langle\bm{P}_{\alpha}\bm{\lambda}_{j},\bm{P}_{\beta}
\bm{\lambda}_{j}\rangle=0$ for all $j=1,\ldots,R$.
Therefore, the formulation in Eq. Cp-Nn is equivalent to the following:

$$ $\displaystyle\min_{\bm{\lambda}_{j}\in\mathbb{R}^{p}}\,~{}\Big{\|}\sum_{j=1}^{ R}\bm{P}_{\alpha}\bm{\Lambda}_{j}\bm{P}_{v}^{\top}-\bm{Y}\Big{\|}_{F}^{2}+ \frac{\gamma}{2}\sum_{j=1}^{R}\langle\bm{P}_{u}^{\top}\bm{P}_{u}+\bm{P}_{v}^{ \top}\bm{P}_{v},\bm{\Lambda}_{j}\rangle$ (17) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}\mathrm{for}~{}j=1,\ldots,R:$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{M}\bm {\lambda}_{j}=\bm{0}_{n\times 1}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\mathrm{ trace}(\bm{P}_{\alpha}\bm{\Lambda}_{j}\bm{P}_{\beta}^{\top})=0$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{ \Lambda}_{j}=\bm{\lambda}_{j}\bm{\lambda}_{j}^{\top}\quad\text{and}\quad\bm{P} _{\alpha\beta}\bm{\lambda}_{j}\geq\bm{0}_{2n\times 1}.$ $$

This concludes the proof.
∎

###### Lemma 2 .

###### Fact 1 .

###### Proof of 1 .

###### Proof of Lemma 2 .

###### Lemma 3 .

###### Proof of Lemma 3 .

###### Proof of Theorem 1 .

## Appendix B Proof of Theorem 2

We first introduce the following lemmas that will be required for stating the proof of Theorem 2.

###### Lemma 4 .

Let $\bm{W}\in\mathbb{R}^{p\times p}$ be a positive semidefinite matrix of rank $k$, with a completely positive principal submatrix $\widehat{\bm{W}}\in\mathbb{R}^{q\times q}$ of CP-rank $K$.
Without loss of generality, we assume $\widehat{\bm{W}}$ is the leading principal submatrix of $\bm{W}$ given by $\widehat{\bm{W}}=\bm{P}\bm{W}\bm{P}^{\top}$ with $\bm{P}=\left[\bm{I}_{q\times q}~{}~{}\bm{0}_{q\times(p-q)}\right]$.
Then, there exists a probability measure $\nu$ and a corresponding random variable $\bm{w}$ such that $\bm{W}$ can be decomposed as

$$ $\bm{W}=\mathbb{E}_{\nu}[\bm{w}\bm{w}^{\top}]=\!\int_{\mathbb{R}^{p}}\!\!\bm{w} \bm{w}^{\top}d\nu(\bm{w}),$ (18) $$

such that $\bm{w}\in\mathbb{R}^{p}$ and $\bm{P}\bm{w}\geq\bm{0}$ a.s. with respect to $\nu$.

###### Proof of Lemma 4 .

Let us define the set $\Theta$ as follows:

$$ $\Theta:=\{\bm{W}\in\mathbb{R}^{p\times p}:\bm{W}\in\mathcal{PSD},~{}~{}\bm{P} \bm{W}\bm{P}^{\top}\!\in\mathcal{CP}\}.$ (19) $$

By Lemma 2, we can find a positive integer $R$ such that any matrix $\bm{W}\in\Theta$ can be factorized as

$$ $\bm{W}=\sum_{j=1}^{R}\bm{w}_{j}\bm{w}_{j}^{\top}~{}~{}~{}\text{for some $\bm{w }_{j}\in\mathbb{R}^{p}~{}~{}~{}$ such that $~{}~{}~{}\bm{P}\bm{w}_{j}\geq\bm{0 }$}.$ (20) $$

Thus, by performing a scaling: $\bm{w}_{j}\rightarrow\sigma_{j}\bm{w}_{j}$ for some $\sigma_{j}\in[0,1]$ such that $\sum_{j=1}^{R}\sigma_{j}^{2}=1$, we can equivalently express $\Theta$ as follows:

$$ $\Theta:=\Bigg{\{}\sum_{j=1}^{R}\sigma_{j}^{2}\bm{w}_{j}\bm{w}_{j}^{\top}:~{} \bm{w}_{j}\in\mathbb{R}^{p},~{}\bm{P}\bm{w}_{j}\geq\bm{0},~{}\sigma_{j}\in[0,1 ],~{}\forall j=1,\ldots,R,~{}~{}\text{and}~{}~{}\sum_{j=1}^{R}\sigma_{j}^{2}=1 \Bigg{\}}.$ (21) $$

We also define the set $\Gamma$ as follows:

$$ $\Gamma=\left\{\int_{\mathbb{R}^{p}}\!\!\bm{w}\bm{w}^{\top}d\nu(\bm{w})~{}:~{} \nu\mbox{ is a probability measure over }\bm{w},~{}\bm{P}\bm{w}\geq\bm{0}~{} \text{a.s. with respect to $\nu$}\right\}.$ (22) $$

It is clear that $\Theta\subseteq\Gamma$, because any element of $\Theta$, expressed as $\sum_{j=1}^{R}\sigma_{j}^{2}\bm{w}_{j}\bm{w}_{j}^{\top}$, can be represented by the measure $\sum_{j=1}^{R}\sigma_{j}^{2}\delta_{\bm{w}_{j}\bm{w}_{j}^{\top}}$, where $\delta_{\bm{w}_{j}\bm{w}_{j}^{\top}}$ is the Dirac delta measure centered at $\bm{w}_{j}\bm{w}_{j}^{\top}$.
The Dirac delta measure $\delta_{x}$ has support at the point $x$, assigning probability $1$ to sets containing $x$ and 0 0 otherwise.

To prove the opposite direction, that $\Gamma\subseteq\Theta$, we define an indicator function $I_{\Theta}(\cdot)$ given by

$$ $I_{\Theta}(\bm{W})=\begin{cases}0&\text{if }\bm{W}\in\Theta\\ +\infty&\text{if }\bm{W}\not\in\Theta\end{cases}$ (23) $$

Since $\Theta$ is a closed convex set, it is clear that $I_{\Theta}$ is a closed convex function.

Let $\bm{W}=\int_{\mathbb{R}_{p}}\bm{w}\bm{w}^{\top}d\nu(\bm{w})\in\Gamma$ for some probability measure $\nu$.
Applying Jensen’s inequality to gives the following:

$$ $I_{\Theta}(\bm{W})=I_{\Theta}\left(\int_{\mathbb{R}_{p}}\bm{w}\bm{w}^{\top}d \nu(\bm{w})\right)\leq\int_{\mathbb{R}^{p}}I_{\Theta}(\bm{w}\bm{w}^{\top})d\nu (\bm{w})=0,$ (24) $$

so $\bm{W}\in\Theta$. This completes the proof.
∎

We recall the following standard result from measure theory without proof, as it will be used in the subsequent analysis:

###### Fact 2 .

Let $f$ be an absolutely integrable function and $\nu$ be a measure defined on the measure space $(\mathbb{R}^{p},\mathcal{B}_{p})$.
If $f\geq 0$ almost surely, then $\int_{\mathbb{R}^{p}}f(\bm{x})d\nu(\bm{x})=0$ if and only if $f=0$ almost surely.

###### Proof of Theorem 2 .

We begin by rewriting Eq. Nn$\!\smallint$-Train in constrained form:

$$ $\displaystyle\min_{\nu:\mathcal{B}_{d}\times\mathcal{B}_{c}\rightarrow[0,1]} \left\|\int_{\mathbb{R}^{d}\times\mathbb{R}^{c}}\bm{s}\bm{v}^{\top}d\nu(\bm{u} ,\bm{v})-\bm{Y}\right\|_{F}^{2}+\frac{\gamma}{2}\int_{\mathbb{R}^{d}\times \mathbb{R}^{c}}(\|\bm{u}\|_{2}^{2}+\|\bm{v}\|_{2}^{2})d\nu(\bm{u},\bm{v})$ (25) $\displaystyle~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}\hat{\bm{s}}(\bm{u}, \bm{v})=\bm{X}\bm{u}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\bm{s}(\bm{u},\bm{v}) =(\hat{\bm{s}}(\bm{u},\bm{v}))_{+}.$ $$

We introduce non-negative variables $\bm{\alpha}(\bm{u},\bm{v}),\bm{\beta}(\bm{u},\bm{v})\in\mathbb{R}_{+}^{n}$ that splits $\bm{s}$ and $\hat{\bm{s}}$ as follows:

$$ $\bm{\alpha}(\bm{u},\bm{v})=\bm{s}(\bm{u},\bm{v})\quad\text{and}\quad\bm{\beta} (\bm{u},\bm{v})=\bm{s}(\bm{u},\bm{v})-\hat{\bm{s}}(\bm{u},\bm{v}).$ (26) $$

For the sake of brevity, we will use the following shorthand notation: $\bm{\alpha}:=\bm{\alpha}(\bm{u},\bm{v})$ and $\bm{\beta}:=\bm{\beta}(\bm{u},\bm{v})$.
Using the equivalence described in Lemma 3, we can reformulate Eq. 25 as follows:

$$ $\displaystyle\min_{\nu:\mathcal{B}_{d}\times\mathcal{B}_{c}\rightarrow[0,1]} \left\|\int_{\mathbb{R}^{d}\times\mathbb{R}^{c}}\bm{\alpha}\bm{v}^{\top}d\nu( \bm{u},\bm{v})-\bm{Y}\right\|_{F}^{2}+\frac{\gamma}{2}\int_{\mathbb{R}^{d} \times\mathbb{R}^{c}}(\|\bm{u}\|_{2}^{2}+\|\bm{v}\|_{2}^{2})d\nu(\bm{u},\bm{v})$ (27) $\displaystyle~{}~{}~{}~{}\mathrm{s.t.}~{}~{}~{}~{}~{}\bm{X}\bm{u}-\bm{\alpha}+ \bm{\beta}=\bm{0}_{n\times 1}$ $\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\,\mathrm{diag}({\bm{ \alpha}}\bm{\beta}^{\top})=\bm{0}_{n\times 1}.$ $$

Let $\nu(\bm{u},\bm{v})$ be a probability measure defined on the support of the variables $\bm{u}$ and $\bm{v}$. Then $\nu(\bm{u},\bm{v})$ represents a feasible solution to Eq. Nn$\!\smallint$-Train.
Define the random vector $\bm{\lambda}:=\bm{\lambda}(\bm{u},\bm{v})\in\mathbb{R}^{p}$ with $p=2n+c+d$ as $\bm{\lambda}:=\bm{\lambda}(\bm{u},\bm{v}):=[\bm{\alpha}^{\top}~{}~{}\bm{\beta}
^{\top}~{}~{}\bm{u}^{\top}~{}~{}\bm{v}^{\top}]^{\top}$.
Also, let $\bm{\Lambda}:=\mathbb{E}_{\nu}[\bm{\lambda}\bm{\lambda}^{\top}]$, where the expectation is with respect to the probability measure $\nu(\bm{u},\bm{v})$.
We can represent the objective function in Eq. 27 as follows:

$$ $\left\|\mathbb{E}_{\nu}[\bm{\alpha}\bm{v}^{\top}]-\bm{Y}\right\|_{F}^{2}+\frac {\gamma}{2}\mathbb{E}_{\nu}[\mathrm{Tr}(\bm{u}\bm{u}^{\top}+\bm{v}\bm{v}^{\top })],$ (28) $$

or equivalently, in terms of $\bm{\Lambda}$, as:

$$ $\left\|\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{v}^{\top}-\bm{Y}\right\|_{F}^{2}+ \frac{\gamma}{2}\langle\bm{P}_{u}^{\top}\bm{P}_{u}+\bm{P}_{v}^{\top}\bm{P}_{v} ,\bm{\Lambda}\rangle,$ (29) $$

which is exactly the objective function in Eq. Cp-Nn.

Then, similar to the proof of Theorem 1, we can reformulate the constraints in Eq. 27 as

$$ $\bm{P}_{\alpha\beta}\bm{\lambda}\geq\bm{0}_{2n\times 1},\quad\langle\bm{P}_{ \alpha}\bm{\lambda},\bm{P}_{\beta}\bm{\lambda}\rangle=0,\quad\text{and}\quad \bm{M}\bm{\lambda}=\bm{0}_{n\times 1}.$ (30) $$

Next, we will show that we can replace these with

$$ $\mathrm{trace}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=0,~{}~{}~{}\mathrm{trace}(\bm{ P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})=0,~{}~{}~{}\bm{\Lambda}\in \mathcal{PSD},~{}~{}~{}\text{and}~{}~{}~{}\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{ P}_{\alpha\beta}^{\top}\in\mathcal{CP},$ (31) $$

which is exactly the feasible set of Eq. Cp-Nn.

By Lemma 4, we can find a probability measure $\nu$ and a corresponding random variable $\bm{\lambda}$ such that any matrix $\bm{\Lambda}\in\mathcal{PSD}$ such that $\bm{P}_{\alpha\beta}\bm{\Lambda}\bm{P}_{\alpha\beta}^{\top}\in\mathcal{CP}$ can be factorized as

$$ $\bm{\Lambda}=\mathbb{E}_{\nu}[\bm{\lambda}\bm{\lambda}^{\top}]=\int_{\mathbb{R }^{p}}\bm{\lambda}\bm{\lambda}^{\top}d\nu(\bm{\lambda}),\quad\text{such that $ \quad\bm{P}_{\alpha\beta}\bm{\lambda}\geq\bm{0}_{2n\times 1}~{}$ a.s. with respect to $\nu$.}$ (32) $$

Given this decomposition, we can show that $\mathrm{trace}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=0$ if and only if $\bm{M}\bm{\lambda}=\bm{0}_{n\times 1}$ a.s., by using 2, since

$$ $\mathrm{trace}(\bm{M}\bm{\Lambda}\bm{M}^{\top})=\int_{\mathbb{R}^{p}}\mathrm{ trace}(\bm{M}\bm{\lambda}\bm{\lambda}^{\top}\bm{M}^{\top})\,d\nu(\bm{\lambda}) =\int_{\mathbb{R}^{p}}\|\bm{M}\bm{\lambda}\|^{2}\,d\nu(\bm{\lambda}).$ (33) $$

Similarly, we can show that $\mathrm{trace}(\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})=0$ if and only if $\langle\bm{P}_{\alpha}\bm{\lambda},\bm{P}_{\beta}\bm{\lambda}\rangle=0$ a.s., since

$$ $\mathrm{trace}(\bm{P}_{\alpha}\bm{\Lambda}\bm{P}_{\beta}^{\top})=\int_{\mathbb {R}^{p}}\mathrm{trace}(\bm{P}_{\alpha}\bm{\lambda}\bm{\lambda}^{\top}\bm{P}_{ \beta}^{\top})\,d\nu(\bm{\lambda})=\int_{\mathbb{R}^{p}}\langle\bm{P}_{\alpha} \bm{\lambda},\bm{P}_{\beta}\bm{\lambda}\rangle\,d\nu(\bm{\lambda}),$ (34) $$

and $\langle\bm{P}_{\alpha}\bm{\lambda},\bm{P}_{\beta}\bm{\lambda}\rangle\geq 0$ a.s. since $\bm{P}_{\alpha}\bm{\lambda}\geq\bm{0}_{n\times 1}$ and $\bm{P}_{\beta}\bm{\lambda}\geq\bm{0}_{n\times 1}$ a.s. by Eq. 32.
Therefore, the constraints in Eq. 30 imply those in Eq. 31, and the constraints in Eq. 31 imply Eq. 30 almost surely.
However, since the objective in Eq. 29 depends only on $\bm{\Lambda}=\mathbb{E}_{\nu}[\bm{\lambda}\bm{\lambda}^{\top}]$, and sets of measure zero do not affect $\bm{\Lambda}$, we can omit the almost surely condition on the constraints without affecting the solution.
∎

###### Lemma 4 .

###### Proof of Lemma 4 .

###### Fact 2 .

###### Proof of Theorem 2 .

## Appendix C Interpretation of the Solution to CP-NN

The solution to Cp-Nn can be interpreted as a moment matrix over the weights and activations with respect to the representing measure for the IWNN.
Based on the factorization in 32, the solution to Cp-Nn is $\bm{\Lambda}=\mathbb{E}_{\nu}[\bm{\lambda}\bm{\lambda}^{\top}]$ for some measure $\nu$ and corresponding random variable $\bm{\lambda}$:

$$ $\bm{\Lambda}=\mathbb{E}_{\nu}\begin{bmatrix}\bm{\alpha}\bm{\alpha}^{\top}&\bm{ \alpha}\bm{\beta}^{\top}&\bm{\alpha}\bm{u}^{\top}&\bm{\alpha}\bm{v}^{\top}\[2 .5pt] \bm{\beta}\bm{\alpha}^{\top}&\bm{\beta}\bm{\beta}^{\top}&\bm{\beta}\bm{u}^{ \top}&\bm{\beta}\bm{v}^{\top}\[2.5pt] \bm{u}\bm{\alpha}^{\top}&\bm{u}\bm{\beta}^{\top}&\bm{u}\bm{u}^{\top}&\bm{u}\bm {v}^{\top}\[2.5pt] \bm{v}\bm{\alpha}^{\top}&\bm{v}\bm{\beta}^{\top}&\bm{v}\bm{u}^{\top}&\bm{v}\bm {v}^{\top}\end{bmatrix}$ (35) $$

Suppose we solve Eq. Cp-Nn to obtain the solution $\bm{\Lambda}$ and consider the submatrix $\mathbb{E}_{\nu}[\bm{u}\bm{u}^{\top}]$.
Using this, we can compute the expected kernel matrix $\mathbf{K}[\bm{X}_{0},\bm{X}_{1}]$ over the measure $\nu$, where $\bm{X}_{0}\in\mathbb{R}^{n_{0}\times d}$ and $\mathbf{X}_{1}\in\mathbb{R}^{n_{1}\times d}$ are two input datasets of size $n_{0}$ and $n_{1}$, as $\mathbf{K}[\bm{X}_{0},\bm{X}_{1}]=\bm{X}_{0}\mathbb{E}_{\nu}[\bm{u}\bm{u}^{
\top}]\bm{X}_{1}^{\top}=\mathbb{E}_{\nu}[(\bm{X}_{0}\mathbf{u})(\bm{X}_{1}\bm{
u})^{\top}]$.
It is thus possible to extract the expected kernel matrix of the infinite-width (or wide) neural network from the solution of Eq. Cp-Nn, which stores the learned features.

## Appendix D Additional Details on Numerical Experiments

We used SGD with small step sizes and trained until convergence. We observed that this approach provides a better approximation of the global solution. Based on this observation, we tuned the learning rate (LR) and the number of iterations for each dataset as:

- •
Random: Initial LR = $10^{-5}$; # of iterations = 500K.
- •
Spiral: Initial LR = $10^{-3}$; # of iterations = 8K.
- •
Iris: Initial LR = $10^{-6}$; # of iterations = 2M.
- •
Ionosphere: Initial LR = $10^{-6}$; # of iterations = 2M for $\gamma=0.1$ and 5M for $\gamma=0.01$.
- •
Pima Indians Diabetes: Initial LR = $10^{-8}$; # of iterations = 5M for $\gamma=0.1$ and 6M for $\gamma=0.01$.
- •
Bank Notes Authentication: Initial LR = $10^{-6}$; # of iterations = 5M.

The SGD training curves for Random and Spiral with varying number of hidden neurons, ranging from 5 to 300, are shown in Figure 2.

Figure: (a) Random Dataset, $\gamma=0.1$.
Refer to caption: x2.png

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: - can the authors comment on the runtime of solving the SDP?
- can the authors provide additional discussion about what kinds of insights can be gained by applying the proposed framework?
soundness: 3
strengths: This is an interesting paper that derives an equivalence between infinite-width RELU network training and solving a certain convex copositive program. This work contributes to the growing literature relating neural network training and convex optimization. The empirical results are also promising. Overall, I think this is an interesting work that provides valuable insight.
summary: In this paper, the authors propose a novel model for infinite-width neural network training based on copositive programming. In theorems 1 and 2, the authors carefully demonstrate that the training of an RELU network of sufficient width using the MSE loss is equivalent to a certain convex copositive program. Next, the authors derive a semidefinite relaxation of the proposed copositive program and an associated rounding scheme to recover weights from the solution to the SDP. On a variety of synthetic and small real datasets, the authors demonstrate that their method recovers neural network weights that achieve smaller MSE (table 2) compared to networks trained via SGD and competitive predictive performance compared to NTK-based kernel methods.
weaknesses: The empirical results seem somewhat weak to me. The authors acknowledge the similarity of their work to earlier work relating copositive programming and RELU network training. Although the existing methods make additional assumptions on the data distribution, a numerical comparison to prior work (e.g. the approximation ratio) would be beneficial in understanding how the proposed framework compares to earlier work in practice or further discussion on the applications of their technique.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: I have some concerns on the empirical evaluation. MSE is not a suitable loss function for classification problems, a selection of regression tasks would have been more representative. Regardless, the hyperparameters chosen for the SGD baseline are quite strange (especially for PIMA Indians). Finally, the dataset is divided in train-test split, this penalizes models that overfit on the training data. However, if the goal is to evaluate how well a model is solving the given optimization problem, it should be evaluated on the training set.

In most cases, the addition of a bias term in the linear projections significantly improves performance, especially for 2-layer ReLU networks. In fact, running the proposed SGD baselines with the addition of a bias term improves performance across the board (gaining up to +15% accuracy). Given this, the lack of a bias term and its implications should probably be noted in the text.

I don't understand the results in Table 2. The metric is the loss, lower mean better, i.e. SDP-NN significantly improves on SDG. However the way the approximation ratio is computed and is presented in the main text implies that the SDP-NN performs worse than SGD.

Some equation are labelled instead of being numbered. I personally find this beaks the flow of the main text as it is not immediately clear that it is a reference to an equation (as most people are used to see numbers). If you want to refer to some equation by a meaningful name, you can use the name in the definition, i.e. instead of "... in (CP-NN)" consider using something like "... in the CP-formulation (equation 2)".
If you decide to keep the labels, at least try to use them to refer to the equations only, not to the concepts they are tied to. For instance in "While the case of countably infinite number of hidden neurons – represented by an infinite sum in (NN) – is captured as a special case of (NN∫) with a discrete probability measure, the (NN∫) model ..." you use (NN∫) to denote both the infinite width model concept as a whole, and its equation.
soundness: 3
strengths: The paper is clear, and reads well.
The supplementary material provides the code to reproduce the results.
Section 2 provides a complete yet brief background on the relevant optimization topics and concepts.
summary: The paper investigates the training of shallow NN using convex optimization techinques.
It provides a brief background on Semi-Definite programming and Copositive Programming. Then it provides a formulation for the training of a NN as a CP problem (for a sufficiently wide 2-layer net with ReLU activations), which is then relaxed to make it computationally feasible.
The tightness of the relaxation is evaluated empirically.
weaknesses: The paper does not consider bias terms in the linear layers.
The tightness of the proposed relaxation is evaluated only empirically.
There is no convergence guarantee for the TOS rounding step.
The rounding step is performed using the critical width, however in practice this is unfeasible.
Time complexity is not considered in the evaluations.
The empirical evaluation is performed on classification tasks, using L2 loss.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: I suggest that the authors should address the weaknesses point out above.
soundness: 3
strengths: - The idea of the proposed method based on convex completely positive program and semidefinite relaxation is interesting. 

- The presentation of the paper is fairly clear.
summary: This paper proposes a novel training method for two-layer neural networks based on a formulation of a convex completely positive program. The authors also propose a semidefinite relaxation which ensures that the problem can be solved in polynomial time. Experiments on several simple synthetic and real datasets are conducted to demonstrate the performance of the proposed method in comparison to the Neural Network Gaussian Process and Neural Tangent Kernel methods.
weaknesses: - The proposed method is only for training wide two-layer neural networks with ReLU activation function. It seems that extending the method to deeper networks is non-trivial, which limits the practical value of the proposed method. It would be beneficial if the authors could discuss the possibility to implement the proposed methods to train modern deep neural networks.

- As the authors have commented, the problem (Cp-Nn) is NP-hard due to the complete positivity constraint, and as far as I can see, the corresponding (Sdp-Nn) is not guaranteed to find the actual solution to (Cp-Nn). It would be helpful if the authors can add discussions about this point.

- To my knowledge, the Neural Network Gaussian Process and Neural Tangent Kernel methods are also settings (mainly about parameter/initialization scale and network width) under which SGD training of neural networks is almost equivalent to a convex optimization procedure. It is widely agreed that NNGP and NTK cannot fully explain the success of deep learning. These models only characterize scenarios where the obtained neural networks have weights very close to their random initialization. As a result, they fail to explain how deep networks learn effective feature representations. Therefore, even if the proposed method achieves competitive results compared with NNGP and NTK, it does not necessarily mean that the proposed method can compete with actual neural network training beyond the “NTK regime”.

- The experiments are not representative enough. All datasets considered in this paper are very small datasets, and the performance of the proposed method on these datasets are insufficient to demonstrate its advantage. I think the authors should at least include results on datasets such as MNIST and CIFAR-10. 

- The experiment setup requires some clarification. As mentioned above, NNGP and NTK represent specific regimes in which SGD training of neural networks may fall [1,2,3,4]. Therefore, it is not reasonable to directly list “SGD”, “NNGP”, “NTK” as different training methods. Discussion how the SGD setting differs from NNGP and NTK and why they are not equivalent according to [1,2,3,4] is necessary to explain the experiment setting clearly.

- The authors should also discuss the relationship between the proposed method and the study in [5].


[1] Simon Du, Xiyu Zhai, Barnabas Poczos, and Aarti Singh. Gradient descent provably optimizes over-parameterized neural networks. ICLR 2019.

[2] Zeyuan Allen-Zhu, Yuanzhi Li, and Zhao Song. A convergence theory for deep learning via over-parameterization. ICML 2019.

[3] Simon Du, Jason Lee, Haochuan Li, Liwei Wang, and Xiyu Zhai. Gradient descent finds global minima of deep neural networks. ICML 2019.

[4] Difan Zou, Yuan Cao, Dongruo Zhou, and Quanquan Gu. Gradient descent optimizes over-parameterized deep ReLU networks." Machine learning 2020

[5] Cong Fang, Yihong Gu, Weizhong Zhang, and Tong Zhang. Convex formulation of overparameterized deep neural networks. IEEE Transactions on Information Theory 2022

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: - Q1: Do you expect Theorem 1 to be applicable in theoretical analyses of 2-layer NNs in different contexts? 

- Q2:  Assume an oracle would give us a solution to the CP-NN problem directly. What is the computational complexity of an algorithm that would recover a solution to the NN-TRAIN problem? 

- Q3: Regarding Theorem 2, are there solutions to the NN-integral-train that are not discrete? 

- Q4: What would be the training loss of the recovered NN-TRAIN solution in Table 2? I am asking to see how far off would the projection back into feasible set be from the actual solution to NN-TRAIN. 

- Q5: What makes the bank notes dataset so different from the others that your method is so significantly worse than all the others, while being pretty good on all the other datasets? 

**Summary:** Despite several flaws outlined above, I still think this paper meets acceptance threshold due to the novel, compact and parameter-light convex formulation of the 2-layer NN training which does not rely on enumerating all the activation patterns of the intermediate layer (as in the previous literature). This result is of theoretical interest and, though simple in nature, it is novel and insightful. Therefore, I recommend acceptance. However, I condition this opinion as well as the score on the satisfying discussion period.

**After the rebuttal:** The authors satisfactorily answered my questions and discussed my concerns. I am keeping my original score.
soundness: 4
strengths: - S1: The strongest part of the paper is, in my opinion, the theoretical contribution (Theorem 1). To the best of my knowledge, this convex formulation significantly improves the number of parameters (as a function of the number of samples and input/output dimension) in the convex program from high-degree polynomial [1], to quadratic. The price the authors pay is the minimal width required to achieve this equivalence, which is also quadratic in the number of samples, as opposed to linear in previous works. Nevertheless, this neat and parameter-light reformulation is of independent theoretical interest and I believe it can be used as a tool in theoretical analyses of 2-layer NNs outside the scope of this paper. I checked the proof of the result and everything seems to be all right with only a small issue that I will discuss in typos. 

- S2: The paper is overall well written, clear and interesting. 

[1] Sahiner, Arda, et al. "Vector-output relu neural network problems are copositive programs: Convex analysis of two layer networks and polynomial-time algorithms." arXiv preprint arXiv:2012.13329 (2020).
summary: The paper studies 2-layer ReLU neural networks trained with MSE loss and $l_2$ regularization. The main contribution of the paper is in showing that beyond certain critical width (of at most quadratic dependence on the number of samples), the standard optimization problem is equivalent to a convex problem over a matrix of quadratic number of parameters in number of samples and input/output dimensions. The optimal values of this problem match.

The equivalence further applies to infinite-width 2-layer ReLU networks with analogical loss function. The optimal solutions of the convex problem can be used to construct optimal solutions of the infinite-width problem. 

Finally, the paper proposes a computationally efficient algorithm that approximates the convex, yet copositive program. This new, semidefinite program is solvable in polynomial time and the solution of the original problem can be recovered from the solution of the SDP problem by a rounding scheme described in the paper, using existing approaches such as the three-operator splitting method. 

The authors then evaluate this method as well as the tightness of the SDP relaxation on some tabular datasets of small size, demonstrating (within these datasets) that the SDP relaxation is rather tight and the practical optimization algorithm works comparably to known optimization methods, despite using number of parameters that is not theoretically guaranteed to work.
weaknesses: - W1: The practical contribution of this paper is limited. Even if we ignored the fact that the practical algorithm introduced in this paper is only suitable to train 2-layer neural networks, the main issue remains -- there already are algorithms in the literature that can obtain approximate solutions of the convex re-formulations using quadratic time per-iteration in the number of samples (see for instance [1]). Moreover, the experimental evaluation is not very convincing. The tested datasets are very small so the results are not very convincing. It is known (both theoretically and empirically), that gradient-based methods work well also with widths linear in the number of samples (see for instance [2], which considers only marginally different setting). Whether your method works well with R constant (the linear dependence on $n$ already comes from lifting into the higher-dimensional space) or at least sub-quadratic in $n$ is questionable (and not addressed in the paper). Moreover, the performance is comparable with already existing and well-established optimization methods (as demonstrated in the paper). 

- W2: While it is written in the paper that the algorithm scales quadratically with the number of samples (which is true in the plain CP-NN formulation), recovering a solution to the original problem actually requires cubic number of parameters and thus also time. This can be seen in equation (3), where it is clear that already the dimension of the $\lambda_j$ vectors scales with $n$ and $R$ scales quadratically with $n$, thus in total the dependence is cubic. Even if an oracle gave us *directly* a solution to CP-NN problem, recovering a solution to the original problem would still require $n^3$ number of parameters. 

- W3 (minor): The computational complexities of all the involved procedures should be explicitly discussed in the paper.

**Small issues:**

- Lemma 3 is wrong, the $w$ can have positive values where $\hat W$ has negative values. A correct version of this lemma removes $\lambda$ completely from the statement and replaces the $\lambda$ in the if and only if statement with $w$. This would also be consistent with the later use of this lemma in the actual proof of Theorem 1. 

- Typo: line 176 you refer to Lemma 4 but should be Lemma 2. 
- Typo: In line 167 index should be also added to the $u, v$ variables. 
- Typo: across the entire proof of theorem 1 including the lemmas that precede it, you use $W$ and $\Lambda$ supposedly to denote the same object. Please only use one of these (preferably $\Lambda$).

[1] Bai, Yatong, Tanmay Gautam, and Somayeh Sojoudi. "Efficient global optimization of two-layer relu networks: Quadratic-time algorithms and adversarial training." SIAM Journal on Mathematics of Data Science 5.2 (2023): 446-474.

[2] Bombari, Simone, Mohammad Hossein Amani, and Marco Mondelli. "Memorization and optimization in deep neural networks with minimum over-parameterization." Advances in Neural Information Processing Systems 35 (2022): 7628-7640.

Meta-review: The manuscript proposes a convex formulation for training two-layer neural networks. While the methodology seems to be only limited to simple architectures, the reviewers overall have a positive evaluation of the manuscript and find the approach interesting. Therefore, the meta-reviewer recommends acceptance of the manuscript.

Decision: Accept (Poster)