Active Learning of Deep Neural Networks via Gradient-Free Cutting Planes

## Abstract

Abstract Active learning methods aim to improve sample complexity in machine learning. In this work, we investigate an active learning scheme via a novel gradient-free cutting-plane training method for ReLU networks of arbitrary depth.
We demonstrate, for the first time, that cutting-plane algorithms, traditionally used in linear models, can be extended to deep neural networks despite their nonconvexity and nonlinear decision boundaries. Our results demonstrate that these methods provide a promising alternative to the commonly employed gradient-based optimization techniques in large-scale neural networks.
Moreover, this training method induces the first deep active learning scheme known to achieve convergence guarantees. We exemplify the effectiveness of our proposed active learning method against popular deep active learning baselines via both synthetic data experiments and sentimental classification task on real datasets.

### 1 Introduction

Large neural network models are now core to artificial intelligence systems. After years of development, current large NN training is still dominated by gradient-based methods, which range from basic gradient descent method to more advanced online stochastic methods such as Adam (Kingma & Ba, ) and AdamW (Loshchilov & Hutter, ). Recent empirical effort has focused on cutting down storage requirement for such optimizers, see (Griewank & Walther, ; Zhao et al., ); accelerating convergence by adding momentum, see (Xie et al., ); designing better step size search algorithms, see (Defazio & Mishchenko, ). Despite its popularity, gradient-based methods suffer from sensitivity to hyperparameters and slow convergence. Therefore, researchers are persistently seeking for alternative training schemes for large NN models, including involve zero-order and second-order algorithms.

Cutting-plane method is a classic optimization algorithm and is known for its fast convergence rate. Research on cutting-plane type methods dates back to 1950s when Ralph Gomory (Gomory, ) first studied it for integer programming and mixed-integer programming problems. Since then, this method has also been heavily investigated for solving nonlinear problems. Different variations of cutting-plane method emerge, including but not limited to center of gravity cutting-plane method, maximum volume ellipsoid cutting-plane method, and analytic center cutting-plane method, which mainly differ in their center-finding strategy.

Historically, deep NN training and cutting plane methods have developed independently, each with its own audience. In this work, we bridge them for the first time by providing a viable cutting-plane-based deep NN training and active learning scheme. Our method finds optimal neural network weights and actively queries additional training points via a gradient-free cutting plane approach. We show that an active learning scheme based on our newly proposed cutting-plane-based strategy naturally inherits classic convergence guarantees of cutting-plane methods. We present synthetic and real data experiments to demonstrate the effectiveness of our proposed methods.

### 2 Notation

We denote the set of integers from 1 to $n$ as $[n]$. We use $\mathcal{B}_{p}:=\{u\in\mathbb{R}^{d}:\|u\|_{p}\leq 1\}$ to denote the unit $\ell_{p}$-norm ball and $\langle\cdot,\cdot\rangle$ to denote the dot product between two vectors. Given a hyperplane $\mathcal{H}:=\{x:x^{T}w=0\}$, we use $\mathcal{H}^{+}$ ($\mathcal{H}^{-}$) to denote the positive (negative) half-space: $\mathcal{H}^{+(-)}:=\{x:x^{T}w\geq(\leq)0\}$. For ReLU, we use notation $(x)_{+}=\max(x,0)$.
We take ${\mathbbm{1}}\{x\in S\}$ as the $1/0$-valued indicator function with respect to the set $S$, evaluated at $x$.

### 3 Preliminaries

Classic cutting-plane method’s usage in different optimization problems
has been heavily studied. To better demonstrate the problem and
offer a more self-contained background, we start with describing basic cutting-plane method’s workflow below.

###### Cutting-Plane Method.

Consider any minimization problem with an objective function $f(\theta)$, where the solution set, denoted as $\Theta$, is a convex set. Cutting-plane method typically assumes the existence of an oracle that, given any input $\theta_{0}$, either confirms that $\theta_{0}\in\Theta$, thereby terminating with $\theta_{0}$ as a satisfactory solution, or returns a pair $(x,y)$ such that $x^{T}\theta_{0}\leq y$ while $x^{T}\theta>y$ for all $\theta\in\Theta$. If the cut is “good enough,” it allows for the elimination of a large portion of the search space, enabling rapid progress toward the true solution set $\Theta$.

Figure: Figure 1: Illustration of a single iteration of general cutting-plane method (Boyd & Vandenberghe, ).
Refer to caption: extracted/5926809/img/cutcut-01.png

The classic convergence results of the cutting-plane method are highly dependent on the quality of the cut in each iteration. For instance, if the center of gravity of the current volume is removed, it guarantees a volume reduction of approximately 63%. Similar results hold for the analytic center and the center of the maximum volume ellipsoid. Figure 1 illustrates a single step of the cutting-plane method, showing how a “good” cut near the center of the current volume induces a much larger volume reduction compared to a “bad” cut near the edge.

###### Cutting-Plane-Based Active Learning (AL) with Linear Models (Louche & Ralaivola, 2015 ) .

Cutting-plane method provides a natural active learning framework to localize a set of deep NN
classifiers with certain classification margins. Prior work (Louche & Ralaivola, ) has studied the use of cutting-plane method in the context of active learning with linear models for binary classification. The setup is the following. One is given a set of unlabeled data $\{x_{1},\cdots,x_{n}\}\in\mathcal{D}$. The authors consider a linear binary classifier $f(x;\theta):=\langle\theta,x\rangle$ with prediction ${\mathrm{sign}}(f(x;\theta))$ for any input $x$. Define the set of model parameters that correctly classify our dataset as $\mathcal{T}_{\mathcal{D}}$, which is a set of linear inequalities. The size of parameter set reflects the level of uncertainty in the classifier. Starting at an initial set $\mathcal{T}^{0}$, the goal is to query additional data points and acquire their labels to reduce the size of $\mathcal{T}^{0}$ to approach $\mathcal{T}_{\mathcal{D}}$, which hopefully would have high test accuracy when the generalization error is small.
The cutting-plane-based active learning framework developed by Louche & Ralaivola () (Algorithm 3) starts with a localized convex set $\mathcal{T}^{0}$. The algorithm is presented with a set of unlabeled data. At each step $t$, it performs the following steps: (i) computes the center of the current parameter space $\theta_{c}^{t}:=\text{center}(\mathcal{T}^{t})$; (ii) queries the label $y_{t}$ for $x_{t}$ from the unlabeled dataset $\mathcal{D}$ which has minimal prediction margin with respect to $\theta_{c}^{t}$; (iii) reduces the parameter space via a cutting-plane in the case of mis-classification: $\mathcal{T}^{t+1}=\mathcal{T}^{t}\cap\{\theta\,|\,y_{t}\langle x_{t},\theta
\rangle>0\}.$ The algorithm terminates when set $\mathcal{T}^{t}$ is small enough or maximum iterations or data budget have been reached. This active learning scheme has strong convergence result inherited from classic cutting-plane method, which has been investigated in (Louche & Ralaivola, ).

###### Disentangling Model Training and Active Learning.

Although Louche & Ralaivola () primarily focuses on the active learning setting, their method implicitly suggests a cutting-plane-based training workflow for linear binary classifiers. To illustrate this, consider a set of training samples $\{(x_{1},y_{1}),\dots,(x_{n},y_{n})\}$. Each pair $(x_{i},y_{i})$ induces a cut on the parameter space $\{\theta\mid y_{i}\langle\theta,x_{i}\rangle>0\}$. The final center, i.e., $\theta_{c}=\text{center}(\{\theta\mid y_{i}\langle\theta,x_{i}\rangle>0,i\in[n
]\})$, accounts for all such cuts while maintaining the desired property ${\mathrm{sign}}(f(x;\theta))=y$ for all training samples. This choice makes $\theta_{c}$ not only an optimal solution, but also a robust choice, as it remains stable under data perturbations.

By disentangling the model training process from the active learning query strategy described in (Louche & Ralaivola, ), we derive a gradient-free cutting-plane-based workflow for training linear binary classifiers. However, this approach has several key limitations: (1) it is restricted to linear models, (2) it requires the data to be linearly separable to ensure an optimal parameter set, and (3) it only supports binary classification tasks. Our current work addresses all three limitations by extending the cutting-plane method to train deep nonlinear neural networks for both classification and regression tasks, without requiring linear separability of the data. Similar to the linear binary classification case, our vanilla training scheme lacks desirable convergence guarantees, as the cuts may occur at the edge of the parameter set. To overcome this, we focus on a cutting-plane-based active learning scheme that enables cuts to be near the center of the parameter set. We provide convergence results for this approach, which, to the best of our knowledge, is the first convergence guarantee for active learning algorithms applied to deep neural networks.

###### Outline.

The paper is organized as follows: in Section 4, we adapt the cutting-plane method for nonlinear model training by transforming the nonlinear training process into a linear programming problem. Section 5 introduces the general gradient-free cutting-plane-based training algorithm for deep NNs. In Section 6, we explore the resulting cutting-plane AL framework and prove its convergence. Section 7 demonstrates the practical effectiveness of our proposed training and active learning methods through extensive experiments. Finally, limitations and conclusions are discussed in Section 8.

### 4 Key Observation: Training ReLU NNs for Binary Classification is Linear Programming

The cutting-plane AL scheme proposed by Louche & Ralaivola () (summarized in Section 3) is designed for linear models like $f(x;\theta)=\langle\theta,x\rangle$. Extending this method to nonlinear models, such as a two-layer ReLU network $f(x;\theta)=(x^{T}\theta_{1})_{+}\theta_{2}$, with $\theta_{1}\in\mathbb{R}^{d\times m}$ and $\theta_{2}\in\mathbb{R}^{m}$, presents additional challenges. Specifically, for a mispredicted data pair $(x_{i},y_{i})$, determining how to cut the parameter space $(\theta_{1},\theta_{2})$ is far more complex, whereas in the linear case, the cut is simply $y_{i}\langle\theta,x_{i}\rangle>0$.

To break this bottleneck and extend the cutting-plane-based learning method to nonlinear models, we observe that training a ReLU network for binary classification can be formulated as a linear programming problem. This insight is crucial for extending the learning scheme in (Louche & Ralaivola, ) to more complex models. We now develop our core idea of reframing binary classification with ReLU models as linear programs. For clarity, we present our results in two theorems: one for two-layer ReLU networks and another for ReLU networks of arbitrary depth. We focus on the two-layer case in the main paper for detailed discussion, deferring the more abstract general case to Appendix E.2. Since the general case is an extension of the two-layer model, focusing on the two-layer case should provide a clearer understanding of the core concepts.

We start with writing the linear program corresponding to the linear model for binary classification tasks as below,

$$ $\displaystyle\underset{}{\text{find}}$ $\displaystyle\theta$ (1) s.t. $\displaystyle y_{i}\langle\theta,x_{i}\rangle\geq 1\quad\forall i.$ $$

Note that intuitively, we want $y_{i}\langle\theta,x_{i}\rangle>0$ to be satisfied for our sign prediction. However, the set $\{\theta|y_{i}\langle\theta,x_{i}\rangle>0\}$ is not compact and is thus not compliant with forms of standard linear programs. This may raise technical issues. We observe that our training data is finite, and thus we can always scale $\theta$ to achieve $y_{i}\langle\theta,x_{i}\rangle\geq c$ for any positive constant $c$ once $y_{i}\langle\theta,x_{i}\rangle>0$ holds, and $y_{i}\langle\theta,x_{i}\rangle\geq c$ for $c>0$ also guarantees $y_{i}\langle\theta,x_{i}\rangle>0$. We pick $c=1$ in (1). With a two-layer ReLU model, we obtain the following problem:

$$ $\displaystyle\underset{}{\text{find}}$ $\displaystyle W_{1},W_{2}$ (2) s.t. $\displaystyle y_{i}(x_{i}^{T}W_{1})_{+}W_{2}\geq 1\quad\forall i.$ $$

Before showing that solving (2) is indeed equivalent to solving a linear program, we first introduce the core concept of activation patterns which we will draw on heavily later. For data matrix $X\in\mathbb{R}^{n\times d}$ and any arbitrary vector $u\in\mathbb{R}^{d}$, we consider the set of diagonal matrices

$$ $\mathcal{D}:=\{\text{diag}({\mathbbm{1}}\{Xu\geq 0\})\}.$ $$

We denote the carnality of set $\mathcal{D}$ as $P$, i.e., $P=|\mathcal{D}|.$ Thus $\{D_{i}\in\mathcal{D},i\in[P]\}$ iterates over all possible activation patterns of ReLU function induced by data matrix $X$. See Definition 3 for more details. With this concept of activation patterns, we can reframe the training of two-layer ReLU model for binary classification as the following linear program:

###### Theorem 4.1 .

When $m\geq 2P$, Problem (2) is equivalent to

$$ $\displaystyle\underset{}{\text{find}}$ $\displaystyle u_{i},u_{i}^{\prime}$ (3) s.t. $\displaystyle y\left(\sum_{i=1}^{P}D_{i}X(u_{i}-u_{i}^{\prime})\right)\geq 1,$ $\displaystyle(2D_{i}-I)Xu_{i}\geq 0,$ $\displaystyle(2D_{i}-I)Xu_{i}^{\prime}\geq 0.$ $$

###### Proof.

See Appendix E.1.
∎

The high-level rationality behind Theorem 4.1 is the observation that

$$ $\displaystyle(XW_{1})_{+}W_{2}$ $\displaystyle=\sum_{i=1}^{m}\left(XW_{1i}\right)_{+}W_{2i}=\sum_{i=1}^{m}\text {diag}({\mathbbm{1}}\{XW_{1i}\geq 0\})XW_{1i}W_{2i}.$ (4) $$

By defining ${K}_{i}:=\text{diag}({\mathbbm{1}}\{XW_{1i}\geq 0\})$ and a set of $\{v_{i},v_{i}^{\prime}\}$ vectors by setting $v_{i}:=W_{1i}W_{2i}$ when $W_{2i}\geq 0$ and 0 0 otherwise, $v_{i}^{\prime}:=-W_{1i}W_{2i}$ when $W_{2i}<0$ and 0 0 otherwise, we have $\text{diag}({\mathbbm{1}}\{XW_{1i}\geq 0\})XW_{1i}W_{2i}={K}_{i}X(v_{i}-v_{i}^
{\prime})$. The expression in (4) thus writes

$$ $(XW_{1})_{+}W_{2}=\sum_{i=1}^{m}{K}_{i}X(v_{i}-v_{i}^{\prime})\text{ where }(2 {K}_{i}-I)Xv_{i}\geq 0,(2{K}_{i}-I)Xv_{i}^{\prime}\geq 0.$ $$

Therefore, if we iterate over all $D_{i}\in\mathcal{D}$, we are guaranteed to reach each ${K}_{i}$. The equivalence in Theorem 4.1 holds in the sense that, which we prove rigorously in Appendix E.1, whenever there is solution $W_{1},W_{2}$ to Problem (2), there is always solution $\{u_{i},u_{i}^{\prime}\}$ to problem (3) and vice versa. Moreover, when an optimal solution $\{u_{i}^{\star},u_{i}^{\prime\star}\}$ to problem (3) has been found, we can explicitly create an optimal solution $\{W_{1}^{\star},W_{2}^{\star}\}$ to Problem (2) from value of $\{u_{i}^{\star},u_{i}^{\prime\star}\}$, see Appendix E.1 for details. Thereafter, for any test point $\tilde{x}$, our sign prediction is simply ${\mathrm{sign}}((\tilde{x}W_{1}^{\star})_{+}W_{2}^{\star})$, which will have the same value as ${\mathrm{sign}}(\sum_{i}(\tilde{x}^{T}u_{i}^{\star})_{+}-(\tilde{x}^{T}u_{i}^{
\prime\star})_{+}).$

We emphasize that our reframing of training a ReLU network for binary classification as linear programming does not eliminate the nonlinearity of the ReLU activation. Instead, this approach works because the ReLU activation patterns for a given training dataset are finite. By looping through these activation patterns, we can explicitly enumerate them. At test time, the ReLU nonlinearity is preserved, as our prediction ${\mathrm{sign}}(\sum_{i}(\tilde{x}^{T}u_{i})_{+}-(\tilde{x}^{T}u_{i}^{\prime})
_{+})$ depends on the sign of $\tilde{x}^{T}u_{i}$ and $\tilde{x}^{T}u_{i}^{\prime}$, ensuring that the expressiveness of the nonlinearity remains intact. A careful reader might note that the number of patterns $P$ increases with the size of the training data, meaning the number of variables in Problem (3) may also grow. Additionally, finding all activation patterns poses a challenge. In Section 7, we demonstrate that subsampling a set of non-duplicate activation patterns performs well in practice. For further grounding, Appendix E.4 outlines an iterative hyperplane filtering method that guarantees the identification of all activation patterns with a reasonable complexity bound.

Now, let us consider ReLU network with $n$ hidden layer for binary classification task

$$ $\displaystyle\underset{}{\text{find}}$ $\displaystyle W_{1},W_{2},\cdots,W_{n+1}$ (5) s.t. $\displaystyle y\odot((\cdots(((XW_{1})_{+}W_{2})_{+}W_{3})_{+}W_{4}\cdots)_{+} W_{n})_{+}W_{n+1}\geq 1.$ $$

We extend the activation patterns involved in Theorem 4.1 to $(n+1)$-layer neural networks. Let $W_{1}\in\mathbb{R}^{d\times m_{1}},W_{2}\in\mathbb{R}^{m_{1}\times m_{2}},W_{3
}\in\mathbb{R}^{m_{2}\times m_{3}},\cdots,W_{n}\in\mathbb{R}^{m_{n-1}\times m_
{n}},W_{n+1}\in\mathbb{R}^{m_{n}}.$ Define $m_{0}:=d$ and the activation pattern in $i$-th layer as

$$ $\mathcal{D}^{(i)}:=\left\{\text{diag}({\mathbbm{1}}\{\left(\cdots\left((Xv_{1} )_{+}v_{2}\right)_{+}v_{3}\cdots\right)_{+}v_{i}\geq 0\})\bigg{|}v_{j}\in \mathbb{R}^{m_{j-1}\times m_{j}}~{}\forall~{}j<i,v_{i}\in\mathbb{R}^{m_{i-1}} \right\}.$ $$

We denote the cardinality of set $\mathcal{D}^{(i)}$ as $P_{i}$, i.e., $P_{i}=|\mathcal{D}^{(i)}|.$ Thus $\{D_{j}^{(i)}\in\mathcal{D}^{(i)},j\in[P_{i}]\}$ iterates over all possible activation patterns at the $i$-th hidden layer. We then reframe Problem (5) as below:

###### Theorem 4.2 .

When $m_{i}\geq\Pi_{i}^{n}2P_{i}$ for each $i\in[n]$, Problem (5) is equivalent to

$$ $\displaystyle\underset{}{\text{find}}$ $\displaystyle u_{j_{n}j_{n-1}j_{n-2}\ldots j_{1}}^{c_{n}c_{n-1}c_{n-2}\cdots c _{1}}$ (6) s.t. $\displaystyle y\odot\sum_{j_{n}=1}^{P_{n}}D_{j_{n}}^{(n)}\left(\mathcal{T}_{1} ^{(n-1)}(D^{(n-1)})-\mathcal{T}_{2}^{(n-1)}(D^{(n-1)})\right)\geq 1$ $\displaystyle(2D_{j_{i}}^{(i)}-I)\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{i-1}}^{( n-1)(n-2)\cdots(i-1)}(D^{(i-1)})\geq 0,2\leq i\leq n{}$ $\displaystyle(2D_{j_{1}}^{(1)}-I)Xu_{j_{n}j_{n-1}\cdots j_{1}}^{c_{n}c_{n-1} \cdots c_{1}}\geq 0,{}$ $$

where $c_{i}\in\{1,2\}$ and

$$ $\displaystyle\mathcal{T}_{c_{n-1}\cdots c_{i}}^{(n-1)\cdots(i)}(D^{(i)})=\sum_ {j_{i}=1}^{P_{i}}D_{j_{i}}^{(i)}\left(\mathcal{T}_{c_{n-1}\cdots c_{i}1}^{(n-1 )\cdots(i)(i-1)}(D^{(i-1)})-\mathcal{T}_{c_{n-1}\cdots c_{i}2}^{(n-1)\cdots(i) (i-1)}(D^{(i-1)})\right),\forall~{}i\leq n-1,$ $\displaystyle\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{1}}^{(n-1)(n-2)\cdots(1)}(D^ {(1)})=\sum_{j_{1}=1}^{P_{1}}D_{j_{1}}^{(1)}X\left(u_{j_{n}j_{n-1}\cdots j_{1} }^{1c_{n-1}c_{n-2}\cdots c_{1}}-{u_{j_{n}j_{n-1}\cdots j_{1}}^{2c_{n-1}c_{n-2} \cdots c_{1}}}\right).$ $$

###### Proof.

See Appendix E.2.
∎

###### Theorem 4.1 .

###### Proof.

###### Theorem 4.2 .

###### Proof.

### 5 Training Deep Neural Networks via Cutting-Planes

Figure: Algorithm 1 Training NN with Cutting Plane Method

With the linear program reframing of training deep ReLU models in place, we now formally introduce our cutting-plane-based NN training scheme for binary classification. We begin with a feasible set of variables in our linear program (6) that contains the optimal solution. For each training sample $(x_{i},y_{i})$, we add the corresponding constraints as cuts. At each iteration, we select the center of the current parameter set, stopping when either the validation loss stabilizes or after a fixed number of iterations, similar to the stopping criteria in gradient-based training of large models.

We present our cutting-plane-based NN training algorithm here in main text. For generalizations, such as (i) relaxing the data distribution to remove the linear separability requirement (Appendix D.1), and (ii) extending from classification to regression (Appendix D.2), we refer readers to Appendix D. We emphasize that both the relaxed data constraint and the ability to handle regression tasks are unique to our method and have not been achieved in prior work.

Algorithm 1 presents the general workflow of how we train deep NNs with gradient-free cutting-plane method, which simply adds a cut corresponding to each training data $(x_{n_{t}},y_{n_{t}})$. When the stopping criterion is satisfied, center of current parameter set $\mathcal{T}^{t}$ is returned. Here we start with the parameter space as the unit $2$-norm ball, which is guaranteed to contain optimal parameters due to scale invariance.

After we get the final $\theta_{c}^{t}$, for any test point $\tilde{x}$, we can directly compute our sign prediction with $\theta_{c}^{t}$. For example, for two-layer ReLU model, $\theta_{c}^{t}$ will be of form $\{u_{i}^{t},u_{i}^{\prime t}\}$, our final sign prediction would simply be ${\mathrm{sign}}(\sum_{i}(\tilde{x}^{T}u_{i}^{t})_{+}-(\tilde{x}^{T}u_{i}^{
\prime t})_{+}).$ For deeper models, the prediction is a bit more complex and is given in equation (17) in Appendix E.2. Notably, the computation of final sign prediction from value of $\theta_{c}^{t}$ always takes a single step, just as one forward pass of original NN model formulation. Moreover, one can also restore optimal NN weights from final $\theta_{c}^{t}$, see Appendix E.1 for two-layer case reconstruction of optimal NN parameters and Appendix E.2 for general deep NN models.

Two key functions in our proposed training scheme are the “center” function and the “cut” function, which we detail below:

- •
Center. The “center” function calculates the center of the convex set $\mathcal{T}^{t}$. There are a couple of notions of centers, such as center of gravity (CG), center of maximum volume ellipsoid (MVE), Chebyshev’s center, and analytic center (Boyd & Vandenberghe ()). Among these, the analytic center is the easiest to compute and is empirically known to be effective (Goffin et al. (), Atlason et al. ()). This is the notion of center that we will adopt to compute the query point in our algorithm. See Appendix G.1 for details.
- •
Cut. The “cut” function determines the cutting planes we get from a specific training data $(x_{n_{t}},y_{n_{t}}).$ For two-layer model, “cut” function would return the constraint set $\{y_{n_{t}}(\sum_{i}D_{in_{t}}x_{n_{t}}^{T}(u_{i}-u_{i}^{\prime}))\geq 1,(2D_{
in_{t}}-1)x_{n_{t}}^{T}u_{i}\geq 0,(2D_{in_{t}}-1)x_{n_{t}}^{T}u_{i}^{\prime}
\geq 0\}$. For deeper NNs, “cut” function would return constraints listed in Problem (6).

Compared to gradient-based NN training scheme, we take a cutting-plane cut for each data point encountered while gradient-based method employs a gradient descent step corresponding to the data query. Moreover, our training scheme is guaranteed to correctly classify all data points we have ever encountered, while gradient-based method has no such guarantees.

### 6 Cutting-Plane-Based Active Learning and Convergence Guarantees

#### 6.1 Algorithm II: Cutting-Plane Localization for Active Learning

Our proposed cutting-plane-based active learning algorithm adapts and extends the generic framework discussed in Section 5 (Algorithm 3). For the sake of simplicity, we present in this section the algorithm specifically for binary classification and with respect to two-layer ReLU NNs. We emphasize that the algorithm can be easily adapted to the case of multi-class and regression, per discussions in Appendix D.2, and for deeper NNs following our reformulation in Theorem 4.2.

Recall the problem formulation for cutting-plane AL with two-layer ReLU NN for binary classification in Equation (3).
Given a training dataset $\mathcal{D}$, we use $X_{\mathcal{D}}$ and $y_{\mathcal{D}}$ to denote the slices of $X$ and $y$ at indices $\mathcal{D}$. Moreover, we succinctly denote the prediction function as:

$$ $\displaystyle\begin{split}f^{\text{two-layer}}(X;\theta):=\sum_{i=1}^{P}(D_{i} X)_{\mathcal{D}}({u^{\prime}}_{i}-u_{i})=\begin{bmatrix}X_{\mathcal{D}}^{1}&-X _{\mathcal{D}}^{1}&\ldots&X_{\mathcal{D}}^{P}&-X_{\mathcal{D}}^{P}\end{bmatrix }\theta,\end{split}$ (7) $$

