DEMYSTIFYING BATCH NORMALIZATION IN RELU NETWORKS: EQUIVALENT CONVEX OPTIMIZATION MODELS
Year: 2021

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.

heading: 1 INTRODUCTION
text: 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.

heading: 1.1 RELATED WORK
text: 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, 2021; 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 (2021) 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;b). 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).

heading: 1.2 OUR CONTRIBUTIONS
text: • 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 ReLU1, vector output networks, and arbitrary convex loss functions.

heading: 1.3 PRELIMINARIES
text: Notation: We denote matrices and vectors as uppercase and lowercase bold letters, respectively, where a subscript indicates a certain element or column. We use 0 (or 1) to denote a vector or matrix of zeros (or ones), where the sizes are appropriately chosen depending on the context. We also use In to denote the identity matrix of size n. To represent Euclidean, Frobenius, and nuclear norms, we use ‖ · ‖2, ‖ · ‖F , and ‖ · ‖∗, respectively. Lastly, we denote the elementwise 0-1 valued indicator function and ReLU activation as 1[x ≥ 0] and (x)+ = max{x, 0}, respectively.
We2 consider an L-layer ReLU network with layer weights W(l) ∈ Rml−1×ml , where m0 = d and mL = C are the input and output dimensions, respectively. Given a training data X ∈ Rn×d and a label matrix Y ∈ Rn×C , we particularly focus on the following regularized training problem
min θ∈Θ L(fθ,L(X),Y) + βR(θ), (1)
where we compactly represent the parameters as θ := {W(l),γl,αl}Ll=1 and the corresponding parameter space as Θ := {{W(l),γ(l),α(l)}Ll=1 : W(l) ∈ Rml−1×ml ,γl ∈ Rml ,αl ∈ Rml ,∀l ∈ [L]}. We note that γl and αl are the parameters of the BN operator, for which we discuss the details below. In addition, L(·, ·) is an arbitrary convex loss function, including squared, hinge, and cross entropy loss, andR(·) is the regularization function for the layer weights with the tuning parameter β > 0. We also compactly define the network output as
fθ,L(X) := A (L−1)W(L), where A(l) := ( BNγ,α ( A(l−1)W(l) )) +
we denote the lth layer activations as A(l) ∈ Rn×ml , and A(0) = X. Here, BNγ,α (·) represents the BN operation introduced in Ioffe & Szegedy (2015) and applies matrices columnwise.
1Presented in Appendix G. 2All the proofs and some extensions are presented in Appendix.
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 W(l) ∈ Rml−1×ml and arbitrary batch of activations denoted as A
(l−1) b ∈ Rs×ml−1 , BN applies to each column j independently as follows
BNγ,α ( A (l−1) b w (l) j ) := (Is − 1s11 T )A (l−1) b w (l) j
‖(Is − 1s11T )A (l−1) b w (l) j ‖2
γ (l) j + α (l) j 1√ n , (2)
where γ(l) and α(l) are trainable parameters that scale and shift the normalized value. In this work, we focus on the full-batch case, i.e, A(l−1)b = A
(l−1) ∈ Rn×ml . 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 weightdecay 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.

heading: 1.4 OVERVIEW OF OUR RESULTS
text: Here, we provide an overview of our main results. To simplify the notation, we consider L-layer ReLU networks with scalar outputs, i.e., mL = C = 1 thus the label vector is y ∈ Rn, and extend the analysis to vector output networks with the label matrix Y ∈ Rn×C in the next sections. The regularized training problems for an L-layer network with scalar output and BN is given by
p∗L := min θ∈Θ
1 2 ‖fθ,L(X)− y‖22 + β 2 L∑ l=1 (∥∥∥γ(l)∥∥∥2 2 + ∥∥∥α(l)∥∥∥2 2 + ∥∥∥W(l)∥∥∥2 F ) , (3)
where we use γ(L) = α(L) = 0 as dummy variables for notational simplicity. Lemma 1.1. The problem in (3) is equivalent to the following optimization problem
min θ∈Θs
1 2 ‖fθ,L(X)− y‖22 + β ∥∥∥w(L)∥∥∥ 1 , (4)
where Θs := {θ ∈ Θ : γ(L−1)j 2 + α (L−1) j 2 = 1, ∀j ∈ [mL−1]}.
Using the equivalence in Lemma 1.1, we now take dual of (4) with respect to the output layer weights w(L) to obtain 3
p∗L ≥ d∗L := max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. max θ∈Θs ∣∣∣∣v> (BNγ,α (A(L−2)w(L−1))) + ∣∣∣∣ ≤ β. (5) Since the original formulation in (3) is a non-convex optimization problem, any solution v in the dual domain yields a lower bound for the primal problem, i.e., p∗L ≥ 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.

