# Demystifying Batch Normalization in ReLU Networks: Equivalent Convex
Optimization Models and Implicit Regularization

Tolga Ergen, Arda Sahiner11footnotemark: 1, Batu Ozturkler, John Pauly,
Morteza Mardani & Mert Pilanci

Department of Electrical Engineering
Stanford University
Stanford, CA 94305, USA
<EMAIL_ADDRESS>
Equal Contribution

###### Abstract

Batch Normalization (BN) is a commonly used technique to accelerate and
stabilize training of deep neural networks. Despite its empirical success, a
full theoretical understanding of BN is yet to be developed. In this work, we
analyze BN through the lens of convex optimization. We introduce an analytic
framework based on convex duality to obtain exact convex representations of
weight-decay regularized ReLU networks with BN, which can be trained in
polynomial-time. Our analyses also show that optimal layer weights can be
obtained as simple closed-form formulas in the high-dimensional and/or
overparameterized regimes. Furthermore, we find that Gradient Descent provides
an algorithmic bias effect on the standard non-convex BN network, and we
design an approach to explicitly encode this implicit regularization into the
convex objective. Experiments with CIFAR image classification highlight the
effectiveness of this explicit regularization for mimicking and substantially
improving the performance of standard BN networks.

###  1 Introduction

Deep neural networks have achieved dramatic progress in the past decade. This
dramatic progress largely hinged on improvements in terms of optimization
techniques. One of the most prominent recent optimization techniques is Batch
Normalization (BN) (Ioffe & Szegedy, 2015). BN is an operation that can be
introduced in between layers to normalize the output of each layer and it has
been shown to be extremely effective in stabilizing and accelerating training
of deep neural networks. Hence, it became standard in numerous state-of-the-
art architectures, e.g., ResNets (He et al., 2016). Despite its empirical
success, it is still theoretically elusive why BN is extremely effective for
training deep neural networks. Therefore, we investigate the mechanisms behind
the success of BN through convex duality.

####  1.1 Related Work

Batch Normalization: One line of research has focused particularly on
alternatives to BN, such as Layer Normalization (Ba et al., 2016), Instance
Normalization (Ulyanov et al., 2016), Weight Normalization (Salimans & Kingma,
2016), and Group Normalization (Wu & He, 2018). Although these techniques
achieved performance competitive with BN, they do not provide any theoretical
insight about its empirical success.

Another line of research studied the effects of BN on neural network training
and identified several benefits. For example, Im et al. (2016) showed that
training deep networks with BN reduces dependence on the parameter
initialization. Wei et al. (2019) analyzed BN via mean-field theory to
quantify its impact on the geometry of the optimization landscape. They
reported that BN flattens the optimization landscape so that it enables the
use of larger learning rates. In addition, Bjorck et al. (2018); Santurkar et
al. (2018); Arora et al. (2018) showed that networks trained with BN achieve
faster convergence and generalize better. Furthermore, Daneshmand et al.
(2020) proved that BN avoids rank collapse so that gradient-based algorithms,
e.g., Stochastic Gradient Descent (SGD), are able to effectively train deep
networks. Even though these studies are important to understand the benefits
of BN, they fail to provide a theoretically complete characterization for
training deep networks with BN.

Convex Neural Networks:  Recently, a series of papers (Pilanci & Ergen, 2020;
Ergen & Pilanci, 2020a; 2021a; 2021b; Sahiner et al., 2021a; b) studied ReLU
networks through the lens of convex optimization theory. Particularly, Pilanci
& Ergen (2020) introduced exact convex representations for two-layer ReLU
networks, which can be trained in polynomial-time via standard convex solvers.
However, this work is restricted to two-layer fully connected networks with
scalar outputs. Later on, Ergen & Pilanci (2021a) first extended this approach
to two-layer scalar output Convolutional Neural Networks (CNNs) with average
and max pooling and provided further improvements on the training complexity.
These results were extended to two-layer fully convolutional networks and two-
layer networks with vector outputs (Sahiner et al., 2021a). However, these
convex approaches are restricted to two-layer ReLU networks without BN, thus,
do not reflect the exact training framework in practice, i.e., regularized
deep ReLU networks with BN. Unlike these prior works, we not only analyze
training problems with batch normalization but also extend the analysis to
deeper architectures in Section 4. Furthermore, we prove that the optimal
layer weight can be obtained as closed-form solutions in the high-dimensional
and/or overparameterized regimes (see Theorem 2.1, 2.4, 4.1, 4.2, G.2).

####  1.2 Our Contributions

  * •

We introduce an exact convex framework to explicitly characterize optimal
solutions to ReLU network training problems with weight-decay regularization
and BN. Thus, we obtain closed-form solutions for the optimal layer weights in
the high-dimensional and overparameterized regime.

  * •

We prove that regularized ReLU network training problems with BN can be
equivalently stated as a finite-dimensional convex problem. As a corollary, we
also show that the equivalent convex problems involve whitened data matrices
unlike the original non-convex training problem. Hence, using convex
optimization, we reveal an implicit whitening effect introduced by BN.

  * •

We demonstrate that GD applied to BN networks provides an implicit
regularization effect by learning high singular value directions of the
training data more aggressively, whereas this regularization is absent for GD
applied to the equivalent whitened data formulation. We propose techniques to
explicitly regularize BN networks to capture this implicit regularization
effect.

  * •

Unlike previous studies, our derivations extend to deep ReLU networks with BN,
CNNs, BN after ReLU111Presented in Appendix G., vector output networks, and
arbitrary convex loss functions.

####  1.3 Preliminaries

Notation: We denote matrices and vectors as uppercase and lowercase bold
letters, respectively, where a subscript indicates a certain element or
column. We use $\mathbf{0}$ (or $\mathbf{1}$) to denote a vector or matrix of
zeros (or ones), where the sizes are appropriately chosen depending on the
context. We also use $\mathbf{I}_{n}$ to denote the identity matrix of size
$n$. To represent Euclidean, Frobenius, and nuclear norms, we use
$\|\cdot\|_{2}$, $\|\cdot\|_{F}$, and $\|\cdot\|_{*}$, respectively. Lastly,
we denote the elementwise 0-1 valued indicator function and ReLU activation as
$\mathbbm{1}[x\geq 0]$ and $\left(x\right)_{+}=\max\\{x,0\\}$, respectively.

We222All the proofs and some extensions are presented in Appendix. consider an
$L$-layer ReLU network with layer weights
$\mathbf{W}^{(l)}\in\mathbb{R}^{m_{l-1}\times m_{l}}$, where $m_{0}=d$ and
$m_{L}=C$ are the input and output dimensions, respectively. Given a training
data $\mathbf{X}\in\mathbb{R}^{n\times d}$ and a label matrix
$\mathbf{Y}\in\mathbb{R}^{n\times C}$, we particularly focus on the following
regularized training problem

$\displaystyle\min_{\theta\in\Theta}\mathcal{L}(f_{\theta,L}(\mathbf{X}),\mathbf{Y})+\beta\mathcal{R}(\theta),$
(1)

where we compactly represent the parameters as
$\theta:=\\{\mathbf{W}^{(l)},\bm{\gamma}_{l},\bm{\alpha}_{l}\\}_{l=1}^{L}$ and
the corresponding parameter space as
$\Theta:=\\{\\{\mathbf{W}^{(l)},\bm{\gamma}^{(l)},\bm{\alpha}^{(l)}\\}_{l=1}^{L}:\mathbf{W}^{(l)}\in\mathbb{R}^{m_{l-1}\times
m_{l}},\bm{\gamma}_{l}\in\mathbb{R}^{m_{l}},\bm{\alpha}_{l}\in\mathbb{R}^{m_{l}},\forall
l\in[L]\\}$. We note that $\bm{\gamma}_{l}$ and $\bm{\alpha}_{l}$ are the
parameters of the BN operator, for which we discuss the details below. In
addition, $\mathcal{L}(\cdot,\cdot)$ is an arbitrary convex loss function,
including squared, hinge, and cross entropy loss, and $\mathcal{R}(\cdot)$ is
the regularization function for the layer weights with the tuning parameter
$\beta>0$. We also compactly define the network output as

$\displaystyle
f_{\theta,L}(\mathbf{X}):=\mathbf{A}^{(L-1)}\mathbf{W}^{(L)},\;\text{ where
}\mathbf{A}^{(l)}:=\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(l-1)}\mathbf{W}^{(l)}\right)\right)_{+}$

we denote the $l^{th}$ layer activations as
$\mathbf{A}^{(l)}\in\mathbb{R}^{n\times m_{l}}$, and
$\mathbf{A}^{(0)}=\mathbf{X}$. Here,
${\mathrm{BN}}_{\gamma,\alpha}\left(\cdot\right)$ represents the BN operation
introduced in Ioffe & Szegedy (2015) and applies matrices columnwise.

######  Remark 1.1.

Note that above we use BN before ReLU activations, which is common practice
and consistent with the way introduced in Ioffe & Szegedy (2015). However, BN
can be placed after ReLU as well, e.g., Chen et al. (2019), and thus in
Section G of Appendix, we will also consider architectures where BN layers are
placed after ReLU.

For a layer with weight matrix $\mathbf{W}^{(l)}\in\mathbb{R}^{m_{l-1}\times
m_{l}}$ and arbitrary batch of activations denoted as
$\mathbf{A}^{(l-1)}_{b}\in\mathbb{R}^{s\times m_{l-1}}$, BN applies to each
column $j$ independently as follows

$\displaystyle{\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(l-1)}_{b}\mathbf{w}^{(l)}_{j}\right):=\frac{(\mathbf{I}_{s}-\frac{1}{s}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(l-1)}_{b}\mathbf{w}^{(l)}_{j}}{\|(\mathbf{I}_{s}-\frac{1}{s}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(l-1)}_{b}\mathbf{w}^{(l)}_{j}\|_{2}}\gamma^{(l)}_{j}+\frac{\alpha^{(l)}_{j}\mathbf{1}}{\sqrt{n}},$
(2)

where $\bm{\gamma}^{(l)}$ and $\bm{\alpha}^{(l)}$ are trainable parameters
that scale and shift the normalized value. In this work, we focus on the full-
batch case, i.e,
$\mathbf{A}^{(l-1)}_{b}=\mathbf{A}^{(l-1)}\in\mathbb{R}^{n\times m_{l}}$. This
corresponds to training the network with GD as opposed to mini-batch SGD. We
note that our empirical findings with GD indicate identical if not better
performance compared to the mini-batch case, which is also consistent with the
previous studies Lian & Liu (2019); Summers & Dinneen (2020).

Throughout the paper, we consider a regression framework with squared loss and
standard weight-decay regularization. Extensions to general convex loss
functions are presented in Appendix B.1. Moreover, below, we first focus on
scalar outputs i.e. $C=1$, and then extend it to vector outputs.

####  1.4 Overview of Our Results

Here, we provide an overview of our main results. To simplify the notation, we
consider $L$-layer ReLU networks with scalar outputs, i.e., $m_{L}=C=1$ thus
the label vector is $\mathbf{y}\in\mathbb{R}^{n}$, and extend the analysis to
vector output networks with the label matrix $\mathbf{Y}\in\mathbb{R}^{n\times
C}$ in the next sections.

The regularized training problems for an $L$-layer network with scalar output
and BN is given by

$\displaystyle
p_{L}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|f_{\theta,L}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\frac{\beta}{2}\sum_{l=1}^{L}\left(\left\|\bm{\gamma}^{(l)}\right\|_{2}^{2}+\left\|\bm{\alpha}^{(l)}\right\|_{2}^{2}+\left\|\mathbf{W}^{(l)}\right\|_{F}^{2}\right),$
(3)

where we use $\bm{\gamma}^{(L)}=\bm{\alpha}^{(L)}=\mathbf{0}$ as dummy
variables for notational simplicity.

######  Lemma 1.1.

The problem in (3) is equivalent to the following optimization problem

$\displaystyle\min_{\theta\in\Theta_{s}}\frac{1}{2}\left\|f_{\theta,L}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\beta\left\|\mathbf{w}^{(L)}\right\|_{1},$
(4)

where
$\Theta_{s}:=\\{\theta\in\Theta:{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}=1,\,\forall
j\in[m_{L-1}]\\}$.

Using the equivalence in Lemma 1.1, we now take dual of (4) with respect to
the output layer weights $\mathbf{w}^{(L)}$ to obtain 333Details regarding the
derivation of the dual problem are presented in Appendix B.3.

$\displaystyle p_{L}^{*}\geq
d_{L}^{*}:=\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{\top}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}\right)\right)_{+}\right|\leq\beta.$
(5)

Since the original formulation in (3) is a non-convex optimization problem,
any solution $\mathbf{v}$ in the dual domain yields a lower bound for the
primal problem, i.e., $p_{L}^{*}\geq d_{L}^{*}$. In this paper, we first show
that strong duality holds in this case, i.e., $p_{L}^{*}=d_{L}^{*}$, and then
derive an exact equivalent convex formulation for the non-convex problem (3).
Furthermore, we even obtain closed-form solutions for the layer weights in
some cases so that there is no need to train a network in an end-to-end
manner.

###  2 Two-layer Networks

In this section, we analyze two-layer networks. Particularly, we first
consider a high-dimensional regime, where $n\leq d$, and then extend the
analysis to arbitrary data matrices by deriving an equivalent convex program.
We also extend the derivations to vector output networks.

We now consider the regularized training problem for a two-layer network with
scalar output and BN. Using the equivalence in Lemma 1.1, the problem can be
stated as

$\displaystyle
p_{2}^{*}=\min_{\theta\in\Theta_{s}}\frac{1}{2}\left\|f_{\theta,2}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\beta\left\|\mathbf{w}^{(2)}\right\|_{1}.$
(6)

We then take dual of (6) with respect to the output weights $\mathbf{w}^{(2)}$
to obtain the following dual problem

$\displaystyle p_{2}^{*}\geq d_{2}^{*}=$
$\displaystyle\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{T}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{X}\mathbf{w}^{(1)}\right)\right)_{+}\right|\leq\beta,$
(7)

where
$\Theta_{s}=\\{\theta\in\Theta:\mathbf{w}^{(1)}\in\mathbb{R}^{d},{\gamma^{(1)}}^{2}+{\alpha^{(1)}}^{2}=1\\}$.
Notice that here we drop the hidden neuron index $j\in[m_{1}]$ since the dual
constraint scans all possible parameters in the continuous set $\Theta_{s}$.

####  2.1 High-dimensional Regime ($n\leq d$)

In the sequel, we prove that the aforementioned dual problem is further
simplified such that one can find an optimal solution (6) in closed-form as
follows.

######  Theorem 2.1.

Suppose $n\leq d$ and $\mathbf{X}$ is full row-rank, then an optimal solution
to (6) is

$\displaystyle\left({\mathbf{w}^{(1)}_{j}}^{*},{{w}^{(2)}_{j}}^{*}\right)=\left(\mathbf{X}^{\dagger}\left((-1)^{j}\mathbf{y}\right)_{+},(-1)^{j}\left(\|\left((-1)^{j}\mathbf{y}\right)_{+}\|_{2}-\beta\right)_{+}\right)$
$\displaystyle\begin{bmatrix}{\gamma^{(1)}_{j}}^{*}\\\
{\alpha^{(1)}_{j}}^{*}\end{bmatrix}=\frac{1}{\|\left((-1)^{j}\mathbf{y}\right)_{+}\|_{2}}\begin{bmatrix}\|\left((-1)^{j}\mathbf{y}\right)_{+}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\left((-1)^{j}\mathbf{y}\right)_{+}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\left((-1)^{j}\mathbf{y}\right)_{+}\end{bmatrix}$

for $j=1,2$. Therefore, strong duality holds, $p_{2}^{*}=d_{2}^{*}$.

####  2.2 Exact Convex Formulation

To obtain an equivalent convex formulation to (6), we first introduce a notion
of hyperplane arrangements; see also Pilanci & Ergen (2020). We first define a
diagonal matrix as
$\mathbf{D}:=\text{diag}(\mathbbm{1}[\mathbf{X}\mathbf{w}\geq 0])$ for an
arbitrary $\mathbf{w}\in\mathbb{R}^{d}$. Therefore, the output of a ReLU
activation can be equivalently written as
$\left(\mathbf{X}\mathbf{w}\right)_{+}=\mathbf{D}\mathbf{X}\mathbf{w}$
provided that $\mathbf{D}\mathbf{X}\mathbf{w}\geq 0$ and
$\left(\mathbf{I}_{n}-\mathbf{D}\right)\mathbf{X}\mathbf{w}\leq 0$ are
satisfied. We can define these two constraints more compactly as
$\left(2\mathbf{D}-\mathbf{I}_{n}\right)\mathbf{X}\mathbf{w}\geq 0$. We now
denote the cardinality of the set of all possible diagonal matrices as $P$,
and obtain the following upperbound $P\leq 2r(e(n-1)/r)^{r}$, where
$r:=\mbox{rank}(\mathbf{X})\leq\min(n,d)$. For the rest of the paper, we
enumerate all possible diagonal matrices (or hyperplane arrangements) in an
arbitrary order and denote them as $\\{\mathbf{D}_{i}\\}_{i=1}^{P}$444We
provide more details in Appendix B.4.

Based on the notion of hyperplane arrangement, we now introduce an equivalent
convex formulation for (6) without an assumption on the data matrix as in the
previous section.

######  Theorem 2.2.

For any data matrix $\mathbf{X}$, the non-convex training problem in (6) can
be equivalently stated as the following finite-dimensional convex program