where $\theta=(u_{1}^{\prime},u_{1},\ldots,u_{P}^{\prime},u_{P})$ with $u_{i},{u^{\prime}}_{i}\in\mathbb{R}^{d}$, and $X_{\mathcal{D}}^{i}$ is a shorthand notation for $(D_{i}X)_{\mathcal{D}}$.
For the further brevity of notation, we denote the ReLU constraints in Equation 3, i.e.
$\left((2D_{i}-I_{n})X\right)_{\mathcal{D}}u_{i}\geq 0,\ \left((2D_{i}-I_{n})X
\right)_{\mathcal{D}}{u^{\prime}}_{i}\geq 0$ for all $i$,
as $C(\mathcal{D}),C^{\prime}(\mathcal{D})$.

With Theorem 4.1 and the linearization of $f^{\text{two-layer}}(X;\theta)$, cutting-plane-based active learning methods become well applicable.
As in Algorithm 3, we restrict the parameter space $\Theta$ to be within the unit $2$-norm ball: $\mathcal{B}_{2}:=\{\theta\,|\,\|\theta\|_{2}\leq 1\}=\{\theta|\sum_{i=1}^{P}(
\|u_{i}\|^{2}_{2}+\|{u^{\prime}}_{i}\|^{2}_{2})\leq 1\}$. For computing the center of the parameter space at each step for queries, we use the analytic center (Definition 1), which is known to be easily computable and has good convergence properties. We refer to Section 6.2 for a more detailed discussion.

###### Definition 1 (Analytic center) .

The analytic center of polyhedron $\mathcal{P}=\{z|a_{i}^{T}z\leq b_{i},i=1,...,m\}$ is given by

$$ $\displaystyle\text{AC}(\mathcal{P}):=\operatorname*{arg\,min}_{z}-\sum_{i=1}^{ m}\log(b_{i}-a_{i}^{T}z)$ (8) $$

We are now ready to present the cutting-plane-based active learning algorithm for deep NNs. For breadth of discussion, we present three versions of the active learning algorithms, each corresponding to the following setups:

- 1.
Cutting-plane AL with query synthesis (Algorithm 2). The cutting-plane oracle gains access to a query synthesis. Therefore, the cut is always active until we have encountered the optimal classifier(s), at which point the algorithm terminates.
- 2.
Cutting-plane AL with limited queries (Algorithm 4). The cutting-plane oracle has access to limited queries. The cut is only performed when the queried candidate mis-classifies the data pair returned by the oracle.
- 3.
Cutting-plane AL with inexact cuts (Algorithm 5). The cutting-plane oracle has access to limited queries. However, the algorithm always performs the cut regardless of whether the queried candidate mis-classifies the data pair returned by the oracle.

Figure: Algorithm 2 Cutting-plane AL for Binary Classification with Query Synthesis

For brevity, we present the algorithm for the first setup here and refer the rest to Appendix C. Algorithm 2 summarizes the proposed algorithm under the first setup, where we have used $\mathcal{D}_{\text{QS}}$ to denote the query synthesis. We note that Algorithm 4 for the second setup is obtained simply by changing $\mathcal{D}_{\text{QS}}$ to $\mathcal{D}_{\text{LQ}}$, which denotes limited query.

The general workflow of Algorithm 2 follows that of Algorithm 3 but on the transformed parameter space via the ReLU networks. We highlight here two key differences in our algorithm: (1) in each iteration, for faster empirical convergence, we query twice for the data point that is classified positively and the one that is classified negatively with highest confidence, respectively.
If one or both of them turn out to be miss-classifications, this informs the active learner well and we expect a large cut.
(2) For two- and three-layer ReLU networks, Ergen & Pilanci () demonstrate that these models can be reformulated as exact convex programs. This allows the option of incorporating a final convex solver into Algorithm 2, applied after the active learning loop with the data collected thus far. This convex reformulation includes regularization, which can improve the performance of our cutting-plane AL in certain tasks. For more details, see Appendix G.2

###### Definition 1 (Analytic center) .

#### 6.2 Convergence Guarantees

We give theoretical examination of the convergence properties of Algorithm 2 and 4 with respect to both the center of gravity (cg) and the center of maximum volume ellipsoid (MVE). Analysis of Algorithm 5 for inexact cuts is given in Appendix F.2. As the analysis of MVE closely parallels that of CG, we refer readers to Appendix F.3 for a detailed discussion, in the interest of brevity. For both centers, we measure the convergence speed with respect to the volume of the localization set $\mathcal{T}^{t}$ and judge the progress in iteration $t$ by the fractional decrease in volume: $\textbf{vol}(\mathcal{T}^{t+1})/\textbf{vol}(\mathcal{T}^{t}).$

To start, we give the definitions of the center of gravity (Boyd & Vandenberghe, ).

###### Definition 2 (Center of gravity (CG)) .

For a given convex body (i.e. a compact convex set with non-empty interior) $C\subseteq\mathbb{R}^{d}$, the centroid, or center of gravity of $C$, denoted $\theta_{G}(C)$, is given by

$$ $\displaystyle\theta_{G}(C)=\frac{1}{\textbf{vol}(C)}\int_{x\in C}xdx.$ $$

Given convex set $\mathcal{C}$, we use the abbreviated notation $\theta_{G}$.

Our analysis on the center of gravity relies on an important proposition (Proposition 1) given by Grünbaum (), which guarantees that in each step of cutting via a hyperplane passing through the centroid of the convex body, a fixed portion of the feasible set is eliminated. A recursive application of this proposition shows that after $t$ steps from the initial step, we obtain the following volume inequality: $\textbf{vol}(\mathcal{T}_{t})\leq(1-1/e)^{t}\textbf{vol}(\mathcal{T}_{0})
\approx(0.63)^{t}\textbf{vol}(\mathcal{T}_{0}).$

Observe that our proposed cutting-plane-based active learning method in Algorithm 2 and 4 uses a modified splitting, where the weight vector in Proposition 1 is substituted by a mapping of the parameters $\theta$ to the feature space via function $f^{\text{two-layer}}$, which depends on point $x_{n_{t}}$ returned by the oracle at the step, along with the associated linear constraints $\mathcal{C}(\{n_{t}\})$, $\mathcal{C}^{\prime}(\{n_{t}\})$. Since $f^{\text{two-layer}}$ is linear in $\theta$ as we recall that

$$ $\displaystyle f^{\text{two-layer}}(x_{n_{t}};\theta)=[x_{n_{t}}^{1}-x_{n_{t}}^ {1}\ ...\ x_{n_{t}}^{P}-x_{n_{t}}^{P}]\theta,$ $$

the set defined by $\{\theta|y_{n_{t}}\cdot f^{\text{two-layer}}(x_{n_{t}};\theta)\geq 0\}$ forms a half-space in the parameter space.
Additionally, since the constraints $\mathcal{C}(\{n_{t}\})$, $\mathcal{C}^{\prime}(\{n_{t}\})$ are linear in $\theta$, the cutting set $\{\theta:y_{n_{t}}\cdot f^{\text{two-layer}}(x_{n_{t}};\theta)\geq 0,\mathcal{
C}(\{n_{t}\}),\mathcal{C}^{\prime}(\{n_{t}\})\}$ in Algorithm 2 and 4 defines a convex polyhedron. This change suggests a non-trivial modification of the results given in Proposition 1. The following theorem is our contribution.

###### Theorem 6.1 (Convergence with Center of Gravity) .

Let $\mathcal{T}\subseteq\mathbb{R}^{d}$ be a convex body and let $\theta_{G}$ denote its center of gravity. The polyhedron cut given in Algorithm 2
and Algorithm 4 (assuming that the cut is active), i.e.,

$$ $\displaystyle\mathcal{T}\cap\{\theta:y_{n}\cdot f^{\text{two-layer}}(x_{n}; \theta)\geq 0,\mathcal{C}(\{n\}),\mathcal{C}^{\prime}(\{n\})\},$ $$

where coupling $(x_{n},y_{n})$ is the data point returned by the cutting-plane oracle after receiving queried point $\theta_{G}$, partitions the convex body $\mathcal{T}$ into two subsets:

$$ $\displaystyle\mathcal{T}_{1}:=\{\theta\in\mathcal{T}:y_{n}\cdot f^{\text{two- layer}}(x_{n};\theta)\geq 0,\mathcal{C}(\{n\}),\mathcal{C}^{\prime}(\{n\})\}$ $\displaystyle\mathcal{T}_{2}:=\{\theta\in\mathcal{T}:y_{n}\cdot f^{\text{two- layer}}(x_{n};\theta)<0,\text{ or }\neg\mathcal{C}(\{n\}),\neg\mathcal{C}^{ \prime}(\{n\})\},$ $$

where $\neg$ denotes the complement of a given set. Then $\mathcal{T}_{1}$ satisfies the following inequality:

$$ $\displaystyle\textbf{vol}(\mathcal{T}_{1})<\left(1-\frac{1}{e}\right)\cdot \textbf{vol}(\mathcal{T}).$ $$

###### Proof.

See Appendix F.1.
∎

###### Definition 2 (Center of gravity (CG)) .

###### Theorem 6.1 (Convergence with Center of Gravity) .

###### Proof.

### 7 Experiments

We validate our proposed training and active learning methods through extensive experiments, comparing them with various popular baselines from scikit-activeml (Kottke et al., ) and DeepAL (Huang, ). Synthetic data experiments are presented in Section 7.1, and real data experiments are presented in Section 7.2. An overview of each baseline is given in Appendix G.3. For implementation details and additional results, refer to Appendix G.

Figure: Figure 2: Decision boundaries for binary classification on the spiral dataset for the cutting-plane AL method using a two-layer ReLU neural network, alongside various deep AL baselines. For compactness, we also include the decision boundaries for the cutting-plane AL method with a three-layer ReLU network in the collage to demonstrate its feasibility. For fairness of comparison, we use the same two-layer ReLU network structure and embedding size of 623 for all methods. We enforce the same hyperparameters for all deep AL baselines and select the best performing number of training epochs at $2000$ and a learning rate at $0.001$ to ensure optimal performance. See Appendix G.3 for details.
Refer to caption: extracted/5926809/cvx_plots/spiral_update2-01.png

#### 7.1 Synthetic Data Experiments

In this section, we present small scale numerical experiments to verify the performance of our algorithm on both classification and regression tasks.

###### Binary Classification on Synthetic Spiral.

We use a synthetic dataset of two intertwined spirals with positive and negative labels, respectively (see generation details in Appendix G.4).
We generate 100 spiral data points with a 4:1 train-test split. Table 3 in Appendix G.4 presents the train and test accuracy of our cutting-plane AL (Algorithm 4) compared to popular deep AL baselines, with all methods evaluated with a query budget of 20 points (25% of train data). Our method achieved perfect accuracy on both sets, outperforming all baselines.
Notably, the strong performance of our cutting-plane AL extends to the 3-layer case, achieving train/test accuracies of 0.71/0.60, while using only a fraction of neurons per layer (57 and 34, resp.) compared to the two-layer cutting-plane AL, which used 623 neurons. This result is illustrated in the corresponding decision boundary plot in Figure 2, where the 3-layer cutting-plane AL is one of the few methods to capture the spiral’s rough shape despite using smaller embeddings, while the two-layer cutting-plane AL, with the same network structure as the baselines, precisely traces the spiral. We emphasize that the superior performance of our cutting-plane AL remains consistent across different random seeds. As shown in the error-bar plot in Figure 15, our approach reliably converges to the optimal classifier faster than all the tested baselines in the number of queries.

###### Quadratic Regression.

We evaluate the performance of our method (Algorithm 6) against the same seven baselines used in the spiral task, along with two additional popular regression AL methods from scikit-activeml: greedy sampling in target space (GreedyT) and KL divergence maximization (kldiv), on a simple quadratic regression task. We generate 100 noise-free data points from the function $y=x^{2}$ and apply the AL methods on a 4:1 train-test split with a query budget of 20 points. Table 5 in Appendix G.5 shows the root mean square error (RMSE) (see Definition 5), and left of Figure 3 provides a visualization of the final predictions of our method compared with a selection of baselines. We refer readers to Figure 18 for the complete version. Our cutting-plane AL achieves the lowest train/test RMSE of 0.01/0.01, representing an over 80%/75% reduction compared to the next best-performing baseline.This superior performance remains consistent across random seeds, as demonstrated in the train/test error-bar plots right of Figure 3, where our method consistently converges faster to the optimal classifier in terms of the number of queries.

Figure: Figure 3: Left: Predictions for the quadratic regression task using the cutting-plane AL method with a two-layer ReLU neural network, alongside various deep AL baselines. For brevity, we present only the most representative examples here and refer to Figure 18 in Appendix G.5 for the full result. Implementation details can be found in Appendix G.5. We note that the linear cutting-plane AL for regression (Algorithm 7) becomes infeasible when solving for the next center after the fourth query. This failure is expected, see Appendix C. Thus the reported prediction for this method is based on 4 queries, while all other methods use 20 queries. Right: Logarithm of mean test/train RMSE across seeds (0-4) versus the number of queries for the two-layer cutting-plane AL and various baselines. The linear cutting-plane method is omitted from the plot due to its infeasibilty.
Refer to caption: extracted/5926809/cvx_plots/reg_all_all.png

#### 7.2 Sentiment Classification for Real Datasets Using LLM Embeddings

To demonstrate the viability of our method when applied to tackle real life tasks, we also explore the concatenation of microsoft Phi-2 (Javaheripi et al., ) model with our two-layer ReLU binary classifier for sentiment classification task on IMDB (Maas et al., ) movie review datasets. Specifically, the dataset contains $50k$ movie reviews collected online with each comment accompanied with binary labels where $+1$ denoting a positive review and $-1$ denoting a negative review. See Appendix G.6 for several training data examples. We test our proposed active learning algorithm conjuncted with our cutting-plane training scheme as described in Section 5 and Section 6 respectively. For our experiment, we first collect all last layer Phi-2 embeddings (corresponding to last token), which is of size $d=2560$, for our training and testing movie reviews as our feature vectors. We follow the implementation details in Appendix G.1 and sample $P=500$ activation patterns to cater for the large dimension feature space.

Figure 4 demonstrates our experimental results involving our method and various baselines. The left most figure compares the classification accuracy between our two-layer ReLU model and the linear model studied in (Louche & Ralaivola, ). As it can be seen, our model achieves higher accuracy within same query budget compared to the linear classifier, which demonstrates the effectiveness of using a nonlinear model in this task. The middle plot compares our active learning method and the random sampling method, both are with two-layer ReLU model. As expected, our active sampling scheme identifies important data points in each iteration, which helps to train a better network within same query budget. The right most plot compares between our newly-introduced cutting-plane-based NN training scheme and classic stochastic gradient descent (SGD) method. For SGD baselines, we take one gradient step corresponding to each data query. For NN trained with SGD, people usually use batched data for gradient computation, thus within as few as $15$ query points, it is reasonable that SGD does not make good progress. On the contrary, our cutting-plane training scheme achieves higher accuracy with this tiny training budget.

Figure: Figure 4: Sentiment analysis on IMDB movie review dataset with two-layer ReLU model. We take Phi-2 embedding as our training features and compare with various baselines. The result shows that the introduction of non-linearity improves upon linear model performance, our active sampling scheme effectively identifies valuable training points compared to random sampling, and our cutting-plane training scheme is more effective than SGD in this setting. See Section 7.2 for details. Linear and our reframed two-layer models are initialized to predict zero while two-layer NN trained with SGD has random weight initialization, thus starting from non-zero prediction.
Refer to caption: extracted/5926809/cvx_plots/imdb-222-01.png

### 8 Conclusion and Limitations

In this work, we introduce a novel cutting-plane-based method for deep neural network training, which, for the first time, enables the application of this approach to nonlinear models. Additionally, our new training scheme removes previous restrictions on the training data distribution and extends the method beyond binary classification to general regression tasks. We also explore an active learning scheme built on our proposed training framework, which inherits convergence guarantees from classic cutting-plane methods. Through both synthetic and real data experiments, we demonstrate the practicality and effectiveness of our training and active learning methods. In summary, our work introduces a novel, gradient-free approach to neural network training, demonstrating for the first time the feasibility of applying the cutting-plane method to neural networks, while also offering the first deep active learning method with convergence guarantees.

Despite its novelty, our current implementation has several key limitations that hinder its competitiveness with large-scale models trained using gradient-based methods. First, although we employ subsampling of activation patterns and propose an iterative filtering scheme (see Appendix E.4), the subsampling process is not exhaustive, which impacts model performance, especially with high-dimensional data. Refining the activation pattern sampling strategy could significantly improve results. Second, we rely on analytic center retrieval during training, which we solve using CVXPY. However, this solution is CPU-bound and becomes inefficient for large-scale problems with many variables. Developing a center-finding algorithm that leverages GPU parallelism is crucial to unlocking the full potential of our training method. Finally, while current large language model (LLM) training often involves cross-entropy loss, our approach has so far been applied only to classification and regression tasks. Extending our method to handle more diverse loss functions presents an exciting avenue for future research.

### 9 Acknowledgments

This work was supported in part by National Science Foundation (NSF) under Grant DMS-2134248; in part by the NSF CAREER Award under Grant CCF-2236829; in part by the U.S. Army Research Office Early Career Award under Grant W911NF-21-1-0242; in part by the Office of Naval Research under Grant N00014-24-1-2164.

In addition, we would like to thank Maximilian Schaller from the Department of Electrical Engineering at Stanford for their valuable contributions to the development of the code used in this paper.

## Appendix

### Appendix A Prior Work

We emphasize the novelty of our work as the first to introduce both cutting-plane-based training schemes and cutting-plane-based active learning methods to deep neural network models. To the best of our knowledge, this work also provides the first convergence guarantees among existing active learning methods for nonlinear neural networks. As discussed in Section 3, the most closely related work is (Louche & Ralaivola, ), which also integrates the concepts of cutting-plane methods, active learning, and machine learning for binary classification tasks. We are not aware of any other studies that closely align with our approach. However, we review additional works that, while less directly connected, share some overlap with our research in key areas.

###### Cutting-Plane Method.

Cutting-plane methods are first introduced by Gomory () for linear programming. Though being considered as ineffective after it was first introduced, it has been later shown by Balas et al. () to be empirically useful when combined with branch-and-bound methods. Cutting-plane method is now heavily used in different commercial MILP solvers. Commonly employed cutting planes for solving convex programs are tightly connected to gradient information, which is also the underline logic for well-known Kelley’s cutting-plane method (Kelley, ). Specifically, consider any minimization problem with objective $f(\theta)$ whose solution set, we denote as $\Theta$, is a convex set. Gradient-based cutting-plane method usually assumes the existence of an oracle such that given any input $\theta_{0}$, it either accepts $\theta_{0}\in\Theta$ - and thus terminates with a satisfactory solution $\theta_{0}$ being found - or it returns a pair $(x,y)$ such that $x^{T}\theta_{0}\leq y$ while $x^{T}\theta>y$ for any $\theta\in\Theta$, i.e., we receive a cutting plane that cuts between the current input $\theta_{0}$ and the desired solution set $\Theta$. Such cut is given by subgradient of $f$ at query point $\theta_{0}$. Consider any subgradient $g\in\nabla f(\theta_{0})$, we have the inequality $f(\theta)\geq f(\theta_{0})+g^{T}(\theta-\theta_{0}),$ which then raises the cut $g^{T}(\theta-\theta_{0})\leq 0$ for minimization problem. Recent development of cutting-plane methods involves those designed for convex-concave games (Jiang et al., ), combinatorial optimization (Lee et al., ), and also application to traditional machine learning tasks such as regularized risk minimization, multiple
kernel learning, and MAP inference in graphical models (Franc et al., ).

###### Active Learning.

With the fast growth of current model sizes, a large amounts of data is need for training an effective deep networks. Compared to just feeding the model with a set of randomly-selected training samples, how to select the most informative data for each training iteration becomes a critical question, especially in neural models (Zhu & Nowak, ; Karzand & Nowak, ). Efficient data sampling is thus increasingly important, especially for current RL-based language model training schemes. With the increasing power of current pretrained models, researchers observed that since most of simple questions can be addressed correctly by the model already, these questions are less important for further boosting the language models’ capacity. Therefore, actively selecting questions that LLM cannot address correctly is key to current LLM training (Lightman et al., ). Based on different learning settings, active learning strategies can be divided into: stream-based selective sampling used when the data is generated continuously (Woodward & Finn, ); pool-based sampling used when a pool of unlabeled data is presented (Gal et al., ); query synthesis methods used when new samples can be generated for labeling. Based on different information measuring schemes, active learning algorithms can be further divided into uncertainty sampling (Lewis & Gale, ) which selects data samples to reduce prediction uncertainty, query-by-committee sampling which involves multiple models for data selection, diversity-weighted method which selects the most diverse data sample, and expected error reduction method (Mussmann et al., ) which selects samples to best reduce models’ expected prediction error.

###### Convex NN.

We note that the idea of introducing hyperplane arrangements to derive equivalent problem formulation for deep NN training task has also been investigated in prior convexification of neural network research. For example, Ergen & Pilanci () has exploited this technique to derive a convex program which is equivalent to two-layer ReLU model training task. More developments involving (Ergen & Pilanci, ) which derives convex programs equivalent to training three-layer CNNs and (Zhang & Pilanci, ) which derives convex programs for diffusion models. However, those work neither derive linear programs as considered in our case, nor did they connect such reframed problems to active learning platform. On the contrast, they mainly focus on solving NN training problem by solving the equivalent convex program (directly) they have derived via convex program solver such as CVXPY (Diamond & Boyd, ).

### Appendix B Key Definitions and Deferred Theorems

#### B.1 Key Definition

A notion central to the linear programming reformulation which enables the feasibility of our proposed cutting-plane based AL method is the notion of hyperplane arrangement. It has been briefly introduced in Section 4. We now give a formal definition.

###### Definition 3 (Hyperplane Arrangement) .

A hyperplane arrangement for a dataset $X\in\mathbb{R}^{n\times d}$, where $X$ contains $x_{i}$ in its rows, is defined as the collection of sign patterns generated by the hyperplanes. Let $\mathcal{A}$ denote the set of all possible hyperplane arrangement patterns:

$$ $\displaystyle\mathcal{A}:=\cup\{\operatorname{sign}(Xw):w\in\mathbb{R}^{d}\},$ (9) $$

where the $\operatorname{sign}$ function is applied elementwise to the product $Xw$.

The number of distinct sign patterns in $\mathcal{A}$ is finite, i.e., $|\mathcal{A}|<\infty$. We define a subset $\mathcal{S}\subseteq\mathcal{A}$, representing the collection of sets corresponding to the positive signs in each element of $\mathcal{A}$:

$$ $\displaystyle\mathcal{S}:=\{\cup_{h_{i}=1}\{i\}:h\in\mathcal{A}\}.$ (10) $$

The cardinality of $\mathcal{S}$, denoted by $P$, represents the number of regions in the partition of $\mathbb{R}^{d}$ created by hyperplanes passing through the origin and orthogonal to the rows of $X$, or more birefly, $P$ is the number of regions formed by the hyperplane arrangement.

As the readers may see in the definition, $P$, the number of regions formed by the hyperplane arrangement, increases with both the dimension and the size of the dataset $X$. To reduce computational costs of our cutting-plane AL algorithm, as we will discuss in Appendix G.1, we sample the number of regions in the partition instead of using the entire $\mathcal{S}$.

In addition to center of gravity (Definition 2) and analytic center (Definition 1), another widely-used notion for center is the center of maximum volume inscribed ellipsoid (MVE), which we also referenced in the main text. We hereby give its definition (Boyd & Vandenberghe, ).

###### Definition 4 (Maximum volume inscribed ellipsoid (MVE)) .

Given a convex body $C\subseteq\mathbb{R}^{d}$, the maximum volume inscribed ellipsoid inside $C$ is found by solving the below optimization problem:

$$ $\displaystyle\begin{split}&\text{maximize \ }\log\det B\\ &\text{subject to \ }\sup_{\|u\|_{2}\leq 1}{\mathbbm{1}}_{C}(\varepsilon(u)) \leq 0,\end{split}$ (11) $$

where we have parametrized the ellipsoid as the image of the unit ball under an affine transformation:

$$ $\displaystyle\varepsilon(u):=\{Bu+c:\|u\|_{2}\leq 1\}$ $$

with $c\in\mathbb{R}^{d}$ and $B\in S^{n}_{++}$.
The optimal value of the variable $u$ in the convex program in 11 gives the center of the maximum volume inscribed ellipsoid of the convex body $C$, which is affine invariant. We denote it as $\theta_{M}(C)$, or $\theta_{M}$ as abbreviation.

A metric important to our evaluation for the regression prediction in Section 7 is the notion of root mean square error, or RMSE.

###### Definition 5 (Root Mean Square Error) .

The root mean square error (RMSE) is a metric used to measure the difference between predicted values and the actual values in a regression model. For a set of predicted values $\hat{y}_{i}$ and true values $y_{i}$, the RMSE is given by:

$$ $\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\hat{y}_{i}-y_{i})^{2}}$ $$

where $n$ is the total number of data points, $\hat{y}_{i}$ is the predicted value, and $y_{i}$ is the actual value.
RMSE provides an estimate of the standard deviation of the prediction errors (residuals), giving an overall measure of the model’s accuracy.

###### Definition 3 (Hyperplane Arrangement) .

###### Definition 4 (Maximum volume inscribed ellipsoid (MVE)) .

###### Definition 5 (Root Mean Square Error) .

#### B.2 Key Theorems

###### Key Theorems in Section 6.2 .

We present the cornerstone theorem to our results in Section 6.2 given by Grünbaum () on the bounds relating to convex body partitioning.

###### Proposition 1 (Grünbaum’s Inequality) .

