Globally Optimal Training of Neural Networks with Threshold Activation
  Functions
Threshold activation functions are highly preferable in neural networks due
to their efficiency in hardware implementations. Moreover, their mode of
operation is more interpretable and resembles that of biological neurons.
However, traditional gradient based algorithms such as Gradient Descent cannot
be used to train the parameters of neural networks with threshold activations
since the activation function has zero gradient except at a single
non-differentiable point. To this end, we study weight decay regularized
training problems of deep neural networks with threshold activations. We first
show that regularized deep threshold network training problems can be
equivalently formulated as a standard convex optimization problem, which
parallels the LASSO method, provided that the last hidden layer width exceeds a
certain threshold. We also derive a simplified convex optimization formulation
when the dataset can be shattered at a certain layer of the network. We
corroborate our theoretical results with various numerical experiments.
Tolga Ergen, Halil Ibrahim Gulluk, Jonathan Lacotte, Mert Pilanci
# Globally Optimal Training of Neural Networks with Threshold Activation Functions

###### Abstract

Threshold activation functions are highly preferable in neural networks due to their efficiency in hardware implementations. Moreover, their mode of operation is more interpretable and resembles that of biological neurons. However, traditional gradient based algorithms such as Gradient Descent cannot be used to train the parameters of neural networks with threshold activations since the activation function has zero gradient except at a single non-differentiable point. To this end, we study weight decay regularized training problems of deep neural networks with threshold activations. We first show that regularized deep threshold network training problems can be equivalently formulated as a standard convex optimization problem, which parallels the LASSO method, provided that the last hidden layer width exceeds a certain threshold. We also derive a simplified convex optimization formulation when the dataset can be shattered at a certain layer of the network. We corroborate our theoretical results with various numerical experiments.

## 1 Introduction

In the past decade, deep neural networks have proven remarkably useful in solving challenging problems and become popular in many applications. The choice of activation plays a crucial role in their performance and practical implementation. In particular, even though neural networks with popular activation functions such as ReLU are successfully employed, they require advanced computational resources in training and evaluation, e.g., Graphical Processing Units (GPUs) (Coates et al., 2013). Consequently, training such deep networks is challenging especially without sophisticated hardware. On the other hand, the threshold activation offers a multitude of advantages: (1) computational efficiency, (2) compression/quantization to binary latent dimension, (3) interpretability. Unfortunately, gradient based optimization methods fail in optimizing threshold activation networks due to the fact that the gradient is zero almost everywhere. To close this gap, we analyze the training problem of deep neural networks with the threshold activation function defined as

\[\sigma_{s}(x)\!:=s\mathbbm{1}\{x\geq 0\}=\begin{cases}s&\text{if }x\geq 0\\ 0&\text{otherwise}\end{cases}, \tag{1}\]

where \(s\in\mathbb{R}\) is a trainable amplitude parameter for the neuron. Our main result is that globally optimal deep threshold networks can be trained by solving a convex optimization problem.

### Why should we care about threshold networks?

Neural networks with threshold activations are highly desirable due to the following reasons:

* Since the threshold activation (1) is restricted to take values in \(\{0,s\}\), threshold neural network models are far more suitable for hardware implementations (Barlett and Downs, 1992; Corwin et al., 1994). Specifically, these networks have significantly lower memory footprint, less computational complexity, and consume less energy (Helwegen et al., 2019).
* Modern neural networks have extremely large number of full precision trainable parameters so that several computational barriers emerge during hardware implementations. One approach to
mitigate these issues is reducing the network size by grouping the parameters via a hash function (Hubara et al., 2017; Chen et al., 2015). However, this still requires full precision training before the application of the hash function and thus fails to remedy the computational issues. On the other hand, neural networks with threshold activations need a minimal amount of bits.
* Another approach to reduce the complexity is to quantize the weights and activations of the network (Hubara et al., 2017) and the threshold activation is inherently in a two level quantized form.
* The threshold activation is a valid model to simulate the behaviour of a biological neuron as detailed in Jain et al. (1996). Therefore, progress in this research field could shed light into the connection between biological and artificial neural networks.

### Related Work

Although threshold networks are essential for several practical applications as detailed in the previous section, training their parameters is a difficult non-differentiable optimization problem due to the discrete nature in (1). For training of deep neural networks with popular activations, the common practice is to use first order gradient based algorithms such as Gradient Descent (GD) since the well known backpropagation algorithm efficiently calculates the gradient with respect to parameters. However, the threshold activation in (1) has zero gradient except at a single non-differentiable point zero, and therefore, one cannot directly use gradient based algorithms to train the parameters of the network. In order to remedy this issue numerous heuristic algorithms have been proposed in the literature as detailed below but they still fail to globally optimize the training objective (see Figure 1).

The Straight-Through Estimator (STE) is a widely used heuristics to train threshold networks (Bengio et al., 2013; Hinton, 2012). Since the gradient is zero almost everywhere, Bengio et al. (2013); Hinton (2012) proposed replacing the threshold activation with the identity function during only the backward pass. Later on, this approach is extended to employ various forms of the ReLU activation function, e.g., clipped ReLU, vanilla ReLU, Leaky ReLU (Yin et al., 2019; Cai et al., 2017; Xiao et al.), during the backward pass. Additionally, clipped versions of the identity function were also used as an alternative to STE (Hubara et al., 2017; Courbariaux et al., 2016; Rastegari et al., 2016).

### Contributions

* We introduce polynomial-time trainable convex formulations of regularized deep threshold network training problems provided that a layer width exceeds a threshold detailed in Table 1.
* In Theorem 2.2, we prove that the original non-convex training problem for two-layer networks is equivalent to standard convex optimization problems.
* We show that deep threshold network training problems are equivalent to standard convex optimization problems in Theorem 3.2. In stark contrast to two-layer networks, deep threshold networks can have a richer set of hyperplane arrangements due to multiple nonlinear layers (see Lemma 3.5).
* In Section 3.1, we characterize the evolution of the set of hyperplane arrangements and consequently hidden layer representation space as a recursive process (see Figure 3) as the network gets deeper.
* We prove that when a certain layer width exceeds \(\mathcal{O}(\sqrt{n}/L)\), the regularized \(L\)-layer threshold network training further simplifies to a problem that can be solved in \(\mathcal{O}(n)\) time.