heading: 2 TWO-LAYER NETWORKS
text: In this section, we analyze two-layer networks. Particularly, we first consider a high-dimensional regime, where n ≤ 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.
3Details regarding the derivation of the dual problem are presented in Appendix B.3.
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
p∗2 = min θ∈Θs
1 2 ‖fθ,2(X)− y‖22 + β ∥∥∥w(2)∥∥∥ 1 . (6)
We then take dual of (6) with respect to the output weights w(2) to obtain the following dual problem
p∗2 ≥ d∗2 = max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. max θ∈Θs ∣∣∣∣vT (BNγ,α (Xw(1))) + ∣∣∣∣ ≤ β, (7) where Θs = {θ ∈ Θ : w(1) ∈ Rd, γ(1) 2 + α(1) 2 = 1}. Notice that here we drop the hidden neuron index j ∈ [m1] since the dual constraint scans all possible parameters in the continuous set Θs.
2.1 HIGH-DIMENSIONAL REGIME (n ≤ 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 ≤ d and X is full row-rank, then an optimal solution to (6) is(
w (1) j
∗ , w
(2) j
∗) = ( X† ( (−1)jy ) + , (−1)j ( ‖ ( (−1)jy ) + ‖2 − β ) + ) [ γ
(1) j
∗
α (1) j
∗ ] =
1
‖ ((−1)jy)+ ‖2
[ ‖ ( (−1)jy ) + − 1 n 11T ( (−1)jy ) + ‖2
1√ n 1T ( (−1)jy ) + ] for j = 1, 2. Therefore, strong duality holds, p∗2 = d ∗ 2.

heading: 2.2 EXACT CONVEX FORMULATION
text: 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 D := diag(1[Xw ≥ 0]) for an arbitrary w ∈ Rd. Therefore, the output of a ReLU activation can be equivalently written as (Xw)+ = DXw provided that DXw ≥ 0 and (In −D) Xw ≤ 0 are satisfied. We can define these two constraints more compactly as (2D− In) Xw ≥ 0. We now denote the cardinality of the set of all possible diagonal matrices as P , and obtain the following upperbound P ≤ 2r(e(n− 1)/r)r, where r := rank(X) ≤ 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 {Di}Pi=14
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 X, the non-convex training problem in (6) can be equivalently stated as the following finite-dimensional convex program
min si,s′i∈Rr+1
1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(si − s′i)− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (‖si‖2 + ‖s′i‖2) s.t. (2Di − In)U′si ≥ 0 (2Di − In)U′s′i ≥ 0 ,∀i,
(8)
where U ∈ Rn×r and U′ ∈ Rn×r+1 are computed via the compact SVD of the zero-mean data matrix, namely (In − 1n11 T )X := UΣVT and U′ := [ U 1√ n 1 ] .
Notice that (8) is a convex program with 2Pr variables and 2Pn constraints, and thus can be globally optimized via interior-point solvers with O(r6 ( n r )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 {U′si}Pi=1 and {U′s′i}Pi=1 multiplied with fixed hyperplane arrangement matrices {Di}Pi=1. 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., U′ T
U′ = Ir+1. 4We provide more details in Appendix B.4.

heading: 2.3 VECTOR OUTPUT NETWORKS
text: Here, we analyze networks with vector outputs, i.e., Y ∈ Rn×C , trained via the following problem
p∗v2 := min θ∈Θ
1 2 ‖fθ,2(X)−Y‖2F + β 2 m1∑ j=1 (∥∥∥w(1)j ∥∥∥2 2 + γ (1) j 2 + α (1) j 2 + ∥∥∥w(2)j ∥∥∥2 2 ) . (9)
Lemma 2.1. The problem in (9) is equivalent to the following optimization problem
p∗2v = min θ∈Θs
1 2 ‖fθ,2(X)−Y‖2F + β m1∑ j=1 ∥∥∥w(2)j ∥∥∥ 2 , (10)
where Θs := {θ ∈ Θ : γ(1)j 2 + α (1) j 2 = 1, ∀j ∈ [m1]}.
Using the equivalence in Lemma 2.1, we then take dual of (10) with respect to the output layer weights W(2) to obtain the following dual problem
p∗2v ≥ d∗2v := max V −1 2 ‖V −Y‖2F + 1 2 ‖Y‖2F s.t. max θ∈Θs ∥∥∥∥VT (BNγ,α (Xw(1))) + ∥∥∥∥ 2 ≤ β. (11)
Theorem 2.3. The non-convex training problem (9) can be cast as the following convex program
p∗v2 = min Si,S′i∈Ki
1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(Si − S′i)−Y ∥∥∥∥∥ 2
F
+ β P∑ i=1 (‖Si‖∗ + ‖S′i‖∗) (12)
for convex sets Ki := conv{zgT : (2Di − In)U′z ≥ 0, ‖g‖2 ≤ 1}, where conv{·} denotes the convex hull of its argument and U′ = [ U 1√ n 1 ] 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 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 Y be a one-hot encoded label matrix, and X be a full row-rank data matrix with n ≤ d, then an optimal solution to (10) admits(
w (1) j
∗ ,w
(2) j ∗) = ( X†yj , (‖yj‖2 − β)+ ej ) ,
[ γ
(1) j
∗
α (1) j
∗ ] = 1
‖yj‖2
[ ‖yj − 1n11
Tyj‖2 1√ n 1Tyj ] ∀j ∈ [C], where ej is the jth ordinary basis vector.

heading: 3 CONVOLUTIONAL NEURAL NETWORKS
text: In this section, we extend our analysis to CNNs. Thus, instead of X, we operate on the patch matrices that are subsets of columns extracted from X and we denote them as Xk ∈ Rn×h, where h is the filter size and k is the patch index. With this notation, a convolution operation between the data matrix X and an arbitrary filter z ∈ Rh can be represented as {Xkz}Kk=1. Therefore, given patch matrices {Xk}Kk=1, regularized training problem with BN can be formulated as follows
p∗2c := min θ∈Θc
1
2 ∥∥∥∥∥∥ m1∑ j=1 K∑ k=1 ( BNCγ,α (Xkzj) ) + wj − y ∥∥∥∥∥∥ 2 + β 2 m1∑ j=1 (‖zj‖22 + γ2j + α2j + w2j ), (13)
where given µj = 1nK ∑K k=1 1 TXkzj , BN is defined as
BNCγ,α (Xkzj) := Xkzj − µj1√∑K
k′=1 ‖Xk′zj − µj1‖22 γ
(1) j + α (1) j 1√ nK .
Using Lemma 1.1, the dual problem with respect to w(2) is
p∗2c ≥ d∗2c := max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. max θ∈Θs ∣∣∣∣∣vT K∑ k=1 ( BNCγ,α (Xkz) ) + ∣∣∣∣∣ ≤ β. (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
min qi,q ′ i∈R h+1
1
2 ∥∥∥∥∥ Pc∑ i=1 K∑ k=1 DikU ′ Mk(qi − q′i) ∥∥∥∥∥+ β Pc∑ i=1 (‖qi‖2 + ∥∥q′i∥∥2) s.t. (2Dik − In)U′Mkqi ≥ 0(2Dik − In)U′Mkq′i ≥ 0, ∀i, k
(15)
where we first define a new matrix M = [X1; X2; . . . XK ] ∈ RnK×h and then use the following notations: the compact SVD of the zero mean form of this matrix is (InK− 1nK11 T )M := UMΣMV T M ,
and U′M := [ UM 1√ nK ] = [U′M1; U ′ M2; . . . U ′ MK ].
Theorem 3.1 proves that (13) can be equivalently stated as a finite-dimensional convex program with 2rcPc variables (see Remark 3.1 for the definitions) and 2nPcK constraints. Since Pc 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 X. Particularly, the number of hyperplane arrangements for M, i.e., denoted as Pc, is upperbounded as Pc ≤ 2rc(e(n − 1)/rc)rc , where rc := rank(M) ≤ h d. Therefore, given a fixed filter size h, Pc is polynomial in n and d even when X is full rank (see Ergen & Pilanci (2021) for more details).

heading: 4 DEEP NETWORKS
text: 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 A(L−2) is Rn. Then an optimal solution for the last two layer weights in (4) admits(
w (L−1) j ∗ , w (L) j ∗) = ( A(L−2) † ( (−1)jy ) + , (−1)j ( ‖ ( (−1)jy ) + ‖2 − β ) + ) [ γ
(L−1) j
∗
α (L−1) j
∗ ] =
1
‖ ((−1)jy)+ ‖2
[ ‖ ( (−1)jy ) + − 1 n 11T ( (−1)jy ) + ‖2
1√ n 1T ( (−1)jy ) + ] for j = 1, 2. Thus, strong duality holds, i.e., p∗L = d ∗ L.

heading: 4.1 VECTOR OUTPUT NETWORKS
text: 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 Y be a one-hot encoded label matrix and the network be overparameterized such that the range of A(L−2) is Rn, then an optimal solution to (3) can be found in closed-form as(
w (L−1) j
∗ ,w
(L) j
∗) = ( A(L−2)
† yj , (‖yj‖2 − β)+ ej
) ,
[ γ
(L−1) j
∗
α (L−1) j
∗ ] = 1
‖yj‖2
[ ‖yj − 1n11
Tyj‖2 1√ n 1Tyj ] ∀j ∈ [C], where ej is the jth ordinary basis vector.

heading: 5 IMPLICIT REGULARIZATION OF GD FOR BN
text: 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.

heading: 5.1 EVIDENCE OF IMPLICIT REGULARIZATION
text: Lemma 5.1. The non-convex problem (9) is equivalent to the following problem (i.e., p∗v2 = p∗v2b)
p∗v2b := min θ∈Θ
1
2 ∥∥∥∥∥ m1∑ j=1 ( Uqj ‖qj‖2 γ (1) j + 1√ n α (1) j ) + w (2) j T −Y ∥∥∥∥∥ 2
F
+ β
2 m1∑ j=1 ( γ (1) j 2 + α (1) j 2 + ∥∥∥w(2)j ∥∥∥2 2 ) .
(16)
(16) is equivalent to a problem with new data matrix U and weight normalization (Salimans & Kingma, 2016). The identity in Lemma 5.1 is achieved by substituting the transformation qj = ΣV>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 fv2(X) and fv2b(U). At iteration k of GD, assume that fkv2(X) = f k v2b(U), i.e. both models have identical weights, besides q (k) j =
ΣV>w (1) j
(k) . Then, the subsequent GD updates to model weights are identical, besides those to
w (1) j and qj . In qj-space, the updates ∆ (k) v2,j and ∆ (k) v2b,j , respectively, admit
∆ (k) v2,j = Σ 2∆ (k) v2b,j .
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.

heading: 5.2 MAKING IMPLICIT REGULARIZATION EXPLICIT
text: Since the convex programs (8), (12), and (15) rely on the whitened data matrix 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(X) := UT (Σ)V>, (17)
with Tk(Σ)ii := σi1{i ≥ k}, which takes the k ≤ r top singular values and removes the rest. We can then solve (1) using fθ,L(g(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 U. We find this low-rank approximation to be effective, which can even improve the generalization of GD applied to BN networks.

heading: 6 NUMERICAL EXPERIMENTS
text: Here5, 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”.
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,m1, β) = (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 m1 =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,m1, β, bs, k, C) = (5x104, 3072, 4096, 10−4, 103, 215, 10) and CIFAR-100 with (n, d,m1, β, bs, k, C) = (5x104, 3072, 1024, 10−4, 103, 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,
5We present additional numerical experiments and details on the experiments in Appendix A.
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-toend 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.

heading: 8 ETHICS AND REPRODUCIBILITY STATEMENTS
text: This paper aims to provide a solid theoretical explanation for a common heuristic approach which improves the training and generalization of deep neural networks. Therefore, we believe that there aren’t any ethical concerns regarding our paper.
For the sake of reproducibility, we provide all the experimental details (including preprocessing, hyperparameter optimization, extensive ablation studies, hardware requirements, and all other implementation details) in Appendix A as well as the source (attached as a supplementary material) to reproduce the experiments in the paper. Similarly, all the proofs and explanations regarding our theoretical analysis can be found in Appendix B, C, D, E, F, G.

heading: Appendix

heading: Table of Contents

heading: A Additional Details and Numerical Experiments 13
text: A.1 Inference in BN models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 A.2 Additional Experiment: Batch Sizes for SGD . . . . . . . . . . . . . . . . . . . 14 A.3 Additional Experiment: Singular Value Truncation Ablation . . . . . . . . . . . 14 A.4 Additional Experiment: Convolutional Architectures . . . . . . . . . . . . . . . 16 A.5 Experiments from Section 5.1–Evidence of Implicit Regularization . . . . . . . 16 A.6 Verification of Closed-Form Solutions . . . . . . . . . . . . . . . . . . . . . . . 17 A.7 Exact Convex Formulation and Implicit Regularization . . . . . . . . . . . . . . 17

heading: B Additional Theoretical Results 21
text: B.1 Extensions to Arbitrary Convex Loss Functions . . . . . . . . . . . . . . . . . . 21 B.2 Proof of Lemma 1.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 B.3 Derivation of the Dual Problem in (5) . . . . . . . . . . . . . . . . . . . . . . . 23 B.4 Hyperplane Arrangements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

heading: C Two-Layer Networks 24
text: C.1 Closed-Form Solution to the Dual Problem in (7) . . . . . . . . . . . . . . . . . 24 C.2 Proof of Theorem 2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 C.3 Proof of Theorem 2.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 C.4 Proof of Lemma 2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 C.5 Proof of Theorem 2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 C.6 Proof of Theorem 2.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

heading: D Convolutional Neural Networks 31
text: D.1 Proof of Theorem 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

heading: E The Implicit Regularization Effect of SGD on BN Networks 31
text: E.1 Proof of Lemma 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 E.2 Proof of Theorem 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

heading: F Deep Networks 33
text: F.1 Proof of Theorem 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 F.2 Proof of Theorem 4.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

heading: G BN after ReLU 34
text: G.1 Proof of Theorem G.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 G.2 Proof of Theorem G.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

heading: A ADDITIONAL DETAILS AND NUMERICAL EXPERIMENTS
text: 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
min si,s′i
1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(si − s′i)− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (‖si‖2 + ‖s′i‖2) + ρ P∑ i=1 1T ( (−(2Di − In)U′si)+ + (−(2Di − In)U ′s′i)+ )
where ρ>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 γ(l)j = 1 and bias α (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(X), for test data X̂, we compute the standard forward pass fθ,L(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.

heading: A.2 ADDITIONAL EXPERIMENT: BATCH SIZES FOR SGD
text: 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× 3 and padding of 1.
We take the first two classes from CIFAR-100 to perform a binary classification task with X∈R1000×3072 and y∈R1000. For (lr, β) = (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.
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.

heading: A.3 ADDITIONAL EXPERIMENT: SINGULAR VALUE TRUNCATION ABLATION
text: 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 X∈R900×856, y∈R900×9.
We evaluated the proposed approaches with (m1, β, bs) = (10000, 10−2, 900), comparing the baseline non-convex SGD approach to the convex truncated formulation for k∈{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.

heading: A.4 ADDITIONAL EXPERIMENT: CONVOLUTIONAL ARCHITECTURES
text: 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 × 3 kernel with padding of 2 and stride of 3, with average pooling to a spatial dimension of 4 × 4. The parameters chosen for this model are given by (n,m1, β, bs, k, C, lr) = (5×104, 1024, 5×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 3a). This phenomenon is essentially due the distribution of the singular values of the data matrix. Particularly, for fully connected layers, we have U∈Rn×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 (σmax/σmin = 5892.5). On the contrary, for CNNs, we operate on the patch matrices so that Umk ∈Rn×h, where h= 27 is the size of each convolutional filter. Since the data matrix is much more well conditioned (σmax/σmin = 271.2) and has a significantly smaller feature dimension, the impact of truncation is less emphasized in our experiment.

heading: A.5 EXPERIMENTS FROM SECTION 5.1–EVIDENCE OF IMPLICIT REGULARIZATION
text: As described in the main paper, for both the standard BN architecture and whitened architecture, we use the parameters (n, d,m1, β, bs, 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 (W,W′), we define the distance measure
dX,i(W,W ′) :=
(v>i W)(v > i W ′)>
‖v>i W‖2‖v>i W′‖2 (18)
where vi are the columns of V from the SVD of the zero-mean data (In − 1n11 >)X. This distance metric thus computes the cosine similarity between W and W′ after being projected to the right singular subspace vi.
For the same networks trained in Figure 2 (in the main paper), we compute {dX,i(W(1) (k) ,W(1) (0) )}di=1 for the BN model (9), and {dU,i(Q(k),Q(0))}di=1 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 threeclass 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.

heading: A.7 EXACT CONVEX FORMULATION AND IMPLICIT REGULARIZATION
text: 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, ConvexTruncated) 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× 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 × 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 × 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× 10−5, 10−4, 10−4, and 5× 10−4, respectively.
For our CIFAR-10 experiments in this section, we always use (m1, β, bs) = (4096, 10−4, 103), while for CIFAR-100, we used (m1, β, bs) = (1024, 10−4, 103), 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.

heading: A.7.1 TRAINING CURVES AND OVERFITTING
text: 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
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 bs∈{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 bs = 1000.
A.7.3 IMPACT OF DIFFERENT LEARNING RATES
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
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 (β1, β2, ) = (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× 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.

heading: A.7.5 COMPARISON OF BN TO NON-BN TWO-LAYER ARCHITECTURES
text: 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,m1, β, 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× 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× 10−10 on CIFAR-10 and 5× 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.

heading: B ADDITIONAL THEORETICAL RESULTS
text: In this section, we present additional theoretical results that are not included in the main paper due to the page limit.

heading: B.1 EXTENSIONS TO ARBITRARY CONVEX LOSS FUNCTIONS
text: 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
min θ∈Θs
L(fθ,L(X),y) + β ∥∥∥w(L)∥∥∥
1 (19)
where L(·,y) is a convex loss function. Then, the dual of (19) is given by
max v −L∗(v) s.t. max
θ∈Θs
∣∣∣vTBNγ,α (A(L−2)w(1))∣∣∣≤β,
where L∗ is the Fenchel conjugate function defined as Boyd & Vandenberghe (2004)
L∗(v) = max z zTv − L(z,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.

heading: B.2 PROOF OF LEMMA 1.1
text: 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 `1 norm regularization penalty on the output layer weights wL.

heading: B.2.1 EFFECTIVE REGULARIZATION IN (3)
text: Here, we prove our claim about the effective regularization. Let us restate the training problem as follows
p∗Lr := min θ∈Θ
1
2 ∥∥∥∥∥∥ mL−1∑ j=1 ( BNγ,α ( A(L−2)w (L−1) j )) + w (L) j − y ∥∥∥∥∥∥ 2
2
+ β
2 L∑ l=1 (∥∥∥γ(l)∥∥∥2 2 + ∥∥∥α(l)∥∥∥2 2 + ∥∥∥W(l)∥∥∥2 F ) ,
(20)
where we use γ(L) =α(L) =0 as dummy variables for notational simplicity. Now, let us rewrite (20) as
p∗Lr = min t≥0 min θ∈Θ
1
2 ∥∥∥∥∥∥ mL−1∑ j=1 ( BNγ,α ( A(L−2)w (L−1) j )) + w (L) j − y ∥∥∥∥∥∥ 2
2
+ β
2 mL−1∑ j=1 ( γ (L−1) j 2 + α (L−1) j 2 + w (L) j 2) + β 2 t
s.t. L−2∑ l=1 (∥∥∥γ(l)∥∥∥2 2 + ∥∥∥α(l)∥∥∥2 2 + ∥∥∥W(l)∥∥∥2 F ) + ∥∥∥W(L−1)∥∥∥2 F ≤ t.
After applying the scaling in Lemma 1.1, we then take the dual with respect to w(L) as in the main paper to obtain the following problem
p∗Lr≥d∗Lr := max t≥0 max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 + β 2 t (21)
s.t. max θ∈Θsr ∣∣∣∣vT (BNγ,α (A(L−2)w(L−1)))+ ∣∣∣∣≤β.
where Θsr :={θ∈Θ :γ(L−1)j 2 +α (L−1) j 2 = 1, ∀j∈ [mL−1], ∑L−2 l=1 (∥∥γ(l)∥∥2 2 + ∥∥α(l)∥∥2 2 + ∥∥W(l)∥∥2 F ) +∥∥W(L−1)∥∥2
F ≤ t}. Since
BNγ,α
( A(L−2)w
(L−1) j
) = (In − 1n11 T )A(L−2)w (L−1) j
‖(In − 1n11T )A(L−2)w (L−1) j ‖2︸ ︷︷ ︸
h(θ′)
γ (L−1) j + 1√ n α (L−1) j ,
where θ′ denotes all the parameters except γ(L−1),α(L−1),w(L). Then, independent of the value t, h(θ′) is always a unit norm vector. Therefore, the maximization constraint in (21) is independent of the norms of the parameters in θ′. Consequently, regularizing the weights in θ′ does not affect the dual characterization in (21).
Based on this observation, (20) can be equivalently stated as
p∗L := min θ∈Θ
1 2 ‖fθ,L(X)− y‖22 + β 2 mL−1∑ j=1 ( γ (L−1) j 2 + α (L−1) j 2 + w (L) j 2) (22)

heading: B.2.2 RESCALING FOR `1 NORM REGULARIZATION IN (4)
text: We first note that similar approaches are also presented in Neyshabur et al. (2014); Savarese et al. (2019); Pilanci & Ergen (2020); Ergen & Pilanci (2020).
For any θ∈Θ, we can rescale the parameters as γ̄(L−1)j =ηjγ (L−1) j , ᾱ (L−1) j =ηjα (L−1) j and w̄ (L) j = w (L) j /ηj , for any ηj>0. Then, the network output becomes
fθ̄,L(X) = mL−1∑ j=1
( (In − 1n11 T )A(L−2)w (L−1) j
‖(In − 1n11T )A(L−2)w (L−1) j ‖2
γ̄ (L−1) j + 1ᾱ (L−1) j ) + w̄(L)
= mL−1∑ j=1
( (In − 1n11 T )A(L−2)w (L−1) j
‖(In − 1n11T )A(L−2)w (L−1) j ‖2
γ (L−1) j + 1α (L−1) j ) + w(L)
=fθ,L(X),
which proves fθ,2(X) =fθ̄,2(X). Moreover, we use the following inequality
1
2 mL−1∑ j=1 ( γ (L−1) j 2 + α (L−1) j 2 + w (L) j 2) ≥ mL−1∑ j=1 ∣∣∣w(L)j ∣∣∣ √ γ (L−1) j 2 + α (L−1) j 2 ,
where the equality is achieved when
ηj =  |w(L)j |√ γ
(L−1) j 2 + α (L−1) j 2  12 is used. Since this scaling does not alter the right-hand side of the inequality, without loss of generality,
we can set √ γ
(L−1) j 2 + α (L−1) j 2 = 1,∀j∈ [mL−1]. Therefore, the right-hand side becomes ‖w(L)‖1
and (22) is equivalent to
min θ∈Θs
1 2 ‖fθ,L(X)− y‖22 + β‖w (L)‖1
where Θs :={θ∈Θ :γ(L−1)j 2 + α (L−1) j 2 = 1, ∀j∈ [mL−1]}.

heading: B.3 DERIVATION OF THE DUAL PROBLEM IN (5)
text: Let us start with reparameterizing the scaled primal problem in (4) as
p∗L= min θ∈Θs,ŷ∈Rn
1 2 ‖ŷ − y‖22 + β ∥∥∥w(L)∥∥∥ 1 s.t. ŷ = mL−1∑ j=1 ( BNγ,α ( A(L−2)w (L−1) j )) + w (L) j .
(23)
where Θs :={θ∈Θ :γ(L−1)j 2 +α (L−1) j 2 = 1, ∀j∈ [mL−1]}. We now form the following Lagrangian
L(v, ŷ,w(L)) = 1
2 ‖ŷ − y‖22 − v
T ŷ + β ∥∥∥w(L)∥∥∥
1 + vT mL−1∑ j=1 ( BNγ,α ( A(L−2)w (L−1) j )) + w (L) j .
Then, the corresponding dual function is
g(v) = min ŷ,w(L)
L(v, ŷ,w(L))
=−1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. ∣∣∣∣v> (BNγ,α (A(L−2)w(L−1)j )) + ∣∣∣∣≤β,∀j∈ [mL−1]. Hence, the dual of (23) with respect to ŷ and w(L) is
p∗L= min θ∈Θs max v g(v) = min θ∈Θs max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. ∣∣∣∣v> (BNγ,α (A(L−2)w(L−1)j )) + ∣∣∣∣≤β,∀j∈ [mL−1].
We then change the order of min-max to obtain the following lower bound
p∗L≥d∗L := max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. max θ∈Θs ∣∣∣∣v> (BNγ,α (A(L−2)w(L−1))) + ∣∣∣∣≤β.

heading: B.4 HYPERPLANE ARRANGEMENTS
text: We first define the set of all hyperplane arrangements for the data matrix X as H := ⋃{ {sign(Xw)} : w∈Rd } ,
where |H|≤2n. We now define a new set to denote the indices with positive signs for each element in the set H as S := { {∪hi=1{i}} : h∈H } . With this definition, we note that given an element S∈S, one can introduce a diagonal matrix D(S)∈Rn×n defined as D(S)ii :=1[i∈S]. Therefore, the output of ReLU activation can be equivalently written as (Xw)+ =D(S)Xw provided that D(S)Xw≥0 and (In −D(S)) Xw≤0 are satisfied. One can define more compactly these two constraints as (2D(S)− In) Xw≥0. We now denote the cardinality of S as P , and obtain the following upperbound
P ≤2 r−1∑ k=0 ( n− 1 k ) ≤2r ( e(n− 1) r )r where r := rank(X)≤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 {Di}Pi=1.

heading: C TWO-LAYER NETWORKS

heading: C.1 CLOSED-FORM SOLUTION TO THE DUAL PROBLEM IN (7)
text: Lemma C.1. Suppose n≤d and X is full row-rank, then an optimal solution to (7) is
v∗=  β (y)+ ‖(y)+‖2 − β (−y)+‖(−y)+‖2 if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 (y)+ − β (−y)+ ‖(−y)+‖2 if β>‖ (y)+ ‖2, β≤ (−y)+ ‖2 β (y)+ ‖(y)+‖2
− (−y)+ if β≤‖ (y)+ ‖2, β >‖ (−y)+ ‖2 y if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2
.
Proof. Since n≤d and X is full row-rank, we have
max θ∈Θs ∣∣∣∣∣vT ( (In − 1n11 T )Xw(1) ‖(In − 1n11T )Xw(1)‖2 γ(1) + 1√ n α(1) ) + ∣∣∣∣∣≤max{‖ (v)+ ‖2, ‖ (−v)+ ‖2} max1T u=0 ‖u‖2=‖s‖2=1 ∥∥∥∥([u 1√n] s) + ∥∥∥∥ 2
≤max { ‖ (v)+ ‖2, ‖ (−v)+ ‖2 } max ‖u‖2=1 ∥∥(u)+∥∥2 ≤max { ‖ (v)+ ‖2, ‖ (−v)+ ‖2 } ,
where s= [ γ(1) α(1) ]T and the equality is achieved when
w(1) =X†(v −min i vi), s=
1
‖v‖2
[ ‖v − 1n11
Tv‖2 1√ n 1Tv
] .
Thus, the dual problem (7) can be reformulated as
max v −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. max
{ ‖ (v)+ ‖2, ‖ (−v)+ ‖2 } ≤β.
The problem above can be maximized using the following optimal dual parameter
v∗=  β (y)+ ‖(y)+‖2 − β (−y)+‖(−y)+‖2 if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 (y)+ − β (−y)+ ‖(−y)+‖2 if β>‖ (y)+ ‖2, β≤ (−y)+ ‖2 β (y)+ ‖(y)+‖2
− (−y)+ if β≤‖ (y)+ ‖2, β >‖ (−y)+ ‖2 y if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2
.

heading: C.2 PROOF OF THEOREM 2.1
text: For the solution in Lemma C.1, the corresponding extreme points for the two-sided constraint are
W(1) ∗ =  ( X† (y)+ ,X † (−y)+ ) if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 X† (−y)+ if β>‖ (y)+ ‖2, β≤ (−y)+ ‖2 X† (y)+ if β≤‖ (y)+ ‖2, β >‖ (y)+ ‖2
0 if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2
and[ γ
(1) 1
∗
α (1) 1
∗ ] =
1
‖ (y)+ ‖2
[ ‖ (y)+ − 1 n11
T (y)+ ‖2 1√ n 1T (y)+
] [ γ
(1) 2
∗
α (1) 2
∗ ] =
1
‖ (−y)+ ‖2
[ ‖ (−y)+ − 1 n11
T (−y)+ ‖2 1√ n 1T (−y)+
] .
We now substitute these solutions into the primal problem (6) and then take the derivative with respect to the output layer weights w(2). Since the primal problem is linear and convex provided that the first layer weights are fixed as W(1) ∗ , we obtain the following closed-form solutions for the output layer weights
w(2) ∗ =  [ ‖ (y)+ ‖2 − β −‖ (−y)+ ‖2 + β ] if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 −‖ (−y)+ ‖2 + β if β>‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 ‖ (y)+ ‖2 − β if β≤‖ (y)+ ‖2, β >‖ (−y)+ ‖2
0 if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2
.
These set of weights, {W(1)∗, γ(1)∗, α(1)∗,w(2)∗}, achieves p∗2 =d∗2, therefore, strong duality holds.

heading: C.3 PROOF OF THEOREM 2.2
text: We first note that if the number of neurons satisfies m1≥m∗ for some m∗∈N, where m∗≤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 m1 =n, and therefore the ReLU network fits any real valued y∈Rn with n neurons. Following the steps in Pilanci & Ergen (2020), we first focus on the dual constraint in (7)
max w(1)
γ(1) 2 +α(1) 2 =1
∣∣∣∣∣vT ( (In − 1n11 T )Xw(1)
‖(In − 1n11T )Xw(1)‖2 γ(1) +
1√ n α(1) ) + ∣∣∣∣∣= max‖q‖2=‖s‖2=1 ∣∣∣∣vT ([Uq 1√n] s) + ∣∣∣∣ = max ‖s′‖
2 =1 ∣∣vT (U′s′)+∣∣ = max i∈[P ] max ‖s′‖
2 =1
(2Di−In)U′s′≥0
∣∣vTDiU′s′∣∣ ,
where q :=ΣVTw(1)/ ∥∥ΣVTw(1)∥∥
2 , s :=
[ γ(1) α(1) ]T , s′ := [ γ(1)qT α(1) ]T , and U′ :=[
U 1√ n
] . Now, let us consider the single side of this constraint and introduce a dual variable
ρ for the constraint as follows
max i∈[P ]
max ‖s′‖
2 =1
(2Di−In)U′s′≥0
vTDiU ′s′= max
i∈[P ] min ρ≥0 max ‖s′‖
2 =1
( vTDiU ′ + ρT (2Di − In)U′ ) s′
= max i∈[P ] min ρ≥0 ∥∥vTDiU′ + ρT (2Di − In)U′∥∥2 . Therefore, we have
max θ∈Θs
vT
( (In − 1n11 T )Xw(1)
‖(In − 1n11T )Xw(1)‖2 γ(1) +
1√ n α(1) ) + ≤β ⇐⇒∀i∈ [P ],min ρ≥0 ∥∥vTDiU′ + ρT (2Di − In)U′∥∥2≤β ∀i∈ [P ], ∃ρi≥0 s.t.
∥∥vTDiU′ + ρTi (2Di − In)U′∥∥2≤β. 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 v,ρi,ρ′i∈Rn,∀i∈ [P ] and 2P constraints as follows
d∗2 = maxv,ρi −1 2 ‖v − y‖22 + 1 2 ‖y‖22 s.t. ∥∥vTDiU′ + ρTi (2Di − In)U′∥∥2≤β∥∥∥−vTDiU′ + ρ′Ti (2Di − In)U′∥∥∥ 2 ≤β , ∀i∈ [P ]. (24)
As noted in Pilanci & Ergen (2020), since the problem in (24) is strictly feasible for the particular choice of parameters v =ρi=ρ′i=0,∀∈ [P ], Slater’s condition and therefore strong duality holds (Boyd & Vandenberghe, 2004). We now write (24) in a Lagrangian form with the parameters λ,λ′∈RP as
d∗2 = min λ,λ′≥0 max v ρi≥0,∀i −1 2 ‖v − y‖22 + 1 2 ‖y‖22 + P∑ i=1 λi(β − ∥∥vTDiU′ + ρTi (2Di − In)U′∥∥2)
+ P∑ i=1 λ′i(β − ∥∥∥−vTDiU′ + ρ′Ti (2Di − In)U′∥∥∥ 2 ).
(25)
We next introduce additional variables ri, r′i∈Rd,∀i∈ [P ] to represent (25) as
d∗2 = min λ,λ′≥0 max v
ρi≥0,∀i min
ri:‖ri‖2≤1 r′i ‖r ′ i‖2≤1
−1 2 ‖v − y‖22 + 1 2 ‖y‖22 + P∑ i=1 λi ( β − ( vTDiU ′ + ρTi (2Di − In)U′ ) ri )
+ P∑ i=1 λ′i ( β + ( −vTDiU′ + ρ′ T i (2Di − In)U′ ) r′i ) .
(26)
Now, we remark that the objective function in (26) is concave in v,ρi,ρ′i and convex in ri, r ′ i,∀i∈ [P ]. In addition, the constraint set for the additional variables, i.e., ‖ri‖2, ‖r′i‖2≤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
d∗2 = min λ,λ′≥0
min ri:‖ri‖2≤1 r′i ‖r ′ i‖2≤1 max v
ρi≥0,∀i −1 2 ‖v − y‖22 + 1 2 ‖y‖22 + P∑ i=1 λi ( β − ( vTDiU ′ + ρTi (2Di − In)U′ ) ri )
+ P∑ i=1 λ′i ( β + ( −vTDiU′ + ρ′ T i (2Di − In)U′ ) r′i ) .
(27)
Thus, we can now perform the maximization over the variables v,ρi,ρ′i to simplify (27) into the following form
p∗2c= min λ,λ′≥0
min ri:‖ri‖2≤1 r′i ‖r ′ i‖2≤1
1
2 ∥∥∥∥∥ P∑ i=1 (λiDiU ′ri − λ′iDiU′r′i)− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (λi + λ ′ i) (28)
s.t. (2Di − In)U′ri≥0, (2Di − In)U′r′i≥0, ∀i∈ [P ].
We then apply a set variable changes as si=λiri and s′i=λ ′ ir ′ i, ∀i∈ [P ] to obtain the following problem
p∗2c= min λ,λ′≥0
min si:‖si‖2≤λi s′i ‖s ′ i‖2≤λ ′ i
1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(si − s′i)− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (λi + λ ′ i) (29)
s.t. (2Di − In)U′si≥0, (2Di − In)U′s′i≥0, ∀i∈ [P ].
Finally, since the optimality conditions require ‖si‖2 =λi and ‖s′i‖2 =λ′i,∀i∈ [P ], we write (29) as a single convex minimization problem as
p∗2c= min si,s′i
1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(si − s′i)− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (‖si‖2 + ‖s′i‖2) (30)
s.t. (2Di − In)U′si≥0, (2Di − In)U′s′i≥0, ∀i∈ [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(
w (1) j1i
∗ , γ
(1) j1i
∗ , α
(1) j1i
∗ , w
(2) j1i
∗) = ( VΣ−1s∗i,1:r,
‖s∗i,1:r‖2 ||s∗i || 1 2 2 , s∗i,r+1 ||s∗i || 1 2 2 , ||s∗i || 1 2 2
) , if s∗i 6= 0
( w
(1) j2i
∗ , γ
(1) j2i
∗ , α
(1) j2i
∗ , w
(2) j2i
∗) = ( VΣ−1s′ ∗
i,1:r, ||s′∗i,1:r||2 ||s′∗i || 1 2 2 , s′ ∗ i,r+1 ||s′∗i || 1 2 2 ,−||s′ ∗ i || 1 2 2
) , if s′ ∗
i 6= 0,
where {s∗i , s′ ∗
i }Pi=1 are the optimal solutions to (8) and and jsi∈ [|Js|] given the definitions J1 := {i : s∗i 6= 0} and J2 :={i : s′ ∗
i 6= 0}. We also note that here the second subscript denotes the selected entries of the corresponding vector, e.g., si,1:r denotes the first r entries of the vector si. Thus, we have m∗= |J1|+ |J2|.

heading: Verifying strong duality:
text: 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 {s∗i , s′ ∗ i }Pi=1 to (30), the convex problem has the following objective
p∗2c= 1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(s∗i − s′ ∗ i )− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (‖s∗i ‖2 + ‖s′ ∗ i ‖2) (31)
and the corresponding non-convex problem objective is as follows
p∗2 = 1
2 ‖fθ,2(X)− y‖22 +
β
2 m1∑ j=1 ( γ (1) j ∗2 + α (1) j ∗2 + w (2) j ∗2 ) . (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
fθ,2(X) = m1∑ j=1  (In − 1n11T )Xw(1)j ∗ ‖(In − 1n11T )Xw (1) j ∗ ‖2 γ (1) j ∗ + 1√ n α (1) j ∗  + w (2) j ∗
= ∑ i∈J1 ( Us∗i,1:r ‖Us∗i,1:r‖2 ||s∗i,1:r||2 ||s∗i || 1 2 2 + 1√ n s∗i,r+1 ||s∗i || 1 2 2 ) + ||s∗i || 1 2 2
− ∑ i∈J2
( Us′ ∗
i,1:r ‖Us′∗i,1:r‖2 ||s′∗i,1:r||2 ||s′∗i || 1 2 2 + 1√ n s′ ∗ i,r+1 ||s′∗i || 1 2 2 ) + ||s′ ∗ i || 1 2 2
= ∑ i∈J1 (U′s∗i )+ − ∑ i∈J2 ( U′s′ ∗ i ) +
= P∑ i=1 DiU ′(s∗i − s′ ∗ i ), (33)
where the last equality follows from the constraints in (30).
Next, we prove that the regularizations terms have the same objective as follows
β
2 m1∑ j=1 ( γ (1) j ∗2 + α (1) j ∗2 + w (2) j ∗2 ) = β 2 ∑ j∈J1 ( ||s∗i,1:r||2 ||s∗i || 1 2 2 )2 + ( s∗i,r+1 ||s∗i || 1 2 2 )2 + ( ||s∗i || 1 2 2 )2 + β
2 ∑ i∈J2 ( ||s′∗i,1:r||2 ||s′∗i || 1 2 2 )2 + ( s′ ∗ i,r+1 ||s′∗i || 1 2 2 )2 + ( −||s′ ∗ i || 1 2 2 )2 = β
2 ∑ i∈J1 ( ||s∗i,1:r||22 ||s∗i ||2 + s∗ 2 i,r+1 ||s∗i ||2 + ||s∗i ||2 )
+ β
2 ∑ i∈J2 ( ||s′∗i,1:r||22 ||s′∗i ||2 + s′ ∗ i,r+1 2 ||s′∗i ||2 + ||s′ ∗ i ||2 ) =β
∑ i∈J1 ||s∗i ||2 + β ∑ i∈J2 ||s′ ∗ i ||2
=β P∑ i=1 (‖s∗i ‖2 + ‖s′ ∗ i ‖2). (34)
Therefore, by (33) and (34), the non-convex objective (32) and the convex objective (31) are the same, i.e., p∗2 =p ∗ 2c.

heading: C.4 PROOF OF LEMMA 2.1
text: We note that similar approaches are also presented in Sahiner et al. (2021a); Ergen & Pilanci (2020).
For any θ∈Θ, we can rescale the parameters as γ̄(1)j =ηjγ (1) j , ᾱ (1) j =ηjα (1) j and w̄ (2) j =w (2) j /ηj , for any ηj>0. Then, the network output becomes
fθ̄,2(X) = m1∑ j=1
( (In − 1n11 T )Xw (1) j
‖(In − 1n11T )Xw (1) j ‖2
γ̄ (1) j + 1ᾱ (1) j ) + w̄ (2)T j
= m1∑ j=1
( (In − 1n11 T )Xw (1) j
‖(In − 1n11T )Xw (1) j ‖2
γ (1) j + 1α (1) j ) + w (2) j T
=fθ,2(X),
which proves fθ,2(X) =fθ̄,2(X). Moreover, we use the following inequality
1
2 m1∑ j=1 ( γ (1) j 2 + α (1) j 2 + ∥∥∥w(2)j ∥∥∥2 2 ) ≥ m1∑ j=1 ∥∥∥w(2)j ∥∥∥ 2 √ γ (1) j 2 + α (1) j 2 ,
where the equality is achieved when
ηj =
 ∥∥∥w(2)j ∥∥∥
2√ γ
(1) j
2 + α
(1) j
2
 1 2
is used. Since this scaling does not alter the right-hand side of the inequality, without loss of generality, we can set √ γ
(1) j
2 + α
(1) j 2 = 1,∀j∈ [m1]. Therefore, the right-hand side becomes ∑m1 j=1 ‖w (2) j ‖2.
We then note that due to the arguments in Appendix B.2.1, we can remove the regularization on w(1)j without loss of generality. Therefore, the final optimization problem is as follows
p∗2v = min θ∈Θs
1 2 ‖fθ,2(X)−Y‖2F + β m1∑ j=1 ∥∥∥w(2)j ∥∥∥ 2 ,
where Θs :={θ∈Θ :γ(1)j 2 + α (1) j 2 = 1, ∀j∈ [m1]}.

heading: C.5 PROOF OF THEOREM 2.3
text: We first note that if the number of neurons satisfies m1≥m∗ for some m∗∈N, where m∗≤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)
max w(1)
γ(1) 2 +α(1) 2 =1
∥∥∥∥∥VT ( (In − 1n11 T )Xw(1)
‖(In − 1n11T )Xw(1)‖2 γ(1) +
1√ n α(1) ) + ∥∥∥∥∥ 2 = max ‖q‖2=‖s‖2=1 ∥∥∥∥VT ([Uq 1√n] s) + ∥∥∥∥ 2
= max ‖s′‖
2 =1 ∥∥VT (U′s′)+∥∥2 = max ‖g‖2=‖s′‖2=1 gTVT (U′s′)+
= max i∈[P ] max ‖g‖2=‖s′‖2=1
(2Di−In)U′s′≥0
gTVTDiU ′s′
= max i∈[P ] max ‖g‖2=‖s′‖2=1
(2Di−In)U′s′≥0
〈 V,DiU ′s′gT 〉
= max i∈[P ] max ‖S‖∗≤1 S∈Ki
〈V,DiU′S〉 ,
where s′= [ γqT α ]T , U′= [ U 1√
n ] , and Ki := conv{s′gT : (2Di − In)U′s′≥0, ‖g‖2≤1}.
Thus, we have
max θ∈Θs ∥∥∥∥∥VT ( (In − 1n11 T )Xw(1) ‖(In − 1n11T )Xw(1)‖2 γ(1) + 1√ n α(1) ) + ∥∥∥∥∥ 2 ≤β ⇐⇒ ∀i∈ [P ], ∃Si∈Ki, ‖Si‖∗≤1 s.t. 〈V,DiU′Si〉≤β.
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
min Si,S′i∈Ki
1
2 ∥∥∥∥∥ P∑ i=1 DiU ′(Si − S′i)−Y ∥∥∥∥∥ 2
F
+ β P∑ i=1 (‖Si‖∗ + ‖S′i‖∗).

heading: C.6 PROOF OF THEOREM 2.4
text: We first restate the dual problem as
max V −1 2 ‖V −Y‖2F + 1 2 ‖Y‖2F s.t. max‖s′‖
2 =1 ∥∥VT (U′s′)+∥∥2≤β, (35) where we use the results from Theorem 2.3 to modify the dual constraint. Since Y is one-hot encoded, (35) can be rewritten as
max V −1 2 ‖V −Y‖2F + 1 2 ‖Y‖2F s.t. max‖u‖2=1 ‖VTu‖2≤β. (36)
This problem has the following closed-form solution
v∗j =
{ β
yj ‖yj‖2 if β≤‖yj‖2 yj otherwise , ∀j∈ [C]. (37)
and the corresponding first layer weights in (36) are
w (1) j
∗ = { X†yj if β≤‖yj‖2 − otherwise , ∀j∈ [o]. (38)
Now let us first denote the set of indices that achieves the extreme point of the dual constraint as E :={j :β≤‖yj‖2, j∈ [C]}. Then the dual objective in (36) using the optimal dual parameter in (37)
d∗2v =− 1
2 ‖V∗ −Y‖2F +
1 2 ‖Y‖2F
=−1 2 ∑ j∈E (β − ‖yj‖2)2 + 1 2 C∑ j=1 ‖yj‖22
=−1 2 β2|E|+ β ∑ j∈E ‖yj‖2 + 1 2 ∑ j /∈E ‖yj‖22. (39)
We next restate the scaled primal problem
p∗2v = min θ∈Θs
1
2 ∥∥∥∥∥∥ m1∑ j=1 ( (In − 1n11 T )Xw (1) j ‖(In − 1n11T )Xw (1) j ‖2 γ (1) j + 1√ n α (1) j ) + w (2) j T −Y ∥∥∥∥∥∥ 2
F
+ β m1∑ j=1 ∥∥∥w(2)j ∥∥∥ 2
(40) and then obtain the optimal primal objective using the hidden layer weights in (38), which yields the following optimal solution(
w (1) j
∗ ,w
(2) j
∗) = {( X†yj , (‖yj‖2 − β) ej ) if β≤‖yj‖2
− otherwise[ γ
(1) j
∗
α (1) j
∗ ] = 1
‖yj‖2
[ ‖yj − 1n11
Tyj‖2 1√ n 1Tyj ] (41) We now evaluate the primal problem objective with the parameters in (41), which yields
p∗2v = 1
2 ∥∥∥∥∥∥ C∑ j=1  (In − 1n11T )Xw(1)j ∗ ‖(In − 1n11T )Xw (1) j ∗ ‖2 γ (1) j ∗ + 1√ n α (1) j ∗  + w (2) j ∗T −Y ∥∥∥∥∥∥ 2
F
+ β C∑ j=1 ∥∥∥w(2)j ∗∥∥∥ 2
= 1
2 ∥∥∥∥∥∥ ∑ j∈E (‖yj‖2 − β) yj ‖yj‖2 eTj −Y ∥∥∥∥∥∥ 2
F
+ β ∑ j∈E (‖yj‖2 − β)
= 1
2 ∑ j∈E ∥∥∥∥β yj‖yj‖2 eTj ∥∥∥∥2 F + 1 2 ∑ j /∈E ‖yjeTj ‖2F + β ∑ j∈E ‖yj‖2 − β2|E|
=−1 2 β2|E|+ 1 2 ∑ j /∈E ‖yj‖22 + β ∑ j∈E ‖yj‖2, (42)
which is the same with (39). Therefore, strong duality holds, i.e., p∗2v =d ∗ 2v , and the layer weights in (41) are optimal for the primal problem (40).

heading: D CONVOLUTIONAL NEURAL NETWORKS

heading: D.1 PROOF OF THEOREM 3.1
text: Following the arguments in Theorem 2.2 and Ergen & Pilanci (2021), 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
max w(1)
γ(1) 2 +α(1) 2 =1
∣∣∣∣∣∣∣ K∑ k=1 vT  Xkz− µ1√∑K k′=1 ‖Xk′z− µ1‖22 γ(1) + 1√ nK 1α(1)  + ∣∣∣∣∣∣∣ = max θ∈Θs ∣∣∣∣∣vTC ( (InK − 1nK11 T )Mz ‖(InK − 1nK11T )Mz‖2 γ(1) + 1√ nK α(1) ) +
∣∣∣∣∣ = max ‖q‖2=‖s‖2=1 ∣∣∣∣vTC([UMq 1√nK ] s)+ ∣∣∣∣
= max ‖s′‖
2 =1 ∣∣vTC (U′Ms′)+∣∣ = max ‖s′‖
2 =1 ∣∣∣∣∣ K∑ k=1 vT (U′Mks ′)+ ∣∣∣∣∣ = max i∈[P ] max ‖s′‖
2 =1
(2Dik−In)U′Mks ′≥0
∣∣∣∣∣ K∑ k=1 vTDikU ′ Mks ′ ∣∣∣∣∣ , where C :=1TK ⊗ In, M= [X1; X2; . . . XK ], and we use the following variable changes q :=ΣMV T Mz, s := [ γ(1) α(1) ]T , s′ := [ γ(1)qT α(1) ]T , and U′M := [ UM 1√ nK ] = [U′M1; U ′ M2; . . . U ′ MK ].
Then, following the steps in Theorem 2.2 and Ergen & Pilanci (2021), the equivalent convex program for (13) is
min qi,q′i
1
2 ∥∥∥∥∥ Pc∑ i=1 K∑ k=1 DikU ′ Mk(qi − q′i) ∥∥∥∥∥+ β Pc∑ i=1 (‖qi‖2 + ‖q ′ i‖2)
s.t. (2Dik − In)U′Mkqi≥0, (2Dik − In)U′Mkq′i≥0, ∀i∈ [Pc],∀k∈ [K].

heading: E THE IMPLICIT REGULARIZATION EFFECT OF SGD ON BN NETWORKS

heading: E.1 PROOF OF LEMMA 5.1
text: We start with the following non-convex training problem
p∗v2 := min θ∈Θ
1
2 ∥∥∥∥∥∥ m1∑ j=1 ( BNγ,α ( Xw (1) j )) + w (2) j T −Y ∥∥∥∥∥∥ 2
F
+ β
2 m1∑ j=1 ( γ (1) j 2 + α (1) j 2 + ∥∥∥w(2)j ∥∥∥2 2 ) ,
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
(In − 1n11 >)X=UΣV>. We can thus equivalently write the problem as
p∗v2 = min θ∈Θ
1
2 ∥∥∥∥∥∥ m1∑ j=1 ( UΣV>w (1) j ‖UΣV>w(1)j ‖2 γ (1) j + 1√ n α (1) j ) + w (2) j T −Y ∥∥∥∥∥∥ 2
F
+ β
2 m1∑ j=1 ( γ (1) j 2 + α (1) j 2 + ∥∥∥w(2)j ∥∥∥2 2 ) ,
Now, simply let qj =ΣV>w (1) j . Noting that ‖Uqj‖2 =‖qj‖2, we arrive at
p∗v2 =p ∗ v2b := min
θ∈Θ
1
2 ∥∥∥∥∥∥ m1∑ j=1 ( Uqj ‖qj‖2 γ (1) j + 1√ n α (1) j ) + w (2) j T −Y ∥∥∥∥∥∥ 2
F
+ β
2 m1∑ j=1 ( γ (1) j 2 + α (1) j 2 + ∥∥∥w(2)j ∥∥∥2 2 ) ,

heading: E.2 PROOF OF THEOREM 5.1
text: When solving (9) with GD, we have the updates
w (1) j
(k+1) =w
(1) j (k) − ηk∇w(1)j L(f k v2(X),Y)
α (1) j
(k+1) =α
(1) j (k) − ηk ( ∇ α (1) j L(fkv2(X),Y) + βα (1) j ) γ
(1) j
(k+1) =γ
(1) j (k) − ηk ( ∇ γ (1) j L(fkv2(X),Y) + βγ (1) j ) w
(2) j
(k+1) =w
(2) j (k) − ηk ( ∇
w (2) j L(fkv2(X),Y) + βw (2) j ) Now, let D(k)j = diag{1[Xw (1) j (k) ≥0]}, R(k) = ∑m1 j=1 ( BNγ,α ( Xw (1) j (k) ))
+
w (2) j
(k)T
−Y be
the residual at step k, and X0 = (In − 1n11 >)X . Then, with squared loss L, we have the following:
∇ w (1) j L(fkv2(X),Y) =
γ (1) j
(k)
‖X0w(1)j (k) ‖2
(D (k) j X 0)>R(k)w (2) j
(k)
− γ
(1) j
(k)
‖X0w(1)j (k) ‖32
( w
(1) j
(k)> (D
(k) j X 0)>R(k)w (2) j
(k)) X0 >
X0w (1) j
(k)
∇ α (1) j L(fkv2(X),Y) =w (2) j
(k)T
R(k) >
D (k) j 1
∇ γ (1) j L(fkv2(X),Y) =
1
‖X0w(1)j (k) ‖2
w (2) j
(k)T
R(k) >
D (k) j X 0w (1) j
(k)
∇ w (2) j L(fkv2(X),Y) =R(k)
> D
(k) j ( X0w(1)j (k) ‖X0w(1)j (k) ‖2 γ (1) j (k) + 1√ n α (1) j ) To express the updates in qj-space for w (1) j , we can compute
∆ (k) v2,j =ΣV >∇ w (1) j L(fkv2(X),Y)
= γ
(1) j
(k)
‖X0w(1)j (k) ‖2
Σ2U>D (k) j R (k)w (2) j
(k) −
γ (1) j
(k)
‖X0w(1)j (k) ‖32
( w
(1) j
(k)> X0 >
D (k) j R (k)w (2) j
(k)) Σ2U>X0w
(1) j
(k)
Now, we consider solving (16) with GD. We assume that all weights are equivalent to those in (9), i.e.
fkv2b(X) =f k v2(X), and q (k) j =ΣV >w (1) j
(k) . In this case, the gradients are given by:
∇qjL(fkv2b(X),Y) = γ
(1) j
(k)
‖qj(k)‖2 (D
(k) j U) >R(k)w (2) j
(k) −
γ (1) j
(k)
‖qj(k)‖32
( qj (k)>(D (k) j U) >R(k)w (2) j (k)) qj (k)
∇ α (1) j L(fkv2b(X),Y) =w (2) j
(k)T
R(k) >
D (k) j 1
∇ γ (1) j L(fkv2b(X),Y) =
1
‖qj(k)‖2 w
(2) j
(k)T
R(k) >
D (k) j Uqj (k)
∇ w (2) j L(fkv2b(X),Y) =R(k)
> D
(k) j ( Uqj(k) ‖qj(k)‖2 γ (1) j (k) + 1√ n α (1) j ) Simplifying, we have
∆ (k) v2b,j =
γ (1) j
(k)
‖qj(k)‖2 U>D (k) j R (k)w (2) j
(k) −
γ (1) j
(k)
‖qj(k)‖32
( qj (k)>U>D (k) j R (k)w (2) j (k)) qj (k)
Now, noting that X0w(1)j (k) =Uqj (k), and thus that ‖X0w(1)j (k) ‖2 =‖qj(k)‖2, we can find that, as desired,
∆ (k) v2b,j =Σ 2∆ (k) v2,j
∇ γ (1) j L(fkv2b(X),Y) =∇γ(1)j L(f k v2(X),Y)
∇ α (1) j L(fkv2b(X),Y) =∇α(1)j L(f k v2(X),Y)
∇ w (2) j L(fkv2b(X),Y) =∇w(2)j L(f k v2(X),Y)

heading: F DEEP NETWORKS

heading: F.1 PROOF OF THEOREM 4.1
text: Applying the same steps in Theorem 2.1 for the data matrix A(L−2), i.e., replacing X with A(L−2), yields the following optimal parameters
W(1) ∗ =  ( A(L−2) † (y)+ ,A (L−2)† (−y)+ ) if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 A(L−2) † (−y)+ if β>‖ (y)+ ‖2, β≤ (−y)+ ‖2 A(L−2) †
(y)+ if β≤‖ (y)+ ‖2, β >‖ (y)+ ‖2 0 if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2[
γ (1) 1
∗
α (1) 1
∗ ] =
1
‖ (y)+ ‖2
[ ‖ (y)+ − 1 n11
T (y)+ ‖2 1√ n 1T (y)+
] [ γ
(1) 2
∗
α (1) 2
∗ ] =
1
‖ (−y)+ ‖2
[ ‖ (−y)+ − 1 n11
T (−y)+ ‖2 1√ n 1T (−y)+
]
w(2) ∗ =  [ ‖ (y)+ ‖2 − β −‖ (−y)+ ‖2 + β ] if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 −‖ (−y)+ ‖2 + β if β>‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 ‖ (y)+ ‖2 − β if β≤‖ (y)+ ‖2, β >‖ (−y)+ ‖2
0 if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2
v∗=  β (y)+ ‖(y)+‖2 − β (−y)+‖(−y)+‖2 if β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 (y)+ − β (−y)+ ‖(−y)+‖2 if β>‖ (y)+ ‖2, β≤ (−y)+ ‖2 β (y)+ ‖(y)+‖2
− (−y)+ if β≤‖ (y)+ ‖2, β >‖ (−y)+ ‖2 y if β>‖ (y)+ ‖2, β >‖ (−y)+ ‖2
.
We now substitute these parameters into the primal (3) and the dual (5) to prove that strong duality holds. For presentation simplicity, we only consider the first case, i.e., β≤‖ (y)+ ‖2, β≤‖ (−y)+ ‖2 and note that the same derivations apply to the other cases as well. We first evaluate the primal objective with these parameters as follows
p∗L= 1
2 ∥∥∥∥∥∥ 2∑ j=1  (In − 1n11T )A(L−2)w(L−1)j ∗ ‖(Is − 1n11T )A(L−2)w (L−1) j ∗ ‖2 γ (l) j + 1√ n α (l) j ∗  + w (L) j ∗ − y ∥∥∥∥∥∥ 2
2
+ β ∥∥∥w(L)∗∥∥∥
1
= 1
2 ∥∥∥∥∥−β (y)+∥∥(y)+∥∥2 + β (−y)+∥∥(−y)+∥∥2 ∥∥∥∥∥ 2
2
+ β( ∥∥(y)+∥∥2 + ∥∥(−y)+∥∥2 − 2β)
=β( ∥∥(y)+∥∥2 + ∥∥(−y)+∥∥2)− β2. (43)
We then evaluate the dual objective as follows
d∗L=− 1
2 ‖v∗ − y‖22 +
1 2 ‖y‖22
=− 1 2 ∥∥∥∥β (y)+‖ (y)+ ‖2 − β (−y)+‖ (−y)+ ‖2 − y ∥∥∥∥2
2
+ 1
2 ‖y‖22
=−β2 + β( ∥∥(y)+∥∥2 + ∥∥(−y)+∥∥2). (44)
By (43) and (44), the set of parameters {W(L−1)∗, γ(L−1)∗, α(L−1)∗,w(L)∗} achieves p∗L=d∗L, therefore, strong duality holds.

heading: F.2 PROOF OF THEOREM 4.2
text: Applying the same steps in Theorem 2.4 for the data matrix A(L−2), i.e., replacing X with A(L−2), yields the following set of parameters(
w (L−1) j
∗ ,w
(L) j
∗) =
{( A(L−2) † yj , (‖yj‖2 − β) ej ) if β≤‖yj‖2
− otherwise[ γ
(L−1) j
∗
α (L−1) j
∗ ] = 1
‖yj‖2
[ ‖yj − 1n11
Tyj‖2 1√ n 1Tyj
]
v∗j =
{ β
yj ‖yj‖2 if β≤‖yj‖2 yj otherwise
, ∀j∈ [C].

heading: G BN AFTER RELU
text: Despite the common practice for BN, i.e., placing BN layers before activation functions, some recent studies show that applying BN after activation function might be more effective (Chen et al., 2019). Thus, we analyze the effects of BN when it is employed after ReLU activations.
After applying the scaling in Lemma 1.1, the primal training problem is as follows
min θ∈Θs
1
2 ∥∥∥∥∥∥ m1∑ j=1 BNγ,α (( Xw (1) j ) + ) w (2) j − y ∥∥∥∥∥∥ 2
2
+ β‖w(2)‖1. (45)
In the sequel, we repeat the same analysis in Section 2. Theorem G.1. For any data matrix X, the non-convex training problem in (45) can be equivalently stated as the following finite-dimensional convex program
min si,s′i
1
2 ∥∥∥∥∥ P∑ i=1 U′i(si − s′i)− y ∥∥∥∥∥ 2
2
+ β P∑ i=1 (‖si‖2 + ‖s′i‖2) s.t. (2Di − In)XViΣ−1i si,1:ri−1≥0 (2Di − In)XViΣ−1i s ′ i,1:ri−1≥0 ,∀i
(46)
provided that (In − 1n11 T )DiX :=UiΣiV T i , U ′ i :=
[ Ui 1√ n 1 ] , ri := rank((In − 1n11
T )DiX), and si,1:ri−1 denotes all the entries of si except the last one.
Theorem G.2. Suppose n≤d and X is full row-rank, then an optimal solution to (45) is given as( w(1) ∗ , w(2) ∗) = ( X†(y −min
i yi), (‖y‖2 − β)+
) , [ γ(1) ∗
α(1) ∗
] = 1
‖y‖2
[ ‖y − 1n11
Ty‖2 1√ n 1Ty
] .
These results indicate that applying BN after ReLU leads to a completely different phenomena. Specifically, Theorem G.1 shows that the hyperplane arrangement matrices DiX are de-meaned and whitened via SVD, namely (In − 1n11 T )DiX=UiΣiV T i . This renders whitened features
U′i= [ Ui 1/ √ n ]
appearing in the convex model. This is in contrast with applying BN before the ReLU layer, which directly operates on the data matrix X.
Theorem G.2 shows that an optimal solution can be found in closed-form in the high-dimensional regime n≤d when X is full row-rank, which is typically satisfied in sample-limited scenarios. Surprisingly, a single neuron optimizes the training cost, which corresponds to a single non-zero block in the group `1 penalized program (46).

heading: G.1 PROOF OF THEOREM G.2
text: Applying the approach in Lemma C.1 to the following constraint
max θ∈Θs ∣∣∣∣∣vT (In − 1 n11 T ) ( Xw(1) ) + ‖(In − 1n11T ) ( Xw(1) ) + ‖2 γ(1) + α(1)vT 1√ n ∣∣∣∣∣ yields
v∗=β y
‖y‖2 Then, the rest of proof is to verify the solutions by substituting them into the primal and dual problems as in Theorem 2.1.

heading: G.2 PROOF OF THEOREM G.1
text: We first focus on the dual constraint as follows
max w(1)
γ(1) 2 +α(1) 2 =1
∣∣∣∣∣ v T (In − 1n11 T ) ( Xw(1) ) + ‖(In − 1n11T ) ( Xw(1) ) + ‖2 γ(1) + α(1) vT1√ n ∣∣∣∣∣ = max i∈[P ] max θ∈Θs
(2Di−In)Xw(1)≥0
∣∣∣∣∣ vT (In − 1n11T )DiXw(1)‖(In − 1n11T )DiXw(1)‖2 γ(1) + α(1) v T1√ n ∣∣∣∣∣ = max i∈[P ] max qi,α,γ
α2+γ2=1
(2Di−In)XViΣ−1i qi≥0
∣∣∣∣vT Uiqi‖qi‖2 γ + αvT 1√n ∣∣∣∣
= max i∈[P ]
max si
‖si‖2=1 (2Di−In)XViΣ−1i s1:ri−1≥0
∣∣∣∣∣∣∣∣v T [ Ui 1√ n ] ︸ ︷︷ ︸
U′i
s ∣∣∣∣∣∣∣∣ . Then, the rest of the proof directly follow from Theorem 2.2.

Review summaries:
The claim of this paper is casting training neural networks with batch normalization to a convex program solvable in polynomial time. The convex reduction sparks an implicit regularization of batch normalization. Taking inspiration from the convex program and the implicit regularization, the authors improve BN. 

The paper studies batch normalization in deep neural networks. For a two-layer network with scalar output and batch normalization, a dual of the problem is derived. It is then shown that in the high-dimensional regime, the dual can be further simplified so that an optimal solution can be computed in closed form. In the general case, the concept of hyperplane arrangements can be used to formulate an equivalent finite-dimensional convex program which can be solved in polynomial time.   The analysis is then extended for vector output networks and different architectures including L-layer neural networks and CNNs. It is then shown that the convex problems tend to fit low singular value directions which will lead to poorer generalization. As a remedy, the authors propose a truncated variant of the problem by obtaining a low-rank approximation of the data.  In the experimental section, it is shown on the CIFAR-100 dataset that the closed form solution leads to superior results compared to the solution obtained by gradient descent. Moreover it is shown that the truncated variant of the approach is needed to achieve a good generalization performance.