Let $\mathcal{T}\in\mathbb{R}^{d}$ be a convex body (i.e. a compact convex set) and let $\theta_{G}$ denote its center of gravity. Let $\mathcal{H}=\{x\in\mathbb{R}^{d}:w^{T}(\theta-\theta_{G})=0\}$ be an arbitrary hyperplane passing through $\theta_{G}$. This plane divides the convex body $\mathcal{T}$ in the two subsets:

$$ $\displaystyle\mathcal{T}_{1}:=\{\theta\in\mathcal{T}:w^{T}\theta\geq w^{T} \theta_{G}\},$ $\displaystyle\mathcal{T}_{2}:=\{\theta\in\mathcal{T}:w^{T}\theta<w^{T}\theta_{ G}\}.$ $$

Then the following relations hold for $i=1,2$:

$$ $\displaystyle\textbf{vol}(\mathcal{T}_{i})\leq(1-(\frac{d}{d+1})^{1/d})\textbf {vol}(\mathcal{T})\leq(1-1/e)\textbf{vol}(\mathcal{T}).$ (12) $$

###### Proof.

See proof of Theorem 2 in Grünbaum ().
∎

###### Proposition 1 (Grünbaum’s Inequality) .

###### Proof.

### Appendix C Deferred Algorithms

#### C.1 More on Cutting-Plane AL for Binary Classification

##### C.1.1 Linear Model

Algorithm 3 describes the original cutting-plane-based learning algorithm proposed by Louche & Ralaivola (). This algorithm, despite having pioneered in bridging the classic cutting-plane optimization algorithm with active learning for the first time, remains limited to linear decision boundary classification and can only be applied to shallow machine learning models. We have demonstrated its inability to handle nonlinear decision boundaries and simple regression tasks in Section 7, where the linearity of the final decision boundary and prediction returned by the cutting-plane AL algorithm is evident. Nevertheless, Algorithm 3 establishes a crucial foundation for the development of our proposed cutting-plane active learning algorithms (Algorithm 2, 5, 4, and 6).

Figure: Algorithm 3 Generic (Linear) Cutting-Plane AL

##### C.1.2 NN Model with Limited Queries

We begin with Algorithm 4, which summarizes our proposed cutting-plane active learning algorithm under the second setup discussed in Section 6.1. In this case, the cutting-plane oracle has access to limited queries provided by the user and, in contrast to its query synthesis alternative (Algorithm 2), only makes the cut if the queried center mis-classifies the returned data point from the oracle. Hence, the convergence speed associated with Algorithm 4 hinges on how often the algorithm queries a center which incorrectly classifies the returned queried points before it reaches the optimal classifier. Therefore, the performance of this algorithm in terms of convergence speed depends not only on the geometry of the parameter version space but also on the effectiveness of the Query function to identify points from the limited dataset which gives the most informative evaluation of the queried center.

While Algorithm 4 still maintains similar rate and convergence guarantees as Algorithm 2, to optimize the empirical performance of Algorithm 4 (see discussions in Section 6.2), we therefore modify the Query function to query twice, once for minimal margin and once for maximal margin, to maximize the chances of the oracle returned points in correctly identifying mis-classification of the queried center. This modification greatly aids the performance of Algorithm 4, allowing the algorithm to make effective classification given very limited data. This is demonstrated in our experiment results in Section 7.

Figure: Algorithm 4 Cutting-plane AL for Binary Classification with Limited Queries

##### C.1.3 NN Model with Inexact Cut

Algorithm 5 summarizes our proposed cutting-plane active learning algorithm under the third setup metioned in Section 6.1. Under this scenario, the cutting-plane oracle has access to limited queries. However, in contrast to Algorithm 4, this cutting-plane AL algorithm always performs the cut regardless of whether the queried center mis-classifies the data point returned by the oracle. This is an interesting extension to consider. On the one hand, it can possibly speed up Algorithm 4 in making a decision boundary as the cut is effective in every iteration. On the other hand, however, the algorithm’s lack of discern for the correctness of the queried candidate presents a non-trivial challenge to evaluate its convergence rate and whether it still maintains convergence guarantees. It turns out that we can still ensure convergence in the case of Algorithm 5, and the convergence rate can be quantified by measuring the “inexactness” of the cut in relation to a cut which directly passes through the queried center. We refer the readers to a detailed discussion on this matter in Appendix F.2. We would also like to emphasize that Algorithm 2, 5, and 4, although written for binary classification, can be easily extended to multi-class by using, for instance, the “one-versus-all” strategy. See Appendix D.2 for details.

While it may be intuitive for the optimal performance of the cutting-plane AL algorithms to translate from binary classification to the multi-class case, it is not entirely evident for us to expect similar performance of the algorithms for regression tasks, where the number of classes $K\rightarrow\infty$. What is surprising is that our cutting-plane AL still maintains its optimal performance on regression tasks, as evidenced by the synthetic toy example using quadratic regression in Section 7. Nevertheless, one can argue that the result is not so surprising after all as it is to be expected in theory due to intuition explained in Appendix D.2.

Figure: Algorithm 5 Cutting-plane AL for Binary Classification with Inexact Cutting

#### C.2 More on Cutting-Plane AL for Regression

##### C.2.1 NN Model with Limited Queries

To generalize our classification cutting-plane AL algorithm to regression tasks, we need to make some adaptations to the cutting criterion and to how the cuts are being made. As the main body of the algorithm is the same across different setups (e.g. query synthesis, limited query, and inexact cuts) except for minor changes as the reader can see in Algorithm 2, 4, and 5, we only present the cutting-plane AL algorithm for regression under limited query. Algorithm 6 summarizes our proposed algorithm for training regression models via cutting-plane active learning. Here $\epsilon>0$ is a threshold value for the $L_{2}-$norm error chosen by the user. Observe that the new $L_{2}$-norm cut of step

$$ $\mathcal{T}^{t+1}\leftarrow\mathcal{T}^{t}\cap\{\theta:\|y_{n_{t}}-f^{\text{ two-layer}}(x_{n_{t}};\theta)\|_{2}\leq\epsilon,\mathcal{C}(\{n_{t}\}), \mathcal{C}^{\prime}(\{n_{t}\})\}$ $$

consists simply of two linear cuts:

$$ $-\epsilon\leq y_{n_{t}}-f^{\text{two-layer}}(x_{n_{t}};\theta)\leq\epsilon.$ $$

We can hence still ensure that the version space remains convex after each cut.

Figure: Algorithm 6 Cutting-plane AL for Regression with Limited Queries

##### C.2.2 Linear Model

Following the adaptation of our cutting-plane AL from classification to regression tasks, we similarly attempt to adapt the original linear cutting-plane AL (Algorithm 3) for regression. Algorithm 7 introduces an $\epsilon>0$ threshold to account for the $L_{2}$-norm error between the predicted value and the actual target. This threshold controls both when a cut is made and the size of the cut. We applied this version of the algorithm (Algorithm 7) in the quadratic regression experiment detailed in Section 7. However, for nonlinear regression data—such as quadratic regression—this approach will necessarily fail. The prediction model $\langle\theta,x\rangle$ in Algorithm 7 is linear, whereas the underlying data distribution follows a nonlinear relationship, i.e., $y=x^{2}$. As a result, once the query budget starts accumulating nonlinearly distributed data points, no linear predictor can satisfy the regression task’s requirements for small values of $\epsilon$.

In particular, after a certain number of iterations $t$, the cut

$$ $\mathcal{T}^{t}\cap\{\theta:\|y_{n_{t}}-\langle\theta,x_{n_{t}}\rangle\|_{2} \leq\epsilon\}$ $$

will eliminate the entire version space (i.e., $\mathcal{T}^{t+1}=\emptyset$). This occurs because the error between the linear prediction and the nonlinear true values (e.g., $y=x^{2}$) cannot be reduced sufficiently, disqualifying all linear predictors. This is precisely what we observed in our quadratic regression example in Section 7. After four queries, the algorithm reported infeasibility in solving for the version space center, as the remaining version space had been reduced to an empty set. This infeasibility is a direct consequence of the mismatch between the linear model’s capacity and the data’s nonlinear nature.

Figure: Algorithm 7 Linear Cutting-Plane AL for Regression

#### C.3 Minimal Margin Query Strategy

Here we note that in our main algorithm 2, we always query points with highest prediction confidence by setting

$$ $(x,y)\leftarrow\operatorname*{arg\,min}_{(x_{i},y_{i})\in\mathcal{D}_{\text{QS }}}sf^{\text{two-layer}}(x_{n_{t}};\theta).$ $$

However, we indeed allow more custom implementation of query selection. For example, an alternative approach is to select data with minimal prediction margin, i.e.

$$ $(x,y)\leftarrow\operatorname*{arg\,min}_{(x_{i},y_{i})\in\mathcal{D}_{\text{QS }}}|f^{\text{two-layer}}(x_{n_{t}};\theta)|.$ $$

We experiment with this query strategy in our real dataset experiments.

### Appendix D Key Generalization to Cutting-Plane AL

In this section, we discuss two important generalizations of our cutting-plane AL method: (i). the relaxation in data distribution requirement and (ii). the extension from classification to regression.

#### D.1 Relaxed Data Distribution Requirement

For linear classifier $f(x;\theta)=x^{T}\theta$, the training data is expected to be linearly separable for an optimal $\theta^{\star}$ to exist. However, this constraint on training data distribution is too restrictive in real scenarios. With ReLU model, due to its uniform approximation capacity, the training data is not required to be linearly separable as long as there is a continuous function $h$ such that ${\mathrm{sign}}(h(x))=y$ for all $(x,y)$ pairs. Due to the discrete nature of sampled data points, this is always satisfiable, thus we can totally remove the prerequisite on training data. For sake of completeness, we provide a version of uniform approximation capacity of ReLU below, with an extended discussion about a sample compression perspective of our cutting-plane based model training scheme thereafter. Given the fact that two-layer model has weaker approximation capacity compared to deeper models, we here consider only two-layer model without loss of generality.
Uniform approximation capacity of single hidden layer NN has been heavily studied (Chen & Chen, ; Chui & Li, ; Costarelli et al., ; Cotter, ; Cybenko, ; ichi Funahashi, ; Gallant & White, ; Hornik, ; Mhaskar & Micchelli, ; Leshno et al., ; Pinkus, ), here we present the version of ReLU network which we find most explicit and close to our setting.

###### Theorem D.1 .

(Theorem 1.1 in (Shen et al., ))
For any continuous function $f\in C([0,1]^{d})$, there exists a two-layer ReLU network $\phi=(xW_{1}+b_{1})_{+}W_{2}+b_{2}$ such that

$$ $\|f-\phi\|_{L^{p}([0,1]^{d})}\in\mathcal{O}\left(\sqrt{d}\omega_{f}\left((m^{2 }\log m)^{-1/d}\right)\right)$ $$

where $\omega_{f}(\cdot)$ is the modulus of continuity of $f$.

The above result can be extended to any $f\in C([-R,R]^{d})$, see Theorem 2.5 in Shen et al. () for discussion. Therefore, for any fixed dimension $d$, as we increase number of hidden neurons, we are guaranteed to approximate $f$ a.e. (under proper condition on $\omega_{f}$). Thus the assumption on training data can be relaxed to existence of some $\{W_{1},b_{1},W_{2},b_{2}\}$ such that ${\mathrm{sign}}(\phi(x))=y$, which happens almost surely. A minor discrepancy here is that the $\phi$ being considered in above Theorem incorporates the bias terms $b_{1},b_{2}$ while our Theorem 4.1 considers two-layer NN of form $(xW_{1})_{+}W_{2}$. We note here that $(xW_{1})_{+}W_{2}$ is essentially the same as the one with layer-wise bias. To see this, we can append the data $x$ by a $1$-value entry to accommodate for $b_{1}$. After this modification, we can always separate a neuron of form $(XW_{1i})_{+}W_{2i}$, set $W_{1(d+1)}=1$ and zeros elsewhere, and $W_{2i}$ would then accommodate for bias $b_{2}.$ Therefore, our relaxed requirement on training data still persists. For deeper ReLU NNs, similar uniform approximation capability has been established priorly and furthermore the general NN form we considered in Theorem 4.2 incorporates the form with layer-wise bias and thus we still have such relaxation on training data distribution, see Appendix E.3 for more explanation.

###### Theorem D.1 .

#### D.2 From Classification to Regression

Algorithm 2 shows how we train deep NN for binary classification task with cutting-plane method. Nevertheless, it can be easily extended to multi-class using, for instance, the so-called “one-versus-all” classification strategy. To illustrate, to extend the binary classifier to handle $K$ classes under this approach, we decompose the multi-class problem into $K$ binary subproblems. Specifically, for each class $\mathcal{C}_{k}$, we define a binary classification task as:

$$ $\text{Classify between }\mathcal{C}_{k}\text{ and }\bigcup_{i\neq k}\mathcal{C }_{i},$ $$

where we classify $\mathcal{C}_{k}$ against all other classes combined as a single class. This creates $K$ binary classification problems, each corresponding to distinguishing one class from the rest.

In fact, our cutting-plane AL can be even applied to the case of regression, where the number of classes $K\rightarrow\infty$. The core intuition behind this is still the uniform approximation capability of nonlinear ReLU model. Given any training data $x$ with its label $y$, we want a model $f(x;\theta)$ to be able to predict $y$ exactly. Here the data label is no longer limited to plus or minus one and can be any continuous real number. For linear model $f(x;\theta)=x^{T}\theta$ considered in (Louche & Ralaivola, ), train such a predictor is almost impossible since it is highly unlikely there exists such a $\theta$ for real dataset. However, with our ReLU model, we are guaranteed there is a set of NN weights that would serve as a desired predictor due to its uniform approximation capacity.

Therefore, for regression task, the training algorithm will be exactly the same as Algorithm 3 instead that the original classification cut $y_{n_{t}}f(x_{n_{t}};\theta)\geq 1$ will be replaced with $f(x_{n_{t}};\theta)=y_{n_{t}}.$ All other activation pattern constraints are leaved unchanged. A minor discrepancy here is that though theoretically sounding, the strict inequality $f(x_{n_{t}};\theta)=y_{n_{t}}$ may raise numerical issues in practical implementation. Thus, we always include a trust region as $y_{n_{t}}-\epsilon\leq f(x_{n_{t}};\theta)\leq y_{n_{t}}+\epsilon$ with some small $\epsilon$ for our experiments. See Algorithm 6 for our implementation details.

### Appendix E Deferred Proofs and Extensions in Section 4

#### E.1 Proof of Theorem 4.1

We prove the equivalence in two directions, we first show that if there exists $W_{1},W_{2}$ to $y\odot((XW_{1})_{+}W_{2})\geq 1$, then we can find solution $\{u_{i},u_{i}^{\prime}\}$ to Problem (3); we then show that when there is solution $\{u_{i},u_{i}^{\prime}\}$ to Problem (3) and the number of hidden neurons $m\geq 2P$, then there exists $W_{1},W_{2}$ such that $y\odot((XW_{1})_{+}W_{2})\geq 1$. We last show that given solution $\{u_{i},u_{i}^{\prime}\}$ to Problem (3), after finding correspondent $W_{1},W_{2}$, the prediction for any input $\tilde{x}$ can be simply computed by $\sum_{i=1}^{P}(\tilde{x}u_{i})_{+}-(\tilde{x}u_{i}^{\prime})_{+}$. We now start with our first part and assume the existence of $W_{1},W_{2}$ to $y\odot((XW_{1})_{+}W_{2})\geq 1$, i.e.,

$$ $\displaystyle y\odot\sum_{j=1}^{m}(XW_{1j})_{+}W_{2j}\geq 1,$ $$

from which we can derive

$$ $y\odot\sum_{j=1}^{m}D^{1}_{j}XW_{1j}W_{2j}\geq 1,$ (13) $$

where $D^{1}_{j}=\text{diag}(\mathbbm{1}\{XW_{1j}\geq 0\})\in\mathbb{R}^{n\times n}$. Now consider set of pairs of $\{u_{j},u_{j}^{\prime}\}$ given by $u_{j}=W_{1j}W_{2j},u_{j}^{\prime}=0$ for $j\in\{j|W_{2j}\geq 0\}$ and $u_{j}=0,u_{j}^{\prime}=-W_{1j}W_{2j}$ for $j\in\{j|W_{2j}<0\}.$ We thus have by (13)

$$ $y\odot\sum_{j=1}^{m}D^{1}_{j}X(u_{j}-u_{j}^{\prime})\geq 1,(2D^{1}_{j}-I)Xu_{j }\geq 0,(2D^{1}_{j}-I)Xu_{j}^{\prime}\geq 0.$ $$

The only discrepancy between our set of $\{u_{j},u_{j}^{\prime}\}$ pairs and our desired solution to Problem (3) is we want to match $D^{1}_{j}$’s in equation (13) to $D_{i}$’s in Problem (3). This is achieved by observing that whenever we have $D^{1}_{a}=D^{1}_{b}=D_{(a,b)}$ for some $a,b\in[m]$, we can merge them as

$$ $\displaystyle D^{1}_{a}X(u_{a}-u_{a}^{\prime})+D^{1}_{b}X(u_{b}-u_{b}^{\prime})$ $\displaystyle=D_{(a,b)}X((u_{a}+u_{b})-(u_{a}^{\prime}+u_{b}^{\prime}))$ $\displaystyle=D_{(a,b)}X(u_{a+b}-u_{a+b}^{\prime}),$ $$

with $(2D_{(a+b)}-I)Xu_{a+b}\geq 0$ and $(2D_{(a+b)}-I)Xu_{a+b}^{\prime}\geq 0$ still hold. We can keep this merging for all activation patterns $\{D_{j}^{1}|j\in[m]\}.$ We are guaranteed to get

$$ $y\odot\sum_{j=1}^{\tilde{m}}\tilde{D}_{j}X(\tilde{u}_{j}-\tilde{u}_{j}^{\prime })\geq 1,(2\tilde{D}_{j}-I)X\tilde{u}_{j}\geq 0,(2\tilde{D}_{j}-I)X\tilde{u}_{ j}^{\prime},\geq 0$ $$

where all $\tilde{D}_{j},j\in[\tilde{m}]$ are different. Note since our $\{D_{i}|i\in[P]\}$ in Problem (3) loop over all possible activation patterns corresponding to $X$, it is always the case that $\tilde{m}\leq P$ and $\tilde{D}_{j}=D_{i}$ for some $i\in[P].$ Thus we get a solution $\{u_{i},u_{i}^{\prime}\}$ to Problem (3) by setting $u_{i}=\tilde{u}_{k},u_{i}^{\prime}=\tilde{u}_{k}^{\prime}$ when $D_{i}=\tilde{D}_{k}$ for some $k\in[\tilde{m}].$ Otherwise we simply set $u_{i}=u_{i}^{\prime}=0.$ This completes our proof of the first direction, we now turn to prove the second direction and assume that there is solution $\{u_{i},u_{i}^{\prime}\}$ to Problem (3) as well as the number of hidden neurons $m\geq 2P$. We aim to show that there exists $W_{1},W_{2}$ such that $y\odot((XW_{1})_{+}W_{2})\geq 1.$ Since we have

$$ $y\odot\sum_{i=1}^{P}\left(D_{i}X(u_{i}-u_{i}^{\prime})\right)\geq 1,(2D_{i}-I) Xu_{i}\geq 0,(2D_{i}-I)Xu_{i}^{\prime}\geq 0,$ $$

we are able to derive the below inequality by setting $v_{i}=u_{i},\alpha_{i}=1$ for $i\in[P]$ and $v_{i}=u_{i-P}^{\prime},\alpha_{i}=-1$ for $i\in[P+1,2P],$

$$ $y\odot\sum_{i=1}^{2P}(Xv_{i})_{+}\alpha_{i}\geq 1.$ $$

Therefore, consider $W_{1}\in\mathbb{R}^{d\times m}$ defined by $W_{1j}=v_{j}$ for $j\in[2P]$ and $W_{1j}=0$ for any $j>2P$, $W_{2}\in\mathbb{R}^{m}$ defined by $W_{2j}=\alpha_{j}$ for $j\in[2P]$ and $W_{2j}=0$ for any $j>2P.$ Then we achieve

$$ $y\odot((XW_{1})_{+}W_{2})=y\odot\sum_{i=1}^{2P}(Xv_{i})_{+}\alpha_{i}\geq 1,$ $$

as desired. Lastly, once a solution $\{u_{i},u_{i}^{\prime}\}$ to Problem (3) is given, we can find corresponding $W_{1},W_{2}$ according to our analysis of second direction above. Then for any input $\tilde{x}$, the prediction given by $(\tilde{x}W_{1})_{+}W_{2}$ is simply $\sum_{i=1}^{P}(\tilde{x}u_{i})_{+}-(\tilde{x}u_{i}^{\prime})_{+}.$

#### E.2 Proof of Theorem 4.2

Similar to proof of Theorem 4.1, we carry out the proof in two directions. We prove first that if there exists solution to Problem (5), then there is also solution to Problem (6). We then show that whenever there exists solution to Problem (6), there is also solution to Problem (5). Concerned with that the notation complexity of $(n+1)$-layer NN might introduce difficulty to follow the proof, we start with showing three-layer case as a concrete example, we then move on to $(n+1)$-layer proof which is more arbitrary.

Proof for three-layer. The three-layer ReLU model being considered is of the form $((XW_{1})_{+}W_{2})_{+}W_{3}.$ The corresponding linear program is given by

$$ find $\displaystyle u_{ij},u_{ij}^{\prime},v_{ij},v_{ij}^{\prime}$ s.t. $\displaystyle y\odot\sum_{j=1}^{P_{2}}D^{(2)}_{j}\left(\sum_{i=1}^{P_{1}}D_{i} ^{(1)}X(u_{ij}-u^{\prime}_{ij})-\sum_{i=1}^{P_{1}}D_{i}^{(1)}X(v_{ij}-v^{ \prime}_{ij})\right)\geq 1,$ $\displaystyle(2D_{i}^{(1)}-I)Xu_{ij}\geq 0,(2D_{i}^{(1)}-I)Xu^{\prime}_{ij} \geq 0,(2D_{i}^{(1)}-I)Xv_{ij}\geq 0,(2D_{i}^{(1)}-I)Xv^{\prime}_{ij}\geq 0,{}$ $\displaystyle(2D_{j}^{(2)}-I)\left(\sum_{i=1}^{P_{1}}D_{i}^{(1)}X(u_{ij}-u^{ \prime}_{ij})\right)\geq 0,(2D_{j}^{(2)}-I)\left(\sum_{i=1}^{P_{1}}D_{i}^{(1)} X(v_{ij}-v^{\prime}_{ij})\right)\geq 0,{}{}{}$ $$

which is a rewrite of (6) with $n=2.$ Firstly, assume there exists solution $\{W_{1},W_{2},W_{3}\}$ to the problem $y\odot((XW_{1})_{+}W_{2})_{+}W_{3}\geq 1$. We want to show there exists $\{u_{ij},u_{ij}^{\prime},v_{ij},v_{ij}^{\prime}\}$ solves the above problem. Note that by $y\odot((XW_{1})_{+}W_{2})_{+}W_{3}\geq 1$, we get

$$ $y\odot\left(\sum_{i=1}^{m_{1}}(XW_{1i})_{+}W_{2i}\right)_{+}W_{3}\geq 1.$ $$

Let $K_{i}$ denote ${\mathrm{sign}}(XW_{1i})$, i.e, $(2K_{i}-I)XW_{1i}\geq 0$, we can write

$$ $y\odot\left(\sum_{i=1}^{m_{1}}K_{i}XW_{1i}W_{2i}\right)_{+}W_{3}\geq 1.$ $$

Expand on the outer layer neurons, we get

$$ $y\odot\sum_{j=1}^{m_{2}}\left(\sum_{i=1}^{m_{1}}K_{i}XW_{1i}W_{2ij}\right)_{+} W_{3j}\geq 1.$ $$

Construct $c_{ij}=W_{1i}W_{2ij}$ whenever $W_{2ij}\geq 0$ and 0 0 otherwise, $c^{\prime}_{ij}=-W_{1i}W_{2ij}$ whenever $W_{2ij}<0$ and 0 0 otherwise, we can write

$$ $y\odot\sum_{j=1}^{m_{2}}\left(\sum_{i=1}^{m_{1}}K_{i}X(c_{ij}-c_{ij}^{\prime}) \right)_{+}W_{3j}\geq 1,(2K_{i}-I)Xc_{ij}\geq 0,(2K_{i}-I)Xc_{ij}^{\prime}\geq 0.$ $$

Denote ${\mathrm{sign}}(\sum_{i=1}^{m_{1}}K_{i}X(c_{ij}-c_{ij}^{\prime}))$ as $K_{j}^{(2)}$, we thus have

$$ $y\odot\sum_{j=1}^{m_{2}}K_{j}^{(2)}\left(\sum_{i=1}^{m_{1}}K_{i}X(c_{ij}-c_{ij }^{\prime})\right)W_{3j}\geq 1,$ $$

with $(2K_{i}-I)Xc_{ij}\geq 0,(2K_{i}-I)Xc_{ij}^{\prime}\geq 0,(2K_{j}^{(2)}-I)(\sum
_{i=1}^{m_{1}}K_{i}X(c_{ij}-c_{ij}^{\prime}))\geq 0.$ We construct $\{d_{ij},d_{ij}^{\prime},e_{ij},e_{ij}^{\prime}\}$ by setting $d_{ij}=c_{ij}W_{3j},d_{ij}^{\prime}=c_{ij}^{\prime}W_{3j}$ when $W_{3j}\geq 0$ and 0 0 otherwise, setting $e_{ij}=-c_{ij}W_{3j},e_{ij}^{\prime}=-c_{ij}^{\prime}W_{3j}$ when $W_{3j}<0$ and 0 0 otherwise, and we will arrive at