$\displaystyle\min_{\mathbf{s}_{i},\mathbf{s}_{i}^{\prime}\in\mathbb{R}^{r+1}}\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}(\mathbf{s}_{i}-\mathbf{s}_{i}^{\prime})-\mathbf{y}\right\|_{2}^{2}+\beta\sum_{i=1}^{P}(\|\mathbf{s}_{i}\|_{2}+\|\mathbf{s}_{i}^{\prime}\|_{2})\text{
s.t.}\begin{split}&(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}\geq
0\\\
&(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}^{\prime}\geq
0\end{split},\forall i,$ (8)

where $\mathbf{U}\in\mathbb{R}^{n\times r}$ and
$\mathbf{U}^{\prime}\in\mathbb{R}^{n\times r+1}$ are computed via the compact
SVD of the zero-mean data matrix, namely
$(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}:=\mathbf{U}\bm{\Sigma}\mathbf{V}^{T}$
and
$\mathbf{U}^{\prime}:=\begin{bmatrix}\mathbf{U}&\frac{1}{\sqrt{n}}\mathbf{1}\end{bmatrix}$.

Notice that (8) is a convex program with $2Pr$ variables and $2Pn$
constraints, and thus can be globally optimized via interior-point solvers
with $\mathcal{O}(r^{6}\left(\frac{n}{r}\right)^{3r})$ complexity.

######  Remark 2.1.

Theorem 2.2 shows that a classical ReLU network with BN can be equivalently
described as a convex combination of linear models
$\\{\mathbf{U}^{\prime}\mathbf{s}_{i}\\}_{i=1}^{P}$ and
$\\{\mathbf{U}^{\prime}\mathbf{s}_{i}^{\prime}\\}_{i=1}^{P}$ multiplied with
fixed hyperplane arrangement matrices $\\{\mathbf{D}_{i}\\}_{i=1}^{P}$.
Therefore, this result shows that optimal ReLU networks with BN are sparse
piecewise linear functions, where sparsity is enforced via the group lasso
regularization in (8). The convex program also reveals an implicit mechanism
behind BN. Particularly, comparing (8) with the one without BN (see eqn. (8)
in Pilanci & Ergen (2020)), we observe that the data matrix for the convex
program is whitened, i.e.,
$\mathbf{U}^{\prime^{T}}\mathbf{U}^{\prime}=\mathbf{I}_{r+1}$.

####  2.3 Vector Output Networks

Here, we analyze networks with vector outputs, i.e.,
$\mathbf{Y}\in\mathbb{R}^{n\times C}$, trained via the following problem

$\displaystyle
p_{v2}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|f_{\theta,2}(\mathbf{X})-\mathbf{Y}\right\|_{F}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left(\left\|\mathbf{w}^{(1)}_{j}\right\|_{2}^{2}+{\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}+\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}^{2}\right).$
(9)

######  Lemma 2.1.

The problem in (9) is equivalent to the following optimization problem

$\displaystyle
p_{2v}^{*}=\min_{\theta\in\Theta_{s}}\frac{1}{2}\left\|f_{\theta,2}(\mathbf{X})-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{j=1}^{m_{1}}\left\|\mathbf{w}^{(2)}_{j}\right\|_{2},$
(10)

where
$\Theta_{s}:=\\{\theta\in\Theta:{\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}=1,\,\forall
j\in[m_{1}]\\}$.

Using the equivalence in Lemma 2.1, we then take dual of (10) with respect to
the output layer weights $\mathbf{W}^{(2)}$ to obtain the following dual
problem

$\displaystyle p_{2v}^{*}\geq d_{2v}^{*}:=$
$\displaystyle\max_{\mathbf{V}}-\frac{1}{2}\left\|\mathbf{V}-\mathbf{Y}\right\|_{F}^{2}+\frac{1}{2}\left\|\mathbf{Y}\right\|_{F}^{2}\text{
s.t.
}\max_{\theta\in\Theta_{s}}\left\|\mathbf{V}^{T}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{X}\mathbf{w}^{(1)}\right)\right)_{+}\right\|_{2}\leq\beta.$
(11)

######  Theorem 2.3.

The non-convex training problem (9) can be cast as the following convex
program

$\displaystyle
p_{v2}^{*}=\min_{\mathbf{S}_{i}}\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{S}_{i}-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{i=1}^{P}\|\mathbf{S}_{i}\|_{c_{i},*}$
(12)

for the norm $\|\cdot\|_{c_{i},*}$ defined as

$\|\mathbf{S}\|_{c_{i},*}:=\min_{t\geq 0}t\text{ s.t. }\mathbf{S}\in
t\mathrm{conv}\\{\mathbf{Z}=\mathbf{h}\mathbf{g}^{T}:(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{h}\geq
0,\,\|\mathbf{Z}\|_{*}\leq 1\\},$

where $\mathrm{conv}\\{\cdot\\}$ denotes the convex hull of its argument and
$\mathbf{U}^{\prime}=\begin{bmatrix}\mathbf{U}&\frac{1}{\sqrt{n}}\mathbf{1}\end{bmatrix}$
as in Theorem 2.2.

######  Remark 2.2.

We note that similar to the scalar output case, vector output BN networks
training problem involves a whitened data matrix. Furthermore, unlike (8), the
nuclear norm regularization in (12) encourages a low rank piecewise linear
function at the network output.

The above problem can be solved in $\mathcal{O}((n/d)^{d})$ time in the worst
case (Sahiner et al., 2021a), which is polynomial for fixed $d$. We note that
the problem in (12) is similar to convex semi non-negative matrix
factorizations (Ding et al., 2008). Moreover, in certain practically relevant
cases, (12) can be significantly simplified so that one can even obtain
closed-form solutions as demonstrated below.

######  Theorem 2.4.

Let $\mathbf{Y}$ be a one-hot encoded label matrix, and $\mathbf{X}$ be a full
row-rank data matrix with $n\leq d$, then an optimal solution to (10) admits

$\displaystyle\begin{split}&\left({\mathbf{w}^{(1)}_{j}}^{*},{\mathbf{w}^{(2)}_{j}}^{*}\right)=\left(\mathbf{X}^{\dagger}\mathbf{y}_{j},\left(\|\mathbf{y}_{j}\|_{2}-\beta\right)_{+}\mathbf{e}_{j}\right),\;\begin{bmatrix}{\gamma^{(1)}_{j}}^{*}\\\
{\alpha^{(1)}_{j}}^{*}\end{bmatrix}=\frac{1}{\|\mathbf{y}_{j}\|_{2}}\begin{bmatrix}\|\mathbf{y}_{j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{y}_{j}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\mathbf{y}_{j}\end{bmatrix}\end{split}$

$\forall j\in[C]$, where $\mathbf{e}_{j}$ is the $j^{th}$ ordinary basis
vector.

###  3 Convolutional Neural Networks

In this section, we extend our analysis to CNNs. Thus, instead of
$\mathbf{X}$, we operate on the patch matrices that are subsets of columns
extracted from $\mathbf{X}$ and we denote them as
$\mathbf{X}_{k}\in\mathbb{R}^{n\times h}$, where $h$ is the filter size and
$k$ is the patch index. With this notation, a convolution operation between
the data matrix $\mathbf{X}$ and an arbitrary filter
$\mathbf{z}\in\mathbb{R}^{h}$ can be represented as
$\\{\mathbf{X}_{k}\mathbf{z}\\}_{k=1}^{K}$. Therefore, given patch matrices
$\\{\mathbf{X}_{k}\\}_{k=1}^{K}$, regularized training problem with BN can be
formulated as follows

$\displaystyle
p_{2c}^{*}:=\min_{\theta\in\Theta_{c}}\frac{1}{2}\left\|\sum_{j=1}^{m_{1}}\sum_{k=1}^{K}\left({\text{BN}}^{C}_{\gamma,\alpha}\left(\mathbf{X}_{k}\mathbf{z}_{j}\right)\right)_{+}w_{j}-\mathbf{y}\right\|_{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}(\|\mathbf{z}_{j}\|_{2}^{2}+\gamma_{j}^{2}+\alpha_{j}^{2}+w_{j}^{2}),$
(13)

where given
$\mu_{j}=\frac{1}{nK}\sum_{k=1}^{K}\mathbf{1}^{T}\mathbf{X}_{k}\mathbf{z}_{j}$,
BN is defined as

$\displaystyle{\text{BN}}^{C}_{\gamma,\alpha}\left(\mathbf{X}_{k}\mathbf{z}_{j}\right)$
$\displaystyle:=\frac{\mathbf{X}_{k}\mathbf{z}_{j}-\mu_{j}\mathbf{1}}{\sqrt{\sum_{k^{\prime}=1}^{K}\|\mathbf{X}_{k^{\prime}}\mathbf{z}_{j}-\mu_{j}\mathbf{1}\|_{2}^{2}}}\gamma^{(1)}_{j}+\frac{\alpha^{(1)}_{j}\mathbf{1}}{\sqrt{nK}}.$

Using Lemma 1.1, the dual problem with respect to $\mathbf{w}^{(2)}$ is

$\displaystyle p_{2c}^{*}\geq d_{2c}^{*}:=$
$\displaystyle\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{T}\sum_{k=1}^{K}\left({\text{BN}}^{C}_{\gamma,\alpha}\left(\mathbf{X}_{k}\mathbf{z}\right)\right)_{+}\right|\leq\beta.$
(14)

Next, we state the equivalent convex program for CNNs.

######  Theorem 3.1.

The non-convex training problem in (13) can be cast as a convex program

$\displaystyle\min_{\mathbf{q}_{i},\mathbf{q}_{i}^{\prime}\in\mathbb{R}^{h+1}}\frac{1}{2}\left\|\sum_{i=1}^{P_{c}}\sum_{k=1}^{K}\mathbf{D}_{ik}\mathbf{U}_{Mk}^{\prime}(\mathbf{q}_{i}-\mathbf{q}_{i}^{\prime})\right\|+\beta\sum_{i=1}^{P_{c}}(\left\|\mathbf{q}_{i}\right\|_{2}+\left\|\mathbf{q}_{i}^{\prime}\right\|_{2})\text{
s.t.}\begin{split}&(2\mathbf{D}_{ik}-\mathbf{I}_{n})\mathbf{U}_{Mk}^{\prime}\mathbf{q}_{i}\geq
0\\\
&(2\mathbf{D}_{ik}-\mathbf{I}_{n})\mathbf{U}_{Mk}^{\prime}\mathbf{q}_{i}^{\prime}\geq
0\end{split},\forall i,k$ (15)

where we first define a new matrix
$\mathbf{M}=[\mathbf{X}_{1};\,\mathbf{X}_{2};\,\ldots\,\mathbf{X}_{K}]\in\mathbb{R}^{nK\times
h}$ and then use the following notations: the compact SVD of the zero mean
form of this matrix is
$(\mathbf{I}_{nK}-\frac{1}{nK}\mathbf{1}\mathbf{1}^{T})\mathbf{M}:=\mathbf{U}_{M}\bm{\Sigma}_{M}\mathbf{V}_{M}^{T}$,
and
$\mathbf{U}_{M}^{\prime}:=\begin{bmatrix}\mathbf{U}_{M}&\frac{\mathbf{1}}{\sqrt{nK}}\end{bmatrix}=\begin{bmatrix}\mathbf{U}_{M1}^{\prime};\mathbf{U}_{M2}^{\prime};\ldots&\mathbf{U}_{MK}^{\prime}\end{bmatrix}$.

Theorem 3.1 proves that (13) can be equivalently stated as a finite-
dimensional convex program with $2r_{c}P_{c}$ variables (see Remark 3.1 for
the definitions) and $2nP_{c}K$ constraints. Since $P_{c}$ is polynomial in
$n$ and $d$ as detailed in the remark below, we can solve (15) in polynomial-
time.

######  Remark 3.1.

The number of hyperplane arrangements for CNNs rely on the rank of the patch
matrices instead of $\mathbf{X}$. Particularly, the number of hyperplane
arrangements for $\mathbf{M}$, i.e., denoted as $P_{c}$, is upperbounded as
$P_{c}\leq 2r_{c}(e(n-1)/r_{c})^{r_{c}}$, where
$r_{c}:=\mathrm{rank}(\mathbf{M})\leq h\ll d$. Therefore, given a fixed filter
size $h$, $P_{c}$ is polynomial in $n$ and $d$ even when $\mathbf{X}$ is full
rank (see Ergen & Pilanci (2021a) for more details).

###  4 Deep Networks

We now analyze $L$-layer neural networks trained via the non-convex
optimization problem in (3). Below, we provide an explicit formulation for the
last two layers’ weights.

######  Theorem 4.1.

Suppose the network is overparameterized such that the range of
$\mathbf{A}^{(L-2)}$ is $\mathbb{R}^{n}$. Then an optimal solution for the
last two layer weights in (4) admits

$\displaystyle\left({\mathbf{w}^{(L-1)}_{j}}^{*},{{w}^{(L)}_{j}}^{*}\right)=\left({\mathbf{A}^{(L-2)}}^{\dagger}\left((-1)^{j}\mathbf{y}\right)_{+},(-1)^{j}\left(\|\left((-1)^{j}\mathbf{y}\right)_{+}\|_{2}-\beta\right)_{+}\right)$
$\displaystyle\begin{bmatrix}{\gamma^{(L-1)}_{j}}^{*}\\\
{\alpha^{(L-1)}_{j}}^{*}\end{bmatrix}=\frac{1}{\|\left((-1)^{j}\mathbf{y}\right)_{+}\|_{2}}\begin{bmatrix}\|\left((-1)^{j}\mathbf{y}\right)_{+}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\left((-1)^{j}\mathbf{y}\right)_{+}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\left((-1)^{j}\mathbf{y}\right)_{+}\end{bmatrix}$

for $j=1,2$. Thus, strong duality holds, i.e., $p_{L}^{*}=d_{L}^{*}$.

####  4.1 Vector Output Networks

Vector output deep networks are commonly used for multi-class classification
problems, where one-hot encoding is a prevalent strategy to convert
categorical variables into a binary representation that can be processed by
neural networks. Although empirical evidence, e.g., Papyan et al. (2020),
shows that deep neural networks trained with one-hot encoded labels exhibit
some common patterns among classes, it is still theoretically elusive how
these patterns emerge. Therefore, in this section, we analyze vector output
deep networks trained using one-hot labels via our advocated convex duality.

The following theorem presents a complete characterization for the last two
layers’ weights.

######  Theorem 4.2.

Let $\mathbf{Y}$ be a one-hot encoded label matrix and the network be
overparameterized such that the range of $\mathbf{A}^{(L-2)}$ is
$\mathbb{R}^{n}$, then an optimal solution to (3) can be found in closed-form
as

$\displaystyle\begin{split}&\left({\mathbf{w}^{(L-1)}_{j}}^{*},{\mathbf{w}^{(L)}_{j}}^{*}\right)=\left({\mathbf{A}^{(L-2)}}^{\dagger}\mathbf{y}_{j},\left(\|\mathbf{y}_{j}\|_{2}-\beta\right)_{+}\mathbf{e}_{j}\right),\;\begin{bmatrix}{\gamma^{(L-1)}_{j}}^{*}\\\
{\alpha^{(L-1)}_{j}}^{*}\end{bmatrix}=\frac{1}{\|\mathbf{y}_{j}\|_{2}}\begin{bmatrix}\|\mathbf{y}_{j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{y}_{j}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\mathbf{y}_{j}\end{bmatrix}\end{split}$

$\forall j\in[C]$, where $\mathbf{e}_{j}$ is the $j^{th}$ ordinary basis
vector.

###  5 Implicit Regularization of GD for BN

(a) BN model (9)

(b) Whitened model (16)

Figure 1: Cosine similarity of network weights relative to initialization for
each singular (in descending order) value direction of the data. Here, we use
the same setting in Figure 2. We compare training the BN model (9) to the
equivalent whitened model (16), noting that the BN model fits high singular
value directions of the data first, while the whitened model fits all
directions equally, thereby overfitting to non-robust features. This shows
implicit regularization effect of SGD.

We first note that proposed convex formulations, e.g., (8), (12) and (15),
utilize whitened data matrices. Thus, when we train them via a standard first
order optimizers such as GD/SGD, the optimizer might follow an unfavorable
optimization path that leads to slower convergence and worse generalization
properties throughout the training. Based on this, we first analyze the impact
of utilizing a whitened data matrix for the non-convex problem in (9) trained
with GD. Then, we characterize this impact as an implicit regularization due
to the choice of optimizer. Finally, we propose an approach to explicitly
incorporate this benign implicit regularization into both convex (8) and non-
convex (9) training objectives.

####  5.1 Evidence of Implicit Regularization

######  Lemma 5.1.

The non-convex problem (9) is equivalent to the following problem (i.e.,
$p_{v2}^{*}=p_{v2b}^{*}$)

$\displaystyle
p_{v2b}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|\sum_{j=1}^{m_{1}}\left(\frac{\mathbf{U}\mathbf{q}_{j}}{\|\mathbf{q}_{j}\|_{2}}{\gamma^{(1)}_{j}}+\frac{\mathbf{1}}{\sqrt{n}}{\alpha^{(1)}_{j}}\right)_{+}{\mathbf{w}^{(2)}_{j}}^{T}-\mathbf{Y}\right\|_{F}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}+\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}^{2}\right).$
(16)

(16) is equivalent to a problem with new data matrix $\mathbf{U}$ and weight
normalization (Salimans & Kingma, 2016). The identity in Lemma 5.1 is achieved
by substituting the transformation
${\mathbf{q}_{j}}=\mathbf{\Sigma}\mathbf{V}^{\top}{\mathbf{w}^{(1)}_{j}}$.

While (9) and (16) have identical global optima, in practice, GD behaves quite
differently when applied to the two problems. Furthermore, these differences
extend to the convex formulations in (8), (12), and (15). In particular, while
these programs are adept at finding the global optimum, they generalize quite
poorly when trained with local search algorithms such as GD. This suggests an
implicit regularization effect of GD, the form of which is demonstrated in the
following theorem.