Figure 1: Training comparison of our convex program in (7) with the non-convex training heuristic STE. We also indicate the time taken to solve the convex programs with markers. For non-convex STE, we repeat the training with \(5\) different initializations. In each case, our convex training algorithms achieve lower objective than all the non-convex heuristics (see Appendix B.5 for details).
## 2 Two-layer threshold networks

We first consider the following two-layer threshold network

\[f_{\theta,2}(\mathbf{X}){=}\sigma_{\mathbf{s}}(\mathbf{X}\mathbf{W}^{(1)}) \mathbf{w}^{(2)}{=}\sum_{j=1}^{m}s_{j}\mathbb{1}\left\{\mathbf{X}\mathbf{w}_{j}^{ (1)}{\geq}0\right\}w_{j}^{(2)}, \tag{2}\]

where the set of the trainable parameters are \(\mathbf{W}^{(1)}\in\mathbb{R}^{d\times m},\mathbf{s}\in\mathbb{R}^{m},\mathbf{ w}^{(2)}\in\mathbb{R}^{m}\) and \(\theta\) is a compact representation for the parameters, i.e., \(\theta:=\{\mathbf{W}^{(1)},\mathbf{s},\mathbf{w}^{(2)}\}\). Note that we include bias terms by concatenating a vector of ones to \(\mathbf{X}\). Next, consider the weight decay regularized training objective

\[\mathcal{P}_{2}^{\mathrm{noncvx}}\,{:=}\,\min_{\mathbf{W}^{(1)},\mathbf{s}, \mathbf{w}^{(2)}}\frac{1}{2}\left\|f_{\theta,2}(\mathbf{X})-\mathbf{y}\right\| _{2}^{2}+\frac{\beta}{2}\sum_{j=1}^{m}\left(\|\mathbf{w}_{j}^{(1)}\|_{2}^{2}+ |s_{j}|^{2}+|w_{j}^{(2)}|^{2}\right)\,. \tag{3}\]

Now, we apply a scaling between variables \(s_{j}\) and \(w_{j}^{(2)}\) to reach an equivalent optimization problem.

**Lemma 2.1** (Optimal scaling).: _The training problem in (3) can be equivalently stated as_

\[\mathcal{P}_{2}^{\mathrm{noncvx}}=\min_{\theta\in\Theta,}\frac{1}{2}\left\|f_ {\theta,2}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\beta\|\mathbf{w}^{(2)}\|_{ 1}\,, \tag{4}\]

_where \(\Theta_{s}{:=}\left\{\theta:|s_{j}|=1,\forall j\in[m]\right\}\)._

We next define the set of hyperplane arrangement patterns of the data matrix \(\mathbf{X}\) as

\[\mathcal{H}(\mathbf{X}){:=}\left\{\mathbb{1}\{\mathbf{X}\mathbf{w}\geq 0\}: \mathbf{w}\in\mathbb{R}^{d}\right\}\subset\{0,1\}^{n}. \tag{5}\]

We denote the distinct elements of the set \(\mathcal{H}(\mathbf{X})\) by \(\mathbf{d}_{1},\ldots,\mathbf{d}_{P}\in\{0,1\}^{n}\), where \(P{:=}\left|\mathcal{H}(\mathbf{X})\right|\) is the number of hyperplane arrangements. Using this fixed set of hyperplane arrangements \(\{\mathbf{d}_{i}\}_{i=1}^{P}\), we next prove that (4) is equivalent to the standard Lasso method (Tibshirani, 1996).

**Theorem 2.2**.: _Let \(m\geq m^{*}\), then the non-convex regularized training problem (4) is equivalent to_

\[\mathcal{P}_{2}^{\mathrm{cvx}}=\min_{\mathbf{w}\in\mathbb{R}^{P}}\frac{1}{2} \left\|\mathbf{D}\mathbf{w}-\mathbf{y}\right\|_{2}^{2}+\beta\|\mathbf{w}\|_{ 1}, \tag{6}\]

_where \(\mathbf{D}=[\mathbf{d}_{1}\quad\mathbf{d}_{2}\quad\ldots\quad\mathbf{d}_{P}]\) is a fixed \(n\times P\) matrix. Here, \(m^{*}\) is the cardinality of the optimal solution, which satisfies \(m^{*}\leq n+1\). Also, it holds that \(\mathcal{P}_{2}^{\mathrm{noncvx}}=\mathcal{P}_{2}^{\mathrm{cvx}}\)._

Theorem 2.2 proves that the original non-convex training problem in (3) can be equivalently solved as a standard convex Lasso problem using the hyperplane arrangement patterns as features. Surprisingly, the non-zero support of the convex optimizer in (6) matches to that of the optimal weight-decay regularized threshold activation neurons in the non-convex problem (3). This brings us two major advantages over the standard non-convex training:

* Since (6) is a standard convex problem, it can be globally optimized without resorting to non-convex optimization heuristics, e.g., initialization schemes, learning rate schedules etc.
* Since (6) is a convex Lasso problem, there exists many efficient solvers (Efron et al., 2004).