$$ $y\odot\sum_{j=1}^{m_{2}}K_{j}^{(2)}\left(\sum_{i=1}^{m_{1}}K_{i}X(d_{ij}-d_{ij }^{\prime})-\sum_{i=1}^{m_{1}}K_{i}X(e_{ij}-e_{ij}^{\prime})\right)\geq 1,$ $$

where

$$ $\displaystyle(2K_{i}-I)Xd_{ij}\geq 0,(2K_{i}-I)Xd_{ij}^{\prime}\geq 0,(2K_{i}- I)Xe_{ij}\geq 0,(2K_{i}-I)Xe_{ij}^{\prime}\geq 0$ $\displaystyle(2K_{j}^{(2)}-I)\sum_{i=1}^{m_{1}}K_{i}X(d_{ij}-d_{ij}^{\prime}) \geq 0,(2K_{j}^{(2)}-I)\sum_{i=1}^{m_{1}}K_{i}X(e_{ij}-e_{ij}^{\prime})\geq 0,$ $$

which is already of form we want. We are left with matching $m_{i}$ with $P_{i}$ for $i\in\{1,2\}$ and matching $D_{j}^{(2)}$ to $K_{j}^{(2)}$, $D_{i}^{(1)}$ to $K_{i}$, $\{u_{ij},u_{ij}^{\prime},v_{ij},v_{ij}^{\prime}\}$ to $\{d_{ij},d_{ij}^{\prime},e_{ij},e_{ij}^{\prime}\}.$ We achieve this by observing that whenever there is duplicate $K_{j_{1}}^{(2)}=K_{j_{2}}^{(2)},$ we can merge the corresponding terms as

$$ $\displaystyle K_{j_{1}}^{(2)}\left(\sum_{i=1}^{m_{1}}K_{i}X(d_{ij_{1}}-d_{ij_{ 1}}^{\prime})-\sum_{i=1}^{m_{1}}K_{i}X(e_{ij_{1}}-e_{ij_{1}}^{\prime})\right)+ K_{j_{2}}^{(2)}\left(\sum_{i=1}^{m_{1}}K_{i}X(d_{ij_{2}}-d_{ij_{2}}^{\prime})- \sum_{i=1}^{m_{1}}K_{i}X(e_{ij_{2}}-e_{ij_{2}}^{\prime})\right)$ $\displaystyle=K^{(2)}_{(j_{1},j_{2})}\left(\sum_{i=1}^{m_{1}}K_{i}X\left((d_{ ij_{1}}+d_{ij_{2}})-(d_{ij_{1}}^{\prime}+d_{ij_{2}}^{\prime})\right)-\sum_{i=1 }^{m_{1}}K_{i}X\left((e_{ij_{1}}+e_{ij_{2}})-(e_{ij_{1}}^{\prime}+e_{ij_{2}}^{ \prime})\right)\right)$ $\displaystyle=K^{(2)}_{(j_{1},j_{2})}\left(\sum_{i=1}^{m_{1}}K_{i}X(d_{i(j_{1} +j_{2})}-d^{\prime}_{i(j_{1}+j_{2})})-\sum_{i=1}^{m_{1}}K_{i}X(e_{i(j_{1}+j_{2 })}-e^{\prime}_{i(j_{1}+j_{2})})\right),$ $$

where $K^{(2)}_{(j_{1},j_{2})}:=K_{j_{1}}^{(2)},d_{i(j_{1}+j_{2})}:=d_{ij_{1}}+d_{ij_
{2}},d^{\prime}_{i(j_{1}+j_{2})}=d_{ij_{1}}^{\prime}+d_{ij_{2}}^{\prime},e_{i(
j_{1}+j_{2})}=e_{ij_{1}}+e_{ij_{2}},$ $e^{\prime}_{i(j_{1}+j_{2})}=e_{ij_{1}}^{\prime}+e_{ij_{2}}^{\prime}$. We have as constraints

$$ $\displaystyle(2K^{(2)}_{(j_{1},j_{2})}-I)\left(\sum_{i=1}^{m_{1}}K_{i}X(d_{i(j _{1}+j_{2})}-d^{\prime}_{i(j_{1}+j_{2})})\right)\geq 0,$ $\displaystyle(2K^{(2)}_{(j_{1},j_{2})}-I)\left(\sum_{i=1}^{m_{1}}K_{i}X(e_{i(j _{1}+j_{2})}-e^{\prime}_{i(j_{1}+j_{2})})\right)\geq 0,$ $\displaystyle(2K_{i}-I)Xd_{i(j_{1}+j_{2})}\geq 0,(2K_{i}-I)Xd^{\prime}_{i(j_{1 }+j_{2})}\geq 0,$ $\displaystyle(2K_{i}-I)Xe_{i(j_{1}+j_{2})}\geq 0,(2K_{i}-I)Xe^{\prime}_{i(j_{1 }+j_{2})}\geq 0.$ $$

Keep such merging until all $K_{j}^{(2)}$ are different, we arrive at

$$ $y\odot\sum_{j=1}^{\tilde{m}_{2}}\overline{K}_{j}^{(2)}\left(\sum_{i=1}^{m_{1}} K_{i}X(\overline{d}_{ij}-\overline{d^{\prime}_{ij}})-\sum_{i=1}^{m_{1}}K_{i}X( \overline{e}_{ij}-\overline{e}^{\prime}_{ij})\right)\geq 1,$ $$

with

$$ $\displaystyle(2K_{i}-I)X\overline{d}_{ij}\geq 0,(2K_{i}-I)X\overline{d}^{ \prime}_{ij}\geq 0,(2K_{i}-I)X\overline{e}_{ij}\geq 0,(2K_{i}-I)X\overline{e}^ {\prime}_{ij}\geq 0,$ $\displaystyle(2\overline{K}_{j}^{(2)}-I)\left(\sum_{i=1}^{m_{1}}K_{i}X( \overline{d}_{ij}-\overline{d}^{\prime}_{ij})\right)\geq 0,(2\overline{K}_{j}^ {(2)}-I)\left(\sum_{i=1}^{m_{1}}K_{i}X(\overline{e}_{ij}-\overline{e}^{\prime} _{ij})\right)\geq 0,$ $$

where $\tilde{m}_{2}\leq P_{2},\overline{K}_{j}^{(2)}\in\{D^{(2)}\}$ and $\overline{K}_{j}^{(2)}$ all different. We now proceed to match $K_{i}$ and $D_{i}^{(1)}.$ Consider $\sum_{i=1}^{m_{1}}K_{i}X(\overline{d}_{ij}-\overline{d}_{ij}^{\prime})$, if $K_{v}=K_{q}$ for some $v,q$, we can merge them as

$$ $K_{v}X(\overline{d}_{vj}-\overline{d}_{vj}^{\prime})+K_{q}X(\overline{d}_{qj}- \overline{d}_{qj}^{\prime})=K_{(v,q)}X(\overline{d}_{(v+q)j}-\overline{d}^{ \prime}_{(v+q)j}),$ $$

where $K_{(v,q)}:=K_{v},\overline{d}_{(v+q)j}:=\overline{d}_{vj}+\overline{d}_{qj},
\overline{d}^{\prime}_{(v+q)j}:=\overline{d}_{vj}^{\prime}+\overline{d}^{
\prime}_{(v+q)j}.$ The constraints are $(2K_{(v,q)}-I)X\overline{d}_{(v+q)j}\geq 0,(2K_{v,q}-I)X\overline{d}_{(v+q)j}^
{\prime}\geq 0$, and we still have

$$ $(2\overline{K}_{j}^{(2)}-I)\left(\sum_{i\in[m_{1}],i\neq v,i\neq q}K_{i}X( \overline{d}_{ij}-\overline{d}_{ij^{\prime}})+K_{(v,q)}X(\overline{d}_{(v+q)j} -\overline{d}_{(v+q)j}^{\prime})\right)\geq 0.$ $$

Continue such merging and also for $\{\overline{e}_{ij},\overline{e}_{ij}^{\prime}\}$, we finally arrive at

$$ $y\odot\sum_{j=1}^{\tilde{m}_{2}}\overline{K}_{j}^{(2)}\left(\sum_{i=1}^{\tilde {m}_{1}}\hat{K}_{i}X(\hat{\overline{d}}_{ij}-\hat{\overline{d}}^{\prime}_{ij}) -\sum_{i=1}^{\tilde{m}_{1}}\hat{K}_{i}X(\hat{\overline{e}}_{ij}-\hat{\overline {e}}^{\prime}_{ij})\right)\geq 1,$ $$

with

$$ $(2\hat{K}_{i}-I)X\hat{\overline{d}}_{ij}\geq 0,(2\hat{K}_{i}-I)X\hat{\overline {d}}^{\prime}_{ij}\geq 0,(2\hat{K}_{i}-I)X\hat{\overline{e}}_{ij}\geq 0,(2\hat {K}_{i}-I)X\hat{\overline{e}}^{\prime}_{ij}\geq 0,$ $$

$$ $(2\overline{K}_{j}^{(2)}-I)\left(\sum_{i=1}^{\tilde{m}_{1}}\hat{K}_{i}X(\hat{ \overline{d}}_{ij}-\hat{\overline{d}}^{\prime}_{ij})\right)\geq 0,(2\overline{ K}_{j}^{(2)}-I)\left(\sum_{i=1}^{\tilde{m}_{1}}\hat{K}_{i}X(\hat{\overline{e}} _{ij}-\hat{\overline{e}}^{\prime}_{ij})\right)\geq 0.$ $$

Now, for any $D_{j}^{(2)}\not\in\{\overline{K}^{(2)}\},$ we set all $u_{ij}=u_{ij}^{\prime}=v_{ij}=v_{ij}^{\prime}=0.$ For $D_{j}^{(2)}=\overline{K}_{j^{\prime}}^{(2)},$ for any $D_{i}^{(1)}\not\in\{\hat{K}\}$, we set all $u_{ij}=u_{ij}^{\prime}=v_{ij}=v_{ij}^{\prime}=0$. For $D_{j}^{(2)}=\overline{K}_{j^{\prime}}^{(2)},D_{i}^{(1)}=\hat{K}_{i^{\prime}},$ we set $u_{ij}=\hat{\overline{d}}_{i^{\prime}j^{\prime}},u^{\prime}_{ij}=\hat{
\overline{d}}^{\prime}_{i^{\prime}j^{\prime}},v_{ij}=\hat{\overline{e}}_{i^{
\prime}j^{\prime}},v^{\prime}_{ij}=\hat{\overline{e}}^{\prime}_{i^{\prime}j^{
\prime}}$. Then we get exactly

$$ $y\odot\left(\sum_{j=1}^{P_{2}}D_{j}^{(2)}\left(\sum_{i=1}^{P_{1}}D_{i}^{(1)}X( u_{ij}-u^{\prime}_{ij})-\sum_{i=1}^{P_{1}}D_{i}^{(1)}X(v_{ij}-v_{ij}^{\prime}) \right)\right)\geq 1,$ $$

and

$$ $(2D_{i}^{(1)}-I)Xu_{ij}\geq 0,(2D_{i}^{(1)}-I)Xu^{\prime}_{ij}\geq 0,(2D_{i}^{ (1)}-I)Xv_{ij}\geq 0,(2D_{i}^{(1)}-I)Xv^{\prime}_{ij}\geq 0,$ $$

$$ $(2D_{j}^{(2)}-I)\left(\sum_{i=1}^{P_{1}}D_{i}^{(1)}X(u_{ij}-u_{ij}^{\prime}) \right)\geq 0,(2D_{j}^{(2)}-I)\left(\sum_{i=1}^{P_{1}}D_{i}^{(1)}X(v_{ij}-v_{ ij}^{\prime})\right)\geq 0,$ $$

which completes the proof of our first direction. We now turn on to prove the second direction, assume there exists $u_{ij},u_{ij}^{\prime},v_{ij},v_{ij}^{\prime}$ such that

$$ $y\odot\left(\sum_{j=1}^{P_{2}}D_{j}^{(2)}\left(\sum_{i=1}^{P_{1}}D_{i}^{(1)}X( u_{ij}-u_{ij}^{\prime})-\sum_{i=1}^{P_{1}}D_{i}^{(1)}X(v_{ij}-v^{\prime}_{ij}) \right)\right)\geq 1,$ $$

with

$$ $(2D_{i}^{(1)}-I)Xu_{ij}\geq 0,(2D_{i}^{(1)}-I)Xu^{\prime}_{ij}\geq 0,(2D_{i}^{ (1)}-I)Xv_{ij}\geq 0,(2D_{i}^{(1)}-I)Xv^{\prime}_{ij}\geq 0,$ $$

$$ $(2D_{j}^{(2)}-I)\left(\sum_{i=1}^{P_{1}}D_{i}X(u_{ij}-u_{ij}^{\prime})\right) \geq 0,(2D_{j}^{(2)}-I)\left(\sum_{i=1}^{P_{1}}D_{i}X(v_{ij}-v_{ij}^{\prime}) \right)\geq 0.$ $$

We want to show that there exists $W_{1},W_{2},W_{3}$ such that

$$ $y\odot((XW_{1})_{+}W_{2})_{+}W_{3}\geq 1.$ $$

We are able to derive

$$ $y\odot\left(\sum_{j=1}^{P_{2}}\left(\sum_{i=1}^{P_{1}}D_{i}X(u_{ij}-u_{ij}^{ \prime})\right)_{+}-\left(\sum_{i=1}^{P_{1}}D_{i}X(v_{ij}-v_{ij}^{\prime}) \right)_{+}\right)\geq 1,$ $$

and furthermore

$$ $y\odot\left(\sum_{j=1}^{P_{2}}\left(\sum_{i=1}^{P_{1}}(Xu_{ij})_{+}-(Xu_{ij}^{ \prime})_{+}\right)_{+}-\left(\sum_{i=1}^{P_{1}}(Xv_{ij})_{+}-(Xv_{ij}^{\prime })_{+}\right)_{+}\right)\geq 1.$ $$

We thus construct $\{\kappa_{ij},\alpha_{ij}\}$ by setting $\kappa_{ij}=u_{ij},\alpha_{ij}=1$ for $i\in[P_{1}],\kappa_{ij}=u_{(i-P_{1})j}^{\prime},\alpha_{ij}=-1$ for $i\in[P_{1}+1,2P_{1}].$ We similarly construct $\{\kappa_{ij}^{\prime},\alpha_{ij}^{\prime}\}$ with $\{v_{ij},v_{ij}^{\prime}\}$, we thus get

$$ $y\odot\left(\sum_{j=1}^{P_{2}}\left(\sum_{i=1}^{2P_{1}}(X\kappa_{ij})_{+} \alpha_{ij}\right)_{+}-\left(\sum_{i=1}^{2P_{1}}(X\kappa_{ij}^{\prime})_{+} \alpha_{ij}^{\prime}\right)_{+}\right)\geq 1.$ $$

We construct $\{u_{j}\}$ by setting $u_{j}=1$ for $j\in[P_{2}]$ and $u_{j}=-1$ for $j\in[P_{2}+1,2P_{2}]$. We construct also $\{\tilde{\kappa}_{ij},\tilde{\alpha}_{ij}\}$ such that $\tilde{\kappa}_{ij}=\kappa_{ij},\tilde{\alpha}_{ij}=\alpha_{ij}$ for $i\in[2P_{1}],j\in[P_{2}]$, $\tilde{\kappa}_{ij}=\kappa^{\prime}_{i(j-P_{2})},\tilde{\alpha}_{ij}=\alpha^{
\prime}_{i(j-P_{2})}$ for $i\in[2P_{1}],j\in[P_{2}+1,2P_{2}]$, we thus get

$$ $y\odot\left(\sum_{j=1}^{2P_{2}}\left(\sum_{i=1}^{2P_{1}}(X\tilde{\kappa}_{ij}) _{+}\tilde{\alpha}_{ij}\right)_{+}u_{j}\right)\geq 1.$ $$

Therefore, we arrive at $y\odot((XW_{1})_{+}W_{2})_{+}W_{3}\geq 1$ by focusing on $j\leq 2P_{2},i\leq 4P_{1}P_{2}$, i.e., setting parameters with $(i,j)$ indices exceeding these thresholds to be all zero. We then set $W_{1i}=\tilde{\kappa}_{ab}$ for $a=\lfloor\frac{i-1}{2P_{1}}\rfloor+1,b=(i-1)\%2P_{1}+1,$ $W_{2ij}=\tilde{\alpha}_{(i-1)\%2P_{1}+1}$ for $i\in[(j-1)*2P_{1}+1,j*2P_{1}]$ and 0 0 otherwise, $W_{3j}=u_{j}.$

Proof for $n+1$-layer. We now provide proof for ReLU model of arbitrary depth. The logic follows the three-layer case proof above, despite the notation now represents $n+1$-layer model for arbitrary $n$. We first assume that there exists $W_{1},W_{2},\cdots,W_{n+1}$ that satisfies problem (5). We want to show that there exists $\{a_{j_{n}\cdots j_{1}}^{c_{n}\cdots c_{1}}\}$ which satisfies problem (6). We first span the inner most neuron to get

$$ $y\odot\left(\cdots\left(\left(\left(\sum_{i_{1}=1}^{m_{1}}(XW_{1i_{1}})_{+}W_{ 2i_{1}}\right)_{+}W_{3}\right)_{+}W_{4}\right)_{+}W_{5}\cdots\right)_{+}W_{n+1 }\geq 1.$ $$

Denote ${\mathrm{sign}}(XW_{1i})$ as $K_{i_{1}}^{(1)}$, i.e., $(2K_{i_{1}}^{(1)}-I)XW_{1i_{1}}\geq 0$. We can then rewrite the above inequality as

$$ $y\odot\left(\cdots\left(\left(\left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}XW_{1 i_{1}}W_{2i_{1}}\right)_{+}W_{3}\right)_{+}W_{4}\right)_{+}W_{5}\cdots\right)_ {+}W_{n+1}\geq 1.$ $$

We then expand the second last inner layer as

$$ $y\odot\left(\cdots\left(\left(\sum_{i_{2}=1}^{m_{2}}\left(\sum_{i_{1}=1}^{m_{1 }}K_{i_{1}}^{(1)}XW_{1i_{1}}W_{2i_{1}i_{2}}\right)_{+}W_{3i_{2}}\right)_{+}W_{ 4}\right)_{+}W_{5}\cdots\right)_{+}W_{n+1}\geq 1.$ $$

We construct $\{b_{i_{1}i_{2}}^{(1)},{b^{\prime}}_{i_{1}i_{2}}^{(1)}\}$ such that $b_{i_{1}i_{2}}^{(1)}=W_{1i_{1}}W_{2i_{1}i_{2}}$ when $W_{2i_{1}i_{2}}\geq 0$ and 0 0 otherwise, ${b^{\prime}}_{i_{1}i_{2}}^{(1)}=-W_{1i_{1}}W_{2i_{1}i_{2}}$ when $W_{2i_{1}i_{2}}<0$ and 0 0 otherwise. Therefore we get

$$ $y\odot\left(\cdots\left(\left(\sum_{i_{2}=1}^{m_{2}}\left(\sum_{i_{1}=1}^{m_{1 }}K_{i_{1}}^{(1)}X\left(b_{i_{1}i_{2}}^{(1)}-{b^{\prime}}_{i_{1}i_{2}}^{(1)} \right)\right)_{+}W_{3i_{2}}\right)_{+}W_{4}\right)_{+}W_{5}\cdots\right)_{+}W _{n+1}\geq 1,$ $$

with constraints $(2K_{i_{1}}^{(1)}-I)Xb_{i_{1}i_{2}}^{(1)}\geq 0,(2K_{i_{1}}^{(1)}-I)X{b^{
\prime}}_{i_{1}i_{2}}^{(1)}\geq 0.$ Let $K_{i_{2}}^{(2)}$ denote ${\mathrm{sign}}(\sum_{i=1}^{m_{1}}K_{i_{1}}^{(1)}X(b_{i_{1}i_{2}}^{(1)}-{b^{
\prime}}_{i_{1}i_{2}}^{(1)}))$, i.e., $(2K_{i_{2}}^{(2)}-I)(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X(b_{i_{1}i_{2}}^{(1
)}-{b^{\prime}}_{i_{1}i_{2}}^{(1)}))\geq 0$, we thus have

$$ $y\odot\left(\cdots\left(\left(\sum_{i_{2}=1}^{m_{2}}K_{i_{2}}^{(2)}\left(\sum_ {i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X\left(b_{i_{1}i_{2}}^{(1)}-{b^{\prime}}_{i_{1 }i_{2}}^{(1)}\right)\right)W_{3i_{2}}\right)_{+}W_{4}\right)_{+}W_{5}\cdots \right)_{+}W_{n+1}\geq 1,$ $$

with constraints

$$ $(2K_{i_{1}}^{(1)}-I)Xb_{i_{1}i_{2}}^{(1)}\geq 0,(2K_{i_{1}}^{(1)}-I)X{b^{ \prime}}_{i_{1}i_{2}}^{(1)}\geq 0,(2K_{i_{2}}^{(2)}-I)\left(\sum_{i_{1}=1}^{m_ {1}}K_{i_{1}}^{(1)}X(b_{i_{1}i_{2}}^{(1)}-{b^{\prime}}_{i_{1}i_{2}}^{(1)}) \right)\geq 0.$ $$

Expand one more hidden layer

$$ $y\odot\left(\cdots\left(\sum_{i_{3}=1}^{m_{3}}\left(\sum_{i_{2}=1}^{m_{2}}K_{i _{2}}^{(2)}\left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X\left(b_{i_{1}i_{2}}^{( 1)}-{b^{\prime}}_{i_{1}i_{2}}^{(1)}\right)\right)W_{3i_{2}i_{3}}\right)_{+}W_{ 4i_{3}}\right)_{+}W_{5}\cdots\right)_{+}W_{n+1}\geq 1.$ $$

Construct $\{b_{i_{1}i_{2}i_{3}}^{(11)},{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(11)}\}$ by setting $b_{i_{1}i_{2}i_{3}}^{(11)}=b_{i_{1}i_{2}}^{(1)}W_{3i_{2}i_{3}}$ and ${b^{\prime}}_{i_{1}i_{2}i_{3}}^{(11)}={b^{\prime}}_{i_{1}i_{2}}^{(1)}W_{3i_{2}
i_{3}}$ when $W_{3i_{2}i_{3}}\geq 0$ and 0 0 otherwise. Construct $\{b_{i_{1}i_{2}i_{3}}^{(12)},{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(12)}\}$ by setting $b_{i_{1}i_{2}i_{3}}^{(12)}=-b_{i_{1}i_{2}}^{(1)}W_{3i_{2}i_{3}}$ and ${b^{\prime}}_{i_{1}i_{2}i_{3}}^{(12)}=-{b^{\prime}}_{i_{1}i_{2}}^{(1)}W_{3i_{2
}i_{3}}$ when $W_{3i_{2}i_{3}}<0$ and 0 0 otherwise. Let $K_{i_{3}}^{(3)}$ denotes ${\mathrm{sign}}(\sum_{i_{2}=1}^{m_{2}}K_{i_{2}}^{(2)}(\sum_{i_{1}=1}^{m_{1}}K_
{i_{1}}^{(1)}X(b_{i_{1}i_{2}}^{(1)}-{b^{\prime}}_{i_{1}i_{2}}^{(1)}))W_{3i_{2}
i_{3}}).$ Thus we have

$$ $\displaystyle y\odot\left(\cdots\left(\sum_{i_{3}=1}^{m_{3}}K_{i_{3}}^{(3)} \left(\sum_{i_{2}=1}^{m_{2}}K_{i_{2}}^{(2)}\left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1 }}^{(1)}X\left(b_{i_{1}i_{2}i_{3}}^{(11)}-{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(11) }\right)\right)\right.\right.\right.$ $\displaystyle\left.\left.\left.-\sum_{i_{2}=1}^{m_{2}}K_{i_{2}}^{(2)}\left( \sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X\left(b_{i_{1}i_{2}i_{3}}^{(12)}-{b^{ \prime}}_{i_{1}i_{2}i_{3}}^{(12)}\right)\right)\right)W_{4i_{3}}\right)_{+}W_{ 5}\cdots\right)_{+}W_{n+1}\geq 1,$ $$

with constraints

$$ $\displaystyle(2K_{i_{1}}^{(1)}-I)Xb_{i_{1}i_{2}i_{3}}^{(11)}\geq 0,(2K_{i_{1}} ^{(1)}-I)X{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(11)}\geq 0$ $\displaystyle(2K_{i_{1}}^{(1)}-I)Xb_{i_{1}i_{2}i_{3}}^{(12)}\geq 0,(2K_{i_{1}} ^{(1)}-I)X{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(12)}\geq 0,$ $\displaystyle(2K_{i_{2}}^{(2)}-I)\left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X \left(b_{i_{1}i_{2}i_{3}}^{(11)}-{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(11)}\right) \right)\geq 0,(2K_{i_{2}}^{(2)}-I)\left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X \left(b_{i_{1}i_{2}i_{3}}^{(12)}-{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(12)}\right) \right)\geq 0,$ $\displaystyle(2K_{i_{3}}^{(3)}-I)\left(\sum_{i_{2}=1}^{m_{2}}K_{i_{2}}^{(2)} \left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X\left(b_{i_{1}i_{2}i_{3}}^{(11)}-{ b^{\prime}}_{i_{1}i_{2}i_{3}}^{(11)}\right)\right)-\sum_{i_{2}=1}^{m_{2}}K_{i_ {2}}^{(2)}\left(\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X\left(b_{i_{1}i_{2}i_{3} }^{(12)}-{b^{\prime}}_{i_{1}i_{2}i_{3}}^{(12)}\right)\right)\right)\geq 0.$ $$

For cleanness, we introduce the following notation, for any $c_{i}\in\{1,2\}$ and $2\leq s\leq n-1,$