######  Theorem 5.1.

Consider (9) and (16) with corresponding models $f_{v2}(\mathbf{X})$ and
$f_{v2b}(\mathbf{U})$. At iteration $k$ of GD, assume that
$f_{v2}^{k}(\mathbf{X})=f_{v2b}^{k}(\mathbf{U})$, i.e. both models have
identical weights, besides
$\mathbf{q}_{j}^{(k)}=\mathbf{\Sigma}\mathbf{V}^{\top}{\mathbf{w}^{(1)}_{j}}^{(k)}$.
Then, the subsequent GD updates to model weights are identical, besides those
to $\mathbf{w}^{(1)}_{j}$ and $\mathbf{q}_{j}$. In $\mathbf{q}_{j}$-space, the
updates $\Delta_{v2,j}^{(k)}$ and $\Delta_{v2b,j}^{(k)}$, respectively, admit

$\displaystyle\Delta_{v2,j}^{(k)}=\mathbf{\Sigma}^{2}\Delta_{v2b,j}^{(k)}.$

######  Remark 5.1.

Theorem 5.1 indicates that when solving (16), GD will fit directions
corresponding to higher singular values of data much slower relative to other
directions compared to (9). By fitting these directions more aggressively, GD
applied to (9) provides an implicit regularization effect that fits to more
robust features. We also note that although our analysis holds only for GD, we
empirically observe and validate the emergence of this phenomenon for networks
trained with SGD as well.

To verify this intuition, in Figure 2, we compare the training and test loss
curves of solving both (9) and (16) with mini-batch SGD applied to CIFAR-10
(Krizhevsky et al., 2014), demonstrating that SGD applied to (9) generalizes
much better than (16). In Figure 1, we plot the cosine similarity of the
weights of both networks compared to their initialization in each singular
value direction.

Our results confirm that the BN model fits high singular value directions much
faster than low ones, and in fact the lower-singular value directions almost
do not move at all from their initialization. In contrast, the whitened model
fits all singular value directions roughly equally, meaning that it overfits
to less robust components of the data. These experiments confirm our
theoretical intuition that SGD applied to BN networks has an implicit
regularization effect.

Figure 2: CIFAR-10 classification accuracies for a two-layer ReLU network in
(9) and the equivalent whitened model (16) with SGD and
$(n,d,m_{1},\beta,\text{bs})=(50000,3072,1000,10^{-4},1000)$, where bs denotes
batch size. The whitened model overfits much more than the standard BN model,
indicating implicit regularization.

####  5.2 Making Implicit Regularization Explicit

Since the convex programs (8), (12), and (15) rely on the whitened data matrix
$\mathbf{U}$, they lack the regularization necessary for generalization when
optimized with GD. Thus, we aim to impose an explicit regularizer which does
not hamper the performance of standard BN networks, but can improve the
generalization of the equivalent programs trained using first-order optimizers
such as GD. One simple idea to prevent the whitened model from fitting low
singular value directions is to obtain a low-rank approximation of the data by
filtering its singular values:

$g(\mathbf{X}):=\mathbf{U}\mathcal{T}(\mathbf{\Sigma})\mathbf{V}^{\top},$ (17)

with $\mathcal{T}_{k}(\mathbf{\Sigma})_{ii}:=\sigma_{i}\mathbbm{1}\\{i\geq
k\\}$, which takes the $k\leq r$ top singular values and removes the rest. We
can then solve (1) using $f_{\theta,L}(g(\mathbf{X}))$ as the network output,
to form what we denote as the truncated problem. This can also reduce the
problem dimension in the whitened-data case since one can simply omit the
columns of $\mathbf{U}$. We find this low-rank approximation to be effective,
which can even improve the generalization of GD applied to BN networks.

###  6 Numerical Experiments

Here555We present additional numerical experiments and details on the
experiments in Appendix A., we present experiments to verify our theory.
Throughout this section, we denote the baseline model (3) and our non-convex
model with truncation as “SGD” and “SGD-Truncated”, respectively. Likewise,
the proposed convex models are denoted as “Convex” and “Convex-Truncated”.

(a) CIFAR-10-$(m_{1},\beta,\text{bs})=(4096,10^{-4},10^{3})$

(b) CIFAR-100-$(m_{1},\beta,\text{bs})=(1024,10^{-4},10^{3})$

Figure 3: Test accuracy of two-layer BN ReLU networks including the standard
non-convex (baseline), its convex equivalent (8), along with their truncated
versions with $k=(215,200)$ respectively. All problems are solved with SGD.
The learning rates are chosen based on test performance and we provide
learning rate and optimizer comparisons in Appendix A. Figure 4: Comparing the
objective values of the closed-form solution in Theorem 2.4, (“Theory”), and
GD (5 trials) applied to the non-convex problem (9) for three-class CIFAR
classification with $(n,d,m_{1},\beta)=(1500,3072,1000,1)$. Here, “x" denotes
the computation time for the closed-form solution.

Closed-Form Solutions: We first verify Theorem 2.4. We consider a three-class
classification task with the CIFAR-100 image classification dataset
(Krizhevsky et al., 2014) trained with squared loss and one-hot encoded
labels, while we set $(n,d,m_{1},\beta)=(1500,3072,1000,1)$, where bs denotes
the batch size. In Figure 4, we compare the objective values achieved by the
closed-form solution and the classical network (6). We observe that our
formula is computationally light-weight and can quickly find a global optimum
that cannot be achieved with GD even after convergence.

Convex Formulation and Implicit Regularization: We now demonstrate the
effectiveness of our convex formulation to solve (9) on CIFAR-10 and CIFAR-100
(Krizhevsky et al., 2014). Here, instead of enumerating all hyperplane
arrangements, we randomly sample via randomly generated Gaussian vectors such
that $m_{1}=P$. We compare the results of training with and without the
explicit regularizer (17) for both cases. We particularly choose $k$ such that
singular values explain 95% of the variance in the training data. Then, we
consider two-layer networks for image classification on CIFAR-10 with
$(n,d,m_{1},\beta,\text{bs},k,C)=(5\text{x}10^{4},3072,4096,10^{-4},10^{3},215,10)$
and CIFAR-100 with
$(n,d,m_{1},\beta,\text{bs},k,C)=(5\text{x}10^{4},3072,1024,10^{-4},10^{3},200,100)$.
In Figure 3, we observe that Convex generalizes poorly (as expected with our
results in Section 5.1). After the singular value truncation, SGD, SGD-
Truncated, and Convex-Truncated generalize equally well for CIFAR-10, and
Convex-Truncated in fact generalizes better for CIFAR-100. This shows that the
optimizer benefits from the convex landscape of the proposed formulation (8).
Furthermore, truncation actually improves the performance of SGD in the case
of CIFAR-100, suggesting this is a good regularizer to use.

###  7 Concluding Remarks

In this paper, we studied training of deep weight-decay regularized ReLU
networks with BN, and introduced a convex analytic framework to characterize
an optimal solution to the training problem. Using this framework, we proved
that the non-convex training problem can be equivalently cast as a convex
optimization problem that can be trained in polynomial-time. Moreover, the
convex equivalents involve whitened data matrix in contrast to the original
non-convex problem using the raw data, which reveals a whitening effect
induced by BN. We also extended this approach to analyze different
architectures, e.g., CNNs and networks with BN placed after ReLU. In the high
dimensional and/or overparameterized regime, our characterization further
simplifies the equivalent convex representations such that closed-form
solutions for optimal parameters can be obtained. Therefore, in these regimes,
there is no need to train a network in an end-to-end manner. In the light of
our results, we also unveiled an implicit regularization effect of GD on non-
convex ReLU networks with BN, which prioritizes learning high singular-value
directions of the training data. We also proposed a technique to incorporate
this regularization to the training problem explicitly, which is numerically
verified to be effective for image classification problems. Finally, we remark
that after this work, similar convex duality arguments were applied to deeper
networks Ergen & Pilanci (2021c; d; e) and Wasserstein Generative Adversarial
Networks (WGANs) Sahiner et al. (2022). Extensions of our convex BN analysis
to these architectures are left for future work.

##### Acknowledgments

This work was partially supported by the National Science Foundation under
grants ECCS- 2037304, DMS-2134248, the Army Research Office.

### References

  * Agrawal et al. (2018) Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd.  A rewriting system for convex optimization problems.  _Journal of Control and Decision_ , 5(1):42–60, 2018.
  * Arora et al. (2018) Sanjeev Arora, Zhiyuan Li, and Kaifeng Lyu.  Theoretical analysis of auto rate-tuning by batch normalization.  2018\.
  * Ba et al. (2016) Jimmy Lei Ba, Jamie Ryan Kiros, and Geoffrey E Hinton.  Layer normalization.  _arXiv preprint arXiv:1607.06450_ , 2016.
  * Bjorck et al. (2018) Nils Bjorck, Carla P Gomes, Bart Selman, and Kilian Q Weinberger.  Understanding batch normalization.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 31, pp. 7694–7705. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper/2018/file/36072923bfc3cf47745d704feb489480-Paper.pdf.
  * Boyd & Vandenberghe (2004) Stephen Boyd and Lieven Vandenberghe.  _Convex optimization_.  Cambridge university press, 2004.
  * Chen et al. (2019) Guangyong Chen, Pengfei Chen, Yujun Shi, Chang-Yu Hsieh, Benben Liao, and Shengyu Zhang.  Rethinking the usage of batch normalization and dropout in the training of deep neural networks.  _arXiv preprint arXiv:1905.05928_ , 2019.
  * Ciarelli & Oliveira (2009) Patrick Marques Ciarelli and Elias Oliveira.  Agglomeration and elimination of terms for dimensionality reduction.  In _2009 Ninth International Conference on Intelligent Systems Design and Applications_ , pp. 547–552. IEEE, 2009.
  * Cover (1965) Thomas M Cover.  Geometrical and statistical properties of systems of linear inequalities with applications in pattern recognition.  _IEEE transactions on electronic computers_ , (3):326–334, 1965.
  * Daneshmand et al. (2020) Hadi Daneshmand, Jonas Kohler, Francis Bach, Thomas Hofmann, and Aurelien Lucchi.  Batch normalization provably avoids rank collapse for randomly initialised deep networks, 2020.
  * Diamond & Boyd (2016) Steven Diamond and Stephen Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  _Journal of Machine Learning Research_ , 17(83):1–5, 2016.
  * Ding et al. (2008) Chris HQ Ding, Tao Li, and Michael I Jordan.  Convex and semi-nonnegative matrix factorizations.  _IEEE transactions on pattern analysis and machine intelligence_ , 32(1):45–55, 2008.
  * Ergen & Pilanci (2020a) Tolga Ergen and Mert Pilanci.  Convex geometry of two-layer relu networks: Implicit autoencoding and interpretable models.  In Silvia Chiappa and Roberto Calandra (eds.), _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , volume 108 of _Proceedings of Machine Learning Research_ , pp. 4024–4033, Online, 26–28 Aug 2020a. PMLR.  URL http://proceedings.mlr.press/v108/ergen20a.html.
  * Ergen & Pilanci (2020b) Tolga Ergen and Mert Pilanci.  Revealing the structure of deep neural networks via convex duality.  2020b.
  * Ergen & Pilanci (2021a) Tolga Ergen and Mert Pilanci.  Implicit convex regularizers of cnn architectures: Convex optimization of two- and three-layer networks in polynomial time.  In _International Conference on Learning Representations_ , 2021a.  URL https://openreview.net/forum?id=0N8jUH4JMv6.
  * Ergen & Pilanci (2021b) Tolga Ergen and Mert Pilanci.  Convex geometry and duality of over-parameterized neural networks.  _Journal of Machine Learning Research_ , 22(212):1–63, 2021b.  URL http://jmlr.org/papers/v22/20-1447.html.
  * Ergen & Pilanci (2021c) Tolga Ergen and Mert Pilanci.  Global optimality beyond two layers: Training deep relu networks via convex programs.  In Marina Meila and Tong Zhang (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 2993–3003. PMLR, 18–24 Jul 2021c.  URL http://proceedings.mlr.press/v139/ergen21a.html.
  * Ergen & Pilanci (2021d) Tolga Ergen and Mert Pilanci.  Path regularization: A convexity and sparsity inducing regularization for parallel relu networks.  _arXiv preprint arXiv:2110.09548_ , 2021d.
  * Ergen & Pilanci (2021e) Tolga Ergen and Mert Pilanci.  Revealing the structure of deep neural networks via convex duality.  In Marina Meila and Tong Zhang (eds.), _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 3004–3014. PMLR, 18–24 Jul 2021e.  URL http://proceedings.mlr.press/v139/ergen21b.html.
  * Grant & Boyd (2014) Michael Grant and Stephen Boyd.  CVX: Matlab software for disciplined convex programming, version 2.1.  http://cvxr.com/cvx, March 2014.
  * He et al. (2015) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 1026–1034, 2015.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Im et al. (2016) Daniel Jiwoong Im, Michael Tao, and Kristin Branson.  An empirical analysis of the optimization of deep network loss surfaces.  _arXiv preprint arXiv:1612.04010_ , 2016.
  * Ioffe & Szegedy (2015) Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In Francis Bach and David Blei (eds.), _Proceedings of the 32nd International Conference on Machine Learning_ , volume 37 of _Proceedings of Machine Learning Research_ , pp. 448–456, Lille, France, 07–09 Jul 2015\. PMLR.  URL http://proceedings.mlr.press/v37/ioffe15.html.
  * Kingma & Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Krizhevsky et al. (2014) Alex Krizhevsky, Vinod Nair, and Geoffrey Hinton.  The CIFAR-10 and -100 datasets.  http://www.cs.toronto.edu/kriz/cifar.html, 2014.
  * Lian & Liu (2019) Xiangru Lian and Ji Liu.  Revisit batch normalization: New understanding and refinement via composition optimization.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pp. 3254–3263. PMLR, 2019.
  * Neyshabur et al. (2014) Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  In search of the real inductive bias: On the role of implicit regularization in deep learning.  _arXiv preprint arXiv:1412.6614_ , 2014.
  * Ojha (2000) Piyush C Ojha.  Enumeration of linear threshold functions from the lattice of hyperplane intersections.  _IEEE Transactions on Neural Networks_ , 11(4):839–850, 2000.
  * Papyan et al. (2020) Vardan Papyan, XY Han, and David L Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  _Proceedings of the National Academy of Sciences_ , 117(40):24652–24663, 2020.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Pilanci & Ergen (2020) Mert Pilanci and Tolga Ergen.  Neural networks are convex regularizers: Exact polynomial-time convex optimization formulations for two-layer networks.  In Hal Daumé III and Aarti Singh (eds.), _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 7695–7705. PMLR, 13–18 Jul 2020.  URL http://proceedings.mlr.press/v119/pilanci20a.html.
  * Sahiner et al. (2021a) Arda Sahiner, Tolga Ergen, John M. Pauly, and Mert Pilanci.  Vector-output relu neural network problems are copositive programs: Convex analysis of two layer networks and polynomial-time algorithms.  In _International Conference on Learning Representations_ , 2021a.  URL https://openreview.net/forum?id=fGF8qAqpXXG.
  * Sahiner et al. (2021b) Arda Sahiner, Morteza Mardani, Batu Ozturkler, Mert Pilanci, and John M. Pauly.  Convex regularization behind neural reconstruction.  In _International Conference on Learning Representations_ , 2021b.  URL https://openreview.net/forum?id=VErQxgyrbfn.
  * Sahiner et al. (2022) Arda Sahiner, Tolga Ergen, Batu Ozturkler, Burak Bartan, John M. Pauly, Morteza Mardani, and Mert Pilanci.  Hidden convexity of wasserstein GANs: Interpretable generative models with closed-form solutions.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=e2Lle5cij9D.
  * Salimans & Kingma (2016) Tim Salimans and Diederik P Kingma.  Weight normalization: A simple reparameterization to accelerate training of deep neural networks.  In _NIPS_ , 2016.
  * Santurkar et al. (2018) Shibani Santurkar, Dimitris Tsipras, Andrew Ilyas, and Aleksander Madry.  How does batch normalization help optimization?  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 31, pp. 2483–2493. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper/2018/file/905056c1ac1dad141560467e0a99e1cf-Paper.pdf.
  * Savarese et al. (2019) Pedro Savarese, Itay Evron, Daniel Soudry, and Nathan Srebro.  How do infinite width bounded norm networks look in function space?  _CoRR_ , abs/1902.05040, 2019.  URL http://arxiv.org/abs/1902.05040.
  * Sion (1958) Maurice Sion.  On general minimax theorems.  _Pacific J. Math._ , 8(1):171–176, 1958.  URL https://projecteuclid.org:443/euclid.pjm/1103040253.
  * Stanley et al. (2004) Richard P Stanley et al.  An introduction to hyperplane arrangements.  _Geometric combinatorics_ , 13:389–496, 2004.
  * Summers & Dinneen (2020) Cecilia Summers and Michael J. Dinneen.  Four things everyone should know to improve batch normalization.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=HJx8HANFDH.
  * Tütüncü et al. (2001) RH Tütüncü, KC Toh, and MJ Todd.  Sdpt3—a matlab software package for semidefinite-quadratic-linear programming, version 3.0.  _Web page http://www. math. nus. edu. sg/mattohkc/sdpt3. html_ , 2001\.
  * Ulyanov et al. (2016) Dmitry Ulyanov, Andrea Vedaldi, and Victor Lempitsky.  Instance normalization: The missing ingredient for fast stylization.  _arXiv preprint arXiv:1607.08022_ , 2016.
  * Wei et al. (2019) Mingwei Wei, James Stokes, and David J Schwab.  Mean-field analysis of batch normalization, 2019.
  * Winder (1966) RO Winder.  Partitions of n-space by hyperplanes.  _SIAM Journal on Applied Mathematics_ , 14(4):811–818, 1966.
  * Wu & He (2018) Yuxin Wu and Kaiming He.  Group normalization.  In _Proceedings of the European conference on computer vision (ECCV)_ , pp. 3–19, 2018.
  * Zhang et al. (2021) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning (still) requires rethinking generalization.  _Communications of the ACM_ , 64(3):107–115, 2021\.

## Appendix

###  Appendix A Additional Details and Numerical Experiments

In this section, we present the details about our experimental setup in the
main paper and provide new numerical experiments.

We first note that since we derive exact convex formulations in Theorem 2.2,
2.3, 3.1, we can use solver such as CVX (Grant & Boyd, 2014) and CVXPY
(Diamond & Boyd, 2016; Agrawal et al., 2018) with the SDPT3 solver (Tütüncü et
al., 2001) to train regularized ReLU networks training problems with BN.
Although these solvers ensure quite fast convergence rates to a global
minimum, they do not scale to moderately large datasets, e.g., CIFAR-10.
Therefore, in our main experiments, we use an equivalent unconstrained form
for the constrained convex optimization problems such that it can be simply
trained via SGD. Particularly, let us consider the constrained convex
optimization problem in (8), for which the constraints can be incorporated
into the objective as follows

$\displaystyle\min_{\mathbf{s}_{i},\mathbf{s}_{i}^{\prime}}\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}(\mathbf{s}_{i}-\mathbf{s}_{i}^{\prime})-\mathbf{y}\right\|_{2}^{2}+\beta\sum_{i=1}^{P}(\|\mathbf{s}_{i}\|_{2}+\|\mathbf{s}_{i}^{\prime}\|_{2})$
$\displaystyle\hskip
99.58464pt+\rho\sum_{i=1}^{P}\mathbf{1}^{T}\left(\left(-(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}\right)_{+}+\left(-(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}^{\prime}\right)_{+}\right)$