\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline
**Result** & **Depth** & **Complexity** & **Minimum width** & **Globally optimal** \\ \hline Theorem 2.2 & \(2\) & \(\mathcal{O}(n^{\alpha})\) & \(m\geq m^{*}\) & \(\mathcal{J}(\text{convex opt})\) \\ \hline Theorem 2.3 & \(2\) & \(\mathcal{O}(n)\) & \(m\geq n+2\) & \(\mathcal{J}(\text{convex opt})\) \\ \hline Theorem 3.2 & \(L\) & \(\mathcal{O}(n^{\gamma}\mathbb{I}_{\mathbb{I}_{\mathbb{I}_{\mathbb{I}_{\mathbb{I} _{\mathbb{I}_{\mathbb{I}_{\mathbb{I}_{\mathbb{I}_{\mathbb{I}_{\mathbb{I}}}}}}}}}}}^ {m_{1}})\) & \(m_{L-1}\geq m^{*}\) & \(\mathcal{J}(\text{convex opt})\) \\ \hline Corollary 3.4 & \(L\) & \(\mathcal{O}(n)\) & \(\exists l:m\geq\mathcal{O}(\sqrt{n}/L)\) & \(\mathcal{J}(\text{convex opt})\) \\ \hline \end{tabular}
\end{table}
Table 1: Summary of our results for the optimization of weight decay regularized threshold network training problems (\(n\): \(\#\) of data samples, \(d\): feature dimension, \(m_{l}\): \(\#\) of hidden neurons in layer \(l\), \(r\): rank of the training data matrix, \(m^{*}\): critical width, i.e., \(\#\) of neurons that obeys \(0\leq m^{*}\leq n+1\))
### Simplified convex formulation for complete arrangements

We now show that if the set of hyperplane arrangements of \(\mathbf{X}\) is complete, i.e., \(\mathcal{H}=\{0,1\}^{n}\) contains all boolean sequences of length \(n\), then the non-convex optimization problem in (3) can be simplified. We call these instances _complete arrangements_. In the case of two-layer threshold networks, complete arrangements emerge when the width of the network exceeds a threshold, specifically \(m\geq n+2\).

We note that the \(m\geq n\) regime, also known as memorization, has been extensively studied in the recent literature (Bubeck et al., 2020; de Dios and Bruna, 2020; Pilanci and Ergen, 2020; Rosset et al., 2007). Particularly, these studies showed that as long as the width exceeds the number of samples, there exists a neural network model that can exactly fit an arbitrary dataset. Vershynin (2020); Bartlett et al. (2019) further improved the condition on the width by utilizing the expressive power of deeper networks and developed more sophisticated weight construction algorithms to fit the data.

**Theorem 2.3**.: _We assume that the set of hyperplane arrangements of \(\mathbf{X}\) is complete, i.e., equal to the set of all length-\(n\) Boolean sequences" \(\mathcal{H}=\{0,1\}^{n}\). Suppose \(m\geq n+2\), then (4) is equivalent to_

\[\mathcal{P}^{\mathrm{cvx}}_{v2}\!:=\min_{\mathbf{\delta}\in\mathbb{R}^{n}}\frac{ 1}{2}\|\mathbf{\delta}-\mathbf{y}\|_{2}^{2}+\beta(\|(\mathbf{\delta})_{+}\|_{\infty}+ \|(-\mathbf{\delta})_{+}\|_{\infty})\,. \tag{7}\]

_and it holds that \(\mathcal{P}^{\mathrm{convx}}_{2}=\mathcal{P}^{\mathrm{cvx}}_{2}\). Also, one can construct an optimal network with \(n+2\) neurons in time \(\mathcal{O}(n)\) based on the optimal solution to the convex problem (7)._

Based on Theorem 2.3, when the data matrix can be shattered, i.e., all \(2^{n}\) possible \(\mathbf{y}\in\{0,1\}^{n}\) labelings of the data points can be separated via a linear classifier, it follows that the set of hyperplane arrangements is complete. Consequently, the non-convex problem in (3) further simplifies to (7).

### Training complexity

We first briefly summarize our complexity results for solving (6) and (7), and then provide the details of the derivations below. Our analysis reveals two interesting regimes:

* (incomplete arrangements) When \(n+1\geq m\geq m^{*}\), we can solve (6) in \(\mathcal{O}(n^{3r})\) complexity, where \(r:=\text{rank}(\mathbf{X})\). Notice this is **polynomial-time** whenever the rank \(r\) is fixed.
* (complete arrangements) When \(m\geq n+2\), we can solve (7) in closed-form, and the reconstruction of the non-convex parameters \(\mathbf{W}^{(1)}\) and \(\mathbf{w}^{(2)}\) only \(\mathcal{O}(n)\) time independent of \(d\).

**Computational complexity of (6):** To solve the optimization problem in (6), we first enumerate all possible hyperplane arrangements \(\{\mathbf{d}_{i}\}_{i=1}^{P}\). It is well known that given a rank-\(r\) data matrix the number of hyperplane arrangements \(P\) is upper bounded by (Stanley et al., 2004; Cover, 1965)

\[P\leq 2\sum_{k=0}^{r-1}\binom{n-1}{k}\leq 2r\left(\frac{e(n-1)}{r}\right)^{r}, \tag{8}\]

where \(r=\text{rank}(\mathbf{X})\leq\min(n,d)\). Furthermore, these can be enumerated in \(\mathcal{O}(n^{r})\)(Edelsbrunner et al., 1986). Then, the complexity for solving (6) is \(\mathcal{O}(P^{3})\approx\mathcal{O}(n^{3r})\)(Efron et al., 2004).

**Computational complexity of (7):** The problem in (7) is the proximal operator of a polyhedral norm. Since the problem is separable over the positive and negative parts of the parameter vector \(\mathbf{\delta}\), the optimal solution can be obtained by applying two proximal steps (Parikh and Boyd, 2014). As noted in Theorem 2.3, the reconstruction of the non-convex parameters \(\mathbf{W}^{(1)}\) and \(\mathbf{w}^{(2)}\) requires \(\mathcal{O}(n)\) time.

### A geometric interpretation

To provide a geometric interpretation, we consider the weakly regularized case where \(\beta\to 0\). In this case, (6) reduces to the following minimum norm interpolation problem

\[\min_{\mathbf{w}\in\mathbb{R}^{P}}\|\mathbf{w}\|_{1}\quad\mathrm{s.t.}\quad \mathbf{D}\mathbf{w}=\mathbf{y}. \tag{9}\]

**Proposition 2.4**.: _The minimum \(\ell_{1}\) norm interpolation problem in (9) can be equivalently stated as_

\[\min_{t\geq 0}t\quad\mathrm{s.t.}\quad\mathbf{y}\in t\mathrm{Conv}\{\pm\mathbf{ d}_{j},\forall j\in[P]\},\]
_where \(\mathrm{Conv}(\mathcal{A})\) denotes the convex hull of a set \(\mathcal{A}\). This corresponds to the gauge function (see Rockafellar (2015)) of the hyperplane arrangement patterns and their negatives. We provide visualizations of the convex set \(\mathrm{Conv}\{\pm\mathbf{d}_{j},\forall j\in[P]\}\) in Figures 3 and 4 (see Appendix) for Example 3.1._

Proposition 2.4 implies that the non-convex threshold network training problem in (3) implicitly represents the label vector \(\mathbf{y}\) as the convex combination of the hyperplane arrangements determined by the data matrix \(\mathbf{X}\). Therefore, we explicitly characterize the representation space of threshold networks. We also remark that this interpretation extends to arbitrary depth as shown in Theorem 3.2.

### How to optimize the hidden layer weights?

After training via the proposed convex programs in (6) and (7), we need to reconstruct the layer weights of the non-convex model in (2). We first construct the optimal hyperplane arrangements \(\mathbf{d}_{i}\) for (7) as detailed in Appendix A.4. Then, we have the following prediction model (2). Notice changing \(\mathbf{w}_{j}^{(1)}\) to any other vector \(\mathbf{w}_{j}^{\prime}\) with same norm and such that \(\mathbb{I}\left\{\mathbf{X}\mathbf{w}_{j}^{(1)}\geq 0\right\}=\mathbb{I} \left\{\mathbf{X}\mathbf{w}_{j}^{\prime}\geq 0\right\}\) does not change the optimal training objective. Therefore, there are multiple global optima and the one we choose might impact the generalization performance as discussed in Section 5.

## 3 Deep threshold networks

We now analyze \(L\)-layer parallel deep threshold networks model with \(m_{L-1}\) subnetworks defined as

\[f_{\theta,L}(\mathbf{X})=\sum_{k=1}^{m_{L-1}}\sigma_{\mathbf{s}^{(L-1)}}( \mathbf{X}_{k}^{(L-2)}\mathbf{w}_{k}^{(L-1)})w_{k}^{(L)}, \tag{10}\]

where \(\theta\!:=\ \{\{\mathbf{W}_{k}^{(l)},\mathbf{s}_{k}^{(l)}\}_{l=1}^{L}\}_{k=1}^{m_{L-1}}\), \(\theta\in\Theta\!:=\{\theta:\mathbf{W}_{k}^{(l)}\in\mathbb{R}^{m_{l-1}\times m _{l}},\mathbf{s}_{k}^{(l)}\in\mathbb{R}^{m_{l}},\forall l,k\}\),

\[\mathbf{X}_{k}^{(0)}\!:=\mathbf{X}, \mathbf{X}_{k}^{(l)}\!:=\sigma_{\mathbf{s}_{k}^{(l)}}(\mathbf{X}_ {k}^{(l-1)}\mathbf{W}_{k}^{(l)}),\forall l\in[L-1],\]

and the subscript \(k\) is the index for the subnetworks (see Ergen and Pilanci (2021a,b); Wang et al. (2023) for more details about parallel networks). We next show that the standard weight decay regularized training problem can be cast as an \(\ell_{1}\)-norm minimization problem as in Lemma 2.1.

**Lemma 3.1**.: _The following \(L\)-layer regularized threshold network training problem_

\[\mathcal{P}_{L}^{\mathrm{noncvx}}\!\!=\!\!\!\min_{\theta\in\Theta}\!\frac{1}{ 2}\|f_{\theta,L}(\mathbf{X})\!-\!\mathbf{y}\|_{2}^{2}\!+\!\frac{\beta}{2}\sum _{k=1}^{m_{L-1}}\!\sum_{l=1}^{L}(\|\mathbf{W}_{k}^{(l)}\|_{F}^{2}\!+\!\| \mathbf{s}_{k}^{(l)}\|_{2}^{2}) \tag{11}\]

_can be reformulated as_

\[\mathcal{P}_{L}^{\mathrm{noncvx}}=\min_{\theta\in\Theta}\frac{1}{2}\left\|f _{\theta,L}(\mathbf{X})-\mathbf{y}\right\|_{2}^{2}+\beta\|\mathbf{w}^{(L)}\|_ {1}, \tag{12}\]

_where \(\Theta_{s}\!:=\{\theta\in\Theta:|s_{k}^{(L-1)}|=1\}\)._

### Characterizing the set of hyperplane arrangements for deep networks

We first define hyperplane arrangements for \(L\)-layer networks with a single subnetwork (i.e., \(m_{L-1}=1\), thus, we drop the index \(k\)). We denote the set of hyperplane arrangements as

\[\mathcal{H}_{L}(\mathbf{X})\!:=\{\mathbb{I}\left\{\mathbf{X}^{(L-2)}\mathbf{w }^{(L-1)}\geq 0\right\}:\theta\in\Theta\}.\]

We also denote the elements of the set \(\mathcal{H}_{L}(\mathbf{X})\) by \(\mathbf{d}_{1},\ldots,\mathbf{d}_{P_{L-1}}\in\{0,1\}^{n}\) and \(P_{L-1}=|\mathcal{H}_{L}(\mathbf{X})|\) is the number of hyperplane arrangements in the layer \(L-1\).

To construct the hyperplane arrangement matrix \(\mathbf{D}^{(l)}\), we define a matrix valued operator as follows

\[\mathbf{D}^{(1)}\!:=\mathcal{A}\left(\mathbf{X}\right)=[\mathbf{d}_{1}\quad \mathbf{d}_{2}\quad\ldots\quad\mathbf{d}_{P_{1}}]\,\ \mathbf{D}^{(l+1)}\!:=\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
Here, the operator \(\mathcal{A}\left(\cdot\right)\) outputs a matrix whose columns contain all possible hyperplane arrangements corresponding to its input matrix as in (5). In particular, \(\mathbf{D}^{(1)}\) denotes the arrangements for the first layer given the input matrix \(\mathbf{X}\). The notation \(\mathbf{D}^{(l)}_{\mathcal{S}}\in\{0,1\}^{n\times m_{l}}\) denotes the submatrix of \(\mathbf{D}^{(l)}\in\{0,1\}^{n\times P_{l}}\) indexed by the subset \(\mathcal{S}\) of its columns, where the index \(\mathcal{S}\) runs over all subsets of size \(m_{l}\). Finally, \(\sqcup\) is an operator that takes a union of these column vectors and outputs a matrix of size \(n\times P_{l+1}\) containing these as columns. Note that we may omit repeated columns and denote the total number of unique columns as \(P_{l+1}\), since this does not change the value of our convex program.

We next provide an analytical example describing the construction of the matrix \(\mathbf{D}^{(l)}\).

**Example 3.1**.: _We illustrate an example with the training data \(\mathbf{X}=[-1\ 1;0\ 1;1\ 1]\in\mathbb{R}^{3\times 2}\). Inspecting the data samples (rows of \(\mathbf{X}\)), we observe that all possible arrangement patterns are_

\[\mathbf{D}^{(1)}=\mathcal{A}\left(\mathbf{X}\right)=\begin{bmatrix}0&0&0&1&1& 1\\ 0&0&1&1&1&0\\ 0&1&1&1&0&0\end{bmatrix}\implies P_{1}=6. \tag{14}\]

_For the second layer, we first specify the number of neurons in the first layer as \(m_{1}=2\). Thus, we need to consider all possible column pairs in (14). We have_

\[\mathbf{D}^{(1)}_{\{1,2\}}=\begin{bmatrix}0&0\\ 0&0\\ 0&1\end{bmatrix} \implies\mathcal{A}\left(\begin{bmatrix}0&0\\ 0&0\\ 0&1\end{bmatrix}\right)=\begin{bmatrix}0&0&1&1\\ 0&0&1&1\\ 0&1&1&0\end{bmatrix}\] \[\mathbf{D}^{(1)}_{\{1,3\}}=\begin{bmatrix}0&0\\ 0&1\\ 0&1\end{bmatrix} \implies\mathcal{A}\left(\begin{bmatrix}0&0\\ 0&1\\ 0&1\end{bmatrix}\right)=\begin{bmatrix}0&0&1&1\\ 0&1&1&0\\ 0&1&1&0\end{bmatrix}\] \[\vdots\]

_We then construct the hyperplane arrangement matrix as_

\[\mathbf{D}^{(2)}=\bigsqcup_{|\mathcal{S}|=2}\mathcal{A}\left(\mathbf{D}^{(1)}_ {\mathcal{S}}\right)=\begin{bmatrix}0&0&0&0&1&1&1&1\\ 0&0&1&1&0&0&1&1\\ 0&1&0&1&0&1&0&1\end{bmatrix},\]

_which shows that \(P_{2}=8\). Consequently, we obtain the maximum possible arrangement patterns, i.e., \(\{0,1\}^{3}\), in the second layer even though we are not able to obtain some of these patterns in the first layer in (14). We also provide a three dimensional visualization of this example in Figure 3._

### Polynomial-time trainable convex formulation

Based on the procedure described in Section 3.1 to compute the arrangement matrix \(\mathbf{D}^{(l)}\), we now derive an exact formulation for the non-convex training problem in (12).

**Theorem 3.2**.: _Suppose that \(m_{L-1}\geq m^{*}\), then the non-convex training problem (12) is equivalent to_

\[\mathcal{P}^{\mathrm{cvx}}_{L}=\min_{\mathbf{w}\in\mathbb{R}^{P_{L-1}}}\frac{1 }{2}\left\|\mathbf{D}^{(L-1)}\mathbf{w}-\mathbf{y}\right\|_{2}^{2}+\beta\| \mathbf{w}\|_{1}, \tag{15}\]

_where \(\mathbf{D}^{(L-1)}\in\{0,1\}^{n\times P_{L-1}}\) is a fixed matrix constructed via (13) and \(m^{*}\) denotes the cardinality of the optimal solution, which satisfies where \(m^{*}\leq n+1\). Also, it holds that \(\mathcal{P}^{\mathrm{noncvx}}_{L}=\mathcal{P}^{\mathrm{cvx}}_{L}\)._

Theorem 3.2 shows that two-layer and deep networks simplify to very similar convex Lasso problems (i.e. (6) and (15)). However, the set of hyperplane arrangements is larger for deep networks as analyzed in Section 3.1. Thus, the structure of the diagonal matrix \(\mathbf{D}\) and the problem dimensionality are significantly different for these problems.

### Simplified convex formulation

Here, we show that the data \(\mathbf{X}\) can be shattered at a certain layer, i.e., \(\mathcal{H}_{l}(\mathbf{X})=\{0,1\}^{n}\) for a certain \(l\in[L]\), if the number of hidden neurons in a certain layer \(m_{l}\) satisfies \(m_{l}\geq C\sqrt{n}/L\). Then we can alternatively formulate (11) as a simpler convex problem. Therefore, compared to the two-layer networks in Section 2.1, we substantially improve the condition on layer width by benefiting from the depth \(L\), which also confirms the benign impact of the depth on the optimization.
**Lemma 3.3**.: _If \(\exists l,C\) such that \(m_{l}\geq C\sqrt{n}/L\), then the set of hyperplane arrangements is complete, i.e., \(\mathcal{H}_{L}(\mathbf{X})=\{0,1\}^{n}\)._

We next use Lemma 3.3 to derive a simpler form of (11).

**Corollary 3.4**.: _As a direct consequence of Theorem 2.3 and Lemma 3.3, the non-convex deep threshold network training problem in (12) can be cast as the following convex program_

\[\mathcal{P}_{L}^{\mathrm{convex}}=\mathcal{P}_{v2}^{\mathrm{cvx}}=\min_{ \boldsymbol{\delta}\in\mathbb{R}^{n}}\frac{1}{2}\|\boldsymbol{\delta}- \mathbf{y}\|_{2}^{2}+\beta(\|(\boldsymbol{\delta})_{+}\|_{\infty}+\|(- \boldsymbol{\delta})_{+}\|_{\infty})\,.\]

Surprisingly, both two-layer and deep networks share the same convex formulation in this case. However, notice that two-layer networks require a condition on the data matrix in Theorem 2.3 whereas the result in Corollary 3.4 requires a milder condition on the layer widths.

### Training complexity

Here, we first briefly summarize our complexity results for the convex training of deep networks. Based on the convex problems in (15) and Corollary 3.4, we have two regimes:

* When \(n+1\geq m_{L-1}\geq m^{*}\), we solve (15) with \(\mathcal{O}(n^{3r}\prod_{k=1}^{L-2}m_{k})\) complexity, where \(r:=\text{rank}(\mathbf{X})\). Note that this is **polynomial-time** when \(r\) and the number of neurons in each layer \(\{m_{l}\}_{l=1}^{L-2}\) are constants.
* When \(\exists l,C:m_{l}\geq C\sqrt{n}/L\), we solve (7) in closed-form, and the reconstruction of the non-convex parameters requires \(\mathcal{O}(n)\) time as proven in Appendix A.9.

**Computational complexity for (15):** We first need to obtain an upperbound on the problem dimensionality \(P_{L-1}\), which is stated in the next result.

**Lemma 3.5**.: _The cardinality of the hyperplane arrangement set for an \(L\)-layer network \(\mathcal{H}_{L}(\mathbf{X})\) can be bounded as \(|\mathcal{H}_{L}(\mathbf{X})|=P_{L-1}\lesssim\mathcal{O}(n^{r\prod_{k=1}^{L-2 }m_{k}})\), where \(r=\mathrm{rank}(\mathbf{X})\) and \(m_{l}\) denotes the number of hidden neurons in the \(l^{th}\) layer._

Lemma 3.5 shows that the set of hyperplane arrangements gets significantly larger as the depth of the network increases. However, the cardinality of this set is still a polynomial term since that \(r<\min\{n,d\}\) and \(\{m_{l}\}_{l=1}^{L-2}\) are fixed constants.

To solve (15), we first enumerate all possible arrangements \(\{\mathbf{d}_{i}\}_{i=1}^{P_{L-1}}\) to construct the matrix \(\mathbf{D}^{(L-1)}\). Then, we solve a standard convex Lasso problem, which requires \(\mathcal{O}(P_{L-1}^{3})\) complexity (Efron et al., 2004). Thus, based on Lemma 3.5, the overall complexity is \(\mathcal{O}(P_{L-2}^{3})\approx\mathcal{O}(n^{3r\prod_{k=1}^{L-2}m_{k}})\).

**Computational complexity for (7):** Since Corollary 3.4 yields (7), the complexity is \(\mathcal{O}(n)\) time.

## 4 Extensions to arbitrary loss functions

In the previous sections, we considered squared error as the loss function to give a clear description of our approach. However, all the derivations extend to arbitrary convex loss. Now, we consider the regularized training problem with a convex loss function \(\mathcal{L}(\cdot,\mathbf{y})\), e.g., hinge loss, cross entropy,

\[\min_{\boldsymbol{\theta}\in\Theta}\mathcal{L}(f_{\theta,L}(\mathbf{X}), \mathbf{y})+\frac{\beta}{2}\sum_{k=1}^{m_{L-1}}\sum_{l=1}^{L}(\|\mathbf{W}_{k} ^{(l)}\|_{F}^{2}+\|\mathbf{s}_{k}^{(l)}\|_{2}^{2}). \tag{16}\]

Then, we have the following generic loss results.

**Corollary 4.1**.: _Theorem 3.2 implies that when \(m_{L-1}\geq m^{*}\), (16) can be equivalently stated as_

\[\mathcal{P}_{L}^{\mathrm{cvx}}=\min_{\mathbf{w}\in\mathbb{R}^{P_{L-1}}} \mathcal{L}\left(\mathbf{D}^{(L-1)}\mathbf{w},\mathbf{y}\right)+\beta\| \mathbf{w}\|_{1}. \tag{17}\]

_Alternatively when \(\mathcal{H}_{L}(\mathbf{X})=\{0,1\}^{n}\), based on Corollary 3.4, the equivalent convex problem is_

\[\min_{\boldsymbol{\delta}\in\mathbb{R}^{n}}\mathcal{L}(\boldsymbol{\delta}, \mathbf{y})+\beta(\|(\boldsymbol{\delta})_{+}\|_{\infty}+\|(-\boldsymbol{ \delta})_{+}\|_{\infty}). \tag{18}\]

Corollary 4.1 shows that (17) and (18) are equivalent to the non-convex training problem in (16). More importantly, they can be globally optimized via efficient convex optimization solvers.
## 5 Experiments

In this section1, we present numerical experiments verifying our theoretical results in the previous sections. As discussed in Section 2.4, after solving the proposed convex problems in (15) and (7), there exist multiple set of weight matrices yielding the same optimal objective value. Therefore, to have a good generalization performance on test data, we use some heuristic methods for the construction of the non-convex parameters \(\{\mathbf{W}^{(l)}\}_{l=1}^{L}\). Below, we provide details regarding the weight construction and review some baseline non-convex training methods.

Footnote 1: We provide additional experiments and details in Appendix B.

**Convex-Lasso:** To solve the problem (15), we first approximate arrangement patterns of the data matrix \(\mathbf{X}\) by generating i.i.d. Gaussian weights \(\mathbf{G}\in\mathbb{R}^{d\times\hat{P}}\) and subsample the arrangement patterns via \(\mathbb{I}\left[\mathbf{X}\mathbf{G}\geq 0\right]\). Then, we use \(\mathbf{G}\) as the hidden layer weights to construct the network. We repeat this process for every layer. Notice that here we sample a fixed subset of arrangements instead of enumerating all possible \(P\) arrangements. Thus, this approximately solves (15) by subsampling its decision variables, however, it still performs significantly better than standard non-convex training.

**Convex-PI:** After solving (7), to recover the hidden layer weights of (2), we solve \(\mathbb{I}\left\{\mathbf{X}\mathbf{w}_{j}^{(1)}\geq 0\right\}=\mathbf{d}_{j}\) as \(\mathbf{w}_{j}^{(1)}=\mathbf{X}^{\dagger}\mathbf{d}_{j}\), where \(\mathbf{X}^{\dagger}\) denotes the pseudo-inverse of \(\mathbf{X}\). The resulting weights \(\mathbf{w}_{j}^{(1)}\) enforce the preactivations \(\mathbf{X}\mathbf{w}_{j}^{(1)}\) to be zero or one. Thus, if an entry is slightly higher or less than zero due to precision issues during the pseudo-inverse, it might give wrong output after the threshold activation. To avoid such cases, we use \(0.5\) threshold in the test phase, i.e., \(\mathbb{I}\left\{\mathbf{X}_{\mathrm{test}}\mathbf{w}_{j}^{(1)}\geq 0.5\right\}\).

**Convex-SVM:** Another approach to solve \(\mathbb{I}\left\{\mathbf{X}\mathbf{w}_{j}^{(1)}\geq 0\right\}=\mathbf{d}_{j}\) is to use Support Vector Machines

Figure 2: In this figure, we compare the classification performance of three-layer threshold networks trained with the setup described in Figure 1 for a single initialization trial. This experiment shows that our convex training approach not only provides the globally optimal training performance but also generalize remarkably well on the test data (see Appendix B.5 for details).
(SVMs), which find the maximum margin vector. Particularly, we set the zero entries of \(\mathbf{d}_{i}\) as \(-1\) and then directly run the SVM to get the maximum margin hidden neurons corresponding to this arrangement. Since the labels are in the form \(\{+1,-1\}\) in this case, we do not need additional thresholding as in the previous approach.

**Nonconvex-STE**(Bengio et al., 2013): This is the standard non-convex training algorithm, where the threshold activations is replaced with the identity function during the backward pass.

**STE Variants:** We also benchmark against variants of STE. Specifically, we replace the threshold activation with ReLU (**Nonconvex-ReLU**(Yin et al., 2019)), Leaky ReLU (**Nonconvex-LReLU**(Xiao et al.)), and clipped ReLU (**Nonconvex-CReLU**(Cai et al., 2017)) during the backward pass.

**Synthetic Datasets:** We compare the performances of **Convex-PI** and **Convex-SVM** trained via (7) with the non-convex heuristic methods mentioned above. We first run each non-convex heuristic for five different initializations and then plot the best performing one in Figure 1. This experiment clearly shows that the non-convex heuristics fail to achieve the globally optimal training performance provided by our convex approaches. For the same setup, we also compare the training and test accuracies for three different regimes, i.e., \(n>d\), \(n=d\), and \(n<d\). As seen in Figure 2, our convex approaches not only globally optimize the training objective but also generalize well on the test data.

**Real Datasets:** In Table 2, we compare the test accuracies of two-layer threshold network trained via our convex formulation in (15), i.e., **Convex-Lasso** and the non-convex heuristics mentioned above. For this experiment, we use CIFAR-10 (Krizhevsky et al., 2014), MNIST (LeCun), and the datasets in the UCI repository (Dua and Graff, 2017) which are preprocessed as in Fernandez-Delgado et al. (2014). Here, our convex training approach achieves the highest test accuracy for most of the datasets while the non-convex heuristics perform well only for a few datasets. Therefore, we also validates the good generalization capabilities of the proposed convex training methods on real datasets.

## 6 Conclusion

We proved that the training problem of regularized deep threshold networks can be equivalently formulated as a standard convex optimization problem with a fixed data matrix consisting of hyperplane arrangements determined by the data matrix and layer weights. Since the proposed formulation parallels the well studied Lasso model, we have two major advantages over the standard non-convex training methods: **1)** We globally optimize the network without resorting to any optimization heuristic or extensive hyperparameter search (e.g., learning rate schedule and initialization scheme); **2)** We efficiently solve the training problem using specialized solvers for Lasso. We also provided a computational complexity analysis and showed that the proposed convex program can be solved in polynomial-time. Moreover, when a layer width exceeds a certain threshold, a simpler alternative convex formulation can be solved in \(\mathcal{O}(n)\). Lastly, as a by product of our analysis, we characterize the recursive process behind the set of hyperplane arrangements for deep networks. Even though this set rapidly grows as the network gets deeper, globally optimizing the resulting Lasso problem still requires polynomial-time complexity for fixed data rank. We also note that the convex analysis proposed in this work is generic in the sense that it can be applied to various architectures including batch normalization (Ergen et al., 2022), vector output networks (Sahiner et al., 2020, 2021), polynomial activations (Bartan and Pilanci, 2021), GANs (Sahiner et al., 2022), autoregressive models (Gupta et al., 2021), and Transformers (Ergen et al., 2022; Sahiner et al., 2022).

\begin{table}
\begin{tabular}{l|c|c|c|c|c|c|c|c|c} \hline \multirow{2}{*}{**Dataset**} & \multicolumn{2}{c|}{**Convex-Lasso (Convex)**} & \multicolumn{2}{c|}{**Nonconvex-LReLU**} & \multicolumn{2}{c|}{**Nonconvex-LReLU**} & \multicolumn{2}{c}{**Nonconvex-LReLU**} \\ \cline{2-10}  & **Accuracy** & **Linear** & **Linear** & **Linear** & **Linear** & **Linear** & **Linear** & **Linear** & **Linear** \\ \hline CIFAR-10 & **0.86\(\pm\)0.008** & **5.9\(\pm\)0.13** & 0.83\(\pm\)0.004 & 83.3\(\pm\)0.004 & 83.5\(\pm\)0.014 & 83.5\(\pm\)0.17 & 0.78\(\pm\)0.004 & 79.1\(\pm\) 0.43 & 0.80\(\pm\)0.004 & 83.1\(\pm\) 0.45 \\ batch & 0.80\(\pm\)0.017 & 7.7\(\pm\)0.02 & 0.80\(\pm\)0.008 & 83.4\(\pm\)0.004 & 83.4\(\pm\)0.004 & 83.6\(\pm\)0.003 & 83.9\(\pm\)0.003 & 83.4\(\pm\)0.21 & 0.80\(\pm\)0.004 & 83.5\(\pm\) 0.11 \\ batch & 0.85\(\pm\)0.001 & 7.7\(\pm\)0.01 & 0.80\(\pm\)0.008 & 83.4\(\pm\)0.004 & 83.4\(\pm\)0.004 & 83.6\(\pm\)0.003 & 83.6\(\pm\)0.003 & 83.4\(\pm\) 0.21 & 0.80\(\pm\) 0.004 & 83.5\(\pm\) 0.11 \\ batch & **0.84\(\pm\)0.001** & 8.3\(\pm\)0.003 & 8.7\(\pm\)0.003 & 8.7\(\pm\)0.003 & 83.5\(\pm\)0.004 & 83.1\(\pm\)0.007 & 83.6\(\pm\)0.003 & 83.7\(\pm\) 0.19 & 0.95\(\pm\) 0.003 & 83.5\(\pm\) 0.11 \\ machine-10 & **0.85\(\pm\)0.001** & 8.4\(\pm\)0.002 & 8.0\(\pm\)0.003 & 8.0\(\pm\)0.001 & 83.4\(\pm\)0.003 & 83.0\(\pm\)0.004 & 83.5\(\pm\) 0.003 & 83.5\(\pm\) 0.16 & 0.87\(\pm\) 0.003 & 52.9\(\pm\) 0.13 \\ deep-10 & **0.85\(\pm\)0.001** & **0.85\(\pm\)0.002** & **0.85\(\pm\)0.001** & 83.0\(\pm\)0.003 & 83.0\(\pm\)0.004 & 83.5\(\pm\)0.004 & 83.5\(\pm\) 0.003 & 83.5\(\pm\) 0.16 & 0.87\(\pm\) 0.003 & 52.9\(\pm\) 0.13 \\ conv-24 & **0.87\(\pm\)0.003** & **2.28\(\pm\)0.003** & **0.87\(\pm\)0.003** & 83.1\(\pm\)0.003 & 83.1\(\pm\)0.003 & 83.0\(\pm\)0.003 & 83.2\(\pm\) 0.12 & 0.72\(\pm\) 0.003 & 83.7\(\pm\) 0.003 & 83.5\(\pm\) 0.003 \\ conv-24 & **0.87\(\pm\)0.003** & **2.28\(\pm\)0.003** & **0.87\(\pm\)0.003** & 83.1\(\pm\)0.003 & 83.1\(\pm\)0.003 & 83.1\(\pm\)0.003 & 83.5\(\pm\) 0.003 & 83.2\(\pm\) 0.003 & 83.5\(\pm\) 0.18 \\ conv-24 & **0.87\(\pm\)0.003** & **2.28\(\pm\)0.003** & **0.75\(\pm\)0.003** & 83.4\(\pm\)0.003 & 83.1\(\pm\)0.003 & 83.1\(\pm\)0.003 & 83.2\(\pm\) 0.003 & 83.5\(\pm\) 0.003 & 83.5\(\pm\) 0.003 & 83.5\(\pm\) 0.10 \\ conv-24 & **0.88\(\pm\)0.003** & **2.20\(\pm\)0.003** & **0.75\(\pm\)0.003** & 83.4\(\pm\)0.003 & 83.0\(\pm\)0.003 & 83.0\(\pm\) 0.003 & 83.2\(\pm\) 0.003 & 83.5\(\pm\) 0.003 & 83.5\(\pm\) 0.003 & 83.5\(\pm\) 0.003 \\ feature & **0.77\(\pm\)0.001** & **0.83\(\pm\)0.003** & **0.56\(\pm\)0.002** & **0.56\(\pm\)0.003** & **0.72\(\pm\)0.003** & **0.56\(\pm\)0.003** & **0.56\(\pm\)0.003** & **0.56\(\pm\)0.004** & **0.56\(\pm\)0.004** & **0.56\(\pm\)0.003** & **0.56\(\pm\)0.003** \\ feature & **0.77\(\pm\)0.001** & **0.83\(\pm\)0.003** & **0.56\(\pm\)0.002** & **0.56\(\pm\)0.001** & **0.72\(\pm\)0.003** & **0.56\(\pm\)0.003** & **0.56\(\pm\)0.003** & **0.56\(\pm\)0.004** & **0.56\(\pm\)0.004** & **0.56\(\pm\)0.003** & **0.56\(\pm\)0.003** \\ \hline Accuracy/Time & **9.13** & **11.13** & **27.13** & **17.13** & **2
## 7 Acknowledgements

This work was partially supported by the National Science Foundation (NSF) under grants ECCS-2037304, DMS-2134248, NSF CAREER award CCF-2236829, the U.S. Army Research Office Early Career Award W911NF-21-1-0242, Stanford Precourt Institute, and the ACCESS - AI Chip Center for Emerging Smart Systems, sponsored by InnoHK funding, Hong Kong SAR.