$$ $\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{s}}^{(n-1)(n-2)\cdots(s)}(K^{(s)})=\sum_{ i_{s}=1}^{m_{s}}K_{i_{s}}^{(s)}\left(\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{s}[c _{s-1}=1]}^{(n-1)(n-2)\cdots(s)(s-1)}(K^{(s-1)})-\mathcal{T}_{c_{n-1}c_{n-2} \cdots c_{s}[c_{s-1}=2]}^{(n-1)(n-2)\cdots(s)(s-1)}(K^{(s-1)})\right).$ $$

When $s=1$,

$$ $\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{2}c_{1}}^{(n-1)(n-2)\cdots(2)(1)}(K^{(1)} )=\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}X\left(b_{i_{n}i_{n-1}\cdots i_{1}}^{(c _{n-1}c_{n-2}\cdots c_{1})}-{b^{\prime}}_{i_{n}i_{n-1}\cdots i_{1}}^{(c_{n-1}c _{n-2}\cdots c_{1})}\right).$ $$

Proceed with the above splitting, under the newly defined notation, we will get

$$ $y\odot\sum_{i_{n}=1}^{m_{n}}K_{i_{n}}^{(n)}\left(\mathcal{T}_{1}^{(n-1)}(K^{(n -1)})-\mathcal{T}_{2}^{(n-1)}(K^{(n-1)})\right)\geq 1$ $$

with constraints

$$ $\displaystyle(2K_{i_{s}}^{(s)}-I)\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{s}c_{s-1 }}^{(n-1)(n-2)\cdots(s)(s-1)}(K^{(s-1)})\geq 0,\forall 2\leq s\leq n-1,$ $\displaystyle(2K_{i_{1}}^{(1)}-I)Xb_{i_{n}i_{n-1}\cdots i_{1}}^{(c_{n-1}c_{n-2 }\cdots c_{1})}\geq 0,(2K_{i_{1}}^{(1)}-I)X{b^{\prime}}_{i_{n}i_{n-1}\cdots i_ {1}}^{(c_{n-1}c_{n-2}\cdots c_{1})}\geq 0.$ $$

Note we already have the form of (6) by combine $\{b_{i_{n}i_{n-1}\cdots i_{1}}^{(c_{n-1}c_{n-2}\cdots c_{1})}\}$ and $\{{b^{\prime}}_{i_{n}i_{n-1}\cdots i_{1}}^{(c_{n-1}c_{n-2}\cdots c_{1})}\}$ into $\{b_{i_{n}i_{n-1}\cdots i_{1}}^{(c_{n}c_{n-1}c_{n-2}\cdots c_{1})}\}$, the only thing left is to match $\{K^{(s)}\}$ with $\{D^{(n)}\}$, we do this by recursion. Consider any layer $l$, assume all $\{K^{(s)}\},s\leq l$ can be matched with $\{D^{(n)}\},n\leq l.$ Now we consider the $(l+1)$-th layer. Note that all $K_{i_{l+1}}^{(l+1)}\in\{D^{(l+1)}\}$ for any $i_{l+1}\in[m_{l+1}].$ If there is duplicate neuron activation patterns, i.e., $K_{a}^{(l+1)}=K_{b}^{(l+1)}$ for some $a\neq b,a,b\in[m_{l+1}].$ Then we merge all lower-level neurons corresponding to $K_{a}^{(l+1)}$ and $K_{b}^{(l+1)}$ by summing up the corresponding (with respect to $i_{1},i_{2},\ldots,i_{l}$ indices) $b$ vectors. Both the layer output and ReLU sign constraints will be preserved for all layers up to $(l+1)$-th layer, and we thus get, after the merging, a new set of $\{\tilde{K}^{(l+1)}\}$ and $\{\tilde{b}\}$ that matches the problem (6) up to layer $(l+1)$, where we just set all parameters to be zero for any $j_{l+1}\in[P_{l+1}]$ such that $D_{j_{l+1}}\not\in\{\tilde{K}^{(l+1)}\}.$ Now, we only need to verify that the last inner layer’s neuron can be matched. By the symmetry between $b_{i_{n}i_{n-1}\cdots i_{1}}^{(c_{n-1}c_{n-2}\cdots c_{1})}$ terms, consider without loss of generality the neuron

$$ $\sum_{i_{1}=1}^{m_{1}}K_{i_{1}}^{(1)}Xb_{i_{n}i_{n-1}\cdots i_{1}}^{[c_{n-1}=1 ][c_{n-2}=1]\cdots[c_{1}=1]},$ (14) $$

which we want to match with

$$ $\sum_{j_{1}=1}^{P_{1}}D_{j_{1}}^{(1)}Xu_{j_{n}j_{n-1}\cdots j_{1}}^{[c_{n}=1][ c_{n-1}=1][c_{n-2}=1]\cdots[c_{1}=1]}.$ (15) $$

Note by construction we know $K_{i_{1}}^{(1)}\in\{D^{(1)}\}$ for any $i_{1}\in[m_{1}].$ If there is any duplicate neurons $K_{c}^{(1)}=K_{d}^{(1)}$ for some $c\neq d,c,d\in[m_{1}].$ We denote $K_{(c,d)}^{(1)}=K_{c}^{(1)}=K_{d}^{(1)}.$ Let $b_{i_{n}i_{n-1}\cdots[i_{1}\leftarrow(c,d)]}^{11\cdots 1}=b_{i_{n}i_{n-1}
\cdots[i_{1}=c]}^{11\cdots 1}+b_{i_{n}i_{n-1}\cdots[i_{1}=d]}^{11\cdots 1}.$ Then we merge $K_{c}^{(1)}$ and $K_{d}^{(1)}$ by replacing them with only one copy of $K_{(c,d)}^{(1)}$, and set the corresponding $b$ vector to be $b_{i_{n}i_{n-1}\cdots[i_{1}=(c,d)]}$. Continue this process until there is no duplicate neurons in the last inner layer. We are guaranteed to get a set of $\{\tilde{K}_{i_{1}}^{(1)}\}$ and corresponding $\{\tilde{b}_{i_{n}i_{n-1}\cdots i_{1}}^{11\cdots 1}\}$ such that $\tilde{K}_{i_{1}}^{(1)}$ belong to $\{D^{(1)}\}$ and are all different. Now assume all outer layers have already been merged by the scheme of summing up corresponding $b$ vectors mentioned above. Using $\{\hat{m}_{s},\hat{K}^{(s)},\hat{b}\}$ to represent the new set of parameters. Then $i_{s}\in[\hat{m}_{s}]$ with $\hat{m}_{s}\leq P_{s}$. The expressions (14) and (15) can be matched by setting $u_{j_{n}j_{n-1}\cdots j_{1}}^{111\cdots 1}=\hat{b}_{i_{n}i_{n-1}\cdots i_{1}}^
{11\cdots 1}$ for $(j_{s},i_{s})$ pairs satisfying $D_{j_{s}}^{(s)}=\hat{K}_{i_{s}}^{(s)}$, and setting to zero if $D_{j_{s}}^{(s)}\not\in[\hat{K}^{(s)}].$ This completes our first direction proof.

We now assume that there exists $\{u_{j_{n}j_{n-1}\cdots j_{1}}^{c_{n}c_{n-1}\cdots c_{1}}\}$ which satisfies problem (6), our goal is to find $(n+1)$-layer NN weights $W_{1},W_{2},\cdots,W_{n+1}$ satisfying (5). Since our $\{u_{j_{n}j_{n-1}\cdots j_{1}}^{c_{n}c_{n-1}\cdots c_{1}}\}$ satisfies all plane arrangement constraints in (6), we thus have, based on the inner most layer’s activation pattern constraint and splits $\{u_{j_{n}j_{n-1}\cdots j_{1}}^{c_{n}c_{n-1}\cdots c_{1}}\}$ into $u_{j_{n}j_{n-1}\cdots j_{1}}^{c_{n-1}c_{n-2}\cdots c_{1}}:=u_{j_{n}j_{n-1}
\cdots j_{1}}^{[c_{n}=1]c_{n-1}c_{n-2}\cdots c_{1}},{u^{\prime}}_{j_{n}j_{n-1}
\cdots j_{1}}^{c_{n-1}c_{n-2}\cdots c_{1}}:=u_{j_{n}j_{n-1}\cdots j_{1}}^{[c_{
n}=2]c_{n-1}c_{n-2}\cdots c_{1}}$,

$$ $\displaystyle y\odot\sum_{j_{n}=1}^{P_{n}}D_{j_{n}}^{(n)}\left(\sum_{j_{n-1}=1 }^{P_{n-1}}D_{j_{n-1}}^{(n-1)}\left(\cdots\sum_{j_{1}=1}^{P_{1}}(Xu_{j_{n} \cdots j_{1}}^{1\cdots 1})_{+}-(X{u^{\prime}}_{j_{n}\cdots j_{1}}^{1\cdots 1}) _{+}\cdots\right)-\right.$ $\displaystyle\left.\sum_{j_{n-1}=1}^{P_{n-1}}D_{j_{n-1}}^{(n-1)}\left(\cdots \sum_{j_{1}=1}^{P_{1}}(Xu_{j_{n}\cdots j_{1}}^{2\cdots 1})_{+}-(X{u^{\prime}}_ {j_{n}\cdots j_{1}}^{2\cdots 1})_{+}\cdots\right)\right)\geq 1.$ $$

Thus we can find some $v_{j_{n}j_{n-1}\cdots j_{2}j_{1}}\in\mathbb{R}^{d},S_{j_{1}}^{(1)},S_{j_{2}}^{
(2)},\cdots,S_{j_{n}}^{(n)}\in\{-1,1\}$ such that

$$ $\displaystyle\sum_{j_{n}=1}^{P_{n}}$ $\displaystyle D_{j_{n}}^{(n)}\left(\mathcal{T}_{1}^{(n-1)}(D^{(n-1)})-\mathcal {T}_{2}^{(n-1)}(D^{(n-1)})\right)\qquad\qquad\qquad\qquad$ (16) $\displaystyle=\sum_{j_{n}=1}^{2P_{n}}\left(\sum_{j_{n-1}=1}^{2P_{n-1}}\left( \sum_{j_{n-2}=1}^{2P_{n-2}}\left(\cdots\sum_{j_{2}=1}^{2P_{2}}\left(\sum_{j_{1 }=1}^{2P_{1}}(Xv_{j_{n}j_{n-1}\cdots j_{1}})_{+}S_{j_{1}}^{(1)}\right)_{+} \right.\right.\right.\qquad\qquad$ $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\left.\left.\left.S_{j_{2}}^{ (2)}\cdots\right)_{+}S_{j_{n-2}}^{(n-2)}\right)_{+}S_{j_{n-1}}^{(n-1)}\right)_ {+}S_{j_{n}}^{(n)}.$ $$

Though the process of finding $\{v_{j_{n}j_{n-1}\cdots j_{1}},S_{j_{i}}^{(i)}\}$ has been outlined in three-layer case proof above and can be extended to $(n+1)$-layer case, we present here an outline of finding such $\{v_{j_{n}j_{n-1}\cdots j_{1}},S_{j_{i}}^{(i)}\}$ for five-layer NN for demonstration. For $n=5$, we have the following

$$ $\displaystyle\sum_{j_{4}=1}^{P_{4}}\left(\left(\sum_{j_{3}=1}^{P_{3}}\left( \sum_{j_{2}=1}^{P_{2}}\left(\sum_{j_{1}=1}^{P_{1}}(Xa_{j_{4}j_{3}j_{2}j_{1}}^{ 111})_{+}-(X{a^{\prime}}_{j_{4}j_{3}j_{2}j_{1}}^{111})_{+}\right)_{+}-\left( \sum_{j_{1}=1}^{P_{1}}(Xa_{j_{4}j_{3}j_{2}j_{1}}^{112})_{+}-(Xa_{j_{4}j_{3}j_{ 2}j_{1}}^{112})_{+}\right)_{+}\right)_{+}\right.\right.$ $\displaystyle\left.\left.-\left(\sum_{j_{2}=1}^{P_{2}}\left(\sum_{j_{1}=1}^{P_ {1}}(Xa_{j_{4}j_{3}j_{2}j_{1}}^{121})_{+}-(X{a^{\prime}}_{j_{4}j_{3}j_{2}j_{1} }^{121})_{+}\right)_{+}-\left(\sum_{j_{1}=1}^{P_{1}}(Xa_{j_{4}j_{3}j_{2}j_{1}} ^{122})_{+}-(X{a^{\prime}}_{j_{4}j_{3}j_{2}j_{1}}^{122})_{+}\right)_{+}\right) _{+}\right)_{+}\right.$ $\displaystyle\left.-\left(\sum_{j_{3}=1}^{P_{3}}\left(\sum_{j_{2}=1}^{P_{2}} \left(\sum_{j_{1}=1}^{P_{1}}(Xa_{j_{4}j_{3}j_{2}j_{1}}^{211})_{+}-(X{a^{\prime }}_{j_{4}j_{3}j_{2}j_{1}}^{211})_{+}\right)_{+}-\left(\sum_{j_{1}=1}^{P_{1}}( Xa_{j_{4}j_{3}j_{2}j_{1}}^{212})_{+}-(X{a^{\prime}}_{j_{4}j_{3}j_{2}j_{1}}^{21 2})_{+}\right)_{+}\right)_{+}\right.\right.$ $\displaystyle\left.\left.-\left(\sum_{j_{2}=1}^{P_{2}}\left(\sum_{j_{1}=1}^{P_ {1}}(Xa_{j_{4}j_{3}j_{2}j_{1}}^{221})_{+}-(X{a^{\prime}}_{j_{4}j_{3}j_{2}j_{1} }^{221})_{+}\right)_{+}-\left(\sum_{j_{1}=1}^{P_{1}}(Xa_{j_{4}j_{3}j_{2}j_{1}} ^{222})_{+}-(X{a^{\prime}}_{j_{4}j_{3}j_{2}j_{1}}^{222})_{+}\right)_{+}\right) _{+}\right)_{+}\right).$ $$

Let $v_{j_{4}j_{3}j_{2}j_{1}}^{111}=\{a_{j_{4}j_{3}j_{2}j_{1}}^{111}\}\cup\{{a^{
\prime}}_{j_{4}j_{3}j_{2}j_{1}}^{111}\}$ and similarly construct $v_{j_{4}j_{3}j_{2}j_{1}}^{112},v_{j_{4}j_{3}j_{2}j_{1}}^{121},v_{j_{4}j_{3}j_{
2}j_{1}}^{122},$ $v_{j_{4}j_{3}j_{2}j_{1}}^{211},v_{j_{4}j_{3}j_{2}j_{1}}^{212},$ $v_{j_{4}j_{3}j_{2}j_{1}}^{221},v_{j_{4}j_{3}j_{2}j_{1}}^{222}.$ Let $S_{j_{1}}^{111},S_{j_{1}}^{112},S_{j_{1}}^{121},S_{j_{1}}^{122},S_{j_{1}}^{211
},S_{j_{1}}^{212},S_{j_{1}}^{221},S_{j_{1}}^{222}$ to be $1$ for $j_{1}\in[P_{1}]$ and to be $-1$ for $j_{1}\in[P_{1}+1,2P_{1}]$. Thus, the above expression is equivalent to

$$ $\displaystyle\sum_{j_{4}=1}^{P_{4}}\left(\left(\sum_{j_{3}=1}^{P_{3}}\left( \sum_{j_{2}=1}^{P_{2}}\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^ {111})_{+}S_{j_{1}}^{111}\right)_{+}-\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_ {3}j_{2}j_{1}}^{112})_{+}S_{j_{1}}^{112}\right)_{+}\right)_{+}\right.\right.$ $\displaystyle\left.\left.-\left(\sum_{j_{2}=1}^{P_{2}}\left(\sum_{j_{1}=1}^{2P _{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{121})_{+}S_{j_{1}}^{121}\right)_{+}-\left( \sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{122})_{+}S_{j_{1}}^{122} \right)_{+}\right)_{+}\right)_{+}\right.$ $\displaystyle\left.-\left(\sum_{j_{3}=1}^{P_{3}}\left(\sum_{j_{2}=1}^{P_{2}} \left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{211})_{+}S_{j_{1}}^{2 11}\right)_{+}-\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{212})_ {+}S_{j_{1}}^{212}\right)_{+}\right)_{+}\right.\right.$ $\displaystyle\left.\left.-\left(\sum_{j_{2}=1}^{P_{2}}\left(\sum_{j_{1}=1}^{2P _{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{221})_{+}S_{j_{1}}^{221}\right)_{+}-\left( \sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{222})_{+}S_{j_{1}}^{222} \right)_{+}\right)_{+}\right)_{+}\right).$ $$

Let $v_{j_{4}j_{3}j_{2}j_{1}}^{11}=\{v_{j_{4}j_{3}j_{2}j_{1}}^{111}\}\cup\{v_{j_{4}
j_{3}j_{2}j_{1}}^{112}\}$ and $\tilde{S}_{j_{1}}^{11}=S_{j_{1}}^{111}.$ Let $S_{j_{2}}^{11}=1$ for $j_{2}\in[P_{2}]$ and $S_{j_{2}}^{11}=-1$ for $j_{2}\in[P_{2}+1,2P_{2}].$ Similarly construct $v_{j_{4}j_{3}j_{2}j_{1}}^{12},\tilde{S}_{j_{1}}^{12},S_{j_{2}}^{12}$, $v_{j_{4}j_{3}j_{2}j_{1}}^{21},\tilde{S}_{j_{1}}^{21},S_{j_{2}}^{21}$, $v_{j_{4}j_{3}j_{2}j_{1}}^{22},\tilde{S}_{j_{1}}^{22},S_{j_{2}}^{22}$. Thus the above expression is equivalent to

$$ $\displaystyle\sum_{j_{4}=1}^{P_{4}}\left(\left(\sum_{j_{3}=1}^{P_{3}}\left( \sum_{j_{2}=1}^{2P_{2}}\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}} ^{11})_{+}\tilde{S}_{j_{1}}^{11}\right)_{+}S_{j_{2}}^{11}\right)_{+}-\left( \sum_{j_{2}=1}^{2P_{2}}\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}} ^{12})_{+}\tilde{S}_{j_{1}}^{12}\right)_{+}S_{j_{2}}^{12}\right)_{+}\right)_{+ }\right.$ $\displaystyle\left.\left(\sum_{j_{3}=1}^{P_{3}}\left(\sum_{j_{2}=1}^{2P_{2}} \left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{21})_{+}\tilde{S}_{j_ {1}}^{21}\right)_{+}S_{j_{2}}^{21}\right)_{+}-\left(\sum_{j_{2}=1}^{2P_{2}} \left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{22})_{+}\tilde{S}_{j_ {1}}^{22}\right)_{+}S_{j_{2}}^{22}\right)_{+}\right)_{+}\right).$ $$

Let $v_{j_{4}j_{3}j_{2}j_{1}}^{1}=\{v_{j_{4}j_{3}j_{2}j_{1}}^{11}\}\cup\{v_{j_{4}j_
{3}j_{2}j_{1}}^{12}\}$, $v_{j_{4}j_{3}j_{2}j_{1}}^{2}=\{v_{j_{4}j_{3}j_{2}j_{1}}^{21}\}\cup\{v_{j_{4}j_
{3}j_{2}j_{1}}^{22}\}.$ Let ${\hat{\tilde{S}}}_{j_{1}}^{(1)}=\tilde{S}_{j_{1}}^{11},$ ${\hat{\tilde{S}}}_{j_{1}}^{(2)}=\tilde{S}_{j_{1}}^{21}.$ Let $\hat{S}_{j_{2}}^{(1)}=S_{j_{2}}^{11},\hat{S}_{j_{2}}^{(2)}=S_{j_{2}}^{21}.$ Let further $S_{j_{3}}^{(1)},S_{j_{3}}^{(2)}$ to take value $1$ for $j_{3}\in[P_{3}]$ and take value $-1$ for $j_{3}\in[P_{3}+1,2P_{3}].$ Therefore the above expression is equivalent to

$$ $\displaystyle\sum_{j_{4}=1}^{P_{4}}\left(\left(\sum_{j_{3}=1}^{2P_{3}}\left( \sum_{j_{2}=1}^{2P_{2}}\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}} ^{1})_{+}{\hat{\tilde{S}}}_{j_{1}}^{(1)}\right)_{+}{\hat{S}_{j_{2}}}^{(1)} \right)_{+}S_{j_{3}}^{(1)}\right)_{+}\right.$ $\displaystyle\left.-\left(\sum_{j_{3}=1}^{2P_{3}}\left(\sum_{j_{2}=1}^{2P_{2}} \left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{2})_{+}{\hat{\tilde{S }}}_{j_{1}}^{(2)}\right)_{+}\hat{S}_{j_{2}}^{(2)}\right)_{+}S_{j_{3}}^{(2)} \right)_{+}\right).$ $$

We once again repeat the above procedure for the outer most layer and we arrive

$$ $\displaystyle\sum_{j_{4}=1}^{P_{4}}\left(\sum_{j_{3}=1}^{2P_{3}}\left(\sum_{j_ {2}=1}^{2P_{2}}\left(\sum_{j_{1}=1}^{2P_{1}}(Xv_{j_{4}j_{3}j_{2}j_{1}}^{1})_{+ }{\hat{\tilde{S}}}_{j_{1}}^{(1)}\right)_{+}{\hat{S}_{j_{2}}}^{(1)}\right)_{+}S _{j_{3}}^{(1)}\right)_{+}S^{\prime}_{j_{4}},$ $$

where $S_{j_{4}}^{\prime}=1$ for $j_{4}\in[P_{4}]$ and $-1$ otherwise. This completes our construction of $\{v_{j_{n}j_{n-1}\cdots j_{1}},S_{j_{i}}^{(i)}\}$ for $n=5.$ Now, we are left with matching the following two expressions

$$ $\displaystyle\sum_{j_{n}=1}^{2P_{n}}\left(\sum_{j_{n-1}=1}^{2P_{n-1}}\left( \sum_{j_{n-2}=1}^{2P_{n-2}}\left(\cdots\sum_{j_{2}=1}^{2P_{2}}\left(\sum_{j_{1 }=1}^{2P_{1}}(Xv_{j_{n}j_{n-1}\cdots j_{1}})_{+}S_{j_{1}}^{(1)}\right)_{+}S_{j _{2}}^{(2)}\cdots\right)_{+}S_{j_{n-2}}^{(n-2)}\right)_{+}S_{j_{n-1}}^{(n-1)} \right)_{+}S_{j_{n}}^{(n)},$ $$

and

$$ $\displaystyle\sum_{i_{n}=1}^{m_{n}}\left(\sum_{i_{n-1}=1}^{m_{n-1}}\left(\sum_ {i_{n-2}=1}^{m_{n-2}}\left(\cdots\sum_{i_{2}=1}^{m_{2}}\left(\sum_{i_{1}=1}^{m _{1}}(XW_{1i_{1}})_{+}W_{2i_{2}i_{1}}\right)_{+}W_{3i_{3}i_{2}}\cdots\right)_{ +}W_{(n-1)i_{n-1}i_{n-2}}\right)_{+}W_{ni_{n}i_{n-1}}\right)_{+}W_{(n+1)i_{n}}.$ $$

This can be done by setting all weights corresponding to indices $\{i_{n}>2P_{n},i_{n-1}>4P_{n}P_{n-1},\cdots,i_{k}>\Pi_{c=k}^{n}2P_{c}\}$ to be all zeros, then for the outer most layer, we set $W_{(n+1)i_{n}}=S_{i_{n}}^{(n)},$ for $2\leq k\leq n$, we set $W_{ki_{k}i_{k-1}}=S_{((i_{k-1}-1)\%2P_{k-1})+1}^{(k-1)}$ for $i_{k-1}\in[(i_{k}-1)*2P_{k-1}+1,i_{k}*2P_{k-1}]$ and 0 0 otherwise. For the inner most layer, we set $W_{1i_{1}}=v_{j_{n}j_{n-1}\cdots j_{1}}$ with

$$ $\displaystyle j_{n}$ $\displaystyle=\lfloor(i_{1}-1)/\Pi_{k=1}^{n-1}2P_{k}\rfloor+1,$ $\displaystyle j_{n-1}$ $\displaystyle=\lfloor((i_{1}-1)\%\Pi_{k=1}^{n-1}2P_{k})/\Pi_{k=1}^{n-2}2P_{k} \rfloor+1,$ $\displaystyle j_{n-2}$ $\displaystyle=\lfloor(((i_{1}-1)\%\Pi_{k=1}^{n-1}2P_{k})\%\Pi_{k=1}^{n-2}2P_{k })/\Pi_{k=1}^{n-3}2P_{k}\rfloor+1,$ $\displaystyle\cdots$ $\displaystyle j_{2}$ $\displaystyle=\lfloor((\cdots)\%\Pi_{k=1}^{2}2P_{k})/2P_{1}\rfloor+1,$ $\displaystyle j_{1}$ $\displaystyle=((\cdots)\%\Pi_{k=1}^{2}2P_{k})\%2P_{1}+1.$ $$

Given any test point $\tilde{x}\in\mathbb{R}^{d}$, the final prediction can be computed by

$$ $\tilde{y}=\sum_{j_{n}=1}^{P_{n}}\left(\mathcal{T}_{1}^{(n-1)}(D^{(n-1)})\right )_{+}-\left(\mathcal{T}_{2}^{(n-1)}(D^{(n-1)})\right)_{+},$ (17) $$

where

$$ $\displaystyle\mathcal{T}_{c_{n-1}\cdots c_{i}}^{(n-1)\cdots(i)}(D^{(i)})=\sum_ {j_{i}=1}^{P_{i}}\left(\mathcal{T}_{c_{n-1}\cdots c_{i}1}^{(n-1)\cdots(i)(i-1) }(D^{(i-1)})\right)_{+}-\left(\mathcal{T}_{c_{n-1}\cdots c_{i}2}^{(n-1)\cdots( i)(i-1)}(D^{(i-1)})\right)_{+},$ $\displaystyle\mathcal{T}_{c_{n-1}c_{n-2}\cdots c_{1}}^{(n-1)(n-2)\cdots(1)}(D^ {(1)})=\sum_{j_{1}=1}^{P_{1}}\left(\tilde{x}^{T}a_{j_{n}j_{n-1}\cdots j_{1}}^{ 1c_{n-1}\cdots c_{1}}\right)_{+}-\left(\tilde{x}^{T}a_{j_{n}j_{n-1}\cdots j_{1 }}^{2c_{n-1}\cdots c_{1}}\right)_{+}.$ $$