where $\rho>0$ is a hyper-parameter to penalize the violating constraints.
Therefore, using the unconstrained form above, we are able to use SGD or other
standard first-order optimizers to optimize the constrained convex
optimization problem in (8).

Unless otherwise stated, all problems are run on a single Titan Xp GPU with 12
GB of RAM, with the Pytorch deep learning library (Paszke et al., 2019). All
SGD models use a momentum parameter of 0.9. All baseline models are trained
with standard weight-decay regularization as outlined in (3). All batch-norm
models are initialized with scale $\gamma^{(l)}_{j}=1$ and bias
$\alpha^{(l)}_{j}=0$, and all other non-convex model weights are initialized
with standard Kaiming uniform initialization (He et al., 2015). For the convex
program, we simply initialize all weights to zero.

####  A.1 Inference in BN models

For the non-convex standard BN network, we perform inference in the
traditional way. All BN problems with SGD take an exponential moving average
of the mini-batch mean and standard deviation with a momentum of 0.1, as is
used as a default in deep learning packages, except in the case of full-batch
GD, in which case a momentum value of 1 is used. In truncated variants, while
the training data is truncated with the explicit regularizer $g(\mathbf{X})$,
for test data $\hat{\mathbf{X}}$, we compute the standard forward pass
$f_{\theta,L}(\hat{\mathbf{X}}$), i.e. without changing the test data in any
fashion, and removing the training-data mean and normalizing by the training-
data standard deviation as found with the exponential moving average method
described above. For the convex program, we can form the weights to the non-
convex standard BN network from the convex program weights, and compute the
test-set predictions with the non-convex standard BN network.

####  A.2 Additional Experiment: Batch Sizes for SGD

The model of BN we present in this paper is with the model with GD without
mini-batching. In this section, we demonstrate that mini-batching does not
change the training or generalization performance of BN networks in any
significant manner. Our observations also align with previous work examining
this subject, such as Lian & Liu (2019); Summers & Dinneen (2020).

In particular, we compare the effect of different batch sizes using a four
layer CNN, where the first three layers are convolutional, followed by BN,
ReLU, and average pooling, and the final layer is a fully-connected layer.
Each CNN layer has 1000 filters with kernel size $3\times 3$ and padding of 1.

We take the first two classes from CIFAR-100 to perform a binary
classification task with $\mathbf{X}\in\mathbb{R}^{1000\times 3072}$ and
$\mathbf{y}\in\mathbb{R}^{1000}$. For $(\text{lr},\beta)=(10^{-6},10^{-4})$,
where lr denotes the learning rate, and weight-decay is applied to all
parameters. We compare the effect of using a batch size of $50,100,500$, and
1000. We run the batch size 1000 case for 501 epochs, and then train the other
batch size cases such that the training time is roughly matched. Due to memory
constraints, the full-batch form of SGD does not fit on a standard GPU–thus,
all models are trained with a standard CPU with 256 GB of RAM. We compare the
train and test accuracy curves in Figure 5, where we see that all batch sizes
perform roughly equivalently in both test and training performance. This
demonstrates that our full-batch model for BN accurately captures the same
dynamics as BN as used in the mini-batch case.

Figure 5: Comparison of different batch sizes (bs) for four-layer CNN with BN
layers, $(m_{l},\text{lr},\beta)=(1000,10^{-6},10^{-4})$. We demonstrate that
different batch sizes perform roughly equivalently in both training and test
convergence.

We note that we also perform an ablation study on batch sizes for the
nonconvex programs in Section 6 in Appendix A.7.2, and find the same
conclusion.

####  A.3 Additional Experiment: Singular Value Truncation Ablation

In order to make clear the effect of the singular value truncation as proposed
in Section 5.2 on generalization performance, we include an additional
ablation study on a new dataset to verify our intuition. In particular, we
evaluated the proposed methods on the CNAE-9 dataset (Ciarelli & Oliveira,
2009), a dataset of 1080 emails with 856 text features over 9 classes.
Following Ciarelli & Oliveira (2009), the first 900 samples were used as the
training set, while the final 180 samples are used as the test set in this
experiment, to obtain $\mathbf{X}\in\mathbb{R}^{900\times 856}$,
$\mathbf{y}\in\mathbb{R}^{900\times 9}$.

We evaluated the proposed approaches with
$(m_{1},\beta,\text{bs})=(10000,10^{-2},900)$, comparing the baseline non-
convex SGD approach to the convex truncated formulation for
$k\in\\{10,100,200,400,856\\}$, where $k$ is the index for which all smaller
singular values are truncated from the training dataset. We trained each model
for roughly an equal amount of time (15 seconds). For all architectures, a
learning rate of $10^{-5}$ was fixed.

We find in Figure 6, that as expected, singular value truncation can
effectively match the training and test performance of the SGD BN baseline,
but only in a certain range of $k$. When $k$ is too low, such as $k=10$,
valuable features of the data are eliminated and the convex program underfits,
whereas when $k$ is too large, such as $k=856$, the convex program overfits
and does not generalize well. Values of $k=100$ and $k=200$ generalize the
best for this problem, and their corresponding truncations preserve $80\%$ and
$91\%$ of the training data respectively. The exact value of $k$ which is
optimal for generalization is a hyper-parameter that can be tuned, though it
is clear that a large range of values can often work well. This ablation lends
additional support to the truncation approach proposed in Section 5.2.

(a) CNAE-9-Training accuracy

(b) CNAE-9-Test accuracy

Figure 6: Comparison of different truncation schemes of convex program
compared to SGD baseline, with
$(m_{1},\beta,\text{bs},\text{lr})=(10000,10^{-2},900,10^{-5})$, and
$k\in\\{10,100,200,400,856\\}$, where $k=856$ is simply the convex program
without truncation. Large values of $k$ overfit, while low values underfit,
and intermediate values match exactly the performance of SGD baseline.

####  A.4 Additional Experiment: Convolutional Architectures

To demonstrate the effectiveness of the convex formulation in convolutional
settings, we provide an additional experiment to show the success of this
approach. In particular, as the baseline architecture, we consider a two-layer
network which consists of a convolutional layer, followed by a channel-wise
BN, followed by a ReLU activation, then subsequently followed by an average
pooling operation, flattening, and a fully-connected layer. To use our convex
formulation (e.g. (15)) for such an architecture, we can split an input image
into spatial blocks (each block corresponds to a region whose outputs will be
averaged), and apply a global-average pooling model as in (15) to each spatial
block separately and combine the predictions from each spatial block together
to obtain one prediction. This suffices as a relaxation of the baseline
convolutional architecture.

To this end, we implement the non-convex and convex architectures on CIFAR-10
with the AdamW optimizer for 50 epochs. The convolution is with a $3\times 3$
kernel with padding of 2 and stride of 3, with average pooling to a spatial
dimension of $4\times 4$. The parameters chosen for this model are given by
$(n,m_{1},\beta,\text{bs},k,C,\text{lr})=(5\times 10^{4},1024,5\times
10^{-4},250,22,10,10^{-2})$.

We compare the results of the baseline, baseline-truncated, convex, and
convex-truncated approaches. In Table 1, we see that the accuracy of this CNN
approach is higher than the fully-connected methods we reported in our paper,
and as expected, there is a strong correspondence between the performance of
the convex and baseline models, again reiterating our strong theoretical
results, even in the case of this convex relaxation.

As illustrated in Table 1, truncation enables the convex approach to obtain
equivalent or higher test accuracy than its non-convex counterparts, while
obtaining similar training accuracy. However, notice that the accuracy
improvement provided by truncation is considerably less than the ones for
fully-connected layers (e.g. Figure 3(a)). This phenomenon is essentially due
the distribution of the singular values of the data matrix. Particularly, for
fully connected layers, we have $\mathbf{U}\in\mathbb{R}^{n\times d}$, and
therefore the dimension we truncate has $d=3072$ features. In this case, the
distribution of singular values follows an exponentially decaying pattern such
that the ratio of maximal and minimum singular values, also known as the
condition number of the data, is quite large
($\sigma_{\max}/\sigma_{\min}=5892.5$). On the contrary, for CNNs, we operate
on the patch matrices so that $\mathbf{U}_{mk}\in\mathbb{R}^{n\times h}$,
where $h=27$ is the size of each convolutional filter. Since the data matrix
is much more well conditioned ($\sigma_{\max}/\sigma_{\min}=271.2$) and has a
significantly smaller feature dimension, the impact of truncation is less
emphasized in our experiment.

Method | Train Accuracy | Test Accuracy
---|---|---
Baseline | 87.33 | 67.06
Baseline Truncated | 78.41 | 66.62
Convex (Ours) | 81.60 | 66.91
Convex Truncated (Ours) | 81.35 | 67.11
Table 1: Results of CIFAR-10 classification with two-layer ReLU CNN
architectures, with a $3\times 3$ kernel size, padding of 2, stride of 3, and
averaging to $4\times 4$ spatial dimensions, with
$(n,m_{1},\beta,\text{bs},k,C,\text{lr})=(5\times 10^{4},1024,5\times
10^{-4},250,22,10,10^{-2})$. Consistent with our results on fully-connected
architectures, the convex architecture matches or improves upon the
performance of the non-convex architecture, verifying our theoretical results.
We also observe the less pronounced impact of truncation in convolutional
experiments.

####  A.5 Experiments from Section 5.1–Evidence of Implicit Regularization

As described in the main paper, for both the standard BN architecture and
whitened architecture, we use the parameters
$(n,d,m_{1},\beta,\text{bs},\text{lr})=(50000,3072,1000,10^{-4},1000,10^{-4})$,
where bs stands for batch size and lr stands for learning rate. Both models
were trained with SGD for 501 epochs, and the learning rate was decayed by a
factor of 2 in the case of a training loss plateau. We now expand on the
definition of cosine similarity to clarify Figure 1 (in the main paper).

For weights $(\mathbf{W},\mathbf{W}^{\prime})$, we define the distance measure

$d_{\mathbf{X},i}(\mathbf{W},\mathbf{W}^{\prime}):=\frac{(\mathbf{v}_{i}^{\top}\mathbf{W})(\mathbf{v}_{i}^{\top}\mathbf{W}^{\prime})^{\top}}{\|\mathbf{v}_{i}^{\top}\mathbf{W}\|_{2}\|\mathbf{v}_{i}^{\top}\mathbf{W}^{\prime}\|_{2}}$
(18)

where $\mathbf{v}_{i}$ are the columns of $\mathbf{V}$ from the SVD of the
zero-mean data
$(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top})\mathbf{X}$. This
distance metric thus computes the cosine similarity between $\mathbf{W}$ and
$\mathbf{W}^{\prime}$ after being projected to the right singular subspace
$\mathbf{v}_{i}$.

For the same networks trained in Figure 2 (in the main paper), we compute
$\\{d_{\mathbf{X},i}({\mathbf{W}^{(1)}}^{(k)},{\mathbf{W}^{(1)}}^{(0)})\\}_{i=1}^{d}$
for the BN model (9), and
$\\{d_{\mathbf{U},i}({\mathbf{Q}}^{(k)},\mathbf{Q}^{(0)})\\}_{i=1}^{d}$ for
the equivalent whitened model (16), for each epoch $k$. We display these in
Figure 1 (in the main paper). These correspond to a measure of how far the
hidden-layer weights move from their initial values over time in the right
singular subspace of the data, with lower values indicating that the weights
have moved further from their initialization.

####  A.6 Verification of Closed-Form Solutions

See Figure 4 for the results of this experiment. The GD baseline model is
trained with a learning rate of $10^{-5}$, and the learning rate was decayed
by a factor of 2 in the case of a training loss plateau. The closed-form
solution takes only 0.578 seconds to compute, whereas we run the GD baseline
model for 40001 epochs, which takes approximately 50 minutes to run per trial.
Each trial behaves slightly differently due to the randomness of the
initialization of each network.

In terms of the generalization performance, the closed-form expression obtains
a test three-class classification accuracy of 47.7%, whereas the baseline GD
networks achieve an average of 60.1% final test accuracy, again indicating the
implicit regularization effect imposed by GD applied to BN networks.

####  A.7 Exact Convex Formulation and Implicit Regularization

In this section, we describe additional details about our main experimental
results of the main paper. In particular, we show the effect of training all
models for a longer duration than displayed in the main paper, the different
learning rates considered, the different batch sizes considered, along with
the effect of using Adam as opposed to SGD, and comparison to a non-BN
architecture.

For these experiments, we trained each model (SGD, SGD-Truncated, Convex,
Convex-Truncated) for roughly the same amount of time, and considered a range
of learning rates for each problem. In particular, for CIFAR-10, we considered
the learning rates $\\{10^{-5},5\times 10^{-5},10^{-4}\\}$, and trained the
SGD, SGD-Truncated, Convex, and Convex-Truncated models for 501, 501, 51, and
201 epochs respectively, and each model was trained for roughly 700 seconds.
In contrast, for CIFAR-100, we considered the learning rates
$\\{10^{-5},5\times 10^{-5},10^{-4}\\}$, and trained the SGD, SGD-Truncated,
Convex, and Convex-Truncated models for 621,1001,41, and 251 epochs
respectively, and each model was trained for roughly 800 seconds.

The range of learning rates was chosen such that all models reached their peak
test accuracy within the allotted time, and the learning rate with the best
final test accuracy was chosen for each model. In particular for CIFAR-10, the
chosen learning rates for the SGD, SGD-Truncated, Convex, and Convex-Truncated
models were $10^{-5}$, $10^{-5}$, $5\times 10^{-5}$, and $10^{-4}$,
respectively. In contrast, for CIFAR-100 the chosen learning rates for the
SGD, SGD-Truncated, Convex, and Convex-Truncated models were $5\times
10^{-5}$, $10^{-4}$, $10^{-4}$, and $5\times 10^{-4}$, respectively.

For our CIFAR-10 experiments in this section, we always use
$(m_{1},\beta,\text{bs})=(4096,10^{-4},10^{3})$, while for CIFAR-100, we used
$(m_{1},\beta,\text{bs})=(1024,10^{-4},10^{3})$, where bs stands for batch
size, as described in the main paper. The only exception is during our batch
size ablation study in Section A.7.2, in which case the batch size is varied.
All models are trained with SGD, except in Section A.7.4.

#####  A.7.1 Training Curves and Overfitting

(a) CIFAR-10-Training accuracy

(b) CIFAR-10-Test accuracy

(c) CIFAR-100-Training accuracy

(d) CIFAR-100-Test accuracy

Figure 7: Here, we provide a longer trained version of the CIFAR-10 and
CIFAR-100 experiments in Figure 3 (in the main paper). We see that training
for longer does not affect the results, as all models begin to overfit, and
they overfit at the same rate. For CIFAR-100, we observe that the baseline
overfits so much that the standard Convex model eventually overtakes it in
generalization performance.

To make more robust our claims in the main paper, we consider training each
network for a longer duration, such that each network can near 100% training
accuracy. We plot the results of both training and test accuracy for CIFAR-10
and CIFAR-100 in Figure 7. We note that extending training for longer only
seems to increase over-fitting, and increasing train accuracy to 100% does not
improve model performance. Therefore, the results of our main paper, which
display the curves for a shorter duration, represent the peak performance of
all algorithms before they begin to overfit.

#####  A.7.2 Impact of Different Batch Sizes

(a) CIFAR-10

(b) CIFAR-100

Figure 8: Batch size (bs) ablation study, considering
$\text{bs}\in\\{100,500,1000\\}$ for SGD (Baseline) and SGD-Truncated for the
experiment in Figure 3 (in the main paper). We see that smaller batch sizes do
not improve the performance of SGD. Therefore, we choose $\text{bs}=1000$ for
all the main experiments.

To further bolster our claim that taking the convex dual form from a full-
batch case of SGD is sufficient for capturing the performance of standard BN
networks, we compare the performance of Convex and Convex-Truncated with SGD
with different batch sizes, namely $\text{bs}\in\\{100,500\\}$. In Figure 8,
we show that different batch sizes do not improve the generalization
performance of BN networks, thereby validating that the convex dual derived
from the full-batch model is accurate at representing BN networks when trained
with mini-batch SGD. Therefore, for all of our other experiments, we use
$\text{bs}=1000$.

#####  A.7.3 Impact of Different Learning Rates

(a) CIFAR-10

(b) CIFAR-100

Figure 9: Learning rate (lr) ablation study for the experiment in Figure 3 (in
the main paper), considering $\text{lr}\in\\{10^{-5},5\times
10^{-5},10^{-4}\\}$ for SGD and SGD-Truncated. We see that larger lrs do not
improve the performance of SGD, because these models with larger lrs reach the
same peak accuracy but overfit quicker. For CIFAR-10, we selected learning
rates $10^{-5}$ for both SGD and SGD-Truncated. For CIFAR-100, we chose
learning rates $5\times 10^{-5}$ for SGD, and $10^{-4}$ for SGD-Truncated.

We plot test curves for all considered learning rates for each program for the
baseline program for CIFAR-10 and CIFAR-100 in Figure 9. We see that the
different learning rates considered do not significantly change the peak
performance, with higher learning rates simply overfitting earlier.

#####  A.7.4 Impact of Different Optimizers

(a) CIFAR-10

(b) CIFAR-100

Figure 10: A counterpart of the experiment in Figure 3 (in the main paper),
where we show test accuracies and all problems are solved with Adam. We see
that these results mirror those with SGD, where Convex-Truncated performs
equally well as or outperforms the baseline.

(a) CIFAR-10-Training accuracy

(b) CIFAR-10-Test accuracy

(c) CIFAR-100-Training accuracy

(d) CIFAR-100-Test accuracy

Figure 11: Here, we provide a longer trained version of the CIFAR-10 and
CIFAR-100 experiments in Figure 10. We see that training for longer does not
change the results, as all models begin to overfit.

(a) CIFAR-10

(b) CIFAR-100

Figure 12: Learning rate (lr) ablation study for the experiment in Figure 10,
where we consider $\text{lr}\in\\{10^{-5},5\times 10^{-5},10^{-4}\\}$ for
Adam. We see that larger learning rates do not improve the performance of Adam
or Convex-Truncated, because these models reach the same peak accuracy but
simply overfit earlier. For all models, we select a learning rate of $10^{-5}$
for the CIFAR-10 experiment in Figure 10(a). We choose a learning rate of
$10^{-5}$ for Adam and Convex-Truncated, and a learning rate of $10^{-4}$ for
Adam-Truncated for the CIFAR-100 experiment in Figure 10(b).

We also consider training both the standard non-convex and our convex
formulations with the Adam optimizer (Kingma & Ba, 2014). Not surprisingly, we
see a similar implicit regularization effect when applied to BN networks. For
the Adam experiments, we consider the same set of learning rates as the SGD
experiments, and fix hyper-parameters
$(\beta_{1},\beta_{2},\epsilon)=(0.9,0.999,10^{-8})$, as are the standard
default values. We run each program for the number of epochs as the SGD
experiments, and consider the same learning rates as the SGD experiments. For
CIFAR-10, the chosen learning rates for all the models was $10^{-5}$. For
CIFAR-100, the chosen learning rates for the Adam, Adam-Truncated, Convex, and
Convex-Truncated models were $10^{-5}$, $10^{-4}$, $5\times 10^{-5}$, and
$10^{-5}$, respectively.

We plot the results of the test curves in Figure 10, and plots of longer
trained curves in Figure 11, as well as the effect of different learning rates
in Figure 12. We find that with Adam, the results mirror that of SGD. For the
best learning rate chosen for each method, the truncated convex formulation
performs the best, followed by the truncated non-convex and standard non-
convex formulations, followed by the standard convex formulation. We can
conclude that our results are robust to the choice of optimizer.

#####  A.7.5 Comparison of BN to non-BN Two-Layer Architectures

(a) CIFAR-10

(b) CIFAR-100

Figure 13: Comparison of two-layer BN networks presented Figure 3 (in the main
paper), to the non-convex two-layer ReLU architecture without BN, and its
equivalent convex dual as presented in Pilanci & Ergen (2020). All problems
are solved with SGD and the learning rates are chosen based on test
performance. We find that both the SGD and Convex-Truncated programs improve
in generalization performance when adding BN, though the margin may vary
depending on the problem. We also observe that the Convex BN program without
truncation does not perform as well as the standard no-BN architecture.

One might speculate what the impact of adding BN to a two-layer network can
have, compared to a simple standard fully-connected architecture. In this
section, we compare the standard two-layer network in both its primal and
convex dual form, as first presented in Pilanci & Ergen (2020), to the BN
networks described in this paper. We let all parameters, i.e.
$(n,d,m_{1},\beta,\text{bs},C)$ be identical to the setting of the main paper,
and again consider a sweep of learning rates as was done for the BN
experiments. For the non-convex primal no BN form, the best performance came
from a learning rate of $5\times 10^{-7}$ on CIFAR-10 and $10^{-5}$ for
CIFAR-100, while for the convex no BN problem, the best performance was with a
learning rate of $5\times 10^{-10}$ on CIFAR-10 and $5\times 10^{-9}$ for
CIFAR-100.

We then compare the best results of each method in Figure 13. We see that for
both CIFAR-10 and CIFAR-100, introducing a BN layer can significantly improve
the performance of the two-layer network architecture, both when compared to
the result found with SGD, and the equivalent convex formulation.

###  Appendix B Additional Theoretical Results

In this section, we present additional theoretical results that are not
included in the main paper due to the page limit.

####  B.1 Extensions to Arbitrary Convex Loss Functions

In this section, we prove that our convex characterization holds for arbitrary
convex loss functions including cross entropy and hinge loss. We first restate
two-layer training problem with arbitrary convex loss functions after the
resealing in Lemma 1.1

$\displaystyle\min_{\theta\in\Theta_{s}}\mathcal{L}(f_{\theta,L}(\mathbf{X}),\mathbf{y})+\beta\left\|\mathbf{w}^{(L)}\right\|_{1}$
(19)

where $\mathcal{L}(\cdot,\mathbf{y})$ is a convex loss function. Then, the
dual of (19) is given by

$\displaystyle\max_{\mathbf{v}}-\mathcal{L}^{*}(\mathbf{v})\mbox{ s.t.
}\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{T}{\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(1)}\right)\right|\leq\beta,$

where $\mathcal{L}^{*}$ is the Fenchel conjugate function defined as Boyd &
Vandenberghe (2004)

$\displaystyle\mathcal{L}^{*}(\mathbf{v})=\max_{\mathbf{z}}\mathbf{z}^{T}\mathbf{v}-\mathcal{L}(\mathbf{z},\mathbf{y})\,.$

The analysis above proves that our convex characterization in the main paper
applies to arbitrary loss function. Therefore, optimal parameters for (3) are
a subset of the set that includes the extreme points of the dual constraint in
(5) independent of loss function.

####  B.2 Proof of Lemma 1.1

Here, we first prove that regularizing the intermediate BN parameters and
hidden layer weights don’t affect the optimization problem therefore the set
of optimal solutions for the problems in (20) and (22) are the same. Then, we
show that (22) can also be cast as an optimization problem with $\ell_{1}$
norm regularization penalty on the output layer weights $\mathbf{w}_{L}$.

#####  B.2.1 Effective Regularization in (3)

Here, we prove our claim about the effective regularization. Let us restate
the training problem as follows

$\displaystyle
p_{Lr}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|\sum_{j=1}^{m_{L-1}}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)\right)_{+}{w}^{(L)}_{j}-\mathbf{y}\right\|_{2}^{2}+\frac{\beta}{2}\sum_{l=1}^{L}\left(\left\|\bm{\gamma}^{(l)}\right\|_{2}^{2}+\left\|\bm{\alpha}^{(l)}\right\|_{2}^{2}+\left\|\mathbf{W}^{(l)}\right\|_{F}^{2}\right),$
(20)

where we use $\bm{\gamma}^{(L)}=\bm{\alpha}^{(L)}=\mathbf{0}$ as dummy
variables for notational simplicity. Now, let us rewrite (20) as

$\displaystyle p_{Lr}^{*}=\min_{t\geq
0}\min_{\theta\in\Theta}\frac{1}{2}\left\|\sum_{j=1}^{m_{L-1}}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)\right)_{+}{w}^{(L)}_{j}-\mathbf{y}\right\|_{2}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{L-1}}\left({\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}+{{w}^{(L)}_{j}}^{2}\right)+\frac{\beta}{2}t$
$\displaystyle\text{s.t.
}\sum_{l=1}^{L-2}\left(\left\|\bm{\gamma}^{(l)}\right\|_{2}^{2}+\left\|\bm{\alpha}^{(l)}\right\|_{2}^{2}+\left\|\mathbf{W}^{(l)}\right\|_{F}^{2}\right)+\left\|\mathbf{W}^{(L-1)}\right\|_{F}^{2}\leq
t.$

After applying the scaling in Lemma 1.1, we then take the dual with respect to
$\mathbf{w}^{(L)}$ as in the main paper to obtain the following problem

$\displaystyle p_{Lr}^{*}\geq d_{Lr}^{*}:=\max_{t\geq
0}\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}+\frac{\beta}{2}t$
(21) $\displaystyle\text{ s.t.
}\max_{\theta\in\Theta_{sr}}\left|\mathbf{v}^{T}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}\right)\right)_{+}\right|\leq\beta.$

where
$\Theta_{sr}:=\\{\theta\in\Theta:{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}=1,\,\forall
j\in[m_{L-1}],\,\sum_{l=1}^{L-2}\left(\left\|\bm{\gamma}^{(l)}\right\|_{2}^{2}+\left\|\bm{\alpha}^{(l)}\right\|_{2}^{2}+\left\|\mathbf{W}^{(l)}\right\|_{F}^{2}\right)+\left\|\mathbf{W}^{(L-1)}\right\|_{F}^{2}\leq
t\\}$. Since

$\displaystyle{\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)=\underbrace{\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\|_{2}}}_{\mathbf{h}(\theta^{\prime})}\gamma^{(L-1)}_{j}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(L-1)}_{j},$

where $\theta^{\prime}$ denotes all the parameters except
${\bm{\gamma}^{(L-1)}},{\bm{\alpha}^{(L-1)}},{\mathbf{w}^{(L)}}$. Then,
independent of the value $t$, $\mathbf{h}(\theta^{\prime})$ is always a unit
norm vector. Therefore, the maximization constraint in (21) is independent of
the norms of the parameters in $\theta^{\prime}$. Consequently, regularizing
the weights in $\theta^{\prime}$ does not affect the dual characterization in
(21).

Based on this observation, (20) can be equivalently stated as

$\displaystyle
p_{L}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|f_{\theta,L}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{L-1}}\left({\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}+{{w}^{(L)}_{j}}^{2}\right)$
(22)

#####  B.2.2 Rescaling for $\ell_{1}$ norm Regularization in (4)

We first note that similar approaches are also presented in Neyshabur et al.
(2014); Savarese et al. (2019); Pilanci & Ergen (2020); Ergen & Pilanci
(2020a; b).

For any $\theta\in\Theta$, we can rescale the parameters as
$\bar{\gamma}^{(L-1)}_{j}=\eta_{j}\gamma^{(L-1)}_{j}$,
$\bar{\alpha}^{(L-1)}_{j}=\eta_{j}\alpha^{(L-1)}_{j}$ and
${\bar{w}}^{(L)}_{j}={w}^{(L)}_{j}/\eta_{j}$, for any $\eta_{j}>0$. Then, the
network output becomes

$\displaystyle f_{\bar{\theta},L}(\mathbf{X})$
$\displaystyle=\sum_{j=1}^{m_{L-1}}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\|_{2}}\bar{\gamma}^{(L-1)}_{j}+\mathbf{1}\bar{\alpha}^{(L-1)}_{j}\right)_{+}{\bar{w}}^{(L)}$
$\displaystyle=\sum_{j=1}^{m_{L-1}}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\|_{2}}\gamma^{(L-1)}_{j}+\mathbf{1}\alpha^{(L-1)}_{j}\right)_{+}{w}^{(L)}$
$\displaystyle=f_{{\theta},L}(\mathbf{X}),$

which proves $f_{\theta,2}(\mathbf{X})=f_{\bar{\theta},2}(\mathbf{X})$.
Moreover, we use the following inequality

$\displaystyle\frac{1}{2}\sum_{j=1}^{m_{L-1}}\left({\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}+{{w}^{(L)}_{j}}^{2}\right)\geq\sum_{j=1}^{m_{L-1}}\left|{{w}^{(L)}_{j}}\right|\sqrt{{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}},$

where the equality is achieved when

$\eta_{j}=\left(\frac{|{w}^{(L)}_{j}|}{\sqrt{{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}}}\right)^{\frac{1}{2}}$

is used. Since this scaling does not alter the right-hand side of the
inequality, without loss of generality, we can set
$\sqrt{{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}}=1,\forall
j\in[m_{L-1}]$. Therefore, the right-hand side becomes
$\|\mathbf{w}^{(L)}\|_{1}$ and (22) is equivalent to

$\displaystyle\min_{\theta\in\Theta_{s}}\frac{1}{2}\left\|f_{\theta,L}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\beta\|\mathbf{w}^{(L)}\|_{1}$

where
$\Theta_{s}:=\\{\theta\in\Theta:{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}=1,\,\forall
j\in[m_{L-1}]\\}$. ∎

####  B.3 Derivation of the Dual Problem in (5)

Let us start with reparameterizing the scaled primal problem in (4) as

$\displaystyle
p_{L}^{*}=\min_{\theta\in\Theta_{s},\hat{\mathbf{y}}\in\mathbb{R}^{n}}\frac{1}{2}\left\|\hat{\mathbf{y}}-\mathbf{y}\right\|_{2}^{2}+\beta\left\|\mathbf{w}^{(L)}\right\|_{1}\text{
s.t.
}\hat{\mathbf{y}}=\sum_{j=1}^{m_{L-1}}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)\right)_{+}{w}^{(L)}_{j}.$
(23)

where
$\Theta_{s}:=\\{\theta\in\Theta:{\gamma^{(L-1)}_{j}}^{2}+{\alpha^{(L-1)}_{j}}^{2}=1,\,\forall
j\in[m_{L-1}]\\}$. We now form the following Lagrangian

$\displaystyle
L(\mathbf{v},\hat{\mathbf{y}},\mathbf{w}^{(L)})=\frac{1}{2}\left\|\hat{\mathbf{y}}-\mathbf{y}\right\|_{2}^{2}-\mathbf{v}^{T}\hat{\mathbf{y}}+\beta\left\|\mathbf{w}^{(L)}\right\|_{1}+\mathbf{v}^{T}\sum_{j=1}^{m_{L-1}}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)\right)_{+}{w}^{(L)}_{j}.$

Then, the corresponding dual function is

$\displaystyle g(\mathbf{v})$
$\displaystyle=\min_{\hat{\mathbf{y}},\mathbf{w}^{(L)}}L(\mathbf{v},\hat{\mathbf{y}},\mathbf{w}^{(L)})$
$\displaystyle=-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\left|\mathbf{v}^{\top}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)\right)_{+}\right|\leq\beta,\forall
j\in[m_{L-1}].$

Hence, the dual of (23) with respect to $\hat{\mathbf{y}}$ and
$\mathbf{w}^{(L)}$ is

$\displaystyle
p_{L}^{*}=\min_{\theta\in\Theta_{s}}\max_{\mathbf{v}}g(\mathbf{v})=\min_{\theta\in\Theta_{s}}\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\left|\mathbf{v}^{\top}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}_{j}\right)\right)_{+}\right|\leq\beta,\forall
j\in[m_{L-1}].$

We then change the order of min-max to obtain the following lower bound

$\displaystyle p_{L}^{*}\geq
d_{L}^{*}:=\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{\top}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{A}^{(L-2)}\mathbf{w}^{(L-1)}\right)\right)_{+}\right|\leq\beta.$

∎

####  B.4 Hyperplane Arrangements

We first define the set of all hyperplane arrangements for the data matrix
$\mathbf{X}$ as

$\displaystyle\mathcal{H}:=\bigcup\big{\\{}\\{\text{sign}(\mathbf{X}\mathbf{w})\\}\,:\,\mathbf{w}\in{}^{d}\big{\\}},$

where $|\mathcal{H}|\leq 2^{n}$. We now define a new set to denote the indices
with positive signs for each element in the set $\mathcal{H}$ as
$\mathcal{S}:=\big{\\{}\\{\cup_{h_{i}=1}\\{i\\}\\}\,:\,\mathbf{h}\in\mathcal{H}\big{\\}}$.
With this definition, we note that given an element $S\in\mathcal{S}$, one can
introduce a diagonal matrix $\mathbf{D}(S)\in{}^{n\times n}$ defined as
$\mathbf{D}(S)_{ii}:=\mathbbm{1}[i\in S]$. Therefore, the output of ReLU
activation can be equivalently written as
$\left(\mathbf{X}\mathbf{w}\right)_{+}=\mathbf{D}(S)\mathbf{X}\mathbf{w}$
provided that $\mathbf{D}(S)\mathbf{X}\mathbf{w}\geq 0$ and
$\left(\mathbf{I}_{n}-\mathbf{D}(S)\right)\mathbf{X}\mathbf{w}\leq 0$ are
satisfied. One can define more compactly these two constraints as
$\left(2\mathbf{D}(S)-\mathbf{I}_{n}\right)\mathbf{X}\mathbf{w}\geq 0$. We now
denote the cardinality of $\mathcal{S}$ as $P$, and obtain the following
upperbound