#### E.3 Explanation of Bias Term for General Case

To show that the $n$-layer NN of form $((\cdots((XW_{1})_{+}W_{2})_{+}W_{3}\cdots)_{+}W_{n-1})_{+}W_{n}$ preserves the same approximation capacity as the biased version $((\cdots((XW_{1}+b_{1})_{+}W_{2}+b_{2})_{+}W_{3}+b_{3}\cdots)_{+}W_{n-1}+b_{n-
1})_{+}W_{n}+b_{n}$, we show that the biased version is incorporated in the form $((\cdots((XW_{1})_{+}W_{2})_{+}W_{3}\cdots)_{+}W_{n-1})_{+}W_{n}$ when the constraint on number of hidden neurons is mild. First note that we can always append the data matrix $X$ with a column of ones to incorporate the inner most bias $b_{1}.$ Then for each outer layer, we can always have an inner neuron to be a pure bias neuron with value one. Then the corresponding outer neuron weight would serve as an outer layer bias.

#### E.4 Activation Pattern Subsampling and Iterative Filtering.

Here we detail more about our hyperplane selection scheme. Take two-layer ReLU model for example, in order to find the activation pattern $\mathcal{D}$ corresponding to the hidden layer, one needs to exhaust the set $\{\text{diag}({\mathbbm{1}}\{Xu\geq 0\})\}$ for all $u\in\mathbb{R}^{d}.$ In our experiments, we adopt a heuristic subsampling procedure, i.e., we usually set a moderate number of hidden neurons $m_{1}$, and we sample a set of Gaussian random vectors $\{u_{1},u_{2},\cdots,u_{n_{1}}\}$ with some random $n_{1}>m_{1}.$ Then we take the set $\{\text{diag}({\mathbbm{1}}\{Xu_{i}\geq 0\})|i\in[n_{1}]\}.$ If $|\{\text{diag}({\mathbbm{1}}\{Xu_{i}\geq 0\})|i\in[n_{1}]\}|>m_{1},$ we take a subset of $m_{1}$ activation patterns there, if not, we increase $n_{1}$ and redo all prior steps until we hit some satisfactory $n_{1}$. This heuristic method always works well in our experiments, see Appendix G.1 for more about our implementation details.

A more rigorous way which exhausts all possible activation patterns can be done via an iterative filtering procedure. We demonstrate here for a toy example. Consider still two-layer model as before, when we are given a data set $X$ of size $n$, a loose upper bound on $|\mathcal{D}|$ is given by $2^{n}$, i.e., each piece of data can take either positive and negative values and they are all independent. Thus one can find all possible $D_{i}$’s by solving

$$ $\displaystyle\underset{}{\text{find}}\qquad\sum_{i}^{2^{n}}\|u_{i}\|_{2}$ $\displaystyle\text{s.t.}\qquad(2D_{i}-I)Xu_{i}\geq 0,\|u_{i}\|_{2}\leq 1$ $$

where $D_{i}$’s loop over all $2^{n}$ possibilities. Then the $D_{i}$’s correspond to non-zero $u_{i}$’s in the solution are feasible plane arrangements. However, this method induces $2^{n}*d$ variables. A more economic way to find all feasible arrangements is to do an iterative filtering with each newly added data. When there is only one non-zero data point $X_{1}$, there always exists $u_{1},u_{2}$ vectors such that $X_{1}^{T}u_{1}>0$ and $X_{1}^{T}u_{2}<0.$ Thus $\mathcal{D}_{1}=\{\begin{bmatrix}1\end{bmatrix},\begin{bmatrix}0\end{bmatrix}\}$ represent all possible sign patterns for this single training data. After the second data point has been added, we know

$$ $\mathcal{D}_{2}\subseteq\left\{\begin{bmatrix}1&0\\ 0&1\end{bmatrix},\begin{bmatrix}1&0\\ 0&0\end{bmatrix},\begin{bmatrix}0&0\\ 0&1\end{bmatrix},\begin{bmatrix}0&0\\ 0&0\end{bmatrix}\right\}.$ $$

Therefore, an upper bound on cardinality of $\mathcal{D}_{2}$ is given by $2^{2}=4.$ However, this upper bound might be pessimistic, for example, if $X_{1}=X_{2}$, then we would expect

$$ $\mathcal{D}_{2}=\left\{\begin{bmatrix}1&0\\ 0&1\end{bmatrix},\begin{bmatrix}0&0\\ 0&0\end{bmatrix}\right\}\subset\left\{\begin{bmatrix}1&0\\ 0&1\end{bmatrix},\begin{bmatrix}1&0\\ 0&0\end{bmatrix},\begin{bmatrix}0&0\\ 0&1\end{bmatrix},\begin{bmatrix}0&0\\ 0&0\end{bmatrix}\right\}.$ $$

Since $X_{1}^{T}u=X_{2}^{T}u$ for any $u$, they will always have the same sign pattern. Things get more complicated with more data points, say for $n$ data points $X_{1},X_{2},\cdots,X_{n}$, it is possible that sign pattern of $X_{n}^{T}u$ can be determined by sign patterns of $X_{1}^{T}u,X_{2}^{T}u,\cdots,X_{n-1}^{T}u$ when there are linear dependency between the data points, which happens more often with larger set of training data.
Thus the true cardinality of $\mathcal{D}_{n}$ might be far smaller than $2^{n}.$ Indeed, one can show that, with $r$ denoting rank of the training data matrix consisting of the first $n$ data points (Pilanci & Ergen, ; Stanley et al., ),

$$ $|\mathcal{D}_{n}|\leq 2r\left(\frac{e(n-1)}{r}\right)^{r},$ $$

which can be much smaller than our pessimistic bound $2^{n}$ especially when training data has small rank. To stay close to the optimal cardinality and avoid solving an optimization problem with $2^{n}*d$ number of variables, what one can do is to find the feasible plane arrangements iteratively at the time when each single data point is added. The underline logic is that plane arrangement patterns which are infeasible to $X_{[1:t-1]}$ will also be infeasible to $X_{[1:t]}$. Therefore, assume one has already found the optimal sign patterns $\mathcal{D}_{t-1}$ for the first $t-1$ training samples which has cardinality $c_{t-1}$, when $X_{t}$ arrives, one needs to solve the following auxiliary problem

$$ $\displaystyle\underset{}{\text{find}}\qquad\sum_{i}^{2c_{t-1}}\|u_{i}\|_{2}$ $\displaystyle\text{s.t.}\qquad\left(2\begin{bmatrix}D_{(t-1)i}&0\\ 0&0\end{bmatrix}-I\right)Xu_{i}\geq 0,i\in[1,\cdots,c_{t-1}]$ $\displaystyle\qquad\quad\left(2\begin{bmatrix}D_{(t-1)i-c_{t-1}}&0\\ 0&1\end{bmatrix}-I\right)Xu_{i}\geq 0,i\in[c_{t-1}+1,\cdots,2c_{t-1}]$ $\displaystyle~{}~{}\qquad\quad\|u_{i}\|_{2}\leq 1,$ $$

which has only $(2*c_{t-1})*d$ number of variables and can be far fewer than $2^{t}*d$. This scheme can also be done lazily each fixed $T$ iterations, one just add all $\{0,1\}$-patterns for the last $T$ data points.

### Appendix F Deferred Proofs and Extensions in Section 6

#### F.1 Deferred Proof of Theorem 6.1

###### Proof.

Notice that the polyhedron cut is simply three consecutive cuts with three hyperplanes:

- •
$\mathcal{H}_{1}:=\{\theta:y_{n}\cdot f^{\text{two-layer}}(x_{n}^{T};\theta)=0\}$;
- •
$\mathcal{H}_{2}:=\{\theta^{e}:(2D_{i}-I)_{n}x_{n}^{T}\theta^{e}=0\}$;
- •
$\mathcal{H}_{3}:=\{\theta^{o}:(2D_{i}-I)_{n}x_{n}^{T}\theta^{o}=0\}$,

where $\theta^{o}$ ($\theta^{e}$) denotes the reduced $\theta=(u^{\prime}_{1},u_{1},...,u^{\prime}_{p},u_{p})$ vector containing only the odd (even) indices.
Since the cuts imposed by the linear inequality constraints $\mathcal{C}_{n}$ and $\mathcal{C}^{\prime}_{n}$ via hyper-planes $\mathcal{H}_{2}$, $\mathcal{H}_{3}$ only reduce the remaining set $\mathcal{T}$, examining the volume remained after just the cut via $\mathcal{H}_{1}$ suffices for the convergence analysis as it bounds $\textbf{vol}(\mathcal{T}_{1})$ from above.

Assuming the cut is active, it follows that $\theta_{G}$ misclassifies $(x_{n},y_{n})$, so $y_{n}\cdot f^{\text{two-layer}}(x_{n};\theta)<0$. Thus, the cut is deep and $\theta_{G}$ is in the interior of $\mathcal{T}_{2}$. By Proposition 1, it follows that

$$ $\displaystyle\textbf{vol}(\mathcal{T}\cap\mathcal{H}^{+}_{1})<\textbf{vol}( \mathcal{T}\cap\mathcal{H}^{+}_{G})\leq(1-1/e)\textbf{vol}(\mathcal{T}),$ $$

where $\mathcal{H}_{G}$ is any hyperplane that goes through $\theta_{G}$ which is parallel to $\mathcal{H}_{1}$. Therefore, at each step $t$ where the cut is active, at least volume of magnitude $\frac{1}{e}\textbf{vol}(\mathcal{T}^{t})$ is cut away. The volume of $\mathcal{T}^{t}$ after $t$ iterations is bounded by:

$$ $\displaystyle\textbf{vol}(\mathcal{T}^{t})<(1-1/e)^{t}\cdot\textbf{vol}( \mathcal{T}^{0}).$ $$

Then as $t\rightarrow\infty$, $(1-1/e)^{t}\rightarrow 0$, it follows that $\textbf{vol}(\mathcal{T}^{t})$ converges to zero and the feasible region shrinks to point(s).

Notice that since $\mathcal{T}$ is a convex body and $\{\mathcal{T}^{t}\}$ is a nested decreasing sequence of convex set with $\mathcal{T}^{t+1}\subseteq\mathcal{T}^{t}$, by the finite intersection property, the intersection of all $\mathcal{T}^{t}$ sets is non-empty and contains the optimal solution $\mathcal{\theta^{*}}$:

$$ $\displaystyle\cap_{t=0}^{\infty}\mathcal{T}^{t}=\{\theta^{*}\}.$ $$

It remains to justify that what the intersection contains is indeed the optimal solution. To see this, observe that the problem that Algorithm 4 is simply a feasibility problem. Since every time the convex set shrinks by intersecting with the constraint set containing feasibility criteria given each new acquired data point $(x_{n_{t}},y_{n_{t}})$, i.e.,

$$ $\{\theta\in\mathcal{T}:y_{n}\cdot f^{\text{two-layer}}(x_{n};\theta)\geq 0, \mathcal{C}(\{n\}),\mathcal{C}^{\prime}(\{n\})\},$ $$

the set shrinks to a finer set that satisfies increasingly more constraints posed by additional acquired data points. This monotonic improvement with the shrinking feasible region implies that the sequence of classifiers $\theta^{t}$ converges to the set of optimal classifiers $\{\theta^{*}\}$.
∎

###### Proof.

#### F.2 Extension to Algorithm 4 : The Case of Inexact Cuts and Convergence

In this section, we discuss convergence results of Algorithm 5 under the third setup described in Section 6.1.

Note that the cut made by both Algorithm 4 and Algorithm 2 are exact in the sense that cuts are only made when the computed center mis-classifies the data points returned by the oracle and is therefore always discarded via the cutting hyperplanes. Algorithm 5, on the other hand, implements inexact cut: cuts are always made regardless of whether the queried center mis-classifies. Nice convergence rate similar to that of Theorem 6.1 can still be guaranteed, as we shall see in the following theorem, as long as the cut is made within a certain neighborhood of the centroid at each step.

First, let us quantify the inexactness of the cut. At each iteration, given a cutting hyperplane $\mathcal{H}_{a}$ of normal vector $a$, which passes through the origin by design, and the computed center $\theta$, the Euclidean distance between $\theta$ and the hyperplane is given by

$$ $\displaystyle h=\|\theta\|_{2}|\cos(\alpha)|,$ (18) $$

where

$$ $\displaystyle\alpha=\arccos(\frac{\theta\cdot a}{\|\theta\|_{2}\|a\|_{2}})$ (19) $$

is the angle between $\theta$ and the normal vector $a$. For simple notation, let us denote the angle between $\theta$ and the cutting plane $\mathcal{H}_{a}$ as $\beta$. Then

$$ $\displaystyle\beta=\frac{\pi}{2}-\alpha.$ (20) $$

Next, we introduce an important extension to Proposition 1 given in Louche & Ralaivola ().

###### Proposition 2 (Generalized Partition of Convex bodies) .

Let $\mathcal{T}\subseteq\mathbb{R}^{d}$ be a convex body. Let $\mathcal{H}_{a}$ be a hyperplane of normal vector $a$. We define the positive (negative) halfspace $\mathcal{T}^{+}$ (resp. $\mathcal{T}^{-}$) of $\mathcal{T}$ with respect to $\mathcal{H}_{a}$ as

$$ $\mathcal{T}^{+}:=\mathcal{T}\cap\{\theta\in\mathbb{R}^{d}:\langle a,\theta \rangle\geq 0\}$ $$

$$ $\mathcal{T}^{-}:=\mathcal{T}\cap\{\theta\in\mathbb{R}^{d}:\langle a,\theta \rangle<0\}$ $$

The following holds true: if $\theta_{G}+\Lambda a\in\mathcal{T}^{+}$ then

$$ $\textbf{vol}(\mathcal{T}^{+})/\textbf{vol}(\mathcal{T})\geq e^{-1}(1-\lambda)^ {d},$ $$

where

$$ $\Lambda=\lambda\Theta_{d}\frac{\textbf{vol}(\mathcal{T})H_{\mathcal{T}^{+}}}{R ^{d}H_{\mathcal{T}^{-}}},$ $$

with $\lambda\in\mathbb{R}$ an arbitrary real such that $\lambda\leq 1$, $\Theta_{d}$ a constant depending only on $d$, $R$ the radius of the $(d-1)$-dimensional ball $\mathcal{B}_{2}$ of volume $\textbf{vol}(\mathcal{B}_{2}):=\textbf{vol}(\mathcal{T}\cap\{\theta\in\mathbb{
R}^{d}:\langle a,\theta\rangle=0\})$ and

$$ $H_{\mathcal{T}^{+}}:=\max_{b\in\mathcal{T}^{+}}b^{T}a\quad(\text{resp. }H_{ \mathcal{T}^{-}}:=\min_{b\in\mathcal{T}^{-}}b^{T}a).$ $$

###### Proof.

See Proof of Theorem 2 in Louche & Ralaivola ().
∎

Note that by symmetry of partitioned convex body with respect to the centroid (an application of Proposition 1), one can similarly establish a bound in the case that

$$ $\displaystyle\theta_{G}+\Lambda a\in\mathcal{T}^{-},$ $$

where $\Lambda$ is defined as in Proposition 2. Then

$$ $\displaystyle\textbf{vol}(\mathcal{T}^{-})/\textbf{vol}(\mathcal{T})\geq e^{-1 }(1-\lambda)^{d}.$ (21) $$

Proposition 2 generalizes Grünbaum’s inequality in Proposition 1 by allowing a cut that is of distance (greater than or equal to) $\Lambda$ along the direction of the normal vector $a$ from the actual center of gravity. To apply Proposition 2, we need to establish an equivalence relation between our quantifiation of the inexactness through the Euclidean distance $h$ (or angle $\alpha$ and or $\beta$) with $\Lambda$.

We are now ready to bound the volume reduction factor for inexact cuts with respect to the center gravity in Algorithm 5.

###### Theorem F.1 (Convergence with Inexact Cuts of Center of Gravity) .

Let $\mathcal{T}\subseteq\mathbb{R}^{d}$ be a convex body and let $\theta_{G}$ denote its center of gravity. Given oracle returned data point $(x_{n},y_{n})$, define hyperplane

$$ $\displaystyle\mathcal{H}_{a}:=\{\theta\in\mathbb{R}^{d}:y_{n}\cdot f^{\text{ two-layer}}(x_{n};\theta)=0\},$ $$

where

$$ $\displaystyle a:=y_{n}\cdot[x_{n}^{1}-x_{n}^{1}\ ...\ x_{n}^{P}-x_{n}^{P}]$ $$

is its associated normal vector.
Then $\mathcal{H}_{a}$ divides $\mathcal{T}$ into a positive half-space $\mathcal{T}^{+}$ and a negative half-space $\mathcal{T}^{-}$,

$$ $\displaystyle\mathcal{T}^{+}:=\mathcal{T}\cap\{\theta\in\mathbb{R}^{d}:\langle a ,\theta\rangle\geq 0\}$ $\displaystyle\mathcal{T}^{-}:=\mathcal{T}\cap\{\theta\in\mathbb{R}^{d}:\langle a ,\theta\rangle<0\},$ $$

and $h:=\frac{|\theta_{G}^{T}a|}{\|a\|_{2}}$ is the Euclidean distance between $\theta_{G}$ and $\mathcal{H}_{a}$.
The inexact polyhedron cut given in Algorithm 5
partitions the convex body $\mathcal{T}$ into two subsets:

$$ $\displaystyle\mathcal{T}_{1}:=\{\theta\in\mathcal{T}:\mathcal{T}^{+},\mathcal{ C}(\{n\}),\mathcal{C}^{\prime}(\{n\})\}$ $\displaystyle\mathcal{T}_{2}:=\{\theta\in\mathcal{T}:\mathcal{T}^{-},\text{ or }\neg\mathcal{C}(\{n\}),\neg\mathcal{C}^{\prime}(\{n\})\},$ $$

where $\neg$ denotes the complement of a given set.
The following holds:

$$ $\displaystyle\textbf{vol}(\mathcal{T}_{1})<(1-e^{-1}(1-\tilde{\lambda})^{d}) \textbf{vol}(\mathcal{T}),$ $$

with

$$ $\displaystyle\tilde{\lambda}=-\frac{R^{d}H_{\mathcal{T}^{-}}}{\Theta_{d} \textbf{vol}(\mathcal{T})H_{\mathcal{T}^{+}}}h>0,$ $$

where
$\Theta_{d}$ is a constant depending only on $d$, $R$ is the radius of the $(d-1)$-dimensional ball $\mathcal{B}_{2}$ of volume $\textbf{vol}(\mathcal{B}_{2}):=\textbf{vol}(\mathcal{T}\cap\mathcal{H}_{a})$ and

$$ $H_{\mathcal{T}^{+}}:=\max_{\theta\in\mathcal{T}^{+}}\theta^{T}a\quad(\text{ resp. }H_{\mathcal{T}^{-}}:=\min_{\theta\in\mathcal{T}^{-}}\theta^{T}a).$ $$

###### Proof.

Observe that at each iteration if $\theta_{G}$ mis-classifies the oracle returned data point $(x_{n},y_{n})$, then we have the same analysis as in Theorem 6.1. In such a case, we have guaranteed accelerated cuts with a strictly better volume reduction rate bound than in the case when cuts are made exactly through the centroid. This is because in this case the centroid is contained in the interior of the eliminated region. On the other hand, when the cut is made such that the centroid is contained in the interior of the positive half-space, we cut away smaller volume than when we make the cut through the centroid. Therefore, to obtain a lower bound on the volume reduction rate for Algorithm 5, it suffices to examine only the latter case, i.e. when $\theta_{G}$ is contained in interior of the positive half-space.

Suppose that $\theta_{G}\in\mathcal{T}^{+}\setminus\mathcal{H}_{a}$. Define real number $\Lambda$ to be such that

$$ $\displaystyle\Lambda<-h:=-\frac{|\theta^{T}a|}{\|a\|_{2}},$ $$

then since $\theta_{G}-ha\in\mathcal{H}_{a}$, it follows that

$$ $\displaystyle\theta_{G}+\Lambda a\in\mathcal{T}^{-}.$ $$

By Proposition 2 and symmetry,
it follows that

$$ $\textbf{vol}(\mathcal{T}^{-})/\textbf{vol}(\mathcal{T})\geq e^{-1}(1-\lambda)^ {d},$ $$

where

$$ $\lambda=\Lambda(\Theta_{d}\frac{\textbf{vol}(\mathcal{T})H_{\mathcal{T}^{+}}}{ R^{d}H_{\mathcal{T}^{-}}})^{-1}$ $$

is a real such that $\lambda\leq 1$, $\Theta_{d}$ is a constant depending only on $d$, $R$ is the radius of the $(d-1)$-dimensional ball $\mathcal{B}_{2}$ of volume $\textbf{vol}(\mathcal{B}_{2}):=\textbf{vol}(\mathcal{T}\cap\mathcal{H}_{a})$ and

$$ $H_{\mathcal{T}^{+}}:=\max_{\theta\in\mathcal{T}^{+}}\theta^{T}a\quad(\text{ resp. }H_{\mathcal{T}^{-}}:=\min_{\theta\in\mathcal{T}^{-}}\theta^{T}a).$ $$

Since

$$ $\displaystyle\Lambda=\lambda\Theta_{d}\frac{\textbf{vol}(\mathcal{T})H_{ \mathcal{T}^{+}}}{R^{d}H_{\mathcal{T}^{-}}}<-h,$ $$

or equivalently

$$ $\displaystyle 0<\lambda<-\frac{R^{d}H_{\mathcal{T}^{-}}}{\Theta_{d}\textbf{vol }(\mathcal{T})H_{\mathcal{T}^{+}}}h,$ $$

it follows that

$$ $\displaystyle\textbf{vol}(\mathcal{T}^{-})\geq e^{-1}(1-\lambda)^{d}\textbf{ vol}(\mathcal{T})>e^{-1}(1-\tilde{\lambda})^{d}\textbf{vol}(\mathcal{T}),$ $$

or equivalently

$$ $\displaystyle\textbf{vol}(\mathcal{T}^{+})<(1-e^{-1}(1-\tilde{\lambda})^{d}) \textbf{vol}(\mathcal{T}),$ $$

where

$$ $\displaystyle\tilde{\lambda}=-\frac{R^{d}H_{\mathcal{T}^{-}}}{\Theta_{d} \textbf{vol}(\mathcal{T})H_{\mathcal{T}^{+}}}h.$ $$

Since

$$ $\displaystyle\mathcal{T}_{1}:=\{\theta\in\mathcal{T}:\mathcal{T}^{+},\mathcal{ C}(\{n\}),\mathcal{C}^{\prime}(\{n\})\}\subseteq\mathcal{T}^{+},$ $$

we have that

$$ $\displaystyle\textbf{vol}(\mathcal{T}_{1})\leq\textbf{vol}(\mathcal{T}^{+})<(1 -e^{-1}(1-\tilde{\lambda})^{d})\textbf{vol}(\mathcal{T}).$ $$

To justify that the derived ratio $\textbf{vol}(\mathcal{T}_{1})/\textbf{vol}(\mathcal{T})$ is valid, i.e. $\textbf{vol}(\mathcal{T}_{1})/\textbf{vol}(\mathcal{T})\in[0,1]$, it suffices to show that $e^{-1}(1-\tilde{\lambda})^{d}\in[0,1]$ for all given $d>0$, which necessarily holds if $0<\tilde{\lambda}\leq 1$.

To see that $0<\tilde{\lambda}\leq 1$, we first define coefficient

$$ $\displaystyle c(\mathcal{H}_{a},\mathcal{T},d):=\Theta_{d}\frac{\textbf{vol}( \mathcal{T})H_{\mathcal{T}^{+}}}{R^{d}H_{\mathcal{T}^{-}}}.$ $$

Since by definition $H_{\mathcal{T}^{+}}>0$, $H_{\mathcal{T}^{-}}<0$, with $\Theta_{d},\textbf{vol}(\mathcal{T})$, $R^{d}>0$, it follows that $c(\mathcal{H}_{a},\mathcal{T},d)<0$.

Since $\tilde{\lambda}:=-c(\mathcal{H}_{a},\mathcal{T},d)^{-1}h$ and $h>0$, we have that $\tilde{\lambda}>0$. Next, to see $\tilde{\lambda}\leq 1$, observe that since we have

$$ $\displaystyle\Lambda=\lambda c(\mathcal{H}_{a},\mathcal{T},d)<-h,$ $$

it follows that

$$ $\displaystyle c(\mathcal{H}_{a},\mathcal{T},d)<-h\leq-\lambda^{-1}h,$ $$

where we have used the fact that $0<\lambda\leq 1$.
It follows that

$$ $\displaystyle\tilde{\lambda}:=-c(\mathcal{H}_{a},\mathcal{T},d)^{-1}h<h^{-1}h=1.$ $$

Hence, $0<\tilde{\lambda}<1$. Then by Theorem 6.1, the volume reduction ratio $1-e^{-1}(1-\tilde{\lambda})^{d}>1-e^{-1}$ offers a strict upper bound for each iteration of Algorithm 5, regardless of whether the centroid $\theta_{G}$ is contained in the positive half-space or the negative half-space.
∎

Given Theorem F.1, at each step $t$ in Algorithm 5, at least a volume of magnitude $e^{-1}(1-\tilde{\lambda})^{d}\textbf{vol}(\mathcal{T}^{t})$ is cut away. The volume of $\mathcal{T}^{t}$ after $t$ iterations is bounded by:

$$ $\displaystyle\textbf{vol}(\mathcal{T}^{t})<(1-e^{-1}(1-\tilde{\lambda})^{d})^{ t}\cdot\textbf{vol}(\mathcal{T}^{0}).$ $$

Since $0<\tilde{\lambda}<1$, we have $0<1-e^{-1}(1-\tilde{\lambda})^{d}\leq 1$. So as $t\rightarrow\infty$, $(1-e^{-1}(1-\tilde{\lambda})^{d})^{t}\rightarrow 0$. It follows that $\textbf{vol}(\mathcal{T}^{t})$ converges to zero as $t\rightarrow\infty$, and the feasible region shrinks to the optimal point(s), as in Theorem 6.1.

Theorem F.1 allows us to quantify the volume reduction rate as a function of the inexactness of cuts through the Euclidean distance between the centroid and the cutting hyperplane, which is monitored by the query sampling method. In particular, in the case of inexact cuts in Algorithm 5 with a minimal margin query sampling scheme, the volume reduction rate converges asymptotically in the number of training data supplied to the cutting-plane oracle to the rate in Proposition 1. This result is stated in the following corollary.

###### Corollary 1 (Asymptotic Convergence of Inexact Cuts with Center of Gravity) .

Assume that data points $x_{i}$’s are uniformly distributed across input space. Let the number of training data supplied to the cutting-pane oracle in Algorithm 5 goes to infinity, i.e. $|\mathcal{D}|\rightarrow\infty$. Then the volume of $\mathcal{T}^{t}$ after $t$ iterations with the polyhedron cut in Algorithm 5 under minimal margin query sampling method converges asymptotically to bound

$$ $\displaystyle\textbf{vol}(\mathcal{T}^{t})<(1-e^{-1})^{t}\cdot\textbf{vol}( \mathcal{T}^{0}).$ $$

###### Proof.

Let us begin with analyzing a single iteration with convex body $\mathcal{T}\subset\mathbb{R}^{d}$ and its centroid $\theta_{G}$.
Notice that as $|\mathcal{D}|$ increases, data points $x_{i}$’s become more densely distributed across the input space. So for any $\epsilon>0$, there exists a sufficiently large $N$ such that for all $|\mathcal{D}|>N$, there exists at least one data point $x_{i}\in\mathcal{D}$, for which $|f^{\text{two-layer}}(x_{i};\theta_{G}):=[x_{i}^{1}\ -x_{i}^{1}\ ...\ x_{i}^{P
}\ -x_{i}^{P}]\cdot\theta_{G}|<\epsilon$. As $|\mathcal{D}|\rightarrow\infty$, $\epsilon$ can be made arbitrarily small, implying that

$$ $\displaystyle\operatorname*{arg\,min}_{x_{i}\in\mathcal{D}_{x}}|f^{\text{two- layer}}(x_{i};\theta_{G})|\rightarrow 0.$ $$

This suggests that the Euclidean distance between $\theta_{G}$ and the oracle returned cutting plane $\mathcal{H}_{a}$ diminishes to 0 0 as $|\mathcal{D}|\rightarrow\infty$.

By Theorem F.1, given an inexact polyhedron cut in Algorithm 5 distance $h$ away from the centroid, the volume of the remaining set is upper bounded by the following:

$$ $\displaystyle\textbf{vol}(\mathcal{T}_{1})\leq(1-e^{-1}(1-\tilde{\lambda})^{d} )\textbf{vol}(\mathcal{T}),$ $$

where $\tilde{\lambda}:=-c(\mathcal{H}_{a},\mathcal{T},d)^{-1}h$ and

$$ $\displaystyle c(\mathcal{H}_{a},\mathcal{T},d):=\Theta_{d}\frac{\textbf{vol}( \mathcal{T})H_{\mathcal{T}^{+}}}{R^{d}H_{\mathcal{T}^{-}}}.$ $$

Then to show that ratio $(1-e^{-1}(1-\tilde{\lambda})^{d})\rightarrow 1-e^{-1}$, or equivalently $\tilde{\lambda}\rightarrow 0$, as $h\rightarrow 0$, it suffices to show that the coefficient term $-c(\mathcal{H}_{a},\mathcal{T},d)=|c(\mathcal{H}_{a},\mathcal{T},d)|$ is lower bounded by a constant $M_{\text{min}}$ so that its inverse is upper bounded. This is easy to show. First observe that dimension dependent constant $\Theta_{d}$ and $\textbf{vol}(\mathcal{T})$ is fixed and finite regardless of the cutting plane $\mathcal{H}_{a}$. Because convex body $\mathcal{T}$ is by definition closed, variables $R,|H_{\mathcal{T}^{-}}|$ are bounded from above. By the same reasoning, $H_{\mathcal{T}^{+}}$ is bounded from below and away from 0 0, for if $H_{\mathcal{T}^{+}}=0$, it follows that $\textbf{vol}(\mathcal{T}^{+})=0$, suggesting the termination and convergence of the algorithm to the optimal solution(s). Hence, there exists a constant $M_{\text{min}}$ such that

$$ $\displaystyle\min_{\mathcal{H}_{a_{i}}}|c(\mathcal{H}_{a_{i}},\mathcal{T},d)| \geq M_{\text{min}}.$ $$

It follows that

$$ $\displaystyle 0\leq\lim_{h\rightarrow 0}\tilde{\lambda}:=\lim_{h\rightarrow 0} -c(\mathcal{H}_{a},\mathcal{T},d)^{-1}h\leq\lim_{h\rightarrow 0}M_{\text{min}} ^{-1}h=0.$ $$

So as $|\mathcal{D}|\rightarrow\infty$, we have that $h\rightarrow 0$, which results in $\tilde{\lambda}\rightarrow 0$. Therefore, the volume of $\mathcal{T}^{t}$ after $t$ iterations in Algorithm 5 follows:

$$ $\displaystyle\lim_{|\mathcal{D}|\rightarrow\infty}\textbf{vol}(\mathcal{T}^{t} )<\lim_{|\mathcal{D}|\rightarrow\infty}(1-e^{-1}(1-\tilde{\lambda})^{d})^{t} \cdot\textbf{vol}(\mathcal{T}^{0})=(1-e^{-1})^{t}\cdot\textbf{vol}(\mathcal{T} ^{0}).$ $$

∎

###### Proposition 2 (Generalized Partition of Convex bodies) .

###### Proof.

###### Theorem F.1 (Convergence with Inexact Cuts of Center of Gravity) .

###### Proof.

###### Corollary 1 (Asymptotic Convergence of Inexact Cuts with Center of Gravity) .

###### Proof.

#### F.3 Convergence w.r.t. Center of the Maximum Volume Ellipsoid

Recall the definition for the center of MVE in Definition 4. We will now show that similar convergence rate as that of the center of gravity can be achieved under center of MVE as well.

###### Theorem F.2 (Convergence with Center of MVE) .

Let $\mathcal{T}\subseteq\mathbb{R}^{d}$ be a convex body and let $\theta_{M}$ denote its center of the maximum volume inscribed ellipsoid. The polyhedron cut given in Algorithm 4 (assuming that the cut is active) and Algorithm 2, i.e.,

$$ $\displaystyle\mathcal{T}\cap\{\theta:y_{n}\cdot f^{\text{two-layer}}(x_{n}; \theta)\geq 0,\mathcal{C}(\{n\}),\mathcal{C}^{\prime}(\{n\})\},$ $$

where coupling $(x_{n},y_{n})$ is the data point returned by the cutting-plane oracle after receiving queried point $\theta_{M}$, partitions the convex body $\mathcal{T}$ into two subsets as in Theorem 6.1. Then $\mathcal{T}_{1}$ satisfies the following inequality:

$$ $\displaystyle\textbf{vol}(\mathcal{T}_{2})<(1-\frac{1}{d})\cdot\textbf{vol}( \mathcal{T}).$ $$

###### Proof.

The proof follows similarly as the proof to Theorem 6.1. For the MVE cutting-plane method, the bound on the volume reduction factor is given by:

$$ $\displaystyle\frac{\textbf{vol}(\mathcal{T}^{k+1})}{\textbf{vol}(\mathcal{T}^{ k})}\leq 1-\frac{1}{d},$ $$

where $k$ denotes the iteration (Boyd & Vandenberghe, ).
Since $f^{\text{two-layer}}(x_{n};\theta)$ is linear in $\theta$ and by the design of the algorithm, the center $\theta_{M}$ is always contained in the interior of the discarded set (or equivalently, the cut is deep), proof for the bound $(1-\frac{1}{d})\cdot\textbf{vol}(\mathcal{T})$ follows similarly as the proof of Theorem 6.1.
∎

###### Theorem F.2 (Convergence with Center of MVE) .

###### Proof.

### Appendix G Experiments Supplementals

#### G.1 Implementation Details

In this section, we provide an overview of the implementation details. For specific aspects, such as baseline implementation and the cutting-plane AL method for regression, we direct readers to the corresponding sections in Appendix G.

In our experiments, we follow exactly the same algorithm workflow as in Algorithm 1 for our model training and Algorithm 4 for active learning with limited queries, which is the version we used conducting all of our experiments in Section 7. In the training of Algorithm 4, we implement the “center” function for analytic center retrieval due to its simple computation formula (see Definition 1). Since the center retrieval problem is of convex minimization form, we solve it with CVXPY (Diamond & Boyd, ) and default to MOSEK (MOSEK ApS, ) as our solver. In our experiments, MOSEK is able to handle all encountered convex optimization programs, given that they are feasible. We involve all training data cuts while we note that some cut dropping method may help alleviate computation overhead, i.e., one may keep the last several cuts only, see Section 2.5 in (Parshakova et al., ) for a discussion on constraint dropping. For activation pattern generation, for two-layer model experiment, say we want $P=m$ patterns, we always generate a set of standard Gaussian vectors $\{u_{i},i\in[n]\}$ for some $n>m$, compute the induced pattern set $\{\mathcal{D}_{i}=\text{diag}({\mathbbm{1}}\{Xu_{i}\geq 0\})\}$, and take $m$ non-duplicate patterns out of it. Furthermore, in our implementation of Algorithm 4, we included a few additional checks for computational efficiency, like skipping the centering step when no cut has been performed before, or book-keeping of discarded data points that should not be re-considered immediately.
For IMDB experiments, we randomly pick $50$ training data points and $20$ test data points for ease of computation, and we use minimal margin query strategy described in Appendix C.3. Through all our data experiments, random seed is fixed to be 0 0 except for error bar plots, where we always take random seeds $\{0,1,2,3,4\}$.

#### G.2 Final Solve Regularization in Cutting-plane AL

In this section, we discuss an optional feature to our cutting-plane AL method via the two-or-three layer ReLU network: the inclusion of a final convex solve which solves for the optimal parameter of the equivalent convex program to the two-or-three layer ReLU network using the data acquired by the AL thus far.

###### Convexifying a Two-layer ReLU Network.

To introduce this final convex solve, we first briefly overview the work of Pilanci & Ergen (). We will focus on the two-layer ReLU network case, as the case for three-layer can be easily extended by referencing the exact convex reformulation of the three-layer ReLU network given in the paper by Ergen & Pilanci ().

Pilanci & Ergen () first introduced a finite dimensional, polynomial-size convex program that globally solves the training problem for two-layer fully connected ReLU networks.
The convexification of ReLU networks can be summarized into a two-step strategy:

- 1.
Project original feature to higher dimensions using convolutional hyperplane arrangements.
- 2.
Convexify using convex regularizers, such as convex variable selection models.

First, as per Step 1), we briefly recall and define a notion of hyperplane arrangements for the neural network in Definition 3 and we can rewrite the ReLU constraint using the partitioned regions by the hyperplanes:

$$ $\displaystyle(2D-I_{n})Xw\geq 0$ (22) $$

Next, per Step 2), we define the primal problem and give its exact finite-dimensional convex formulation. Given data matrix $X\in\mathbb{R}^{n\times d}$, we consider a two-layer network $f(X;\theta):\mathbb{R}^{n\times d}\rightarrow\mathbb{R}^{n}$ with $m$ neurons:

$$ $\displaystyle f(X;\theta)=\sum_{j=1}^{m}(Xw_{j})_{+}\alpha_{j},$ (23) $$

where $w_{j}\in\mathbb{R}^{d}$, $\alpha_{j}\in\mathbb{R}$ are weights for hidden and output layers respectively, and $\theta=\{w_{j},\alpha_{j}\}_{j=1}^{m}$. Supplied additionally with a label vector $y\in\mathbb{R}^{n}$ and a regularization parameter $\beta>0$, the primal problem of the 2-layer fully connected ReLU network is the following:

$$ $\displaystyle p^{*}:=\min_{\{w_{j},\alpha_{j}\}_{j=1}^{m}}\frac{1}{2}\|f(X; \theta)-y\|^{2}_{2}+\frac{\beta}{2}\sum_{j=1}^{m}(\|w_{j}\|^{2}_{2}+\alpha_{j} ^{2}).$ (24) $$

The above non-convex objective is transformed into an equivalent convex program.
For a detailed derivation, please refer to Pilanci & Ergen ().
We will give an overview of the key steps here.
First, re-represent the above optimization problem as an equivalent $\ell_{1}$ penalized minimization,

$$ $\displaystyle p^{*}=\min_{\|w_{j}\|_{2}\leq 1\forall j\in[m]}\min_{\{\alpha_{j }\}_{j=1}^{m}}\frac{1}{2}\|f(X;\theta)-y\|^{2}_{2}+\beta\sum_{j=1}^{m}|\alpha_ {j}|.$ (25) $$

Using strong duality, the exact semi-infinite convex program to objective 24 is obtained:

$$ $\displaystyle p^{*}=\max_{v\in\mathbb{R}^{n}\text{ s.t. }|v^{T}(Xw)_{+}|\leq \beta\forall w\in\mathcal{B}_{2}}-\frac{1}{2}\|y-v\|^{2}_{2}+\frac{1}{2}\|y\|^ {2}_{2}.$ (26) $$

Under certain regularity conditions in Pilanci & Ergen (), the equivalent finite-dimensional convex formulation is the following:

$$ $\displaystyle\begin{split}p^{*}=&\min_{\{u_{i},{u^{\prime}}_{i}\}_{i=1}^{P},u_ {i},{u^{\prime}}_{i}\in\mathbb{R}^{d}\forall i}\frac{1}{2}\|\sum_{i=1}^{P}D_{i }X({u^{\prime}}_{i}-u_{i})-y\|^{2}_{2}\\ &+\beta(\sum_{i=1}^{P}(\|u_{i}\|_{2}+\|u^{\prime}_{i}\|_{2}))\end{split}$ (27) $$

subject to the constraints that

$$ $(2D_{i}-I_{n})Xu_{i}\geq 0,\ (2D_{i}-I_{n})X{u^{\prime}}_{i}\geq 0,\ \forall i.$ $$

Here, $\beta$ is a regularization parameter, and we denote $\{u_{i}^{*},{u^{\prime}}_{i}^{*}\}_{i=1}^{P}$ as the solution to objective 27.

To see that the convex program outlined in 27 is indeed an exact reformulation of the two-layer ReLU network, we run the convex solve on the spiral dataset and the quadratic regression dataset in Section 7 and compare it with standard stochastic gradient descent technique. We also note here that for the sake of brevity, we will be referring to the solving of the exact convex program with respect to the two-or-three layer ReLU networks interchangeably as “convex solve” or “final solve”.

First, using the same spiral dataset of randomly selected $80$ points from $100$ points Spiral $(k_{1}=13,k_{2}=0.5,n_{\text{shape}}=50)$ (detailed in Appendix G.4) and $\beta=0.001$, the convex solve achieves an objective value of 0.0008, while stochastic gradient descent levels out at an objective value greater than 0.001. The final decision boundary on the spiral using the convex solve (left) and SGD (right) is shown in Figure 5, where we have used marker “x” to indicate the train points passed into each solver and triangle marker to indicate test points.

Figure: Figure 5: Decision boundary made by convex solve (left) and stochastic gradient descent (right) on the Spiral Dataset of $n_{\text{shape}}=50$.
Refer to caption: x1.png

While both methods capture the spiral perfectly in this case, the advantage of the convex solve becomes more evident given more complex (e.g. in terms of shape) data. This is demonstrated in Figure 6, where we have used a spiral dataset of randomly selected $80$ points from $100$ points Spiral $(k_{1}=13,k_{2}=0.5,n_{\text{shape}}=100)$.

Figure: Figure 6: Decision boundary made by convex solve (left) and stochastic gradient descent (right) on the Spiral Dataset of $n_{\text{shape}}=100$.
Refer to caption: x2.png

Similar results hold for the regression case. The convex solve for the two-layer ReLU network returns an objective value of $4.142\times 10^{-5}$. This aligns with the result of Figure 7 as the final prediction of the convex solves aligns perfectly with the actual quadratic function.

Figure: Figure 7: Prediction made by convex solve (red) and stochastic gradient descent (cyan) on the regression dataset.
Refer to caption: x3.png

###### Integration of the Final Convex Solve.

Now we talk about how we can incorporate the convex program in 27 into our cutting-plane AL method to potentially aid its performance. We first provide some justifications.

As we have mentioned in Section 6, since the prediction function

$$ $\displaystyle\sum_{i=1}^{P}D_{i}X({u^{\prime}}_{i}-u_{i}):=f^{\text{two-layer} }(X;\theta)=[X^{1}-X^{1}\ ...\ X^{P}-X^{P}]\theta,$ $$

is linear in $\theta:=(u_{1}^{\prime},u_{1},\ldots,u_{P}^{\prime},u_{P})$ with $u_{i},{u^{\prime}}_{i}\in\mathbb{R}^{d}$ along with the ReLU cuts, we preserve the convexity of the parameter space after each cut. And hence, the final solve becomes well applicable. To introduce this step, we use the short-hand notations in Section 6 and also, for the sake of brevity, we denote the exact convex objective in Equation 27 as the following:

$$ $\displaystyle\begin{split}f^{\text{obj}}(\mathcal{D};\theta,\beta)=&\min_{ \theta}\frac{1}{2}\|\sum_{i=1}^{P}(D_{i}X)_{\mathcal{D}}({u^{\prime}}_{i}-u_{i })-y_{\mathcal{D}}\|^{2}_{2}+\beta(\sum_{i=1}^{P}(\|u_{i}\|_{2}+\|u^{\prime}_{ i}\|_{2})),\end{split}$ (28) $$

where $X_{\mathcal{D}}$ and $y_{\mathcal{D}}$ are the slices of $X$ and $y$ at indices $\mathcal{D}$.

Figure: Algorithm 8 Cutting-plane AL for Binary Classification with Limited Queries using Final Solve

Algorithm 8 shows how the final convex solve is incorporated into our cutting-plane AL method via the two-layer ReLU network with limited queries. Other variations such as query synthesis (Algorithm 2), inexact cuts (Algorithm 5), as well as regression (Algorithm 6) can be adapted in the same way. In Algorithm 8, SOLVE() is a convex solve that solves the exact convex formulation in Equation 27 with all the selected training data pairs $(x_{i},y_{i})\mid i\in\mathcal{D}_{\text{AL}}$ from the active learning loop. In our implementation, we use CVXPY and default to solver CLARABEL (Goulart & Chen, ).

Upon examining Algorithm 8 and the objective function of the equivalent convex program $f^{\text{obj}}(\mathcal{D}_{\text{AL}};\theta^{t},\beta)$ in 27, a key difference is the introduction of a regularization term, i.e. $\beta(\sum_{i=1}^{P}(\|u_{i}\|_{2}+\|u^{\prime}_{i}\|_{2}))$. In more complex tasks, the inclusion of the final convex solve, and thus an additional regularization, tends to aid the cutting-plane AL method in achieving faster convergence in the number of queries. This is the case, for example, for the spiral dataset in our experiment (see Appendix G.4). However, for simpler tasks, such as the quadratic regression prediction, adding regularization could slow down the convergence (see Appendix G.5). Therefore, in our experiments, we emphasize that we have used the best cutting-plane AL method, chosen from the version with or without the final convex solve. This is explicitly noted in the supplementary sections, namely Appendix G.4 for the spiral experiment and Appendix G.5 for the quadratic regression task. We also remark that due to the limitation in theory on the convexification of neural networks, such a convex solve is only viable up to a three-layer ReLU network.

#### G.3 Active Learning Baselines

In this section, we introduce the various baselines used in the spiral task and the quadratic regression task in Section 7. We also elaborate on the implementation details of each baseline discussed.

In addition to random sampling, which serves as a baseline for all other AL methods, and the linear cutting-plane AL method (Algorithm 3 for classification and Algorithm 6 for regression), which provides a baseline for our cutting-plane NN algorithm (Algorithm 4), demonstrating our extension from linear to nonlinear decision boundaries in the cutting-plane AL framework, we also survey popular AL algorithms from the scikit-activeml library (Kottke et al., ) and the DeepAL package (Huang, ).

We first introduce the baselines we implemented from the DeepAL package.

- •
Entropy Sampling (Settles, ). This technique selects samples for labeling based on the entropy of the predicted class probabilities. Recall that the entropy for a sample $x$ with predicted class probabilities $p(y|x)$ is given by:
$H(y|x)=-\sum_{c}p(y=c|x)\log p(y=c|x)$
Higher entropy indicates greater uncertainty, making such samples good candidates for active learning.
- •
Bayesian Active Learning Disagreement (BALD) with Dropout (Gal et al., ). BALD aims to choose samples that maximize the mutual information $I[y,\theta|x,\mathcal{D}_{train}]$
between predictions $y$ and model parameters $\theta$, given a sample $x$ and training data $\mathcal{D}_{train}$. Mathematically, this is expressed as:
$I[y,\theta|x,\mathcal{D}_{train}]=H[y|x,\mathcal{D}_{train}]-\mathbb{E}_{p(
\theta|\mathcal{D}_{train})}[H[y|x,\theta]]$
BALD with dropout extends this approach by using dropout during inference to approximate Bayesian uncertainty, allowing for efficient estimation of uncertainty in deep learning models through Monte Carlo dropout.
- •
Least Confidence (Lewis & Gale, ). This strategy selects the samples where the model is least confident in its most likely prediction. For a given sample $x$, it is measured by the confidence of the predicted class $\hat{y}$, as follows:
$\text{LC}(x)=1-\max_{c}p(y=c|x)$
Samples with lower confidence values are considered more uncertain and thus more informative for labeling.

Next, we introduce the baselines surveyed from scikit-activeml.

- •
Query by Committee (QBC) (Seung et al., ). The Query-by-Committee strategy uses an ensemble of estimators (a “committee”) to identify samples on which there is disagreement. The committee members vote on the label of each sample, and the sample with the highest disagreement is selected for labeling. This disagreement is often quantified using measures like vote entropy:
$H_{vote}(x)=-\sum_{c}\frac{v_{c}}{C}\log\frac{v_{c}}{C}$
where $v_{c}$ is the number of votes for class $c$ and $C$ is the total number of committee members. This strategy focuses on reducing uncertainty by selecting instances where committee members are in conflict.
- •
Greedy Sampling on the Feature Space (GreedyX) (Wu et al., ). GreedyX implements greedy sampling on the feature space, aiming to select samples that increase the diversity of the feature space the most. The method iteratively selects samples that maximize the distance between the newly selected sample and the previously chosen ones, ensuring that the selected subset of samples represents diverse regions of the feature space. This is often formulated as:
$\max_{x_{i}\in\mathcal{D}}\min_{x_{j}\in\mathcal{Q}}\|x_{i}-x_{j}\|^{2}$
where $\mathcal{D}$ is the set of all data points and $\mathcal{Q}$ is the set of already selected query points.
- •
Greedy Sampling on the Target Space (GreedyT) (Wu et al., ). This query strategy initially selects samples to maximize diversity in the feature space and then shifts to maximize diversity in both the feature and target spaces. Alternatively, it can focus solely on the target space (denoted as GSy). This approach attempts to increase the representativeness of the selected samples in terms of both input features and target values. The optimization can be formulated as:
$\max_{x_{i}\in\mathcal{D}}\min_{x_{j}\in\mathcal{Q}}\|x_{i}-x_{j}\|^{2}\quad
\text{and/or}\quad\max_{y_{i}\in\mathcal{D}}\min_{y_{j}\in\mathcal{Q}}\|y_{i}-
y_{j}\|^{2}$
where $\mathcal{D}$ is the set of all data points and $\mathcal{Q}$ is the set of already selected query points, applied either to the feature space, the target space, or both.
- •
KL Divergence Maximization (Kldiv) (Elreedy et al., ). This strategy selects samples that maximize the expected Kullback-Leibler (KL) divergence between the predicted and true distributions of the target values. In this method, it is assumed that the target probabilities for different samples are independent. The KL divergence is a measure of how one probability distribution diverges from a second, reference distribution, and is given by:
$D_{KL}(P\parallel Q)=\sum_{i}P(x_{i})\log\frac{P(x_{i})}{Q(x_{i})}$
where $P$ is the true distribution and $Q$ is the predicted distribution. This method balances exploration and exploitation by focusing on the areas where the current model’s uncertainty is greatest.

We used Query-By-Committee and Greedy Sampling on the Feature Space for both the spiral experiment and the quadratic regression experiment, and used Greedy Sampling on the Target and KL Divergence Maximization only for the regression task. As the readers may see in the descriptions provided above, both KL Divergence Maximization and Greedy Sampling on the Target Space are primarily suited for regression tasks due to the way they handle target distributions and diversity in the target space, which are more relevant in regression scenarios where the output values are continuous. Therefore, we have only included them for the regression task to offer a more robust comparison to our AL method.

Figure: Algorithm 9 Deep Active Learning Baseline