$\displaystyle P\leq 2\sum_{k=0}^{r-1}{n-1\choose k}\leq
2r\left(\frac{e(n-1)}{r}\right)^{r}\,$

where $r:=\mbox{rank}(\mathbf{X})\leq\min(n,d)$ (Ojha, 2000; Stanley et al.,
2004; Winder, 1966; Cover, 1965). For the rest of the paper, we enumerate all
possible hyperplane arrangements in an arbitrary order and denote them as
$\\{\mathbf{D}_{i}\\}_{i=1}^{P}$.

###  Appendix C Two-Layer Networks

####  C.1 Closed-Form Solution to the Dual Problem in (7)

######  Lemma C.1.

Suppose $n\leq d$ and $\mathbf{X}$ is full row-rank, then an optimal solution
to (7) is

$\displaystyle\mathbf{v}^{*}=\begin{cases}\beta\frac{\left(\mathbf{y}\right)_{+}}{\|\left(\mathbf{y}\right)_{+}\|_{2}}-\beta\frac{\left(-\mathbf{y}\right)_{+}}{\|\left(-\mathbf{y}\right)_{+}\|_{2}}&\text{
if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\left(\mathbf{y}\right)_{+}-\beta\frac{\left(-\mathbf{y}\right)_{+}}{\|\left(-\mathbf{y}\right)_{+}\|_{2}}&\text{
if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\beta\frac{\left(\mathbf{y}\right)_{+}}{\|\left(\mathbf{y}\right)_{+}\|_{2}}-\left(-\mathbf{y}\right)_{+}&\text{
if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\mathbf{y}&\text{ if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\end{cases}.$

###### Proof.

Since $n\leq d$ and $\mathbf{X}$ is full row-rank, we have

$\displaystyle\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{T}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}\|_{2}}\gamma^{(1)}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}\right)_{+}\right|$
$\displaystyle\leq\max\left\\{\|\left(\mathbf{v}\right)_{+}\|_{2},\|\left(-\mathbf{v}\right)_{+}\|_{2}\right\\}\max_{\begin{subarray}{c}\mathbf{1}^{T}\mathbf{u}=0\\\
\left\|\mathbf{u}\right\|_{2}=\left\|\mathbf{s}\right\|_{2}=1\end{subarray}}\left\|\left(\begin{bmatrix}\mathbf{u}&\frac{\mathbf{1}}{\sqrt{n}}\end{bmatrix}\mathbf{s}\right)_{+}\right\|_{2}$
$\displaystyle\leq\max\left\\{\|\left(\mathbf{v}\right)_{+}\|_{2},\|\left(-\mathbf{v}\right)_{+}\|_{2}\right\\}\max_{\begin{subarray}{c}\left\|\mathbf{u}\right\|_{2}=1\end{subarray}}\left\|\left(\mathbf{u}\right)_{+}\right\|_{2}$
$\displaystyle\leq\max\left\\{\|\left(\mathbf{v}\right)_{+}\|_{2},\|\left(-\mathbf{v}\right)_{+}\|_{2}\right\\},$

where $\mathbf{s}=\begin{bmatrix}\gamma^{(1)}&\alpha^{(1)}\end{bmatrix}^{T}$
and the equality is achieved when

$\displaystyle\mathbf{w}^{(1)}=\mathbf{X}^{\dagger}(\mathbf{v}-\min_{i}v_{i}),\;\mathbf{s}=\frac{1}{\|\mathbf{v}\|_{2}}\begin{bmatrix}\|\mathbf{v}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{v}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\mathbf{v}\end{bmatrix}.$

Thus, the dual problem (7) can be reformulated as

$\displaystyle\max_{\mathbf{v}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}\text{
s.t.
}\max\left\\{\|\left(\mathbf{v}\right)_{+}\|_{2},\|\left(-\mathbf{v}\right)_{+}\|_{2}\right\\}\leq\beta.$

The problem above can be maximized using the following optimal dual parameter

$\displaystyle\mathbf{v}^{*}=\begin{cases}\beta\frac{\left(\mathbf{y}\right)_{+}}{\|\left(\mathbf{y}\right)_{+}\|_{2}}-\beta\frac{\left(-\mathbf{y}\right)_{+}}{\|\left(-\mathbf{y}\right)_{+}\|_{2}}&\text{
if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\left(\mathbf{y}\right)_{+}-\beta\frac{\left(-\mathbf{y}\right)_{+}}{\|\left(-\mathbf{y}\right)_{+}\|_{2}}&\text{
if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\beta\frac{\left(\mathbf{y}\right)_{+}}{\|\left(\mathbf{y}\right)_{+}\|_{2}}-\left(-\mathbf{y}\right)_{+}&\text{
if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\mathbf{y}&\text{ if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\end{cases}.$

∎

####  C.2 Proof of Theorem 2.1

For the solution in Lemma C.1, the corresponding extreme points for the two-
sided constraint are

$\displaystyle{\mathbf{W}^{(1)}}^{*}=\begin{cases}\left(\mathbf{X}^{\dagger}\left(\mathbf{y}\right)_{+},\mathbf{X}^{\dagger}\left(-\mathbf{y}\right)_{+}\right)&\text{
if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\mathbf{X}^{\dagger}\left(-\mathbf{y}\right)_{+}&\text{ if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\mathbf{X}^{\dagger}\left(\mathbf{y}\right)_{+}&\text{ if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\mathbf{0}&\text{ if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\end{cases}$

and

$\displaystyle\begin{bmatrix}{\gamma^{(1)}_{1}}^{*}\\\
{\alpha^{(1)}_{1}}^{*}\end{bmatrix}=\frac{1}{\|\left(\mathbf{y}\right)_{+}\|_{2}}\begin{bmatrix}\|\left(\mathbf{y}\right)_{+}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\left(\mathbf{y}\right)_{+}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\left(\mathbf{y}\right)_{+}\end{bmatrix}\quad\begin{bmatrix}{\gamma^{(1)}_{2}}^{*}\\\
{\alpha^{(1)}_{2}}^{*}\end{bmatrix}=\frac{1}{\|\left(-\mathbf{y}\right)_{+}\|_{2}}\begin{bmatrix}\|\left(-\mathbf{y}\right)_{+}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\left(-\mathbf{y}\right)_{+}\end{bmatrix}.$

We now substitute these solutions into the primal problem (6) and then take
the derivative with respect to the output layer weights $\mathbf{w}^{(2)}$.
Since the primal problem is linear and convex provided that the first layer
weights are fixed as ${\mathbf{W}^{(1)}}^{*}$, we obtain the following closed-
form solutions for the output layer weights

$\displaystyle{\mathbf{w}^{(2)}}^{*}=\begin{cases}\begin{bmatrix}\|\left(\mathbf{y}\right)_{+}\|_{2}-\beta\\\
-\|\left(-\mathbf{y}\right)_{+}\|_{2}+\beta\end{bmatrix}&\text{ if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil-\|\left(-\mathbf{y}\right)_{+}\|_{2}+\beta&\text{ if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta\leq\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil\|\left(\mathbf{y}\right)_{+}\|_{2}-\beta&\text{ if
}\beta\leq\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\\\
\hfil 0&\text{ if
}\beta>\|\left(\mathbf{y}\right)_{+}\|_{2},\;\beta>\|\left(-\mathbf{y}\right)_{+}\|_{2}\end{cases}.$

These set of weights,
$\\{{\mathbf{W}^{(1)}}^{*},{\gamma^{(1)}}^{*},{\alpha^{(1)}}^{*},{\mathbf{w}^{(2)}}^{*}\\}$,
achieves $p_{2}^{*}=d_{2}^{*}$, therefore, strong duality holds.

∎

####  C.3 Proof of Theorem 2.2

We first note that if the number of neurons satisfies $m_{1}\geq m^{*}$ for
some $m^{*}\in\mathbb{N}$, where $m^{*}\leq n+1$, then strong duality holds,
i.e., $p_{2}^{*}=d_{2}^{*}$ (see Section 2 of Pilanci & Ergen (2020)). In the
weakly regularized regime where the network interpolates the training data,
there are even simpler ways to verify this upperbound, e.g., Zhang et al.
(2021) show how to construct a solution that exactly interpolates the training
data with $m_{1}=n$, and therefore the ReLU network fits any real valued
$\mathbf{y}\in\mathbb{R}^{n}$ with $n$ neurons.

Following the steps in Pilanci & Ergen (2020), we first focus on the dual
constraint in (7)

$\displaystyle\max_{\begin{subarray}{c}\mathbf{w}^{(1)}\\\
{\gamma^{(1)}}^{2}+{\alpha^{(1)}}^{2}=1\end{subarray}}\left|\mathbf{v}^{T}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}\|_{2}}\gamma^{(1)}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}\right)_{+}\right|$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{q}\right\|_{2}=\left\|\mathbf{s}\right\|_{2}=1\end{subarray}}\left|\mathbf{v}^{T}\left(\begin{bmatrix}\mathbf{U}\mathbf{q}&\frac{\mathbf{1}}{\sqrt{n}}\end{bmatrix}\mathbf{s}\right)_{+}\right|$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\left|\mathbf{v}^{T}\left(\mathbf{U}^{\prime}\mathbf{s}^{\prime}\right)_{+}\right|$
$\displaystyle=\max_{i\in[P]}\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\\\
(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}^{\prime}\geq
0\end{subarray}}\left|\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{s}^{\prime}\right|,$

where
$\mathbf{q}:=\bm{\Sigma}\mathbf{V}^{T}\mathbf{w}^{(1)}/\left\|\bm{\Sigma}\mathbf{V}^{T}\mathbf{w}^{(1)}\right\|_{2}$,
$\mathbf{s}:=\begin{bmatrix}\gamma^{(1)}&\alpha^{(1)}\end{bmatrix}^{T}$,
$\mathbf{s}^{\prime}:=\begin{bmatrix}\gamma^{(1)}\mathbf{q}^{T}&\alpha^{(1)}\end{bmatrix}^{T}$,
and
$\mathbf{U}^{\prime}:=\begin{bmatrix}\mathbf{U}&\frac{\mathbf{1}}{\sqrt{n}}\end{bmatrix}$.
Now, let us consider the single side of this constraint and introduce a dual
variable $\bm{\rho}$ for the constraint as follows

$\displaystyle\max_{i\in[P]}\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\\\
(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}^{\prime}\geq
0\end{subarray}}\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{s}^{\prime}$
$\displaystyle=\max_{i\in[P]}\min_{\bm{\rho}\geq
0}\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\left(\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right)\mathbf{s}^{\prime}$
$\displaystyle=\max_{i\in[P]}\min_{\bm{\rho}\geq
0}\left\|\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2}.$

Therefore, we have

$\displaystyle\max_{\theta\in\Theta_{s}}\mathbf{v}^{T}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}\|_{2}}\gamma^{(1)}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}\right)_{+}\leq\beta\iff$
$\displaystyle\forall i\in[P],\min_{\bm{\rho}\geq
0}\left\|\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2}\leq\beta$
$\displaystyle\iff$ $\displaystyle\forall
i\in[P],\,\exists\bm{\rho}_{i}\geq\mathbf{0}\text{ s.t.
}\left\|\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2}\leq\beta.$

Similar arguments also hold for the negative side of the absolute value
constraint. Therefore, we can reformulate (7) as a convex optimization problem
with variables
$\mathbf{v},\bm{\rho}_{i},\bm{\rho}_{i}^{\prime}\in\mathbb{R}^{n},\forall
i\in[P]$ and $2P$ constraints as follows

$\displaystyle d_{2}^{*}=$
$\displaystyle\max_{\begin{subarray}{c}\mathbf{v},\bm{\rho}_{i}\end{subarray}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}~{}\text{
s.t.
}\begin{split}&\left\|\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2}\leq\beta\\\
&\left\|-\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{\prime^{T}}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2}\leq\beta\end{split},\;\forall
i\in[P].$ (24)

As noted in Pilanci & Ergen (2020), since the problem in (24) is strictly
feasible for the particular choice of parameters
$\mathbf{v}=\bm{\rho}_{i}=\bm{\rho}_{i}^{\prime}=\mathbf{0},\forall\in[P]$,
Slater’s condition and therefore strong duality holds (Boyd & Vandenberghe,
2004). We now write (24) in a Lagrangian form with the parameters
$\bm{\lambda},\bm{\lambda}^{\prime}\in\mathbb{R}^{P}$ as

$\displaystyle
d_{2}^{*}=\min_{\bm{\lambda},\bm{\lambda}^{\prime}\geq\mathbf{0}}\max_{\begin{subarray}{c}\mathbf{v}\\\
\bm{\rho}_{i}\geq\mathbf{0},\forall
i\end{subarray}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}$
$\displaystyle+\sum_{i=1}^{P}\lambda_{i}(\beta-\left\|\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2})$
$\displaystyle+\sum_{i=1}^{P}\lambda_{i}^{\prime}(\beta-\left\|-\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{\prime^{T}}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right\|_{2}).$
(25)

We next introduce additional variables
$\mathbf{r}_{i},\mathbf{r}_{i}^{\prime}\in\mathbb{R}^{d},\forall i\in[P]$ to
represent (C.3) as

$\displaystyle
d_{2}^{*}=\min_{\bm{\lambda},\bm{\lambda}^{\prime}\geq\mathbf{0}}\max_{\begin{subarray}{c}\mathbf{v}\\\
\bm{\rho}_{i}\geq\mathbf{0},\forall
i\end{subarray}}\min_{\begin{subarray}{c}\mathbf{r}_{i}:\|\mathbf{r}_{i}\|_{2}\leq
1\\\ \mathbf{r}_{i}^{\prime}\>\|\mathbf{r}_{i}^{\prime}\|_{2}\leq
1\end{subarray}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}$
$\displaystyle+\sum_{i=1}^{P}\lambda_{i}\left(\beta-\left(\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right)\mathbf{r}_{i}\right)$
$\displaystyle+\sum_{i=1}^{P}\lambda_{i}^{\prime}\left(\beta+\left(-\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{\prime^{T}}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right)\mathbf{r}_{i}^{\prime}\right).$
(26)

Now, we remark that the objective function in (C.3) is concave in
$\mathbf{v},\bm{\rho}_{i},\bm{\rho}_{i}^{\prime}$ and convex in
$\mathbf{r}_{i},\mathbf{r}_{i}^{\prime},\forall i\in[P]$. In addition, the
constraint set for the additional variables, i.e.,
$\|\mathbf{r}_{i}\|_{2},\|\mathbf{r}_{i}^{\prime}\|_{2}\leq 1$, are convex and
compact. Then, we use Sion’s minimax theorem (Sion, 1958) to change the order
of the inner max-min problems as

$\displaystyle
d_{2}^{*}=\min_{\bm{\lambda},\bm{\lambda}^{\prime}\geq\mathbf{0}}\min_{\begin{subarray}{c}\mathbf{r}_{i}:\|\mathbf{r}_{i}\|_{2}\leq
1\\\ \mathbf{r}_{i}^{\prime}\>\|\mathbf{r}_{i}^{\prime}\|_{2}\leq
1\end{subarray}}\max_{\begin{subarray}{c}\mathbf{v}\\\
\bm{\rho}_{i}\geq\mathbf{0},\forall
i\end{subarray}}-\frac{1}{2}\|\mathbf{v}-\mathbf{y}\|_{2}^{2}+\frac{1}{2}\|\mathbf{y}\|_{2}^{2}$
$\displaystyle+\sum_{i=1}^{P}\lambda_{i}\left(\beta-\left(\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{T}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right)\mathbf{r}_{i}\right)$
$\displaystyle+\sum_{i=1}^{P}\lambda_{i}^{\prime}\left(\beta+\left(-\mathbf{v}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}+\bm{\rho}_{i}^{\prime^{T}}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\right)\mathbf{r}_{i}^{\prime}\right).$
(27)

Thus, we can now perform the maximization over the variables
$\mathbf{v},\bm{\rho}_{i},\bm{\rho}_{i}^{\prime}$ to simplify (C.3) into the
following form

$\displaystyle p_{2c}^{*}=$
$\displaystyle\min_{\bm{\lambda},\bm{\lambda}^{\prime}\geq\mathbf{0}}\min_{\begin{subarray}{c}\mathbf{r}_{i}:\|\mathbf{r}_{i}\|_{2}\leq
1\\\ \mathbf{r}_{i}^{\prime}\>\|\mathbf{r}_{i}^{\prime}\|_{2}\leq
1\end{subarray}}\frac{1}{2}\left\|\sum_{i=1}^{P}\left(\lambda_{i}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{r}_{i}-\lambda_{i}^{\prime}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{r}_{i}^{\prime}\right)-\mathbf{y}\right\|_{2}^{2}+\beta\sum_{i=1}^{P}(\lambda_{i}+\lambda_{i}^{\prime})$
(28) $\displaystyle\text{ s.t.
}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{r}_{i}\geq
0,\;(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{r}_{i}^{\prime}\geq
0,\;\forall i\in[P].$

We then apply a set variable changes as
$\mathbf{s}_{i}=\lambda_{i}\mathbf{r}_{i}$ and
$\mathbf{s}_{i}^{\prime}=\lambda_{i}^{\prime}\mathbf{r}_{i}^{\prime},\,\forall
i\in[P]$ to obtain the following problem

$\displaystyle p_{2c}^{*}=$
$\displaystyle\min_{\bm{\lambda},\bm{\lambda}^{\prime}\geq\mathbf{0}}\min_{\begin{subarray}{c}\mathbf{s}_{i}:\|\mathbf{s}_{i}\|_{2}\leq\lambda_{i}\\\
\mathbf{s}_{i}^{\prime}\>\|\mathbf{s}_{i}^{\prime}\|_{2}\leq\lambda_{i}^{\prime}\end{subarray}}\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}(\mathbf{s}_{i}-\mathbf{s}_{i}^{\prime})-\mathbf{y}\right\|_{2}^{2}+\beta\sum_{i=1}^{P}(\lambda_{i}+\lambda_{i}^{\prime})$
(29) $\displaystyle\text{ s.t.
}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}\geq
0,\;(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}^{\prime}\geq
0,\;\forall i\in[P].$

Finally, since the optimality conditions require
$\|\mathbf{s}_{i}\|_{2}=\lambda_{i}$ and
$\|\mathbf{s}_{i}^{\prime}\|_{2}=\lambda_{i}^{\prime},\forall i\in[P]$, we
write (29) as a single convex minimization problem as

$\displaystyle p_{2c}^{*}=$
$\displaystyle\min_{\mathbf{s}_{i},\mathbf{s}_{i}^{\prime}}\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}(\mathbf{s}_{i}-\mathbf{s}_{i}^{\prime})-\mathbf{y}\right\|_{2}^{2}+\beta\sum_{i=1}^{P}(\|\mathbf{s}_{i}\|_{2}+\|\mathbf{s}_{i}^{\prime}\|_{2})$
(30) $\displaystyle\text{ s.t.
}(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}\geq
0,\;(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}_{i}^{\prime}\geq
0,\;\forall i\in[P].$

Now, we note that since (24) is convex and Slater’s condition holds as shown
above, we have $p_{2c}^{*}=d_{2}^{*}$. In addition, as proven in Pilanci &
Ergen (2020) strong duality also holds for the original non-convex problem,
i.e., $p_{2}^{*}=d_{2}^{*}$, therefore, we get
$p_{2}^{*}=p_{2c}^{*}=d_{2}^{*}$.

Constructing an optimal solution to the primal problem:
One can also construct an optimal solution to the non-convex training problem
in (6) as follows

$\displaystyle\left({\mathbf{w}^{(1)}_{j_{1i}}}^{*},{\gamma^{(1)}_{j_{1i}}}^{*},{\alpha^{(1)}_{j_{1i}}}^{*},{{w}^{(2)}_{j_{1i}}}^{*}\right)=\bigg{(}\mathbf{V}\mathbf{\Sigma}^{-1}\mathbf{s}_{i,1:r}^{*},\frac{\|\mathbf{s}_{i,1:r}^{*}\|_{2}}{||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}},\frac{{s}_{i,r+1}^{*}}{||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}},||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}\bigg{)},\mbox{
if }\mathbf{s}_{i}^{*}\neq 0$
$\displaystyle\left({\mathbf{w}^{(1)}_{j_{2i}}}^{*},{\gamma^{(1)}_{j_{2i}}}^{*},{\alpha^{(1)}_{j_{2i}}}^{*},{{w}^{(2)}_{j_{2i}}}^{*}\right)=\bigg{(}\mathbf{V}\mathbf{\Sigma}^{-1}\mathbf{s}_{i,1:r}^{\prime^{*}},\frac{||\mathbf{s}_{i,1:r}^{\prime^{*}}||_{2}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}},\frac{{s}_{i,r+1}^{\prime^{*}}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}},-||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}\bigg{)},\mbox{
if }\mathbf{s}_{i}^{\prime^{*}}\neq 0,$

where $\\{\mathbf{s}_{i}^{*},\mathbf{s}_{i}^{\prime^{*}}\\}_{i=1}^{P}$ are the
optimal solutions to (8) and and $j_{si}\in[|\mathcal{J}_{s}|]$ given the
definitions $\mathcal{J}_{1}:=\\{i:\mathbf{s}_{i}^{*}\neq 0\\}$ and
$\mathcal{J}_{2}:=\\{i:\mathbf{s}_{i}^{\prime^{*}}\neq 0\\}$. We also note
that here the second subscript denotes the selected entries of the
corresponding vector, e.g., $\mathbf{s}_{i,1:r}$ denotes the first $r$ entries
of the vector $\mathbf{s}_{i}$. Thus, we have
$m^{*}=|\mathcal{J}_{1}|+|\mathcal{J}_{2}|$.

Verifying strong duality:
Next, we verify strong duality, i.e., $p_{2}^{*}=p_{2c}^{*}=d_{2}^{*}$, by
plugging the optimal solution to the non-convex objective in (3) and the
corresponding convex objective in (30).

Given the optimal solutions
$\\{\mathbf{s}_{i}^{*},\mathbf{s}_{i}^{\prime^{*}}\\}_{i=1}^{P}$ to (30), the
convex problem has the following objective

$\displaystyle p_{2c}^{*}=$
$\displaystyle\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}(\mathbf{s}_{i}^{*}-\mathbf{s}_{i}^{\prime^{*}})-\mathbf{y}\right\|_{2}^{2}+\beta\sum_{i=1}^{P}(\|\mathbf{s}_{i}^{*}\|_{2}+\|\mathbf{s}_{i}^{\prime^{*}}\|_{2})$
(31)

and the corresponding non-convex problem objective is as follows

$\displaystyle
p_{2}^{*}=\frac{1}{2}\left\|f_{\theta,2}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{*^{2}}+{\alpha^{(1)}_{j}}^{*^{2}}+{{w}^{(2)}_{j}}^{*^{2}}\right).$
(32)

We first prove that the first terms in both objectives have the same value. To
do that, we compute the output of the non-convex neural network with the
optimal weight construction above

$\displaystyle f_{\theta,2}(\mathbf{X})$
$\displaystyle=\sum_{j=1}^{m_{1}}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}{\mathbf{w}^{(1)}_{j}}^{*}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}{\mathbf{w}^{(1)}_{j}}^{*}\|_{2}}{\gamma^{(1)}_{j}}^{*}+\frac{\mathbf{1}}{\sqrt{n}}{\alpha^{(1)}_{j}}^{*}\right)_{+}{{w}^{(2)}_{j}}^{*}$
$\displaystyle=\sum_{i\in\mathcal{J}_{1}}\left(\frac{\mathbf{U}\mathbf{s}_{i,1:r}^{*}}{\|\mathbf{U}\mathbf{s}_{i,1:r}^{*}\|_{2}}\frac{||\mathbf{s}_{i,1:r}^{*}||_{2}}{||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}}+\frac{\mathbf{1}}{\sqrt{n}}\frac{{s}_{i,r+1}^{*}}{||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}}\right)_{+}||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}$
$\displaystyle-\sum_{i\in\mathcal{J}_{2}}\left(\frac{\mathbf{U}\mathbf{s}_{i,1:r}^{\prime^{*}}}{\|\mathbf{U}\mathbf{s}_{i,1:r}^{\prime^{*}}\|_{2}}\frac{||\mathbf{s}_{i,1:r}^{\prime^{*}}||_{2}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}}+\frac{\mathbf{1}}{\sqrt{n}}\frac{{s}_{i,r+1}^{\prime^{*}}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}}\right)_{+}||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}$
$\displaystyle=\sum_{i\in\mathcal{J}_{1}}\left(\mathbf{U}^{\prime}\mathbf{s}_{i}^{*}\right)_{+}-\sum_{i\in\mathcal{J}_{2}}\left(\mathbf{U}^{\prime}\mathbf{s}_{i}^{\prime^{*}}\right)_{+}$
$\displaystyle=\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}(\mathbf{s}_{i}^{*}-\mathbf{s}_{i}^{\prime^{*}}),$
(33)

where the last equality follows from the constraints in (30).

Next, we prove that the regularizations terms have the same objective as
follows

$\displaystyle\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{*^{2}}+{\alpha^{(1)}_{j}}^{*^{2}}+{{w}^{(2)}_{j}}^{*^{2}}\right)$
$\displaystyle=\frac{\beta}{2}\sum_{j\in\mathcal{J}_{1}}\left(\left(\frac{||\mathbf{s}_{i,1:r}^{*}||_{2}}{||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}}\right)^{2}+\left(\frac{s_{i,r+1}^{*}}{||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}}\right)^{2}+\left(||\mathbf{s}_{i}^{*}||_{2}^{\frac{1}{2}}\right)^{2}\right)$
$\displaystyle+\frac{\beta}{2}\sum_{i\in\mathcal{J}_{2}}\left(\left(\frac{||\mathbf{s}_{i,1:r}^{\prime^{*}}||_{2}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}}\right)^{2}+\left(\frac{s_{i,r+1}^{\prime^{*}}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}}\right)^{2}+\left(-||\mathbf{s}_{i}^{\prime^{*}}||_{2}^{\frac{1}{2}}\right)^{2}\right)$
$\displaystyle=\frac{\beta}{2}\sum_{i\in\mathcal{J}_{1}}\left(\frac{||\mathbf{s}_{i,1:r}^{*}||_{2}^{2}}{||\mathbf{s}_{i}^{*}||_{2}}+\frac{s_{i,r+1}^{*^{2}}}{||\mathbf{s}_{i}^{*}||_{2}}+||\mathbf{s}_{i}^{*}||_{2}\right)$
$\displaystyle+\frac{\beta}{2}\sum_{i\in\mathcal{J}_{2}}\left(\frac{||\mathbf{s}_{i,1:r}^{\prime^{*}}||_{2}^{2}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}}+\frac{{s_{i,r+1}^{\prime^{*}}}^{2}}{||\mathbf{s}_{i}^{\prime^{*}}||_{2}}+||\mathbf{s}_{i}^{\prime^{*}}||_{2}\right)$
$\displaystyle=\beta\sum_{i\in\mathcal{J}_{1}}||\mathbf{s}_{i}^{*}||_{2}+\beta\sum_{i\in\mathcal{J}_{2}}||\mathbf{s}_{i}^{\prime^{*}}||_{2}$
$\displaystyle=\beta\sum_{i=1}^{P}(\|\mathbf{s}_{i}^{*}\|_{2}+\|\mathbf{s}_{i}^{\prime^{*}}\|_{2}).$
(34)

Therefore, by (C.3) and (C.3), the non-convex objective (32) and the convex
objective (31) are the same, i.e., $p_{2}^{*}=p_{2c}^{*}$. ∎

####  C.4 Proof of Lemma 2.1

We note that similar approaches are also presented in Sahiner et al. (2021a);
Ergen & Pilanci (2020b).

For any $\theta\in\Theta$, we can rescale the parameters as
$\bar{\gamma}^{(1)}_{j}=\eta_{j}\gamma^{(1)}_{j}$,
$\bar{\alpha}^{(1)}_{j}=\eta_{j}\alpha^{(1)}_{j}$ and
$\bar{\mathbf{w}}^{(2)}_{j}=\mathbf{w}^{(2)}_{j}/\eta_{j}$, for any
$\eta_{j}>0$. Then, the network output becomes

$\displaystyle f_{\bar{\theta},2}(\mathbf{X})$
$\displaystyle=\sum_{j=1}^{m_{1}}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}_{j}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}_{j}\|_{2}}\bar{\gamma}^{(1)}_{j}+\mathbf{1}\bar{\alpha}^{(1)}_{j}\right)_{+}{\bar{\mathbf{w}}_{j}}^{(2)^{T}}$
$\displaystyle=\sum_{j=1}^{m_{1}}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}_{j}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}_{j}\|_{2}}\gamma^{(1)}_{j}+\mathbf{1}\alpha^{(1)}_{j}\right)_{+}{\mathbf{w}^{(2)}_{j}}^{T}$
$\displaystyle=f_{{\theta},2}(\mathbf{X}),$

which proves $f_{\theta,2}(\mathbf{X})=f_{\bar{\theta},2}(\mathbf{X})$.
Moreover, we use the following inequality

$\displaystyle\frac{1}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}+\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}^{2}\right)\geq\sum_{j=1}^{m_{1}}\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}\sqrt{{\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}},$

where the equality is achieved when

$\eta_{j}=\left(\frac{\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}}{\sqrt{{\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}}}\right)^{\frac{1}{2}}$

is used. Since this scaling does not alter the right-hand side of the
inequality, without loss of generality, we can set
$\sqrt{{\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}}=1,\forall j\in[m_{1}]$.
Therefore, the right-hand side becomes
$\sum_{j=1}^{m_{1}}\|\mathbf{w}^{(2)}_{j}\|_{2}$.

We then note that due to the arguments in Appendix B.2.1, we can remove the
regularization on $\mathbf{w}^{(1)}_{j}$ without loss of generality.
Therefore, the final optimization problem is as follows

$\displaystyle
p_{2v}^{*}=\min_{\theta\in\Theta_{s}}\frac{1}{2}\left\|f_{\theta,2}(\mathbf{X})-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{j=1}^{m_{1}}\left\|\mathbf{w}^{(2)}_{j}\right\|_{2},$

where
$\Theta_{s}:=\\{\theta\in\Theta:{\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}=1,\,\forall
j\in[m_{1}]\\}$. ∎

####  C.5 Proof of Theorem 2.3

We first note that if the number of neurons satisfies $m_{1}\geq m^{*}$ for
some $m^{*}\in\mathbb{N}$, where $m^{*}\leq nC+1$, then strong duality holds,
i.e., $p_{2v}^{*}=d_{2v}^{*}$, (see Section A.3.1 of Sahiner et al. (2021a)
for details). Following the steps in Sahiner et al. (2021a) and Theorem 2.2,
we again focus on the dual constraint in (11)

$\displaystyle\max_{\begin{subarray}{c}\mathbf{w}^{(1)}\\\
{\gamma^{(1)}}^{2}+{\alpha^{(1)}}^{2}=1\end{subarray}}\left\|\mathbf{V}^{T}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}\|_{2}}\gamma^{(1)}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}\right)_{+}\right\|_{2}$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{q}\right\|_{2}=\left\|\mathbf{s}\right\|_{2}=1\end{subarray}}\left\|\mathbf{V}^{T}\left(\begin{bmatrix}\mathbf{U}\mathbf{q}&\frac{\mathbf{1}}{\sqrt{n}}\end{bmatrix}\mathbf{s}\right)_{+}\right\|_{2}$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\left\|\mathbf{V}^{T}\left(\mathbf{U}^{\prime}\mathbf{s}^{\prime}\right)_{+}\right\|_{2}$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{g}\right\|_{2}=\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\mathbf{g}^{T}\mathbf{V}^{T}\left(\mathbf{U}^{\prime}\mathbf{s}^{\prime}\right)_{+}$
$\displaystyle=\max_{i\in[P]}\max_{\begin{subarray}{c}\left\|\mathbf{g}\right\|_{2}=\left\|\mathbf{s}^{\prime}\right\|_{2}=1\\\
(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}^{\prime}\geq
0\end{subarray}}\mathbf{g}^{T}\mathbf{V}^{T}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{s}^{\prime}$
$\displaystyle=\max_{i\in[P]}\max_{\begin{subarray}{c}\left\|\mathbf{g}\right\|_{2}=\left\|\mathbf{s}^{\prime}\right\|_{2}=1\\\
(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}^{\prime}\geq
0\end{subarray}}\left\langle\mathbf{V},\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{s}^{\prime}\mathbf{g}^{T}\right\rangle$
$\displaystyle=\max_{i\in[P]}\max_{\begin{subarray}{c}\mathbf{S}\in\mathcal{K}_{i}\end{subarray}}\left\langle\mathbf{V},\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{S}\right\rangle,$

where
$\mathbf{s}^{\prime}=\begin{bmatrix}\gamma\mathbf{q}^{T}&\alpha\end{bmatrix}^{T}$,
$\mathbf{U}^{\prime}=\begin{bmatrix}\mathbf{U}&\frac{\mathbf{1}}{\sqrt{n}}\end{bmatrix}$,
and
$\mathcal{K}_{i}:=\mathrm{conv}\\{\mathbf{Z}=\mathbf{s}^{\prime}\mathbf{g}^{T}:(2\mathbf{D}_{i}-\mathbf{I}_{n})\mathbf{U}^{\prime}\mathbf{s}^{\prime}\geq
0,\,\|\mathbf{Z}\|_{*}\leq 1\\}$. Thus, we have

$\displaystyle\max_{\theta\in\Theta_{s}}\left\|\mathbf{V}^{T}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}\|_{2}}\gamma^{(1)}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}\right)_{+}\right\|_{2}\leq\beta\iff\forall
i\in[P],\,\exists\mathbf{S}_{i}\in\mathcal{K}_{i}\text{ s.t.
}\left\langle\mathbf{V},\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{S}_{i}\right\rangle\leq\beta.$

Similar arguments also hold for the negative side of the absolute value
constraint. We then directly follow the steps in Sahiner et al. (2021a) and
achieve the following convex program

$\displaystyle\min_{\mathbf{S}_{i}}\frac{1}{2}\left\|\sum_{i=1}^{P}\mathbf{D}_{i}\mathbf{U}^{\prime}\mathbf{S}_{i}-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{i=1}^{P}\|\mathbf{S}_{i}\|_{c_{i},*},$

where

$\displaystyle\|\mathbf{S}\|_{c_{i},*}:=\min_{t\geq 0}t\text{ s.t.
}\mathbf{S}\in t\mathcal{K}_{i}$

denotes a constrained version of the standard nuclear norm. ∎

####  C.6 Proof of Theorem 2.4

We first restate the dual problem as

$\displaystyle\max_{\mathbf{V}}-\frac{1}{2}\|\mathbf{V}-\mathbf{Y}\|_{F}^{2}+\frac{1}{2}\|\mathbf{Y}\|_{F}^{2}\mbox{
s.t.
}\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\left\|\mathbf{V}^{T}\left(\mathbf{U}^{\prime}\mathbf{s}^{\prime}\right)_{+}\right\|_{2}\leq\beta,$
(35)

where we use the results from Theorem 2.3 to modify the dual constraint. Since
$\mathbf{Y}$ is one-hot encoded, (35) can be rewritten as