We now discuss the implementation details for each of the aforementioned baselines.
We implement random sampling and all DeepAL baselines using the DeepAL pipeline, while the scikit-activeml baselines are implemented within their respective framework. We note that the original DeepAL framework given by Huang () only implements active learning for classification tasks. Thus, we modified the pipeline to allow AL methods surveyed from DeepAL to also be able to handle regression tasks. Since the modifications are minor, such as changing the loss criterion from cross-entropy to root mean square error loss, but not structural, we omit this distinction here and refer the readers to our submitted codes for details. Algorithm 9 outlines the general workflow of a deep active learning baseline. Here, $\mathcal{D}_{L}$ and $\mathcal{D}_{U}$ refer to the set of labeled and unlabeled training data respectively, and $n_{\text{init}}$ is the number of initial data to be randomly selected and labeled before the active learning loop.

It remains for us to discuss implementation details for scikit-activeml baselines.
It is noteworthy that the default scikit-activeml active learning method for both classification and regression is implemented using its own classifiers or regressors, such as the mixture model classifiers for classification and the normal inverse Chi kernel regressor, instead of trained using customizable deep neural nets, as in the case of DeepAL. For the sake of breadth in surveying popular active learning baselines, we use the Python package skorch, a scikit-learn compatible wrapper for PyTorch models (Viehmann et al., ), to integrate custom neural networks with scikit-activeml classifiers for the classification tasks, and used popular scikit-activeml active learning methods trained with its default regressor, the normal inverse Chi kernel regressor, for the regression tasks. In the classification case, we enforce the same ReLU architecture on all baselines for fairness. And in the regression case, to ensure the optimal performance of the scikit-activeml baselines for the robustness of comparison, we use the so-called “bagging regressor” in scikit-activeml, which is an ensemble method that fits multiple base regressors (in our case, 4) on random subsets of the original dataset using bootstrapping (i.e., sampling with replacement) to update the base normal inverse Chi kernel regressor. This method is known to improve query learning strategies by reducing variance and improving robustness, ensuring that the baselines perform optimally (Abe & Mamitsuka, ).

Further further implementation details, such as how we implement the QBC method in the classification case using Skorch by creating an ensemble of 5 classifiers using the exact same ReLU architecture with parameters in Table 2 but with different random state of 0-4 respectively, and additional details on the specific implementation, please refer to our submitted codes.

Finally, we would like to again emphasize that towards a fair and robust comparison, we select the best performing number of epochs (for a discussion on this, refer to Appendix G.4) and learning rate for AL baselines using deep neural networks and use enhancements such as bagging regressors for the scikit-activeml baselines trained on its default regressors.

#### G.4 Synthetic Spiral Experiment

In this section, we provide additional implementation details and supplementary results on the synthetic spiral binary classification experiment in Section 7.1.

###### Data Generation.

We generate the two intertwined spiral used in Section 7 according to the following: the coordinates $(x_{1},x_{2})$ for the $i$-th data point of the spiral with label $y$ are generated as

$$ $\begin{array}[]{cc}x_{1}=\frac{r\cos(\phi)y}{k_{1}}+k_{2},\quad x_{2}=\frac{r \sin(\phi)y}{k_{1}}+k_{2},\quad r=k_{3}\left(\frac{k_{4}-i}{k_{4}}\right), \quad\phi=k_{5}i\pi,\end{array}$ $$

where $k_{1},\ldots,k_{5}$ are positive coefficients. In the implementation, index $i$ is normalized to lie within the range from 0 0 to $n_{\text{shape}}-1$, where $n_{\text{shape}}$ controls the shape of the spiral, so that the shape is independent of the total number of data points generated. The bigger $n_{\text{shape}}$ is, the more complex the spiral (see for instance the spirals in Figure 5 versus Figure 6). In our experiment, we use the spiral dataset with $k_{1}=13,k_{2}=0.5,n_{\text{shape}}=50$.

###### AL Implementation Details.

To implement the cutting-plane AL method, we use $1000$ simulations to randomly sample the number of partitions in the hyperplane arrangements (Definition 3). This determines the embedding size (or the number of neurons). Table 1 summarizes the statistics with respect to each seed. In our implementation of the deep AL baselines, we adjust the embedding size accordingly when using a different seed to ensure fair comparisons.

**Table 1: Number of neurons for each seed (0-4) sampled using $1000$ simulations for the spiral dataset.**
| Seed | 0 | 1 | 2 | 3 | 4 |
| --- | --- | --- | --- | --- | --- |
| Neurons | 623 | 607 | 605 | 579 | 627 |

For all deep active learning baselines used in the Spiral case, we have set the hyper-parameters to be according to Table 2. We empirically selected learning rate $0.001$ from the choices $\{0.1,0.01,0.001\}$ as it tends to give the best result among the three for all baselines in both the classification and regression tasks. We also choose the number of epochs to be $2000$ as it also gives the best result among choices $\{20,200,2000\}$ and show significant improvement from both $20$ and $200$.

**Table 2: Hyper-parameters of deep AL baselines’ training networks with the Stochastic Gradient Descent (SGD) optimizer**
| Epochs | Learning Rate | Train Batch Size | Test Batch Size | Momentum | Weight Decay |
| --- | --- | --- | --- | --- | --- |
| 2000 | 0.001 | 16 | 10 | 0.9 | 0.003 |

In fact, the performance of the deep AL baselines hinges much on the number of epochs in the train step (see Algorithm 9. Even when the AL method gains access to the full data, if the number of epochs is small, we would still not obtain satisfactory classification result. In the following, we explore the effect of number of epochs on classification accuracy for the baselines. We fix seed = 0 and use the same spiral data generation used in Section 7.1 with a 4:1 train/test split along with a budget of 20 queried points. Figure 8 to Figure 13 show the classification boundary with varying epochs = 20, 200, and 2000.

Figure: (a) epoch = 20
Refer to caption: x4.png

Figure: (a) epoch = 20
Refer to caption: x6.png

Figure: (a) epoch = 20
Refer to caption: x8.png

Figure: (a) epoch = 20
Refer to caption: x10.png

Figure: (a) epoch = 20
Refer to caption: x12.png

Figure: (a) epoch = 20
Refer to caption: x15.png

For most baselines, we observe a significant improvement in the decision boundary between 200 and 2000 epochs. To ensure a robust comparison across methods, we set the number of epochs to 2000. This mini experiment also underscores the inefficiency of gradient-based training, as the accuracy rate heavily depends on the number of training epochs, highlighting the need for a large number of iterations to achieve satisfactory results.

###### Experiment Results.

We now give the deferred supplementary experiment results in Section 7.1. To start, per discussion in Appendix G.2, we would like to discuss the inclusion of regularization through the final convex solve of the equivalent convex program of the two-layer ReLU networks.

Figure 14 shows the final decision boundary with (left) and without (right) the convex solve with the same data and setup as in the spiral experiment in Section 7.1. We abbreviate “after final solve” as AFS and “before final solve” as BFS. It is evident that the inclusion of regularization improves the performance of our cutting-plane active learning method, increasing the accuracy rate on train/test set from 0.84/0.60 to 1.00/1.00. This is therefore the method we use for our cutting-plane AL algorithm. We will see that regularization with the final convex solve does not always speed up the convergence of our cutting-plane AL method in the number of queries. It tends to work less optimally when used on relatively simple dataset, as we will see in Appendix G.5.

Figure: (a) AFS
Refer to caption: x18.png

Table 3 presents the train and test accuracies on the binary spiral dataset. This corroborates the decision boundaries shown in Figure 2.

**Table 3: Train and test accuracies of binary classification on the Spiral ($k_{1}=12,k_{2}=0.5,n_{\text{shape}}=50$) dataset for cutting-plane AL via the 2-layer ReLU NN and various deep AL baselines using seed = 0.**
| Method | Train Accuracy | Test Accuracy |
| --- | --- | --- |
| Cutting-plane (ours) | 1.00 | 1.00 |
| Linear Cutting-plane | 0.50 | 0.50 |
| Random Sampling | 0.73 | 0.70 |
| Entropy Sampling | 0.76 | 0.70 |
| BALD with Dropout | 0.71 | 0.65 |
| Least Confidence Sampling | 0.64 | 0.55 |
| Greedy Sampling (GreedyX) | 0.59 | 0.60 |
| Query By Committee (qbc) | 0.59 | 0.60 |

Finally, to demonstrate the consistency of optimal performance of our proposed method and to highlight its fast convergence, we plot the mean train/test accuracy rates against the number of queries for our method and the baselines with error-bar generated by running 5 experiments on seeds 0-4. Figure 15 demonstrates the result.

Figure: (a) Test
Refer to caption: x20.png

#### G.5 Quadratic Regression Experiment

###### Data Generation and AL Implementation Details.

We generate the experiment dataset in this section simply according to the quadratic equation $y=x^{2}$ without adding any noise. Since the dimension of the regression dataset goes down from $d=3$ in the spiral dataset (with the third dimension acting as bias for the ReLU networks) to $d=2$ in the regression dataset, we increase the number of simulations to 2000 to randomly sample partitions for hyperplane arrangements to maintain an adequate embedding size. As in the spiral task, Table 4 summarizes the number of neurons sampled with respect to each seed. Once again, in our implementation of the deep AL baselines, we adjust the embedding
size accordingly when using a different seed to ensure fair comparisons.

**Table 4: Number of neurons for each seed (0-4) sampled using $2000$ simulations for the quadratic regression dataset.**
| Seed | 0 | 1 | 2 | 3 | 4 |
| --- | --- | --- | --- | --- | --- |
| Neurons | 160 | 160 | 157 | 159 | 160 |

As in the spiral experiment, we set the hyper-parameters for all deep active learning baselines according to Table 2. In the quadratic regression task, we similarly observe improved performance in deep AL baselines with a higher number of training epochs and a lower learning rate of $1\times 10^{-3}$, ensuring that the baselines provide a robust comparison. For baselines from scikit-activeml, we use the bagging regressor along with the default regressor to enhance their performance (see Appendix G.3). For our cutting-plane AL method (Algorithm 6) and the linear cutting-plane AL method (Algorithm 7), we choose the threshold value $\epsilon=1\times 10^{-3}$.

###### Experiment Results.

We now present the deferred experiment results in the quadratic regression experiment. To start, we discuss the performance of our cutting-plane AL method with and without convex solve. While using the full $80$ training data, the stand-alone convex solve achieves perfect prediction of quadratic regression (see Figure 7), Figure 16 shows that including the regularization in our cutting-plane AL method with a query budget of 20 leads to slightly suboptimal predictions compared to the non-regularized version. This could be because of the reduced complexity in the quadratic regression task, where the relationship between the features and the target is simple and well-behaved and the model is less likely to overfit. Therefore, the incorporation of regularization could slow down the parameter updates, leading to slower convergence.

Figure: (a) AFS
Refer to caption: x22.png

**Table 5: Train and test RMSE for the quadratic regression task ($y=x^{2}$) using the cutting-plane AL with a 2-layer ReLU neural network and various deep AL baselines with seed = 0.**
| Method | Train RMSE | Test RMSE |
| --- | --- | --- |
| Cutting-plane (ours) | 0.0111 | 0.0100 |
| Linear Cutting-plane (infeasible) | 0.7342 | 0.7184 |
| Random Sampling | 0.0824 | 0.0483 |
| Entropy Sampling | 0.0599 | 0.0529 |
| BALD with Dropout | 0.0568 | 0.0408 |
| Least Confidence Sampling | 0.0599 | 0.0529 |
| Greedy Sampling (GreedyX) | 0.0745 | 0.0447 |
| Query By Committee (qbc) | 0.1245 | 0.1366 |
| KL Divergence Maximization (kldiv) | 0.1356 | 0.0811 |
| Greedy Sampling Target (GreedyT) | 0.3296 | 0.3106 |

The trend plot in Figure 17 clearly illustrates the faster convergence of our cutting-plane AL method without regularization (BFS), compared to the regularized version (AFS), and significantly outpacing all other baselines. Therefore, for the quadratic regression task, we use the cutting-plane AL method without final solve (BFS). Nevertheless, it is noteworthy that both AFS and BFS cutting-plane AL method significantly outperforms all baselines and converging to the optimal prediction at a considerably faster rate.

Now we present the deferred results for the quadratic regression experiment in Section 7.1. Table 5 documents the train and test RMSE for our proposed cutting-plane AL against various deep AL baselines. In addition, Figure 18 presents the corresponding prediction made by our proposed method and the complete set of surveyed baselines.

Figure: (a) Test
Refer to caption: x24.png

Figure: Figure 18: Predictions of various AL algorithms for quadratic regression task using the cutting-plane AL with a two-layer ReLU neural network and various deep AL baselines (complete).
Refer to caption: extracted/5926809/cvx_plots/reg_final.png

#### G.6 IMDB Data Examples

**Table 6: Examples of IMDB movie reviews**
| Examples | Content | Label |
| --- | --- | --- |
| Review 1 | “If you like original gut wrenching laughter you will like this movie. | positive |
| If you are young or old then you will love this movie, hell even …” |  |  |
| Review 2 | “An American Werewolf in London had some funny parts, but this | negative |
| one isn’t so good. The computer werewolves are just awful…” |  |  |
| Review 3 | “”Ardh Satya” is one of the finest film ever made in Indian Cinema. | positive |
| Directed by the great director Govind Nihalani…” |  |  |

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 3
presentation: 3
questions: See weakness section.
soundness: 2
strengths: 1. The paper introduces a novel perspective of neural network training by leveraging cutting plane algorithms.
2. The algorithm also has significant theoretical guarantees in its convergence and effectiveness.
3. Experiments are conducted to showcase the effectiveness of the cutting plane approach
summary: The paper proposes a novel neural network training method by using the cutting plane algorithm. Extending from previous cutting plane methods to multi-layer ReLU network, the authors decompose the problem of data fitting to a set of linear programming problems. The authors also introduce a new active learning strategy that queries the most confident samples, hoping it finds misclassified samples to shrink the parameter set significantly.
weaknesses: I have several concerns, especially towards the scalability of the approach.
1. My biggest fear is the method will not work well with larger scale neural networks. This is because neural networks, with moderate capacity, are universal approximaters. For example, with a small ResNet-18, the model can fit arbitrary set of labels for the CIFAR-10 dataset with perfect accuracy. This suggests that in most vision/language applications involving large scale neural networks, the activation pattern $D$ of the last layer can simply be the set of $diag(\{0, 1\}^n)$. In such a case, cutting algorithms are not effective in training neural networks, nor doing active learning.
2. The active learning experiments all have extremely small budgets. This is potentially because the active learning strategy in this paper is somewhat counterintuitive. The algorithm queries the most confident samples, and hope for the model to be misclassifying them. However, for any neural network trained with a moderate amount of labeled samples and with gradient-based approach, the most confident samples are almost always correctly classified. I suspect the reason existing algorithms (with gradient based training) underperform is because, under limited budget, it is very easy for neural networks to make mistakes. However, in almost all deep learning applications, the model is only effective with high accuracies, i.e., the model needs to learn about uncertain samples. I also wonder if the gain from the cutting plane active learning algorithm only comes from the gradient-free training, instead of the active learning strategy.
3. The paper primarily relies on the property of ReLU activation, which is not generalizable to modern transformer architectures, where training the model is definitely not a linear program.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: $\bullet$ 
I think that the one of advantages of cutting-plane can be computational costs. Can you provide any discussion about this issue? In detail,  can you provide empirical runtime comparisons between the proposed method and gradient-based approaches on the datasets tested. This would help quantify any computational advantages.

$\bullet$ 
Is there any issue in the selection of querying (unlabeled) data points when multiple data points are targeted on a query?  How the method handles batch-mode active learning, where multiple samples are selected in each iteration, and whether there are any theoretical guarantees in that scenario.
soundness: 3
strengths: The proposed algorithm is interesting, and it has a solid theoretical background, especially for the separable cases. The algorithm is simple, and intuitive in the application aspects.
summary: This paper utilizes the cutting planes for active learning with two or multiple layers of MLP having ReLU activation functions. The cutting plane approach is free from the gradient and is usually used in the SGD, etc.  Cutting plane is the procedure to make the hyper-plane region smaller for the parameter estimation, along with the observed data points.  The cutting plane has been used in the linear binary classification. However, they show that the MLP with a ReLU in binary classification can be converted into linear programming with a constraint. From this point, they proposed a new algorithm, addressing the application of cutting plane to DNN, and active learning.

Theoretical validation based on linear programming is shown, and various experiments were done. On synthetic and IDBM datasets, the active learning of the proposed algorithm is dominant over Random, Entropy, Least confidence, BALD dropout, QBC, etc.
weaknesses: $\bullet$ The first issue is the applicability in real-world active learning with DNN. Many studies of active learning with DNN consider the CNN layer, ResNet architecture, transformer, etc in the image classification, natural langauge processing, and medical image analsys.  The linear layers with ReLU activation is limited in the use of active learning.

$\bullet$ The second issue is the range of experiments in the compared alg., and datasets. The datasets are too small compared to other active learning studies (usually considering 5~10 datasets such as tiny image-net or CIFAR100). In the aspects of algorithms, recent algorithms such as BADGE, BAIT, LDM-based, etc. are not considered. Especially, the BADGE is based on the gradient. The comparison of BADGE and cutting-plane can reveal the advantages or disadvantages of two algorithms in the overall performance and fast convergence in the early steps w/wo gradients.  It can highlight the proposed algorithm compared to gradient-based query algorithm. 

$\bullet$ Furthermore, the effect of the depth of layers can be an interesting issue. The advanced issues are not thoroughly studied.

---

confidence: 2
contribution: 3
initial_rating: 5
presentation: 4
questions: 1. How does the proposed method's computational complexity scale with the number of neurons and input dimensions? Could the authors provide a detailed analysis?

2. Have the authors considered extending this framework to other neural network architectures beyond ReLU networks, such as transformers or convolutional networks?

3. Could the authors elaborate on potential approaches to improve the scalability of their method, particularly in addressing the CPU-based solver limitation?

4. Is there a potential way to combine the benefits of both gradient-based and cutting-plane methods in a hybrid approach that could leverage the strengths of both paradigms?
soundness: 3
strengths: 1. The paper presents a significant theoretical contribution by providing the first convergence guarantees for deep active learning, addressing a crucial gap in the literature. The mathematical foundations are rigorous and well-explained.

2. The proposed method effectively eliminates common challenges associated with gradient-based methods, such as hyperparameter sensitivity and slow convergence, while maintaining competitive performance.

3. The method shows promising results in both classification and regression tasks when using small-scale datasets, indicating potential for specific use cases.
summary: This paper introduces a novel gradient-free approach for training deep neural networks using cutting-plane methods, along with an innovative active learning scheme. The authors successfully extend traditional cutting-plane methods, previously limited to linear models, to deep neural networks by reformulating ReLU network training as linear programming. The authors validate their method through experiments on both synthetic datasets and real-world sentiment classification tasks, showing superior performance compared to existing active learning baselines. The paper bridges the gap between classical optimization techniques and modern deep learning, offering a theoretically grounded alternative to gradient-based training methods.
weaknesses: 1. The method fundamentally lacks scalability due to its heavy reliance on CPU-based convex solvers and high computational complexity. The computational and memory requirements grow significantly with model size and data dimensionality, making it impractical for most real-world applications. This severely limits its practical utility in modern deep learning applications.

2. The method relies on a random sampling approach to collect ReLU activation patterns, which is not exhaustive and lacks theoretical guarantees for pattern coverage. This sampling process could miss important activation patterns, especially in high-dimensional data, potentially leading to suboptimal solutions. The paper does not provide sufficient analysis of how this limitation affects the model's performance and what patterns might be missed.

3. The method is extremely limited in its applicability, as it can only be used with neural networks using linear layers and ReLU activation functions. This restriction makes it incompatible with modern architectures such as transformers, convolutional neural networks, or networks using other activation functions, significantly limiting its practical impact in the field.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: see above
soundness: 2
strengths: 1.	Novelty in terms of the extension to Relu networks using the linear programming
2.	Provide justification to the linear programming approach.
3.	convergence analysis (mostly relying on the work of L. and R (2005))
summary: The paper describes an extension of the work of Louch and Ralaivola (2005) on using cutting plane methodology for active learning (AL) to Nonlinear cases using Relu networks. The work contributes some ingredients necessary to implement into Relu NNs using linear programming, provides justification for that, and also provide convergence guarantees. Some empirical validation is provided as for the methods capabilities with some standard and basic AL benchmarks.
weaknesses: There are several major weaknesses in the level of contribution of this paper, as well as in the empirical validation of the method. There are also various issues of clarity. Altogether  that makes this paper problematic for acceptance. 

1.	Effectiveness of the proposed methodology: Since the paper of L and R (2005) active learning and in particular deep active learning has made quite a progress. The active learning problem itself is meant to tackle the scalability of data labeling.

a.	 Much of Deep AL today work son massive data sets, in batch mode querying. These doesn’t not exist in this paper’s empirical validation. Instead the authors show their case on a simple 2D spiral binary classification and regressions problem and on some 50K data set of text data. By all means this is not a convincing demonstration of solving the scalability issues that AL tackles. More data sets are needed, and it would be also desired to test again R and L, on the data sets they used, even just as an ablation study (which is missing here too.)

b.	Multiclass extensions and batch mode querying are not touched here

c.	The volumetric stopping criterion is also not even touched upon, in particular its dependence on the dimension of the hypothesis space. How would a user designate it?

d.	Im missing some standard baselines for empirical comparison with active DL, such as Badge by J Ash et al, or Coreset by Sener et. al., other graph based Deep AL that uses graph built in the feature space and can address non-linearity  as such (e.g. Diffusion-based Deep Active Learning Kushnir et. al.)

2.	Contribution wise the convergence analysis is too similar to L and R (2005). It would be at least appropriate to address that Theorem 6.1 is similar to L and R.  One should also address the initial work of Tong and Kohler (JMLR, 2002) on cutting planes type methodology in AL. Also the effectiveness of the proposed method inlight of the the work done so far on deep AL is marginal in my view, and is not proved in this paper on enough examples and in terms of addressing key factors on scalability of AL

3.	Clarity:

a.	 the images, their label axes and their text fonts are reduced in size and sometime also not readable.

b.	Image 1 doesn’t correspond to the notation in the paper

c.	Statements such as line 93: “the size of parameters set reflects level of uncertainty” need to be made clear

d.	It is confusing to address theta as input in the context of DL

e.	Equation in line 163 definition is not clear: is it the space of all such diagonals matrix X? or the space of all vectors activated by the positive product?

f.	What is the benchmark accuracy for the IMDB? lines are quite flat. How does it compare with the R & L (2005) performance?

4.	The reader is redirected for many important details to the appendix, it disrupts the flow of the paper.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 2
questions: Please see weaknesses.
soundness: 2
strengths: The problem of using convex solvers for training neural networks is timely and relevant, and the theoretical results seem reasonable. The paper is detailed with several algorithms for various settings and numerical experiments.
summary: This paper introduces a scheme for simultaneous training and query selection of data for a (ReLU) deep neural network within an active learning context. The paper map the training of a ReLU DNN to a linear program and derive cutting plane algorithms to solve this program. This further leads to an active learning selection scheme. Finally, the paper validates the framework on several synthetic datasets.
weaknesses: This paper does not include much of the rich related literature around MIP formulations of neural networks (See [1]-[5] and all of their related references). For instance, the problem of training a ReLU neural network with 1 hidden layer to global optimality via polynomial complexity was addressed in [5]. Moreover, MIP/cutting planes/bisection search formulations of active learning for deep neural networks has also been well-studied (see [6]-[8]).

Furthermore, the numerical experiments focus primarily on synthetic experiments on the spiral data set & quadratic regression, giving some analysis on IMDB. Most importantly, these experiments all rely on small 2 or 3-layer ReLU, which limits analysis of the scaling of the algorithms to larger, more realistic neural networks.


[1] Anderson, Ross, et al. "Strong mixed-integer programming formulations for trained neural networks." Mathematical Programming 183.1 (2020): 3-39.

[2] Toro Icarte, Rodrigo, et al. "Training binarized neural networks using MIP and CP." Principles and Practice of Constraint Programming: 25th International Conference, CP 2019, Stamford, CT, USA, September 30–October 4, 2019, Proceedings 25. Springer International Publishing, 2019.

[3] Fischetti, Matteo, and Jason Jo. "Deep neural networks and mixed integer linear optimization." Constraints 23.3 (2018): 296-309.

[4] Bienstock, Daniel, Gonzalo Muñoz, and Sebastian Pokutta. "Principled deep neural network training through linear programming." Discrete Optimization 49 (2023): 100795.

[5] Arora, Raman, et al. "Understanding deep neural networks with rectified linear units." arXiv preprint arXiv:1611.01491 (2016).

[6] Sener, Ozan, and Silvio Savarese. "Active learning for convolutional neural networks: A core-set approach." arXiv preprint arXiv:1708.00489 (2017).

[7] Mahmood, Rafid, Sanja Fidler, and Marc T. Law. "Low-Budget Active Learning via Wasserstein Distance: An Integer Programming Approach." International Conference on Learning Representations.

[8] Ducoffe, Mélanie. Active learning and input space analysis for deep networks. Diss. COMUE Université Côte d'Azur (2015-2019), 2018.

Meta-review: a:

This paper proposes a gradient-free cutting-plane training method for ReLU networks, by utilzing the insight that MLP with a ReLU in binary classification can be converted into linear programming. The paper also proposes an active learning scheme. Convergence guarantee of the training method is shown as a Theorem. Experiments with synthetic and real-world datasets show the effectiveness of the proposed approach.

b:

A novel method is proposed which is a gradient-free cutting-plane training method for neural nets. Experiments show the effectiveness of the proposed method. Theoretical convergence analysis is shown. Intuition/motivation is clear.

c:

Reviewers discussed weaknesses in terms of effectiveness/practicality of the proposed methods (such as the problem/task and architecture), lack of significance of the convergence analysis, and scalability issues.

d:

Although the paper works on an interesting problem and it seems that substantial effort went into preparing the paper, a rebuttal to the reviews were not provided during the rebuttal period. This left concerns of the reviewers unaddressed.

Decision: Reject