$\displaystyle\max_{\mathbf{V}}-\frac{1}{2}\|\mathbf{V}-\mathbf{Y}\|_{F}^{2}+\frac{1}{2}\|\mathbf{Y}\|_{F}^{2}\mbox{
s.t.
}\max_{\left\|\mathbf{u}\right\|_{2}=1}\|\mathbf{V}^{T}\mathbf{u}\|_{2}\leq\beta.$
(36)

This problem has the following closed-form solution

$\displaystyle\mathbf{v}_{j}^{*}=\begin{cases}\beta\frac{\mathbf{y}_{j}}{\|\mathbf{y}_{j}\|_{2}}&\text{
if }\beta\leq\|\mathbf{y}_{j}\|_{2}\\\ \ \hfil\mathbf{y}_{j}&\text{ otherwise
}\end{cases},\quad\forall j\in[C].$ (37)

and the corresponding first layer weights in (36) are

$\displaystyle{\mathbf{w}^{(1)}_{j}}^{*}=\begin{cases}\mathbf{X}^{\dagger}\mathbf{y}_{j}&\text{
if }\beta\leq\|\mathbf{y}_{j}\|_{2}\\\ \ \hfil-&\text{ otherwise
}\end{cases},\quad\forall j\in[o].$ (38)

Now let us first denote the set of indices that achieves the extreme point of
the dual constraint as
$\mathcal{E}:=\\{j:\beta\leq\|\mathbf{y}_{j}\|_{2},j\in[C]\\}$. Then the dual
objective in (36) using the optimal dual parameter in (37)

$\displaystyle d_{2v}^{*}$
$\displaystyle=-\frac{1}{2}\|\mathbf{V}^{*}-\mathbf{Y}\|_{F}^{2}+\frac{1}{2}\|\mathbf{Y}\|_{F}^{2}$
$\displaystyle=-\frac{1}{2}\sum_{j\in\mathcal{E}}(\beta-\|\mathbf{y}_{j}\|_{2})^{2}+\frac{1}{2}\sum_{j=1}^{C}\|\mathbf{y}_{j}\|_{2}^{2}$
$\displaystyle=-\frac{1}{2}\beta^{2}|\mathcal{E}|+\beta\sum_{j\in\mathcal{E}}\|\mathbf{y}_{j}\|_{2}+\frac{1}{2}\sum_{j\notin\mathcal{E}}\|\mathbf{y}_{j}\|_{2}^{2}.$
(39)

We next restate the scaled primal problem

$\displaystyle
p_{2v}^{*}=\min_{\theta\in\Theta_{s}}\frac{1}{2}\left\|\sum_{j=1}^{m_{1}}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}_{j}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}\mathbf{w}^{(1)}_{j}\|_{2}}\gamma^{(1)}_{j}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}_{j}\right)_{+}{\mathbf{w}^{(2)}_{j}}^{T}-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{j=1}^{m_{1}}\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}$
(40)

and then obtain the optimal primal objective using the hidden layer weights in
(38), which yields the following optimal solution

$\displaystyle\begin{split}&\left({\mathbf{w}^{(1)}_{j}}^{*},{\mathbf{w}^{(2)}_{j}}^{*}\right)=\begin{cases}\left(\mathbf{X}^{\dagger}\mathbf{y}_{j},\left(\|\mathbf{y}_{j}\|_{2}-\beta\right)\mathbf{e}_{j}\right)&\text{
if }\beta\leq\|\mathbf{y}_{j}\|_{2}\\\ \hfil-&\text{ otherwise }\end{cases}\\\
&\begin{bmatrix}{\gamma^{(1)}_{j}}^{*}\\\
{\alpha^{(1)}_{j}}^{*}\end{bmatrix}=\frac{1}{\|\mathbf{y}_{j}\|_{2}}\begin{bmatrix}\|\mathbf{y}_{j}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T}\mathbf{y}_{j}\|_{2}\\\
\frac{1}{{\sqrt{n}}}\mathbf{1}^{T}\mathbf{y}_{j}\end{bmatrix}\end{split}$ (41)

We now evaluate the primal problem objective with the parameters in (41),
which yields

$\displaystyle p_{2v}^{*}$
$\displaystyle=\frac{1}{2}\left\|\sum_{j=1}^{C}\left(\frac{(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}{\mathbf{w}^{(1)}_{j}}^{*}}{\|(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{T})\mathbf{X}{\mathbf{w}^{(1)}_{j}}^{*}\|_{2}}{\gamma^{(1)}_{j}}^{*}+\frac{\mathbf{1}}{\sqrt{n}}{\alpha^{(1)}_{j}}^{*}\right)_{+}{\mathbf{w}^{(2)}_{j}}^{*^{T}}-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{j=1}^{C}\left\|{\mathbf{w}^{(2)}_{j}}^{*}\right\|_{2}$
$\displaystyle=\frac{1}{2}\left\|\sum_{j\in\mathcal{E}}\left(\|\mathbf{y}_{j}\|_{2}-\beta\right)\frac{\mathbf{y}_{j}}{\|\mathbf{y}_{j}\|_{2}}\mathbf{e}_{j}^{T}-\mathbf{Y}\right\|_{F}^{2}+\beta\sum_{j\in\mathcal{E}}(\|\mathbf{y}_{j}\|_{2}-\beta)$
$\displaystyle=\frac{1}{2}\sum_{j\in\mathcal{E}}\left\|\beta\frac{\mathbf{y}_{j}}{\|\mathbf{y}_{j}\|_{2}}\mathbf{e}_{j}^{T}\
\right\|_{F}^{2}+\frac{1}{2}\sum_{j\notin\mathcal{E}}\|\mathbf{y}_{j}\mathbf{e}_{j}^{T}\|_{F}^{2}+\beta\sum_{j\in\mathcal{E}}\|\mathbf{y}_{j}\|_{2}-\beta^{2}|\mathcal{E}|$
$\displaystyle=-\frac{1}{2}\beta^{2}|\mathcal{E}|+\frac{1}{2}\sum_{j\notin\mathcal{E}}\|\mathbf{y}_{j}\|_{2}^{2}+\beta\sum_{j\in\mathcal{E}}\|\mathbf{y}_{j}\|_{2},$
(42)

which is the same with (C.6). Therefore, strong duality holds, i.e.,
$p_{2v}^{*}=d_{2v}^{*}$, and the layer weights in (41) are optimal for the
primal problem (40). ∎

###  Appendix D Convolutional Neural Networks

####  D.1 Proof of Theorem 3.1

Following the arguments in Theorem 2.2 and Ergen & Pilanci (2021a), strong
duality holds for the CNN training problem in (13), i.e.,
$p_{2c}^{*}=d_{2c}^{*}$, and the dual constraint can be equivalently written
as follows

$\displaystyle\max_{\begin{subarray}{c}\mathbf{w}^{(1)}\\\
{\gamma^{(1)}}^{2}+{\alpha^{(1)}}^{2}=1\end{subarray}}\left|\sum_{k=1}^{K}\mathbf{v}^{T}\left(\frac{\mathbf{X}_{k}\mathbf{z}-\mu\mathbf{1}}{\sqrt{\sum_{k^{\prime}=1}^{K}\|\mathbf{X}_{k^{\prime}}\mathbf{z}-\mu\mathbf{1}\|_{2}^{2}}}\gamma^{(1)}+\frac{1}{\sqrt{nK}}\mathbf{1}\alpha^{(1)}\right)_{+}\right|$
$\displaystyle=\max_{\theta\in\Theta_{s}}\left|\mathbf{v}^{T}\mathbf{C}\left(\frac{(\mathbf{I}_{nK}-\frac{1}{nK}\mathbf{1}\mathbf{1}^{T})\mathbf{M}\mathbf{z}}{\|(\mathbf{I}_{nK}-\frac{1}{nK}\mathbf{1}\mathbf{1}^{T})\mathbf{M}\mathbf{z}\|_{2}}\gamma^{(1)}+\frac{\mathbf{1}}{\sqrt{nK}}\alpha^{(1)}\right)_{+}\right|$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{q}\right\|_{2}=\left\|\mathbf{s}\right\|_{2}=1\end{subarray}}\left|\mathbf{v}^{T}\mathbf{C}\left(\begin{bmatrix}\mathbf{U}_{M}\mathbf{q}&\frac{\mathbf{1}}{\sqrt{nK}}\end{bmatrix}\mathbf{s}\right)_{+}\right|$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\left|\mathbf{v}^{T}\mathbf{C}\left(\mathbf{U}_{M}^{\prime}\mathbf{s}^{\prime}\right)_{+}\right|$
$\displaystyle=\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\end{subarray}}\left|\sum_{k=1}^{K}\mathbf{v}^{T}\left(\mathbf{U}_{Mk}^{\prime}\mathbf{s}^{\prime}\right)_{+}\right|$
$\displaystyle=\max_{i\in[P]}\max_{\begin{subarray}{c}\left\|\mathbf{s}^{\prime}\right\|_{2}=1\\\
(2\mathbf{D}_{ik}-\mathbf{I}_{n})\mathbf{U}_{Mk}^{\prime}\mathbf{s}^{\prime}\geq
0\end{subarray}}\left|\sum_{k=1}^{K}\mathbf{v}^{T}\mathbf{D}_{ik}\mathbf{U}_{Mk}^{\prime}\mathbf{s}^{\prime}\right|,$

where $\mathbf{C}:=\mathbf{1}^{T}_{K}\otimes\mathbf{I}_{n}$,
$\mathbf{M}=[\mathbf{X}_{1};\,\mathbf{X}_{2};\,\ldots\,\mathbf{X}_{K}]$, and
we use the following variable changes
$\mathbf{q}:=\bm{\Sigma}_{M}\mathbf{V}_{M}^{T}\mathbf{z}$,
$\mathbf{s}:=\begin{bmatrix}\gamma^{(1)}&\alpha^{(1)}\end{bmatrix}^{T}$,
$\mathbf{s}^{\prime}:=\begin{bmatrix}\gamma^{(1)}\mathbf{q}^{T}&\alpha^{(1)}\end{bmatrix}^{T}$,
and
$\mathbf{U}_{M}^{\prime}:=\begin{bmatrix}\mathbf{U}_{M}&\frac{\mathbf{1}}{\sqrt{nK}}\end{bmatrix}=\begin{bmatrix}\mathbf{U}_{M1}^{\prime};\mathbf{U}_{M2}^{\prime};\ldots&\mathbf{U}_{MK}^{\prime}\end{bmatrix}$.

Then, following the steps in Theorem 2.2 and Ergen & Pilanci (2021a), the
equivalent convex program for (13) is

$\displaystyle\min_{\mathbf{q}_{i},\mathbf{q}_{i}^{\prime}}\frac{1}{2}\left\|\sum_{i=1}^{P_{c}}\sum_{k=1}^{K}\mathbf{D}_{ik}\mathbf{U}_{Mk}^{\prime}(\mathbf{q}_{i}-\mathbf{q}_{i}^{\prime})\right\|+\beta\sum_{i=1}^{P_{c}}(\left\|\mathbf{q}_{i}\right\|_{2}+\left\|\mathbf{q}_{i}^{\prime}\right\|_{2})$
$\displaystyle\text{ s.t.
}(2\mathbf{D}_{ik}-\mathbf{I}_{n})\mathbf{U}_{Mk}^{\prime}\mathbf{q}_{i}\geq
0,\;(2\mathbf{D}_{ik}-\mathbf{I}_{n})\mathbf{U}_{Mk}^{\prime}\mathbf{q}_{i}^{\prime}\geq
0,\;\forall i\in[P_{c}],\forall k\in[K].$

∎

###  Appendix E The Implicit Regularization Effect of SGD on BN Networks

####  E.1 Proof of Lemma 5.1

We start with the following non-convex training problem

$\displaystyle
p_{v2}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|\sum_{j=1}^{m_{1}}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{X}\mathbf{w}^{(1)}_{j}\right)\right)_{+}{\mathbf{w}^{(2)}_{j}}^{T}-\mathbf{Y}\right\|_{F}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}+\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}^{2}\right),$

where we note that hidden layer weights are not regularized without loss of
generality due to the arguments in Appendix B.2.1. Now, consider the SVD of
the zero mean form of the data matrix
$(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top})\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$.
We can thus equivalently write the problem as

$\displaystyle
p_{v2}^{*}=\min_{\theta\in\Theta}\frac{1}{2}\left\|\sum_{j=1}^{m_{1}}\left(\frac{\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}\mathbf{w}^{(1)}_{j}}{\|\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}\mathbf{w}^{(1)}_{j}\|_{2}}\gamma^{(1)}_{j}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}_{j}\right)_{+}{\mathbf{w}^{(2)}_{j}}^{T}-\mathbf{Y}\right\|_{F}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}+\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}^{2}\right),$

Now, simply let
$\mathbf{q}_{j}=\mathbf{\Sigma}\mathbf{V}^{\top}\mathbf{w}^{(1)}_{j}$. Noting
that $\|\mathbf{U}\mathbf{q}_{j}\|_{2}=\|\mathbf{q}_{j}\|_{2}$, we arrive at

$\displaystyle
p_{v2}^{*}=p_{v2b}^{*}:=\min_{\theta\in\Theta}\frac{1}{2}\left\|\sum_{j=1}^{m_{1}}\left(\frac{\mathbf{U}\mathbf{q}_{j}}{\|\mathbf{q}_{j}\|_{2}}\gamma^{(1)}_{j}+\frac{\mathbf{1}}{\sqrt{n}}\alpha^{(1)}_{j}\right)_{+}{\mathbf{w}^{(2)}_{j}}^{T}-\mathbf{Y}\right\|_{F}^{2}+\frac{\beta}{2}\sum_{j=1}^{m_{1}}\left({\gamma^{(1)}_{j}}^{2}+{\alpha^{(1)}_{j}}^{2}+\left\|\mathbf{w}^{(2)}_{j}\right\|_{2}^{2}\right),$

∎

####  E.2 Proof of Theorem 5.1

When solving (9) with GD, we have the updates

$\displaystyle{\mathbf{w}^{(1)}_{j}}^{(k+1)}$
$\displaystyle={\mathbf{w}^{(1)}_{j}}^{(k)}-\eta_{k}\nabla_{\mathbf{w}^{(1)}_{j}}\mathcal{L}(f_{v2}^{k}(\mathbf{X}),\mathbf{Y})$
$\displaystyle{\alpha^{(1)}_{j}}^{(k+1)}$
$\displaystyle={\alpha^{(1)}_{j}}^{(k)}-\eta_{k}\Big{(}\nabla_{\alpha^{(1)}_{j}}\mathcal{L}(f_{v2}^{k}(\mathbf{X}),\mathbf{Y})+\beta\alpha^{(1)}_{j}\Big{)}$
$\displaystyle{\gamma^{(1)}_{j}}^{(k+1)}$
$\displaystyle={\gamma^{(1)}_{j}}^{(k)}-\eta_{k}\Big{(}\nabla_{\gamma^{(1)}_{j}}\mathcal{L}(f_{v2}^{k}(\mathbf{X}),\mathbf{Y})+\beta\gamma^{(1)}_{j}\Big{)}$
$\displaystyle{\mathbf{w}^{(2)}_{j}}^{(k+1)}$
$\displaystyle={\mathbf{w}^{(2)}_{j}}^{(k)}-\eta_{k}\Big{(}\nabla_{\mathbf{w}^{(2)}_{j}}\mathcal{L}(f_{v2}^{k}(\mathbf{X}),\mathbf{Y})+\beta\mathbf{w}^{(2)}_{j}\Big{)}$

Now, let
$\mathbf{D}_{j}^{(k)}=\mathrm{diag}\\{\mathbbm{1}[\mathbf{X}{\mathbf{w}^{(1)}_{j}}^{(k)}\geq
0]\\}$,
$\mathbf{R}^{(k)}=\sum_{j=1}^{m_{1}}\left({\mathrm{BN}}_{\gamma,\alpha}\left(\mathbf{X}{\mathbf{w}^{(1)}_{j}}^{(k)}\right)\right)_{+}{{\mathbf{w}^{(2)}_{j}}^{(k)}}^{T}-\mathbf{Y}$
be the residual at step $k$, and
$\mathbf{X}^{0}=(\mathbf{I}_{n}-\frac{1}{n}\mathbf{1}\mathbf{1}^{\top})\mathbf{X}$
. Then, with squared loss $\mathcal{L}$, we have the following:

$\displaystyle\nabla_{\mathbf{w}^{(1)}_{j}}\mathcal{L}(f_{v2}^{k}(\mathbf{X}),\mathbf{Y})$
$\displaystyle=\frac{{\gamma^{(1)}_{j}}^{(k)}}{\|\mathbf{X}^{0}{\mathbf{w}^{(1)}_{j}}^{(k)}\|_{2}}(\mathbf{D}_{j}^{(k)}\mathbf{X}^{0})^{\top}\mathbf{R}^{(k)}{\mathbf{w}^{(2)}_{j}}^{(k)}$
$\displaystyle-\frac{{\gamma^{(1)}_{j}}^{(k)}}{\|\mathbf{X}^{0}{\mathbf{w}^{(1)}_{j}}^{(k)}\|_{2}^{3}}\Big{(}{{\mathbf{w}^{(1)}_{j}}^{(k)}}^{\top}(\mathbf{D}_{j}^{(k)}\mathbf{X}^{0})^{\top}\mathbf{R}^{(k)}{\mathbf{w}^{(2)}_{j}}^{(k)}\Big{)}{\mathbf{X}^{0}}^{\top}{\mathbf{X}^{0}}\boldsymbol{\mathbf{w}^{(1)}_{j}}^{(k)}$
$\displaystyle\nabla_{\alpha^{(1)}_{j}}\mathcal{L}(f_{v2}^{k}(\mathbf{X}),\mathbf{Y})$