\begin{document}

\doparttoc 
\faketableofcontents 

\title{A Greedy Algorithm For Building Compact Binary Activated Neural Networks}

\author{
  Benjamin Leblanc \\
  Université Laval\\
  Canada \\
  \texttt{benjamin.leblanc.2@ulaval.ca} \\
  \And
  Pascal Germain \\
  Université Laval\\
  Canada \\
  \texttt{pascal.germain@ift.ulaval.ca} \\
}

\maketitle

\begin{abstract}
We study binary activated neural networks in the context of regression tasks, provide guarantees on the expressiveness of these particular networks and propose a greedy algorithm for building such networks. Aiming for predictors having small resources needs, the greedy approach does not need to fix in advance an architecture for the network: this one is built one layer at a time, one neuron at a time, leading to predictors that aren't needlessly wide and deep for a given task. Similarly to boosting algorithms, our approach guarantees a training loss reduction every time a neuron is added to a layer. This greatly differs from most binary activated neural networks training schemes that rely on stochastic gradient descent (circumventing the 0-almost-everywhere derivative problem of the binary activation function by surrogates such as the straight through estimator or continuous binarization). We show that our method provides compact and sparse predictors while obtaining similar performances to state-of-the-art methods for training binary activated networks.
\end{abstract}

\section{Introduction}
\label{section:introduction}

Deep neural networks (DNNs) have achieved successes in various domains, namely in natural language processing \citep{DBLP:journals/spm/X12a} and computer vision \citep{DBLP:conf/nips/KrizhevskySH12, DBLP:journals/corr/HowardZCKWWAA17}.

These advances allows one with important computational and memory resources to train truly deep and wide neural networks on enormous corpus of text or image datasets containing millions of examples.
It is a widespread thought that the predictive performance of a NN increase with its number of parameters, motivating researchers to always train bigger and bigger NNs \citep{DBLP:journals/corr/SimonyanZ14a, DBLP:journals/corr/Krizhevsky14, DBLP:conf/cvpr/HeZRS16, DBLP:conf/cvpr/SzegedyLJSRAEVR15}. The main drawback of these models is that their high computational and memory costs prohibit their deployment in scenarios where the resources are limited. Various techniques have been proposed to lighten DNNs needs; in addition to being natural regularizer \citep{DBLP:conf/nips/HubaraCSEB16, DBLP:journals/corr/abs-1904-08444}, the use of weights encoded by a small number of bits \citep{DBLP:journals/pr/QinGLBSS20, DBLP:conf/nips/CourbariauxBD15, DBLP:conf/icml/MengBK20} and quantized activation functions, notably binary activation functions \citep{DBLP:conf/nips/SoudryHM14, DBLP:conf/nips/HubaraCSEB16, DBLP:conf/cvpr/WangLT0019}, have shown promising results. Also, methods such as network pruning \citep{DBLP:journals/corr/HanPTD15, DBLP:conf/iclr/RendaFC20, DBLP:journals/corr/HuPTT16, DBLP:journals/corr/IandolaMAHDK16, DBLP:conf/mobisys/GoetschalckxMWV18}, weight sharing \citep{DBLP:conf/icml/ChenWTWC15} or matrix factorisation \citep{DBLP:conf/bmvc/JaderbergVZ14} allow compacting a trained NN, resulting in predictors having less parameters.

The aforementioned processes involve many manipulations~:~first, it is necessary to train several deep and wide NNs (with various combinations of hyperparameters and architectures). Then, they must be pruned using, once again, various techniques, having no guarantee that trained network that obtain the best performances will still be the best pruned network. Plus, as implied by the Lottery Ticket Hypothesis from \citet{DBLP:conf/iclr/FrankleC19}, in order to obtain an efficient compact network from a pruned bigger network, the latter must be big enough so that the chances are high that, within a few different random seeds, it contains a compact but performing sub-network. Thus, it seems like a lot of work for obtaining, in the end, a small and compact predictor. Working backward would be an interesting avenue in order to obtain such networks without needing important computational resources to train wide DNNs first; starting with small NNs and making them grow, just like a decision tree is built. Unfortunately, while it's been shown that greedy approaches can scale to colossal problems such as ImageNet \citep{DBLP:conf/icml/BelilovskyEO19}, the literature concerning greedy approach for training NNs is pretty arid \citep{DBLP:journals/corr/KulkarniK17, DBLP:conf/nips/BengioLPL06, DBLP:journals/corr/abs-1905-10409}, and even more for networks with quantized components.

We propose Binary Greedy Network (BGN), a greedy algorithm for building one layer at a time, one neuron at a time (until a performance criteria is fulfilled) a fully connected binary activated neural networks for regression tasks. Therefore, the architecture of the resulting predictor is not fixed in advance:~it will vary depending on the complexity of the tackled problem. Theoretically well-grounded, this approach uses the fact that a layer's binary activation can be seen as a hyperplane separating the preceding layer's output space. The parameters of hidden layers neurons are obtained by computing a Lasso regression on the output of the preceding layer. The obtained predictors are thus compact (because of the greedy approach), sparse (because of the use of Lasso regression), rapidly trainable and efficient in a testing setting (because of the binary activation functions). The algorithm, inspired by Adaboost \citep{DBLP:conf/icml/FreundS96}, has convergence guarantees on the training mean squared error (MSE). At each layer, a lower bound on the obtainable MSE by the global network is minimized.

The contributions of this work are two folds. On the one hand, we make a comprehensive study of binary activated networks (\autoref{section:BANN}). 
On the other hand, we propose a greedy algorithm for building binary activated neural networks that is based on the study of these particular predictors (\autoref{section:algorithm}).

\section{Related work}
\label{section:related}
\paragraph{Binary activated neural networks (BANNs)} Neural networks learning mostly relies on stochastic gradient descent (SGD), which requires the objective function to be fully differentiable. Because the derivative of a binary function (a function whose output can only take two different values) is zero almost everywhere, this technique cannot directly be used for training BANNs. Most of the literature concerning BANNs propose workarounds allowing using the SGD algorithm. A simple way to estimate the gradient of the binary function is the straight-through estimator \cite{DBLP:journals/corr/BengioLC13, DBLP:conf/nips/HubaraCSEB16}, which uses the identity function as a surrogate of the gradient and which has then been refined over time \cite{DBLP:journals/corr/abs-1812-11800, DBLP:journals/corr/abs-1808-00278}. While such methods are convenient, they lack solid theoretical groundings. It has also been proposed to use continuous binarization \cite{DBLP:conf/iccv/GongLJLHLYY19, DBLP:journals/corr/abs-1911-09464, DBLP:conf/icassp/SakrCWGS18}; continuous activation functions are used, increasingly resembling a binary activation functions over the training iterations. Another branch of BANNs training algorithms consists of assuming a probability distribution on the weights of the network; doing so, one can work with the expectation of each layers output and train BANNs with the SGD algorithm \cite{DBLP:conf/nips/SoudryHM14, DBLP:conf/nips/LetarteGGL19}.

\paragraph{Greedy neural network training} Greedy approaches for training DNNs usually refers to the use of a layer-wise training scheme \cite{DBLP:conf/icml/BelilovskyEO19, DBLP:journals/corr/KulkarniK17, DBLP:conf/nips/BengioLPL06, DBLP:journals/corr/abs-2107-04466, DBLP:conf/evoW/CustodeTBCCV20}.  The usual procedure is the following: first, the full network is trained. Then, the first layer of the network is fixed. The network is retrained, where each layer but the first is fine-tuned. The second layer of the network is then fixed, and so on.

\section{Background and notation}
\label{section:background}
We focus our study on multivariate regression tasks. Each task is characterized by a dataset $S = \{(\mathbf{x}_i,\mathbf{y}_i)\}_{i=1}^m$ containing $m$ examples, each one described by features $\mathbf{x}\in\mathcal{X}$ and labels $\ybf\in\Rbb^d$. From this dataset, we want to train a binary activated neural network (BANN), {\it i.e.}, with activation functions on hidden layers having two possible output values. Without loss of generality, every binary activation function can be described by three parameters $t$, $h_1$, $h_2 \in \mathbb{R}$, where $h_1 < h_2$, as follow :
\begin{align*}
f_{t,h_1,h_2}(x) = \begin{aligned}
    	\begin{dcases}
    	h_1&\textup{if } x < t,\\
    	h_2&\textup{otherwise},
    	\end{dcases}
    	\end{aligned}
\end{align*}
where the value of $f_{t,h_1,h_2}(t)$ is arbitrarily chosen (most of the time either $h_1$ or $h_2$). We will be using the sign activation function, using the following parametrization: $t = 0$, $h_1 = -1$ and $h_2 = +1$.

We consider fully-connected BANNs composed of $l \in \mathbb{N}^*$ layers $L_k$ of size $d_k$, $\forall k \in \{1,\dots,l\}$. We denote $d_0$ the size of the network input and $d_l$ the size of its output (thus, $d_l=d$). We call $l$ the depth of the network, $d_k$ the width of layer $L_k$ and $\max_{k \in \{1, 2, \dots, l\}}{d_k}$ the width of the network. The sequence $\mathbf{d}\eqdots\langle d_k\rangle_{k=0}^{l}$ constitutes the neural network architecture. We denote a specific BANN predictor by $B:\Xcal\to\Rbb^{d_l}$. Every predictor $B$ is characterized by the values of its weights $\{\Wbf_k\}_{k=1}^{l}$ and biases $\{\bbf_k\}_{k=1}^{l}$, where $\Wbf_k \in \mathbb{R}^{d_{k-1} \cdot d_k}$ and $\bbf_k \in \mathbb{R}^{d_k}$. Each layer $L_k$ can be seen as a function 
\begin{equation}\label{eq:layer_function}
L_k(\xbf) = f_k\left(\Wbf_k\xbf+\bbf_k\right) = f_k\left(\begin{bmatrix} \wbf_{k,1}\cdot\xbf+b_{k,1} \\ \vdots \\ \wbf_{k,d_k}\cdot\xbf+b_{k,d_k} \end{bmatrix}\right),
\end{equation}
where $f_k(\cdot)$, for $k\in\{1,\dots,l-1\}$, is a binary function acting element-wise and $f_l(\cdot)$ is the identity function. The whole network $B$ corresponds to the composition of each of its layers : $B(\xbf) = (L_l \circ L_{l-1} \circ \dots \circ L_1)(\xbf)$. We denote $L_{i:j}(\xbf) = (L_j \circ L_{j-1} \circ \dots \circ L_i)(\xbf), i,j \in \{0,\dots,L\}, i<j$. While $\mathcal{X}$ denotes the input space (domain) of the network, $\mathcal{L}_{i:j}$ denotes the image of the function given by the composition of layers $L_i$ to $L_j$ (with $\mathcal{L}_{i}\eqdots\mathcal{L}_{i:i}$) and $\mathcal{B}\eqdots \Lcal_{1:l}$ the image of the whole network.

\section{Binary activated neural networks: a study}
\label{section:BANN}

\subsection{The binary activation function}

When considering a BANN predictor, is the parametrization of the activation function matters? While most previous work on BANNs use of the sign function \cite{DBLP:journals/corr/abs-1911-09464, DBLP:journals/corr/abs-1812-11800} or the threshold function (parametrized as follows : $t = 0$, $h_1 = 0$ and $h_2 = 1$) \cite{DBLP:conf/nips/LetarteGGL19, DBLP:conf/nips/HubaraCSEB16, DBLP:journals/corr/abs-1808-00278}, attempts have been made to directly consider $h_1$, $h_2$ and/or $t$ as hyperparameters to optimize \cite{DBLP:conf/icassp/SakrCWGS18}. 

We emphasis that the choices of $t$, $h_1$ and $h_2$ doesn't matter, as long as $h_1 \neq h_2$ and that no restrictions are made toward weights and biases values.
\renewcommand\theprop{1}
\begin{prop}\label{prop:bin_func_param}
Let $B$ be a BANN parametrized by $\{\Wbf_k,\bbf_k\}_{k=1}^l$ such that its activation functions are $f_{t,h_1,h_2}$. The same goes for $B^*$, having binary functions $f_{t^*,h_1^*,h_2^*}$. For all $t,h_1,h_2,t^*,h_1^*,h_2^*\in\mathbb{R}$, there exist weights and biases $\{\Wbf_k^*,\bbf_k^*\}_{k=1}^l$ parametrizing $B^*$ such that
\begin{enumerate}[label=(\roman*)]
    \item $B^*$ has the same architecture as $B$\,;
    \item $B^*(\xbf) = B(\xbf)\ \forall \xbf \in \mathbb{R}^{d_0}$\,;
    \item $
\Wbf_{k}^* = \begin{aligned}
                \begin{dcases}
                \Wbf_k&\textup{if } k = 1\ ,\\
                \Wbf_k\frac{h_1-h_2}{h_1^*-h_2^*}&\textup{if } k \in \{2,\dots,l\}\, ,\\
                \end{dcases}
                \end{aligned}\quad\ \ 
\bbf_{k}^* = \begin{aligned}
                \begin{dcases}
                \bbf_{k}+(t-t^*)&\textup{if } k = 1\ ,\\
                \bbf_{k}+\mathbf{c}_k&\textup{if } k \in \{2,\dots,l\}\, ,\\
                \end{dcases}
                \end{aligned}\\
$
\end{enumerate}
with $\mathbf{c}_k = (t-t^*)-(h_1-h_2)(h_1^*-h_2^*)\left(\sum_j\wbf_{k,j}\right)$.
\end{prop}

The proof of all presented propositions and theorems can be found in the appendices. Let $\mathcal{P}_{\mathbf{d}; t,h_1,h_2}$ be the space of neural network predictors obtainable while having architecture $\mathbf{d}$ and binary activation functions $f_{t,h_1,h_2}$. \autoref{prop:bin_func_param} informs us that if no restrictions are made toward weights and biases values, then the predictor space $\mathcal{P}_{\mathbf{d}; t,h_1,h_2}$ is independent of $t$, $h_1$ and $h_2$. However, this is not necessarily true if the weights or biases are constrained to specific values (-1 or +1 \cite{DBLP:conf/nips/HubaraCSEB16, DBLP:conf/icml/MengBK20}, ternary values \cite{DBLP:journals/corr/LiL16, DBLP:conf/iclr/ZhuHMD17, DBLP:journals/corr/abs-2011-00580, DBLP:conf/sips/HwangS14}, low-bit values \cite{DBLP:conf/aaai/LengDLZJ18, DBLP:conf/cvpr/ZhouYWC18, DBLP:journals/jmlr/LybrandS21, DBLP:journals/corr/abs-1806-08342}, etc.) In these cases, $t$, $h_1$ and $h_2$ might just be yet some other hyperparameters of the model.

\subsection{The particularity of the regularized BANN}

In machine learning, regularization is known to help prevent overfitting a problem. It is done by adding a \textit{regularization term} to the loss function, oftenly the sum of the squared weights (Ridge regularization) or the absolute weights (Lasso regularization). By constraining the weights to smaller values, the model might just lose in complexity and shouldn't be as prone to overfitting than before. Problem is, binary activation functions such as sign and threshold are not sensible to the magnitude of their input, due to the fact that they're parametrized with $t=0$ : $\sgn(x) = \sgn(c\,x)$, $\forall c \in \mathbb{R}^*_+$. Therefore, when adding a regularization term to the objective function, chances are that weights (and biases too) on layers preceding activation layers will shrink down to 0. 

As in \citet{DBLP:journals/corr/abs-1812-11800}, the use of a modified Manhattan regularization $R(\Wbf) = \sum_{i,j}\big|\alpha-|w_{i,j}|\big|$ or a modified Euclidean regularization $R(\Wbf) = \sum_{i,j} (\alpha-|w_{i,j}|)^2$, with $\alpha\in\mathbb{R}^*_+$, could be useful, even though it introduces a new hyperparameter. Those two methods push the weights toward values $\{-\alpha,+\alpha\}$; it could be use to enforce binary weights by setting $\alpha = 1$. However, the value of $\alpha$ must be chosen carefully: for example, with the modified Euclidean regularization and $\alpha = \textup{mean}(\Wbf)$, as in \citet{DBLP:journals/corr/abs-1812-11800}, we have that $R(\Wbf) > R(c\Wbf), c\in(0,1)$, while the predictions aren't affected by that operation on the parameters of a layer preceding an activation layer.

\subsection{A hyperplane arrangement}
\label{section:BANN_hyperplanes}

Each neuron of layers $L_1$ to $L_{l-1}$ of a BANN can be seen as a hyperplane acting in the input space of the network ($L_1$) or the output space of the preceding layer ($L_2$ to $L_{l-1}$). Indeed, given an input $\xbf$, these neurons apply a binary function on a linear function of $\xbf$ (see \autoref{eq:layer_function}). The result indicates whether $\xbf$ is on one side ($-1$) or the other ($+1$) of a hyperplane parametrized by the coefficients $\{\wbf,b\}$ characterizing the neuron. Therefore, a layer composing a BANN is nothing else than a hyperplane arrangement, and the output $f_k(\xbf) \in \{-1,+1\}^{d_k}$ of layer $L_k$ indicates in which region created by the hyperplane arrangement is the input $\xbf$. \autoref{fig:hyperplanes} illustrates that phenomenon.

\begin{figure}
\caption{Different layers of a BANN with architecture $\mathbf{d}
\label{fig:hyperplanes}
\end{figure}

From the observation that the output space of layers $L_1$ to $L_{l-1}$ is discrete and finite, the following relation can be established for any BANN:~$|\mathcal{B}| \leq |\mathcal{L}_{1:l-1}| \leq \dots \leq |\mathcal{L}_{1}|$. In a few words, this means that the cardinality of the image of the function $(L_k \circ L_{k-1} \circ \dots \circ L_1)(\cdot)$ cannot grow with $k, k<l$. Therefore, when using a BANN predictor, the width of the leading hidden layers shouldn't be neglected, since it has the most impact on the capacity of the network to assign different predictions to different inputs. In fact, it seems that a good rule of thumb for selecting the architecture of a fully-connected BANN is to have $d_{l-1} < d_{l-2} < \dots < d_{1}$.

Having a shallow first hidden layer means that the network will only be able to render a few unique predictions. As already mentioned, layer $L_{k+1}$ will define regions in the output space of layer $L_{k}$ and $|\mathcal{L}_{1:k+1}| \leq |\mathcal{L}_{1:k}|$, meaning that, from the point of view of the network, some of the regions defined by layer $L_{k}$ are grouped together by layer $L_{k+1}$. A binary activated neural network can thus be defined that way: $L_1$ cut the input space of a BANN into a set of regions; $L_2$ to $L_{l-1}$ group several of these regions together; $L_l$ assigns a prediction to each grouped region.

The following proposition uses the logic previously developed in order to establish bounds on the training performance (mean squared error) of the network. Those bounds are obtained by analyzing every regions independently created by a hidden layer of a BANN.

\renewcommand\theprop{2.1}
\begin{prop}
Let $B$ be a BANN of depth $l$. Let $Q(B,S)$ be the mean squared error of predictor $B$ on dataset $S$ and $||\cdot||$ be the Euclidean norm. Then :
\begin{align*}
Q(B,S)
\eqdots 
\sum_{\mathclap{(\xbf, \mathbf{y})\in S}} \frac{\|B(\xbf) - \mathbf{y}\|^2}m
\geq \sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\frac{\left|\ybf_{\pbf,j}^{(l-1)}\right|}{m}\textup{Var}\left(\ybf_{\pbf,j}^{(l-1)}\right) \geq \dots \geq  \sum_{j=1}^{d_1}\sum_{\pbf\in\mathcal{L}_1}\frac{\left|\ybf_{\pbf,j}^{(1)}\right|}{m}\textup{Var}\left(\ybf_{\pbf,j}^{(1)}\right),
\end{align*}
where $m = |S|$ and $\ybf_{\pbf,j}^{(k)} = \{y_{i,j}\ |\  (\mathbf{x}_i,\mathbf{y}_i)\in S,\ (L_k \circ \dots \circ L_1)(\xbf_i) = \pbf\} \forall k\in\{1,\dots,l-1\}$.
\label{prop:bin_bnd}
\end{prop}

Note that a similar bound can be obtained for classification tasks on the 0-1 loss (see appendices). The bound can be interpreted as the weighted sum of the variance of the examples falling in every region created by a layer; the less the points falling thee a same region (at any hidden layer in the network) have dissimilar labels, the less is the obtainable mean squared error.

\subsection{Expressiveness of BANNs over ReLU networks}

Given a regression task, we naturally expect that a BANN needs more parameters that ReLU networks to obtain similar accuracy, but how many more? In order to quantify the expressiveness power of such NNs, we consider the required number of parameters in order to estimate any function from $F_{n,d}$ with an error of at most $\epsilon$. $F_{n,d}$ corresponds to the space of functions from a Sobolev space $\mathcal{W}^{n,\infty}(\textup{[0,1]}^d)$, $n,d \in \mathbb{N}^*$, having a norm upper bounded by 1. This metric is borrowed from previous works studying  the expressiveness power of several NN architectures, from ReLU networks \cite{DBLP:journals/nn/Yarotsky17, DBLP:conf/colt/Yarotsky18, DBLP:journals/corr/abs-1806-01528} to networks with quantized weights \cite{DBLP:conf/iclr/DingLXS19}, and even networks having both ReLU and binary activation functions \cite{DBLP:conf/iclr/DingLXS19, DBLP:conf/iclr/LiangS17}, but never BANNs. For more details on Sobolev functional spaces, see appendices. 

Inspired by the work of \citet{DBLP:journals/nn/Yarotsky17}, we develop in \autoref{theo:lb} lower and upper bounds on that required number of parameters.

\begin{theorem}
For $n, d \in \mathbb{N}^*$, there exists a binary activated network capable of estimating any function from $F_{n,d}$ with maximum absolute error $\epsilon$ (i) for which the width of the first hidden layer (and of the global network) is in $O\big(c\epsilon^{-(d/n+1)}\big)$, where $c$ is a constant depending on $n$ and $d$, (ii) the depth is constant, depending on $n$ and $d$ and (iii) the number of weights is in $\Omega\big(c\epsilon^{-d/n}\textup{ln}^{-1}(\epsilon^{-1})\big)$.
\label{theo:lb}
\end{theorem}

Notice that in \autoref{theo:lb}, since the depth is constant, then the upper bound from \emph{(i)} also holds for the number of weights. Equations (\ref{eq:bnd1}-\ref{eq:bnd4}) below summarize lower and upper bounds for BANNs as well as those holding for ReLU networks, which are taken from \citet{DBLP:journals/nn/Yarotsky17}. For comparison purposes, we also report the ratios between selected bounds by Equations (\ref{eq:ratio1}-\ref{eq:ratio3}).

\begin{minipage}[c]{0.44\textwidth}
\begin{equation}\label{eq:bnd1}
    \hspace{-10mm}\Omega\ \textup{ReLU}: c\epsilon^{-(d/2n)}
\end{equation}

\begin{equation}\label{eq:bnd2}
    \Omega\ \textup{BANN}: c\epsilon^{-d/n}\textup{ln}^{-1}(\epsilon^{-1})
\end{equation}

\begin{equation}\label{eq:bnd3}
    \hspace{-2.5mm}\textup{O ReLU}: c\epsilon^{-d/n}\textup{ln}(\epsilon^{-1})
\end{equation}

\begin{equation}\label{eq:bnd4}
    \hspace{-8.5mm}\textup{O BANN}: c\epsilon^{-(d/n+1)}
\end{equation}
\end{minipage}
\begin{minipage}[c]{0.02\textwidth}
\ 
\end{minipage}
\begin{minipage}[c]{0.54\textwidth}
\begin{equation}\label{eq:ratio1}
    \quad\textup{O ReLU} / \Omega\ \textup{ReLU}: \epsilon^{-(d/2n)}\textup{ln}(\epsilon^{-1})
\end{equation}

\begin{equation}\label{eq:ratio2}
    \hspace{-5.5mm}\textup{O BANN} / \Omega\ \textup{BANN}: \epsilon^{-1}\textup{ln}(\epsilon^{-1})
\end{equation}

\begin{equation}\label{eq:ratio3}
    \hspace{-1.5mm}\textup{O BANN} / \textup{O ReLU}: \epsilon^{-1} /\ \textup{ln}(\epsilon^{-1})
\end{equation}
\end{minipage}

Notice that the ratio O/$\Omega$ only depend on $\epsilon$ for BANNs, while it also depends on $d$ and $n$ for ReLU networks. Therefore, with BANNs, the tightness of the bounds does not depends on the input dimension while it matters for ReLU networks. Since using a large $n$ is restrictive, because $|F_{n+1,d}|\leq|F_{n,d}| \forall n, d \in \mathbb{N}^*$, and that the caracteristics dimension $d$ could be big when it comes to function estimation problems, then the bounds might just be tighter for BANNs than ReLU networks for several function estimation problems. Notice that the ratio O Binary on O ReLU increases less than linearly with $\epsilon^{-1}$. As such, the number of weights needed for a BANN in order to obtain the same performance as a ReLU network on some regression tasks might not be as important as one could have imagined, even when the required maximum absolute error is truly small.

\section{The Binary Greedy Network (BGN) algorithm}
\label{section:algorithm}
From the study of binary activated neural networks of the previous section, and particularly the observations made in \autoref{section:BANN_hyperplanes}, we design an algorithm for greedily building such predictors. We leverage the bounds obtained from \autoref{prop:bin_bnd} as an objective function for parameterizing the hidden layers of the BANNs. We first explain the core algorithm principle in the simple cases of single-hidden layers neural networks and univariate regression tasks. Then, we extend the algorithm to deep and sparse BANNs and multivariate regression tasks.

\subsection{BGN algorithm for single-hidden layer neural networks on univariate regression tasks}

Consider a single-hidden layer BANN (that is, $l=2$) with a real-valued input space  ($\xbf\in\Rbb^{d_0}$) and an univariate output ($y\in\Rbb^{d_2}$, with $d_2 =1$).
In order to learn a compact BANN architecture, the proposed BGN algorithm iteratively construct the hidden layer by adding one neuron at the time.
Thus, we start by building a network with architecture $\mathbf{d}^{(1)} = \langle d_0,1,d_2\rangle$, then adding a neuron on the hidden layer in order to obtain an architecture $\mathbf{d}^{(2)} = \langle d_0,2,d_2\rangle$ and so on. Recall that in a BANN, every neuron of a hidden layer can be seen as a hyperplane. The idea is to iteratively minimize the bound from \autoref{prop:bin_bnd} as the layer is built, so that the final predictor has the lowest lower bound on the obtainable MSE. Using the notation defined in \autoref{prop:bin_bnd}, we aim to parameterize the first neuron of the hidden layer as follows:
\begin{equation}\label{eq:simple_hyperplane}
\underset{\mathbf{w},b}{\argmin}\ \left( \frac{|\ybf_{-1}|}{m}\textup{Var}\left(\ybf_{-1}\right)+\frac{|\ybf_{+1}|}{m}\textup{Var}\left(\ybf_{+1}\right)\right),
\end{equation}
with $m = |S|$ and $\ybf_{\pm1} = \{y_{j}\ |\ (\mathbf{x}_j,y_j)\in S,\ \sgn(\mathbf{w}\cdot\mathbf{x}_j+b) = \pm1\}$. The parameterization of the following neurons to be added to the hidden layer will aim to minimize that same quantity, but with modified labels. Therefore, we want to place a hyperplane in the input space such that the weighted sum of the variance of examples falling on each of its side is minimized. Since it is not a convex problem, we rely on a heuristic, which is based on the fact that the hyperplane location we are seeking aims to separate the points with "large" labels values from points with "small" labels values. Note that the heuristic use an auxiliary linear learning algorithm to find the weights of the added neurons, which is reminiscent of the Boosting algorithms that add a predictor to the learned ensemble at every iteration \citep{DBLP:conf/icml/FreundS96}.

\begin{heuristic}
Use a linear regression algorithm to obtain a weight vector $\mathbf{w}$ (do not consider the obtained bias value). Once $\mathbf{w}$ is fixed, find $b$ minimizing \autoref{eq:simple_hyperplane}.
\label{heuristic:regression}
\end{heuristic}

Note that, in \autoref{heuristic:regression}, finding $b$ minimizing \autoref{eq:simple_hyperplane} can be computed in O($m\ ln\ m$) operations. It is done by sorting the examples by their distance to the hyperplane parametrized by $\mathbf{w}$, without bias, and computing the obtained MSE on the $m+1$ possible separation of the sorted dataset. Once the hyperplane is fixed, that is, when parameters $(\mathbf{w}_{1,t},b_{1,t})$ of the network are fixed, the way parameter $w_{2,t}$ is obtained and $b_{2}$ is updated is presented in \autoref{prop:BGN_min}, where $c_t$ corresponds to $w_{2,t}$ and $d_t$ corresponds to the update on $b_{2}$ during iteration $t$. \autoref{algo:bgn} illustrates and summarizes the process of building the single-hidden layer BANN.
\renewcommand\theprop{3}
\begin{prop}
In \autoref{algo:bgn}: $(c_t, d_t) = \argmin_{c, d}\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)^2$.
\label{prop:BGN_min}
\end{prop}

\begin{algorithm}[t]
	\caption{Binary Greedy Network (BGN) - One hidden layer predictor}\label{algo:bgn}
	\begin{algorithmic}[1]
	\STATE \textbf{Input} : $S = \{(\xbf_1,y_1),\dots,(\xbf_m,y_m)\}$, $y \in \mathbb{R}$, $\xbf\in\mathbb{R}^d$, a dataset
	\STATE $\qquad\quad$ $T$, a number of iteration
	\STATE $\qquad\quad$ \textup{Heuristic(}$\cdot$\textup{)}, a method for placing hyperplanes
	\STATE Set $\mathbf{r}^{(1)} = \ybf = (y_1,\ldots, y_m)$ 
	\STATE \textbf{For} $t = 1,\dots,T$ :
	\STATE $\qquad$ $\wbf_t, b_t = \textup{Heuristic}(S^{(t)})$,
	\quad where $S^{(t)}= \{(\xbf_1,r_1^{(t)}),\dots,(\xbf_m,r_m^{(t)})\}$\label{line:heuristic}
	\STATE $\qquad$ $c_t, d_t = \frac{1}{2}\left(\rho_+^{(t)}-\rho_-^{(t)}\right), \frac{1}{2}\left(\rho_+^{(t)}+\rho_-^{(t)}\right)$, \quad with $\rho_\pm^{(t)} = \frac{\sum_{i:\sgn(\xbf_i\cdot\mathbf{w}_t+b_t)=\pm1} r_i^{(t)}}{\sum_{i:\sgn(\xbf_i\cdot\mathbf{w}_t+b_t)=\pm1} 1}$
	\STATE $\qquad$ $r^{(t+1)}_i=r^{(t)}_i - c_t\ \sgn(\mathbf{w}_t\cdot\mathbf{x}_i+b_t) - d_t$\ \ $\forall i \in \{1,\dots,m\}$ \label{line:update_r}
	\STATE \textbf{Output} : $BGN_T(\xbf)$ = $\sum_{t=1}^T c_t\ \sgn(\wbf_t\cdot\xbf+b_t)+d_t$
	\end{algorithmic} 
\end{algorithm}

The following proposition informs us that the objective we chose for our hyperplanes to minimize not only optimize the minimum obtainable MSE of the global network, but ensure that at every iteration of the algorithm, \textit{i.e}., every time a neuron is added to the network, the MSE obtained by the network on its training dataset diminishes.
\renewcommand\theprop{4}
\begin{prop}
In \autoref{algo:bgn}, we have 
\begin{center}
$Q(BGN_{t-1},S) - Q(BGN_{t},S) = (c_t+d_t)(c_t-d_t) > 0, \ \forall t \in \mathbb{N}\setminus\{0,1\}.$
\end{center}
\label{prop:BGN_dec}
\end{prop}

\subsection{Generalizing the BGN algorithm to sparse and deep predictors and multivariate tasks}

A key motivation of our work is to find the model with the right complexity for the task at hand. To enforce the \emph{compactness} and the \emph{sparsity} of the learned BNN, the following two improvements are proposed to the basic BGN algorithm (\autoref{algo:bgn}).\footnote{The full algorithms implementing the improvements presented in this section are given in the appendices.}

\begin{improvement}
In order to tackle the drawbacks of the greedy building of a layer, at iteration $t>1$, the oldest hyperplane that has been placed in the input space (i.e., a hidden neuron) is removed before applying another iteration of the algorithm. Then, if the training error diminishes, the replaced hyperplane (neuron) is kept. Otherwise, the former is placed back. 

\label{improv:3}
\end{improvement}

\begin{improvement}
Instead of using linear regression as a heuristic for placing the hyperplanes (neurons), Lasso regression \citep{tibshirani96} can be used; the sparse predictor obtained results in a neuron with some zero-valued weights.  

\label{improv:1}
\end{improvement}

The next improvement is used in order to tackle multivariate regression, that is when the output of the network is a real-valued vector.

\begin{improvement}
To tackle multivariate regression problems where the output is in $\Rbb^d$, we consider an example $(\xbf_i,\ybf_i)=(\xbf_i, (y_{i,1},\ldots,y_{i,d}))$ as $d$ examples with a single output value when placing the hyperplanes (neurons). That is, we use the Heuristic of \autoref{algo:bgn}--Line~\ref{line:heuristic} with 
\begin{align*}
S^{(t)} = \textstyle\bigcup_{i=1}^m\{(\xbf_i,r_{i,j}^{(t)})\}_{j=1}^d,\mbox{ with the residuals initialized to } r^{(1)}_{i,j}=y_{i,j},
\end{align*}
and updated coordinate-wise at Line~\ref{line:update_r}. See the appendices for the whole algorithmic procedure, in which $c_t$ and $d_t$ are also replaced by vectors of dimension $d$ to produce a multivariate output.
\label{improv:2}
\end{improvement}

Note that the result from \autoref{prop:BGN_dec} also holds in the multivariate case, but the MSE now diminishes of $(\mathbf{c}_t+\mathbf{d}_t)(\mathbf{c}_t-\mathbf{d}_t)$ at iteration $t$. Finally, the last improvement allows extending our BGN algorithm to more than one hidden layer.

\begin{improvement}
A deep BANN is constructed by building the first layer until the MSE on a validation set does not improve more than a chosen criterion; then, a second hidden layer is built using the same technique, but uses the output of the first hidden layer for features of its training dataset; and so on. 
\label{improv:4}
\end{improvement}

\section{Numerical experiments}
\label{section:experiments}

We conduct two different sets of experiments (all of them being performed on 2 NVIDIA GPUs GeForce RTX 2080 Ti). The first one is to test if BGN can compete with other techniques from the literature for training BANNs in terms of predictive performance on different regression tasks. The second one aims to verify whether BGN yields compact predictors and if pruned networks can obtain similar performances to the small networks built by BGN.

\subsection{Performance experiments}
\label{subsec:perf}

\paragraph{Baselines} The first benchmark we compare BGN to is a BANN trained by SGD, with the help of the straight-through estimator (BNN$^*$). We also compare to BNN+ \cite{DBLP:journals/corr/abs-1812-11800}, which was inspired by Binarized Neural Network \cite{DBLP:conf/nips/HubaraCSEB16} but make use of a refined version of the straight-through estimator. The next two benchmarks use different continuous binarization tricks: Quantization Networks (QN) \cite{DBLP:journals/corr/abs-1911-09464} and Bi-real net$^*$ \cite{DBLP:journals/corr/abs-1808-00278} (adapted to fully-connected architectures). Note that test performances for continuous binarization algorithms have been calculated with the use of binary activation functions. We finally present the results obtained by BinaryConnect \cite{DBLP:conf/nips/CourbariauxBD15}, as oftenly seen in the BANN training algorithm's literature.

\paragraph{Experimental setup} For each benchmarks, three learning rates were tested (0.1, 0.01, 0.001), three different width (100, 500 and 1000 neurons) and depth (1, 2 and 3 hidden layers). For the BNN+ algorithm, $L_1$ or $L_2$ regularization has been tested, with coefficients $1e-6$ and $1e-7$. Both QN and BNN+ have a specific hyperparameter (respectively $T_{\textup{Start}}$ and $\beta$), for which three different values were tested (respectively 5, 10, 20 and 1, 2, 5). The Adam \cite{DBLP:journals/corr/KingmaB14} optimization algorithm was used for every benchmark, along with Kaiming uniform initialization \cite{DBLP:conf/iccv/HeZRS15}. The maximum number of epoch was set to 200. We used the BGN algorithm with Improvements~\ref{improv:1}, \ref{improv:2}, \ref{improv:3} and \ref{improv:4}. Since the BGN algorithm was overfitting too easily many datasets, we used the following procedure for selecting the regularization parameter of the Lasso regression:~the parameter is set to 1e5; if the Lasso regression leads to a solution with all coefficients being 0s (except for the bias), the regularization parameter is divided by 1.5. Parameters $r$ and $q$ (patience) mentioned in Improvements \ref{improv:3} and \ref{improv:4} are set to 10 and 20 respectively.

\paragraph{Datasets} We chose a variety of datasets having different number of examples (315 to 36~736) and input and output cardinality (4 to 20 and 1 to 6). There are 4 different multivariate tasks and 7 univariate tasks. The architecture and hyperparameter set yielding the best average validation MSE over 3 random initialization and train / validation separation are selected, and we report in \autoref{tab:results_binary} the corresponding architecture and mean test MSE. See the appendices for specific values and more details on the experiment setting, as well as extended details on the considered datasets.

\renewcommand{1.5}{1.25}

\begin{table}[t]
    \centering

    \caption{Experiment results over 3 repetitions for selected models: MSE obtained on test dataset (Error$_{T}$) of the selected predictor, width ($l$) and depth ($d$). Lowest MSE per dataset are bolded.}
    \label{tab:results_binary}
    \setlength{\tabcolsep}{2.5pt}
    {\tiny
    \begin{center}
    \begin{tabular}{|l||ccc||ccc|ccc|ccc|ccc|ccc|}
    \hline
    \multirow{2}{*}{Dataset} & \multicolumn{3}{c||}{BGN} & \multicolumn{3}{c|}{BC} & \multicolumn{3}{c|}{BNN$^*$} & \multicolumn{3}{c|}{BNN+} & \multicolumn{3}{c|}{Bi-real net$^*$} & \multicolumn{3}{c|}{QN}\\
    \cline{2-19}
    & Error$_T$ & $l$ & $d$ & Error$_T$ & $l$ & $d$ & Error$_T$ & $l$ & $d$ & Error$_T$ & $l$ & $d$ & Error$_T$ & $l$ & $d$ & Error$_T$ & $l$ & $d$\\
    \hline
    \hline
    \rowcolor{Gray}
    bike hour \cite{bike} & 3573.67 & 259 & 2 & 21881.79 & 1000 & 3 & 1948.73 & 1000 & 3 & \textbf{1712.45} & 1000 & 2 & 1818.36 & 1000 & 3 & 2008.47 & 1000 & 2\\
    carbon \cite{carbon1} & \textbf{0.0048} & 40 & 1 & 5.45 & 100 & 2 & 0.82 & 500 & 1 & 0.95 & 100 & 1 & 1.28 & 1000 & 1 & 0.94 & 500 & 3\\
    \rowcolor{Gray}
    diabete \cite{efron2004least} & \textbf{3755.93} & 6 & 1 & 9134.43 & 500 & 1 & 8919.54 & 1000 & 1 & 8999.32 & 500 & 2 & 8816.83 & 100 & 1 & 8811.23 & 500 & 2\\
    housing \cite{HARRISON197881} & \textbf{0.32} & 166 & 1 & 1.65 & 100 & 2 & 2.19 & 500 & 1 & 2.09 & 1000 & 1 & 2.36 & 500 & 1 & 2.16 & 500 & 1\\
    \rowcolor{Gray}
    hung pox \cite{rozemberczki2021chickenpox} & 223.76 & 3 & 1 & 235.75 & 100 & 2 & 250.07 & 100 & 2 & 229.83 & 100 & 3 & 242.55 & 100 & 3 & \textbf{214.27} & 500 & 3\\
    ist stock usd \cite{Akbilgic2014ANH} & 2.27 & 5 & 1 & 22787.86 & 100 & 1 & 2.98 & 1000 & 3 & 1.95 & 1000 & 1 & 2.42 & 1000 & 2 & \textbf{1.94} & 1000 & 1\\
    \rowcolor{Gray}
    parking \cite{10.1007/978-3-319-59513-9_11} & \textbf{6954.87} & 1000 & 1 & 176286.81 & 1000 & 1 & 32705.94 & 500 & 3 & 11684.12 & 1000 & 3 & 16543.25 & 1000 & 2 & 107347.98 & 1000 & 2\\
    power plant \cite{TUFEKCI2014126, HeysemKayaLocalAG} & 17.18 & 187 & 2 & 13527.9 & 500 & 2 & 18.34 & 1000 & 1 & \textbf{16.66} & 1000 & 2 & 16.77 & 1000 & 2 & 17.14 & 1000 & 2\\
    \rowcolor{Gray}
    solar flare \cite{Dua:2019} & \textbf{0.66} & 13 & 1 & 1.11 & 1000 & 3 & 0.8 & 100 & 3 & 0.71 & 100 & 2 & 0.82 & 100 & 3 & 0.73 & 100 & 1\\
    stock portfolio \cite{portfolio} & \textbf{0.0660} & 32 & 1 & 13.58 & 100 & 3 & 3.78 & 100 & 1 & 119.32 & 1000 & 3 & 2.55 & 100 & 1 & 0.36 & 500 & 3\\
    \rowcolor{Gray}
    turbine \cite{turbine} & 27.91 & 361 & 2 & 3477.78 & 100 & 2 & 21.92 & 500 & 3 & \textbf{19.02} & 1000 & 2 & 21.8 & 500 & 3 & 25.54 & 1000 & 2\\
    \hline
    \end{tabular}
    \end{center}}
\end{table}

\subsection{Pruning experiments}\label{subsec:pruning}

\textbf{Procedure} \ \ We chose two datasets from \autoref{subsec:perf} to test the ability of the considered algorithms to yield compact predictors : Istanbul Stock USD and Power Plant (two datasets where BGN is outperformed by a benchmark). We then took the benchmark algorithm that resulted in the best performances over three different random seed in both of these tasks (respectively BNN+ and QN), with the corresponding set of hyperparameters, and alternately prune the trained network by a factor 0.5 and fine-tune it until only a few parameters were left active (non-zero). The fine-tuning is done by applying the same algorithm as for training, but the learning rate is divided by 10 and the maximum number of epoch is divided by 2. We then compared the reported mean test loss (on 5 different random seeds) over the pruning operation to the test loss BGN obtained while building its predictor. The number of active weights is a common metric for comparing pruning methods \cite{DBLP:conf/mlsys/BlalockOFG20}.

\textbf{Pruning methods} \ \ 
We chose five pruning methods that have been shown to be competitive with more elaborate and complex methods \cite{DBLP:journals/corr/abs-1903-01611, DBLP:journals/corr/abs-1902-09574, DBLP:conf/nips/HanPTD15, DBLP:journals/corr/HanMD15} : Global Magnitude Pruning (GMP), Layerwise Magnitude Pruning (LMP), Global Gradient Magnitude Pruning (GGMP), Layerwise Gradient Magnitude Pruning (LGMP). Respectively, at each iteration of the methods, the pruned weights are those with the lowest absolute~:~value anywhere in the network; value for each layer; value $\times$ associated gradient, evaluated on a batch of inputs, anywhere in the network; value $\times$ associated gradient, for each layer. We also included Random Pruning (each weight is pruned with fixed probability), which is a common straw man in pruning experiments. Since QN uses a continuous binarization trick, the gradient-based pruning methods aren't applicable here.

\begin{figure}
\caption{Pruning experiments on the Power Plant dataset. Dotted lines show where BGN has converged.}
\label{fig:pruning1}

\caption{Pruning experiments on the Istanbul Stock USD dataset. Dotted lines show where BGN has converged.}
\label{fig:pruning2}
\end{figure}

\subsection{Discussion}
\label{section:discussion}

The performance experiments reported in \autoref{tab:results_binary} clearly show that BGN is competitive on various regression tasks when compared to several algorithms from the literature for training BANNs. This includes both univariate and multivariate regression tasks, as well as small and big datasets, datasets with a few features per example, etc. On some datasets, it seems that all of the baseline methods fail to provide a good predictor, despite all of the hyperparameter combination that were tested (BGN was used without the need for testing a set of hyperparameters, while the benchmarks had between 27 and 324 hyperparameters combinations to test - per dataset). The most surprising result occurs on the diabete dataset, where BGN's MSE is more than two times smaller then the best runner-up.

\autoref{tab:results_binary} also highlights that the predictors given by BGN are consistently shallower and narrower than those selected by the baselines. On simple tasks, BGN achieves competitive prediction using very few parameters, whereas the baselines provide predictors containing up to a million parameters. Obtaining, from the baselines, a compact but efficient predictor is not as easy as training a huge one and then pruning it. \autoref{fig:pruning1} and \autoref{fig:pruning2} both are inline with that assertion: on both a big (\autoref{fig:pruning1}) and a small (\autoref{fig:pruning2}) dataset, pruning the best predictor obtained from a baseline until it contains as few active parameters as the predictor given by BGN truly impact the obtained test MSE. A few pruning methods were tested, but even then, it seems that the best pruning method depends on the problem or the fine-tuning algorithm. When the size of the predictor is not a considered factor, baselines and BGN seems to be truly comparable, but when the predictors are forced to be considerably small, BGN definitely shines.

\section{Conclusion}
\label{section:conclusion}
We presented lower and upper bounds on the expressiveness of binary activated neural networks (BANNs) in the context of regression tasks and provided a comprehensive study of these particular networks. Using the observations that were made, we designed a greedy algorithm for building compact and sparse BANNs that is theoretically well-grounded. The approach abandoned the classical stochastic gradient descent scheme to build the predictor one neuron at a time, one layer at a time, without the need of fixing in advance an architecture. A con of the greedy approach was circumvent by letting neurons being removed and replaced. The proposed method also has a training loss-decrease guarantee over each iteration.

We showed that the method is able to tackle both small and big regression problems, obtaining MSE comparable to state-of-the art methods for training BANNs. We provided empirical evidences that our method is able to provide truly sparse and compact predictors capable of achieving good performances on various regression tasks when compared to pruned networks having similar numbers of non-zero parameters. 

At the moment, BGN algorithm is limited to classical fully-connected networks and regression tasks. 
Future work also includes the possibility of building in a similar way ReLU network and convolutional neural networks. Finally, since the method yields sparse and compact predictor, it would be interesting to test the interpretability properties of our approach on real world-datasets. 

\textbf{Societal Impact.} \ 
When machine learning algorithms are deployed in the real world, ethical concerns must be considered. 
In particular, it is hazardous to use them in automated decision systems that might affect the welfare of individuals. The same apply to the proposed algorithm; however, our focus on predictors that aren't needlessly complex for a given task is in line with the objective of attenuating those concerns. 

\section*{Checklist}

\begin{enumerate}

\item For all authors...
\begin{enumerate}
  \item Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
    \answerYes{}
  \item Did you describe the limitations of your work?
    \answerYes{See section \autoref{section:conclusion}}
  \item Did you discuss any potential negative societal impacts of your work?
    \answerYes{See the end of \autoref{section:conclusion}.}
  \item Have you read the ethics review guidelines and ensured that your paper conforms to them?
    \answerYes{}
\end{enumerate}

\item If you are including theoretical results...
\begin{enumerate}
  \item Did you state the full set of assumptions of all theoretical results?
    \answerYes{Part of the appendices and \autoref{section:BANN}.}
        \item Did you include complete proofs of all theoretical results?
    \answerYes{Part of the appendices}
\end{enumerate}

\item If you ran experiments...
\begin{enumerate}
  \item Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)?
    \answerNo{The code will be released on GitHub with the final version of the paper.}
  \item Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?
    \answerYes{Mostly in \autoref{section:experiments}, with extended informations in appendices.}
        \item Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)?
    \answerNo{Will be part of the appendices}
        \item Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)?
    \answerYes{See the beginning of \autoref{section:experiments}.}
\end{enumerate}

\item If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
\begin{enumerate}
  \item If your work uses existing assets, did you cite the creators?
    \answerNA{}
  \item Did you mention the license of the assets?
    \answerNA{}
  \item Did you include any new assets either in the supplemental material or as a URL?
    \answerNA{}
  \item Did you discuss whether and how consent was obtained from people whose data you're using/curating?
    \answerNA{}
  \item Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content?
    \answerNA{}
\end{enumerate}

\item If you used crowdsourcing or conducted research with human subjects...
\begin{enumerate}
  \item Did you include the full text of instructions given to participants and screenshots, if applicable?
    \answerNA{}
  \item Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable?
    \answerNA{}
  \item Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation?
    \answerNA{}
\end{enumerate}

\end{enumerate}

\appendix

\part{Appendix} 
\parttoc

\section{Supplementary material}
\label{section:sup}

\subsection{Extended details on the performance experiments}

\textbf{Tested hyperparameters}
\begin{itemize}
    \item \# Hidden layers : 1,2,3
    \item Width : 100, 500, 1000
    \item Learning rates : 0.1, 0.01, 0.001
    \item Regularization type (BNN+) : $L_1$, $L_2$
    \item Regularization values (BNN+) : $10^{-6}$, $10^{-7}$
    \item $\beta$ values (BNN+) : 1, 2, 5
    \item $T_{\textup{Start}}$ values (QN) : 5, 10, 20 (while $T_{\textup{At epoch n}}$ = $T_{\textup{Start}} * n$)
\end{itemize}

\textbf{Fixed hyperparameters}
\begin{itemize}
    \item Initialization  : Kaiming uniform  \cite{DBLP:conf/iccv/HeZRS15}
    \item Batch size : 512 for big datasets (over 9000 examples) and 64 for small datasets; see \autoref{tab:algo_overview}
    \item Maximum number of epochs : 200
    \item Loss function : Mean squared error (MSE)
    \item Optimization algorithm : Adam \cite{DBLP:journals/corr/KingmaB14}
    \begin{itemize}
        \item $\epsilon$ = 0.001
        \item $\rho_1 = 0.9, \rho_2 = 0.999$
        \item $\delta = 1e-8$
        \item $\lambda$ = 0
    \end{itemize}
    \item Learning rate decay : plateau
    \item Patience : 5
    \item Early stop : 20
\end{itemize}

The test set was always composed of 25\% of the total dataset (if the train-test separation isn't inerant to the problem), while the validation set was 20\% of the remaining data.\\

\renewcommand{1.5}{1.35}

\begin{table}[H]
    \centering
    \caption{Experiment results over 3 repetitions for selected models: mean MSE obtained on test dataset for the selected predictor with standard deviation.}
    \label{tab:results_binary_ext}
    \setlength{\tabcolsep}{9pt}
    {\small
    \begin{center}
    \begin{tabular}{|l||c||c|c|c|c|c|}
    \hline
    Dataset & BGN & BC & BNN$^*$ & BNN+ & Bi-real net$^*$ & QN\\
    \hline
    \hline
    \rowcolor{Gray}
     & 3573.67  & 21881.79 & 1948.73 & 1712.45 & 1818.36 & 2008.47\\
    \rowcolor{Gray}
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \rowcolor{Gray}
    \multirow{-3}{*}{bike hour} & 208.64 & 6954.39 & 134.03 & 144.65 & 137.1 & 120.36\\
     & 0.0048 & 5.45 & 0.82 & 0.95 & 1.28 & 0.94\\
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \multirow{-3}{*}{carbon} & 0.0002 & 2.84 & 0.32 & 0.52 & 0.54 & 0.08\\
    \rowcolor{Gray}
    & 3755.93  & 9134.43 & 8919.54 & 8999.32 & 8816.83 & 8811.23\\
    \rowcolor{Gray}
    & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \rowcolor{Gray}
    \multirow{-3}{*}{diabete} & 178.15 & 432.55 & 941.61 & 918.98 & 1102.82 & 721.51\\
     & 0.32 & 1.65 & 2.19 & 2.09 & 2.36 & 2.16\\
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \multirow{-3}{*}{housing} & 0.02 & 0.01 & 0.19 & 0.1 & 0.19 & 0.13\\
    \rowcolor{Gray}
     & 223.76 & 235.75 & 250.07 & 229.83 & 242.55 & 214.27\\
    \rowcolor{Gray}
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \rowcolor{Gray}
    \multirow{-3}{*}{hung pox} & 14.82 & 28.42 & 44.63 & 34.45 & 106.58 & 9.87\\
     & 2.27 & 22787.86 & 2.98 & 1.95 & 2.42 & 1.94\\
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \multirow{-3}{*}{ist stock usd} & 0.06 & 26955.24 & 0.63 & 0.13 & 0.12 & 0.08\\
    \rowcolor{Gray}
     & 6954.87 & 176286.81 & 32705.94 & 11684.12 & 16543.25 & 107347.98\\
    \rowcolor{Gray}
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \rowcolor{Gray}
    \multirow{-3}{*}{parking} & 2108.30 & 863.17 & 5146.44 & 2849.72 & 2690.26 & 10455.96\\
     & 17.18 & 13527.9 & 18.34 & 16.66 & 16.77 & 17.14\\
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \multirow{-3}{*}{power plant} & 0.42 & 13271.0 & 1.0 & 1.3 & 1.34 & 0.67\\
    \rowcolor{Gray}
     & 0.66 & 1.11 & 0.8 & 0.71 & 0.82 & 0.73\\
    \rowcolor{Gray}
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \rowcolor{Gray}
    \multirow{-3}{*}{solar flare} & 0.01 & 0.38 & 0.13 & 0.07 & 0.12 & 0.11\\
     & 0.0660 & 13.58 & 3.78 & 119.32 & 2.55 & 0.36\\
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \multirow{-3}{*}{stock portfolio} & 0.0332 & 8.41 & 0.57 & 205.89 & 0.28 & 0.2\\
    \rowcolor{Gray}
     & 27.91 & 3477.78 & 21.92 & 19.02 & 21.8 & 25.54\\
    \rowcolor{Gray}
     & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$  & $\pm$ \\
    \rowcolor{Gray}
    \multirow{-3}{*}{turbine} & 2.09 & 2841.67 & 1.01 & 1.34 & 0.66 & 0.68\\
    \hline
    \end{tabular}
    \end{center}}
\end{table}

\begin{table}[H]
    \centering
    \caption{Hyperparameters of selected models over 3 repetitions.}
    \label{tab:results_hyperparams}
    \setlength{\tabcolsep}{10pt}
    {\scriptsize
    \begin{center}
    \begin{tabular}{llccccccc}
    \hline
    Dataset & Algorithm & Width & Depth & lr & reg. type & $\lambda$ & $T_{\textup{start}}$ & $\beta$\\
    \hline\hline
\multirow{5}{*}{bike hour} & BC & 1000 & 3 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 1000 & 3 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 2 & 0.1 & L2 & 1e-07 & --- & 2\\
 & Bi-real net$^*$ & 1000 & 3 & 0.1 & --- & ---  & --- & ---\\
 & QN & 1000 & 2 & 0.01 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{carbon} & BC & 100 & 2 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 500 & 1 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 100 & 1 & 0.1 & L1 & 1e-07 & --- & 5\\
 & Bi-real net$^*$ & 1000 & 1 & 0.1 & --- & ---  & --- & ---\\
 & QN & 500 & 3 & 0.01 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{diabete} & BC & 500 & 1 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 1000 & 1 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 500 & 2 & 0.1 & L1 & 1e-06 & --- & 5\\
 & Bi-real net$^*$ & 100 & 1 & 0.1 & --- & ---  & --- & ---\\
 & QN & 500 & 2 & 0.1 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{housing} & BC & 100 & 2 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 500 & 1 & 0.01 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 1 & 0.01 & L1 & 1e-06 & --- & 2\\
 & Bi-real net$^*$ & 500 & 1 & 0.01 & --- & ---  & --- & ---\\
 & QN & 500 & 1 & 0.01 & --- & ---  & 10 & ---\\
\hline
\multirow{5}{*}{hung pox} & BC & 100 & 2 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 100 & 2 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 100 & 3 & 0.1 & L1 & 1e-07 & --- & 5\\
 & Bi-real net$^*$ & 100 & 3 & 0.1 & --- & ---  & --- & ---\\
 & QN & 500 & 3 & 0.1 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{ist stock usd} & BC & 100 & 1 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 1000 & 3 & 0.001 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 1 & 0.1 & L2 & 1e-07 & --- & 1\\
 & Bi-real net$^*$ & 1000 & 2 & 0.001 & --- & ---  & --- & ---\\
 & QN & 1000 & 1 & 0.1 & --- & ---  & 10 & ---\\
\hline
\multirow{5}{*}{parking} & BC & 1000 & 1 & 0.001 & --- & ---  & --- & ---\\
 & BNN$^*$ & 500 & 3 & 0.01 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 3 & 0.01 & L2 & 1e-06 & --- & 2\\
 & Bi-real net$^*$ & 1000 & 2 & 0.01 & --- & ---  & --- & ---\\
 & QN & 1000 & 2 & 0.01 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{power plant} & BC & 500 & 2 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 1000 & 1 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 2 & 0.1 & L2 & 1e-06 & --- & 5\\
 & Bi-real net$^*$ & 1000 & 2 & 0.1 & --- & ---  & --- & ---\\
 & QN & 1000 & 2 & 0.01 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{solar flare} & BC & 1000 & 3 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 100 & 3 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 100 & 2 & 0.1 & L1 & 1e-07 & --- & 5\\
 & Bi-real net$^*$ & 100 & 3 & 0.01 & --- & ---  & --- & ---\\
 & QN & 100 & 1 & 0.1 & --- & ---  & 5 & ---\\
\hline
\multirow{5}{*}{stock portfolio} & BC & 100 & 3 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 100 & 1 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 3 & 0.01 & L2 & 1e-06 & --- & 1\\
 & Bi-real net$^*$ & 100 & 1 & 0.1 & --- & ---  & --- & ---\\
 & QN & 500 & 3 & 0.01 & --- & ---  & 10 & ---\\
\hline
\multirow{5}{*}{turbine} & BC & 100 & 2 & 0.1 & --- & ---  & --- & ---\\
 & BNN$^*$ & 500 & 3 & 0.1 & --- & ---  & --- & ---\\
 & BNN+ & 1000 & 2 & 0.1 & L1 & 1e-07 & --- & 2\\
 & Bi-real net$^*$ & 500 & 3 & 0.1 & --- & ---  & --- & ---\\
 & QN & 1000 & 2 & 0.01 & --- & ---  & 5 & ---\\
\hline
    \end{tabular}
    \end{center}}
\end{table}

\subsection{Extended details on the pruning experiments}

Figures (\ref{fig:std_1}-\ref{fig:std_2}) present the same results showed on figures (\ref{fig:pruning1}-\ref{fig:pruning2}), but with standard deviation. Overall, the predictors given by the BGN algorithm seem to be more stable in terms of performance then the pruned baselines.

\begin{figure}
\caption{Pruning experiments on the Power Plant dataset with standard deviation over 3 random seeds. Dotted lines show where BGN has converged.}
\label{fig:std_1}
\end{figure}

\begin{figure}
\caption{Pruning experiments on the Istanbul Stock USD dataset with standard deviation over 3 random seeds. Dotted lines show where BGN has converged.}
\label{fig:std_2}
\end{figure}

\subsection{Extended details on the benchmarks and datasets}

\renewcommand{1.5}{1.5}

\begin{table}[H]
    \centering
    \caption{Benchmarks overview}
    \label{tab:methods_overview}
    \setlength{\tabcolsep}{3pt}
    {\small
    \begin{tabular}{l||c|ccccc}
    \cline{1-7}
    Algorithm & BGN & BC & BNN$^*$ & BNN+ & Bi-real net$^*$ & QN\\
    \cline{1-7}
    Weights & $\mathbb{R}$ & $\{-1,+1\}$ & $\mathbb{R}$ & $\mathbb{R}$ & $\mathbb{R}$ & $\mathbb{R}$\\
    Activations output & $\{-1,+1\}$ & $\mathbb{R}$ & $\{-1,+1\}$ & $\{0,1\}$ & $\{-1,+1\}$ & $\{0,1\}$\\
    Uses batch norm$^*$ & False & True & True & True & True & True\\
    Uses regularization & False & False & False & True & False & False\\
    \cline{1-7}
    \end{tabular}
    }
\end{table}

$^*$For every algorithm making use of batch normalization, each activation layer is preceded by a batch norm layer. The intrinsic parameters of batch norm layers aren't considered in \autoref{subsec:pruning}.

\begin{table}[H]
    \centering
    \caption{Datasets overview (C = Categorical, F = Floats, I = Integers)}
    \label{tab:algo_overview}
    \setlength{\tabcolsep}{4pt}
    {\small
    \begin{tabular}{lllllllll}
    \toprule
    \# & Dataset & Full name & Taken from & Source & $|\mathbf{y}|$ & $|\mathbf{x}|$ & Type of $\mathbf{x}$ & Nb. ex.\\
    \midrule
    1 & bike hour & Bike sharing dataset & UCI Repo. & \cite{bike} & 1 & 16 & F/I & 17 389 \\
    2 & carbon & Carbon nanotubes & UCI Repo. & \cite{carbon1} & 3 & 5 & F & 10 721\\
    3 & diabete & Diabetes & SKLearn & \cite{efron2004least} & 1 & 10 & F/I & 442 \\
    4 & housing & California housing & SKLearn & \cite{HARRISON197881} & 1 & 8 & F & 20 640 \\
    5 & hung pox & Hungarian chickenpox cases & UCI Repo. & \cite{rozemberczki2021chickenpox} & 1 & 20 & F & 521 \\
    6 & ist stock usd & Istanbul stock exchange (USD) & UCI Repo. & \cite{Akbilgic2014ANH} & 1 & 8 & F & 536 \\
    7 & parking & Parking Birmingham & UCI Repo. & \cite{10.1007/978-3-319-59513-9_11} & 1 & 4 & F & 35 717 \\
    8 & power plant & Combined cycle power plant & UCI Repo. & \cite{TUFEKCI2014126, HeysemKayaLocalAG} & 1 & 4 & F & 9568 \\
    9 & solar flare & Solar flare & UCI Repo. & \cite{Dua:2019} & 3 & 10 & C & 1389 \\
    10 & stock portfolio & Stock portfolio performance & UCI Repo. & \cite{portfolio} & 6 & 6 & F & 315 \\
    11 & turbine & Gas Turbine CO and NOx & UCI Repo. & \cite{turbine} & 2 & 10 & F & 36 733 \\ \addlinespace[-1.5mm]
    & & emission data set & & & & & & \\
    \bottomrule
    \end{tabular}
    }
\end{table}

\subsection{Algorithms}

Algorithms (\ref{algo:bgn_improv3}-\ref{algo:bgn_improv4}) respectively show the base BGN algorithm enhanced with Improvements (\ref{improv:3}, \ref{improv:2}, \ref{improv:4}). They allow a clear understanding of how the base algorithm changes with each one of them. \autoref{algo:bgn_general} depicts the complete BGN algorithm (base algorithm + every proposed improvement); it is the BGN algorithm that has been used in the experiments of \autoref{section:experiments}.

\begin{algorithm}[H]
	\caption{\autoref{algo:bgn} + \autoref{improv:3} (possibility to replace the hyperplanes)} 
	\begin{algorithmic}[1]\label{algo:bgn_improv3}
	\STATE \textbf{Input} : $S = \{(\xbf_1,y_1),\dots,(\xbf_m,y_m)\}$, $y \in \mathbb{R}$, $\xbf\in\mathbb{R}^d$, a dataset;
	\STATE $\qquad\quad$ $T$, a number of iteration; \textup{Heuristic(}$\cdot$\textup{)}, a method for placing hyperplanes
	\STATE Set $\mathbf{r}^{(1)} = \ybf = (y_1,\ldots, y_m)$ 
	\STATE \textbf{For} $t = 1,\dots,T$ :
	\STATE $\qquad$ $\wbf_t, b_t = \textup{Heuristic}(S^{(t)})$,
	\quad where $S^{(t)}= \{(\xbf_1,r_1^{(t)}),\dots,(\xbf_m,r_m^{(t)})\}$
	\STATE $\qquad$ $c_t, d_t = \frac{1}{2}\left(\rho_+^{(t;t)}-\rho_-^{(t;t)}\right), \frac{1}{2}\left(\rho_+^{(t;t)}+\rho_-^{(t;t)}\right)$, \quad with $\rho_\pm^{(t;k)} = \frac{\sum_{i:\sgn(\xbf_i\cdot\mathbf{w}_k+b_k)=\pm1} r_i^{(t)}}{\sum_{i:\sgn(\xbf_i\cdot\mathbf{w}_k+b_k)=\pm1} 1}$
	\STATE $\qquad$ $r^{(t+1)}_i=r^{(t)}_i - c_t\ \sgn(\mathbf{w}_t\cdot\mathbf{x}_i+b_t) - d_t$\ \ $\forall i \in \{1,\dots,m\}$
	\STATE $\quad$\textbf{If} $t>1$:
	\STATE $\quad\quad$ $k=0$
	\STATE $\quad\quad$ \textbf{While} $Q(BGN_t,S)$ diminishes \textbf{And} $k < t$:
	\STATE $\quad\quad\quad$ $k=k+1$
	\STATE $\quad\quad\quad$ $r^{(t)}_i=r^{(t)}_i + c_k\ \sgn(\mathbf{w}_k\cdot\mathbf{x}_i+b_k) + d_k$\ \ $\forall i \in \{1,\dots,m\}$
	\STATE $\quad\quad\quad$ $\wbf_k, b_k = \textup{Heuristic}(S^{(t)})$
	\STATE $\quad\quad\quad$ $c_k, d_k = \frac{1}{2}\left(\rho_+^{(t;k)}-\rho_-^{(t;k)}\right), \frac{1}{2}\left(\rho_+^{(t;k)}+\rho_-^{(t;k)}\right)$
	\STATE $\quad\quad\quad$ $r^{(t)}_i=r^{(t)}_i - c_k\ \sgn(\mathbf{w}_k\cdot\mathbf{x}_i+b_k) - d_k$\ \ $\forall i \in \{1,\dots,m\}$
	\STATE \textbf{Output} : $BGN_T(\xbf)$ = $\sum_{t=1}^T c_t\ \sgn(\wbf_t\cdot\xbf+b_t)+d_t$
	\end{algorithmic} 
\end{algorithm}

\begin{algorithm}[H]
	\caption{\autoref{algo:bgn} + \autoref{improv:2} (multivariate tasks)} 
	\begin{algorithmic}[1]\label{algo:bgn_improv2}
	\STATE \textbf{Input} : $S = \{(\xbf_1,\ybf_1),\dots,(\xbf_m,\ybf_m)\}$, $\ybf \in \mathbb{R}^k$, $\xbf\in\mathbb{R}^d$, a dataset
	\STATE $\qquad\quad$ $T$, a number of iteration; \textup{Heuristic(}$\cdot$\textup{)}, a method for placing hyperplanes
	\STATE Set $r_{i,j}^{(1)} = (y_{1,j},\ldots, y_{m,j})\forall i \in \{1,\dots,m\},j \in \{1,\dots,k\}$
	\STATE \textbf{For} $t = 1,\dots,T$ :
	\STATE $\qquad$ $S_j^{(t)} = \{(\xbf_1,r^{(t)}_{1,j}),\dots,(\xbf_m,r^{(t)}_{m,j})\} \forall j \in \{1,\dots,k\}$
	\STATE $\qquad$ $S^{*(t)} = \bigcup_{j=1}^k S_j^{(t)}$
	\STATE $\qquad$ $\wbf_t, b_t = \textup{Heuristic}(S^{*(t)})$
	\STATE $\qquad$ $c_{t,j}, d_{t,j} = \frac{1}{2}\left(\rho_{j+}^{(t)}-\rho_{j-}^{(t)}\right), \frac{1}{2}\left(\rho_{j+}^{(t)}+\rho_{j-}^{(t)}\right) \forall j \in \{1,\dots,k\}$, 
	\STATE $\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$ with $\rho_{j\pm}^{(t)} = \frac{\sum_{i:\sgn(\xbf_i\cdot\mathbf{w}_t+b_t)=\pm1} r_{i,j}^{(t)}}{\sum_{i:\sgn(\xbf_i\cdot\mathbf{w}_t+b_t)=\pm1} 1}$
	\STATE $\qquad$ $r^{(t+1)}_{i,j}=r^{(t)}_{i,j} - c_{t,j}\ \sgn(\mathbf{w}_t\cdot\mathbf{x}_i+b_t) - d_{t,j}$\ \ $\forall i \in \{1,\dots,m\}, j\in\{1,\dots,k\}$
	\STATE \textbf{Output} : $BGN_T(\xbf)$ = $\sum_{t=1}^T \mathbf{c}_t\ \sgn(\wbf_t\cdot\xbf+b_t)+\mathbf{d}_t$
	\end{algorithmic} 
\end{algorithm}

\begin{algorithm}[H]
	\caption{\autoref{algo:bgn} + \autoref{improv:4} (deep neural networks building)} 
	\begin{algorithmic}[1]\label{algo:bgn_improv4}
	\STATE \textbf{Input} : $S = \{(\xbf_1,y_1),\dots,(\xbf_m,y_m)\}$, $y \in \mathbb{R}$, $\xbf\in\mathbb{R}^d$, a dataset
	\STATE $\qquad\quad$ $T$, a number of iteration; \textup{Heuristic(}$\cdot$\textup{)}, a method for placing hyperplanes
	\STATE Set $\mathbf{r}^{(1)} = \ybf = (y_1,\ldots, y_m)$ 
	\STATE Set $\{\mathbf{x}_i^{(1)}\}_{i=1}^m = \{\mathbf{x}_i\}_{i=1}^m$ 
	\STATE \textbf{For} $k = 1,\dots,l-1$ :
	\STATE $\quad$ \textbf{For} $t = 1,\dots,T$ :
	\STATE $\quad\quad$ $\wbf_{k,t}, b_{k,t} = \textup{Heuristic}(S^{(k,t)})$,
	\quad where $S^{(k,t)}= \{(\xbf_1^{(k)},r_1^{(t)}),\dots,(\xbf_m^{(k)},r_m^{(t)})\}$
	\STATE $\quad\quad$ $c_t, d_t = \frac{1}{2}\left(\rho_+^{(k,t)}-\rho_-^{(k,t)}\right), \frac{1}{2}\left(\rho_+^{(k,t)}+\rho_-^{(k,t)}\right)$; $\rho_\pm^{(k,t)} = \frac{\sum_{i:\sgn(\xbf_i^{(k)}\cdot\mathbf{w}_{k,t}+b_{k,t})=\pm1} r_i^{(t)}}{\sum_{i:\sgn(\xbf_i^{(k)}\cdot\mathbf{w}_{k,t}+b_{k,t})=\pm1} 1}$
	\STATE $\quad\quad$ $r^{(t+1)}_i=r^{(t)}_i - c_t\ \sgn(\mathbf{w}_{k,t}\cdot\mathbf{x}_i^{(k)}+b_{k,t}) - d_t$\ \ $\forall i \in \{1,\dots,m\}$
	\STATE $\quad$ $\mathbf{x}_i^{(k+1)} = (L_{k}\circ\dots\circ L_1)(\mathbf{x}_i) \forall i \in \{1,\dots,m\}$
	\STATE \textbf{Output} : $BGN_T(\xbf)$ = $\sum_{t=1}^T c_t\ (L_{l-1}\circ\dots\circ L_1)(\mathbf{x})+d_t$
	\end{algorithmic} 
\end{algorithm}

\begin{algorithm}[H]
	\caption{The complete BGN algorithm}
	\begin{algorithmic}[1]\label{algo:bgn_general}
	\STATE \textbf{Input} : $S = \{(\xbf_1,\ybf_1),\dots,(\xbf_m,\ybf_m)\}$, $\ybf \in \mathbb{R}^k$, $\xbf\in\mathbb{R}^d$, a dataset
	\STATE $\qquad\quad$ $T$, a number of iteration; \textup{Heuristic(}$\cdot$\textup{)}, a method for placing hyperplanes
	\STATE Set $r_{i,j}^{(1)} = (y_{1,j},\ldots, y_{m,j})\forall i \in \{1,\dots,m\},j \in \{1,\dots,k\}$
	\STATE Set $\{\mathbf{x}_i^{(1)}\}_{i=1}^m = \{\mathbf{x}_i\}_{i=1}^m$
	\STATE \textbf{For} $z = 1,\dots,l-1$ :
	\STATE $\quad$ \textbf{For} $t = 1,\dots,T$ :
	\STATE $\quad\quad$ $S_j^{(z,t)} = \{(\xbf_1^{(z)},r^{(t)}_{1,j}),\dots,(\xbf_m^{(z)},r^{(t)}_{m,j})\} \forall j \in \{1,\dots,k\}$
	\STATE $\quad\quad$ $S^{*(z,t)} = \bigcup_{j=1}^k S_j^{(z,t)}$
	\STATE $\quad\quad$ $\wbf_{z,t}, b_{z,t} = \textup{Heuristic}(S^{*(z,t)})$
	\STATE $\quad\quad$ $c_{t,j}, d_{t,j} = \frac{1}{2}\left(\rho_{j+}^{(z,t;t)}-\rho_{j-}^{(z,t;t)}\right), \frac{1}{2}\left(\rho_{j+}^{(z,t;t)}+\rho_{j-}^{(z,t;t)}\right) \forall j \in \{1,\dots,k\}$, 
	\STATE $\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$ with $\rho_{j\pm}^{(z,t;p)} = \frac{\sum_{i:\sgn(\xbf_i^{(z)}\cdot\mathbf{w}_{z,p}+b_{z,p})=\pm1} r_{i,j}^{(t)}}{\sum_{i:\sgn(\xbf_i^{(z)}\cdot\mathbf{w}_{z,p}+b_{z,p})=\pm1} 1}$
	\STATE $\quad\quad$ $r^{(t+1)}_{i,j}=r^{(t)}_{i,j} - c_{t,j}\ \sgn(\mathbf{w}_{z,t}\cdot\mathbf{x}_i^{(z)}+b_{z,t}) - d_{t,j}$\ \ $\forall i \in \{1,\dots,m\}, j\in\{1,\dots,k\}$
	\STATE $\quad\quad$ \textbf{If} $t>1$:
	\STATE $\quad\quad\quad$ $p=0$
	\STATE $\quad\quad\quad$ \textbf{While} $Q(BGN_t,S)$ diminishes \textbf{And} $p < t$:
	\STATE $\quad\quad\quad\quad$ $p=p+1$
	\STATE $\quad\quad\quad\quad$ $r^{(t)}_{i,j}=r^{(t)}_{i,j} + c_{p,j}\ \sgn(\mathbf{w}_{z,p}\cdot\mathbf{x}_i^{(z)}+b_{z,p}) + d_{p,j}$\ \ $\forall i \in \{1,\dots,m\},j\in\{1,\dots,k\}$
	\STATE $\quad\quad\quad\quad$ $\wbf_{z,p}, b_{z,p} = \textup{Heuristic}(S^{*(z,t)})$
	\STATE $\quad\quad\quad\quad$ $c_{p,j}, d_{p,j} = \frac{1}{2}\left(\rho_{j+}^{(z,t;p)}-\rho_{j-}^{(z,t;p)}\right), \frac{1}{2}\left(\rho_{j+}^{(z,t;p)}+\rho_{j-}^{(z,t;p)}\right)\forall j\in\{1,\dots,k\}$
	\STATE $\quad\quad\quad\quad$ $r^{(t)}_{i,j}=r^{(t)}_{i,j} - c_{p,j}\ \sgn(\mathbf{w}_{z,p}\cdot\mathbf{x}_i^{(z)}+b_{z,p}) - d_{p,j}$\ \ $\forall i \in \{1,\dots,m\},j\in\{1,\dots,k\}$
	\STATE $\quad$ $\mathbf{x}_i^{(z+1)} = (L_{z}\circ\dots\circ L_1)(\mathbf{x}_i) \forall i \in \{1,\dots,m\}$
	\STATE \textbf{Output} : $BGN_T(\xbf)$ = $\sum_{t=1}^T \mathbf{c}_t\ (L_{l-1}\circ\dots\circ L_1)(\mathbf{x})+\mathbf{d}_t$
	\end{algorithmic} 
\end{algorithm}

\subsection{Sobolev spaces} 

We denote a Sobolev space $\mathcal{W}^{n,p}(\textup{[0,1]}^d)$, where $n,d \in \mathbb{N}^*$, $p \in [1,\infty)$, defined as the space of functions on [0,1]$^d$ lying in $L^p$, along with their weak derivatives up to an order $n$:
$$
\mathcal{W}^{n,p}(\textup{[0,1]}^d) = \{f\in L^{p}(\textup{[0,1]})\ |\ \forall\mathbf{n} \textup{ such that } |\mathbf{n}|\leq n, D^{\mathbf{n}}f \in L^p(\textup{[0,1]})\}\ ,
$$
where $\mathbf{n} \in \mathbb{N}^{*d}$, $|\mathbf{n}|$ = $\sum_i^d n_i$, $D^{\mathbf{n}}f$ is a weak derivative of $f$ and $L^p$ denotes a Lebesgue space. We focus on $\mathcal{W}^{n,\infty}(\textup{[0,1]}^d)$, where the functions from that space have a norm that is defined as follows :
$$
\begin{aligned}
||f||_{\mathcal{W}^{n,\infty}} &= \underset{p\rightarrow+\infty}{\textup{lim}} \left(\sum_{|\mathbf{n}|\leq n}\left|\left|D^\mathbf{n}u\right|\right|^p_{L^p}\right)^{1/p}\\
&= \underset{\mathbf{n}:|\mathbf{n}|\leq n}{\textup{max}}\ \underset{\mathbf{x}\in\textup{[0,1]}^d}{\textup{sup ess}}\left|D^\mathbf{n}f(\mathbf{x})\right|\ .
\end{aligned}
$$
An intuitive way to define a Sobolev space $\mathcal{W}^{n,\infty}(\textup{[0,1]}^d)$ is to see it as the space of functions in $C^{n-1}(\textup{[0,1]}^d)$ such that all their derivatives of order $n-1$ are Lipschitz continuous. In our results, we consider the unit ball $F_{n,d}$ of functions from $\mathcal{W}^{n,\infty}(\textup{[0,1]}^d)$ a norm upper bounded by 1:
\begin{equation}\label{eq:unit_ball}
\begin{aligned}
F_{n,d} = \{f \in \mathcal{W}^{n,\infty}(\textup{[0,1]}^d) : ||f||_{\mathcal{W}^{n,\infty}(\textup{[0,1]}^d)} \leq 1\}.
\end{aligned}
\end{equation}

\subsection{On why it makes sens to use (regularized) linear rergession to parametrize the hidden neurons of BANNs for regression tasks}

One of the claim that was made in the paper is that the BGN algorithm has solid theorical groundings. Therefore, how can we explain the use of (reguarized) linear regression for parametrizing the neurons of the BANNs? In what way does it makes sens and why does it work?

Let's consider the first hidden layer building. Adding a new neuron consist in adding a hyperplane in the input space. A hyperplane is parametrized by a vector to which it is orthogonal, plus a bias. Therefore, using linear regression to parametrize the hyperplane consists in directly using the direction toward which the linear predictor is pointing to parametrize the "orthogonal vector".

We saw that the hyperplane must, in order to obtain good empirical results, separates the input space in a way such that the term of \autoref{eq:simple_hyperplane} is the smallest possible. We saw that this equation corresponds to the weighted sum of the variances of the points in both resulting regions; therefore, the hyperplane has to separates points with "high-valued labels" to those having "small-values labels". 

Using linear regression in order to parametrize the hyperplane makes sens because the resulting vector will indicates a direction in which, in the input space, the labels grow. Using this direction to orientate the hyperplane than surely is a good way to separates points with "high-valued labels" to those having "small-values labels".

\section{On the sparseness of BGN predictors}
\label{section:sparse}
The sparseness of the predictors the BGN yields has been discussed in the paper, but not empirically supported. Two experiments have been conducted in order to demonstrate it. On two different datasets, networks with hidden layers of width respectively 100-100-20 and 100-65 have been built. \autoref{fig:sparse_1} and \autoref{fig:sparse_2} clearly show that the BGN predictors are composed of many 0-valued parameters, meaning that only the relevant parameters are non-zeros. Also, those two figures highlight the fact that BGN predictors are consistent; their number of parameters are quite stable.

\begin{figure}
\caption{Mean and standard deviation over three repetitions of the number of non-zero parameters of the predictor yeilded by BGN over it's training on the Carbon nanotubes dataset versus the maximum number of parameter a NN with a similar architecture could have.}
\label{fig:sparse_1}
\end{figure}

\begin{figure}
\caption{Mean and standard deviation over three repetitions of the number of non-zero parameters of the predictor yeilded by BGN over it's training on the Hungarian chickenpox cases dataset versus the maximum number of parameter a NN with a similar architecture could have.}
\label{fig:sparse_2}
\end{figure}

\subsection{Computational complexity}

Suppose a feedforward neural network with $l$ layers of width $d_k$. There are $m$ examples, $d_0$ features, the labels are in $\mathbb{R}^{d_l}$ and the network is trained for $t$ epochs. The remove-replace iterations can be done $j$ times per neuron. Then:

Time complexity of the BGN algorithm (\autoref{algo:bgn}) : 
$\begin{aligned}
\begin{dcases}
d_k \cdot d_0^3&\textup{if }\ ,\\
d_k \cdot d_0^2 \cdot m&\textup{if }\ ,\\
d_k \cdot m \textup{ ln } m&\textup{if }\ ,\\
\end{dcases}
\end{aligned}$\\
Time complexity of the BGN algorithm (\autoref{algo:bgn_improv3}) : 
$\begin{aligned}
\begin{dcases}
j \cdot d_k \cdot d_0^3&\textup{if }\ ,\\
j \cdot d_k \cdot d_0^2 \cdot m&\textup{if }\ ,\\
j \cdot d_k \cdot m \textup{ ln } m&\textup{if }\ ,\\
\end{dcases}
\end{aligned}$\\
Time complexity of the BGN algorithm (\autoref{algo:bgn_improv2}) : 
$\begin{aligned}
\begin{dcases}
d_k \cdot d_0^3&\textup{if }\ ,\\
d_k \cdot d_0^2 \cdot m \cdot d_l&\textup{if }\ ,\\
d_k \cdot m \cdot d_l \textup{ ln } (m \cdot d_l)&\textup{if }\ ,\\
\end{dcases}
\end{aligned}$\\
Time complexity of the BGN algorithm (\autoref{algo:bgn_improv4}) : 
$\begin{aligned}
\begin{dcases}
d_k^4 \cdot l&\textup{if }\ ,\\
d_k^3 \cdot l \cdot m&\textup{if }\ ,\\
d_k \cdot l \cdot m \textup{ ln } m&\textup{if }\ ,\\
\end{dcases}
\end{aligned}$\\
Time complexity of the BGN algorithm (\autoref{algo:bgn_general}) : 
$\begin{aligned}
\begin{dcases}
j \cdot d_k^4 \cdot l&\textup{if }\ ,\\
j \cdot d_k^3 \cdot l \cdot m \cdot d_l&\textup{if }\ ,\\
j \cdot d_k \cdot l \cdot m \cdot d_l \textup{ ln } (m \cdot d_l)&\textup{if }\ ,\\
\end{dcases}
\end{aligned}$\\

Here is how the computational complexity of the algorithms has been found : 
\begin{itemize}
    \item In \autoref{algo:bgn}, the following independent steps are made $d$ times :
    \begin{itemize}
        \item Using a linear regression algorithm to find the values of $\mathbf{w}_t$ is made in O($k^3 + m\cdot k^2$) operations (if the least squares algorithm is used, then $k^3$ appears because of a matrix inversion and $m\cdot k^2$ because of matrices multiplications; if the least angle regression (LARS\cite{lars}) algorithm is used, with $k<m$, then the conclusion is the same, while allowing the use of Lasso regression without impact the required number of operations\cite{lars});
        \item As mentioned earlier, finding $b$ minimizing \autoref{eq:simple_hyperplane} can be computed in O$(m\ \ln\ m)$ operations, mainly because the $m$ examples have to be sorted;
        \item Finding the optimal $c_t$ and $d_t$ values is done in O($m\cdot k$) operations;
        \item Calculating $r^{(t+1)}$ is done in O($m\cdot k$) operations;
        \item Leading to having \autoref{algo:bgn} in O($d \cdot (k^3 + m\cdot k^2 + m \ln m)$).
    \end{itemize}
    \item In \autoref{algo:bgn_improv3}, the independent steps aren't made $d$ times, but at most $j \cdot d$ times.
    \item In \autoref{algo:bgn_improv2}, the training dataset doesn't have $m$ examples, but $m \cdot d_y$ examples.
    \item In \autoref{algo:bgn_improv4}, the first hidden layer is built with the same number of operations as in \autoref{algo:bgn}: O($d \cdot (k^3 + m\cdot k^2 + m \ln m)$). The second layer, as well as the third one, and so on, are built in O($d \cdot (d^3 + m\cdot d^2 + m \ln m)$). Therefore, \autoref{algo:bgn_improv4} is in O($l \cdot d \cdot (d^3 + m\cdot d^2 + m \ln m)$) = O($l \cdot d^4 + l \cdot m \cdot d^3 + l \cdot d \cdot m \ln m)$).\\
\end{itemize}

Finally, since we study a regression context and that the problems we're tackling only have a few features compared to many examples, we have that for the considered datasets, the complete BGN algorithm is in O$(j \cdot d_k^3 \cdot l \cdot m \cdot d_l)$. Note that the time complexity of the backpropagation algorithm is in O($t \cdot d_k^2 \cdot l \cdot m$). Since, most of the time, $d_l$ is neglectable and $j$ is similar to $d_k$, then we obtain O$(d_k^4 \cdot l \cdot m)$ versus O($t \cdot d_k^2 \cdot l \cdot m$) for the backpropagation algorithm. It seems that the backpropagation algorithm is theoretically faster than BGN. Let's look at empirical evidences.\\

\subsection{Experiments}

\section{Mathematical results}
\label{section:math}
\subsection{Proof of \autoref{prop:bin_func_param}}

 In the main paper, there is a small error in \autoref{prop:bin_func_param}'s statement; at point (iii), the expression of $\mathbf{b}^*$ is inexact. It does not change the conclusions that were drawn from \autoref{prop:bin_func_param}. The proposition presented below is corrected.

\paragraph{Proposition 1.}
\textit{Let $B$ be a BANN parametrized by $\{\Wbf_k,\bbf_k\}_{k=1}^l$ such that its activation functions are $f_{t,h_1,h_2}$ on layers $L_1$ to $L_{l-1}$ and the identity function on layer $L_l$. The same goes for $B^*$, having binary functions $f_{t^*,h_1^*,h_2^*}$. For all $t,h_1,h_2,t^*,h_1^*,h_2^*\in\mathbb{R}$, $h_1>h_2$ and $h_1^* > h_2^*$, there exist weights and biases $\{\Wbf_k^*,\bbf_k^*\}_{k=1}^l$ parametrizing $B^*$ such that}
\begin{enumerate}[label=(\roman*)]
    \item $B^*$ \textit{has the same architecture as} $B$\,;
    \item $B^*(\xbf) = B(\xbf)\ \forall \xbf \in \mathbb{R}^{d_0}$\,;
    \item $
\Wbf_{k}^* = \begin{aligned}
                \begin{dcases}
                \Wbf_k&\textup{if } k = 1\ ,\\
                \Wbf_k\frac{h_1-h_2}{h_1^*-h_2^*}&\textup{if } 1 < k \leq l\, ,\\
                \end{dcases}
                \end{aligned}\qquad
\bbf_{k}^* = \begin{aligned}
                \begin{dcases}
                \bbf_{k}+\delta_t&\textup{if } k = 1\ ,\\
                \bbf_{k}-\mathbf{c}_k+\delta_t&\textup{if } 1 < k < l\ ,\\
                \bbf_{k}-\mathbf{c}_k&\textup{if } k = l ,\\
                \end{dcases}
                \end{aligned}\\
$
\end{enumerate}
\textit{with} $\delta_t = t^*-t$ \textit{and} $\mathbf{c}_k = \frac{(h_1^*-h_1)(h_1-h_2)}{h_1^*-h_2^*}\left(\sum_j\wbf_{k,j}\right)$.

\begin{proof}
The proof relies on finding weights and biases $\{\Wbf_k^*,\bbf_k^*\}_{k=1}^l$ such that $B^*(\xbf) = B(\xbf)\ \forall \xbf \in \mathbb{R}^{d_0}$. To do so, we will ensure that $\{\Wbf_k^*,\bbf_k^*\}_{k=1}^l$ are such that the input and output of every binary activation function is the same in both networks. We will first tackle the single-hidden layer case before extending it to deep binary activated neural networks. 

In order to lighten the notation, we denote $f = f_{t,h_1,h_2}$ and $f^* = f_{t^*,h_1^*,h_2^*}$.

First of all, knowing that $h^*_1 - h^*_2 \neq 0$:
$$
h_1 = \frac{h_1-h_2}{h_1^*-h_2^*}h_1^*-\left(\frac{h_1-h_2}{h_1^*-h_2^*}h_1^*-h_1\right)\\
= ah_1^*-k_1\ ,
$$

$$
h_2 = \frac{h_1-h_2}{h_1^*-h_2^*}h_2^*-\left(\frac{h_1-h_2}{h_1^*-h_2^*}h_2^*-h_2\right)\\
= ah_2^*-k_2\ ,
$$

Let's demonstrate that $k_1 = k_2$:
$$
\begin{aligned}
h_1-h_2 &= h_1-h_2\ ,\\
 &= (h_1^*-h_2^*)\frac{h_1-h_2}{h_1^*-h_2^*}\ ,\\
 &= h_1^*\frac{h_1-h_2}{h_1^*-h_2^*}-h_2^*\frac{h_1-h_2}{h_1^*-h_2^*}\ ,\\
\Rightarrow \frac{h_1-h_2}{h_1^*-h_2^*}h_2^*-h_2 &= \frac{h_1-h_2}{h_1^*-h_2^*}h_1^*-h_1\\
\Rightarrow k_1 &= k_2\ .\\
\end{aligned}
$$

By setting $k=k_1=k_2$ and $\delta = t^* - t$, we have te following relations between the parametrization of the binary functions from $B$ of the one from $B^*$:
$$
t = t^* - \delta\ ,\qquad h_1 = ah_1^*-k\ ,\qquad h_2 = ah_2^*-k\ .
$$
where $\delta$, $a$ and $k$ are known constants. Let us consider the first layer of the networks. In order to obtain $B^*(\xbf) = B(\xbf)\ \forall \xbf \in \mathbb{R}^{d_0}$, we must have the following relation:
\begin{equation}\label{eq:prop1_lay1}
L_1^*(w_{1,i,j}^*x+b_{1,i}^*) = \begin{aligned}
    	\begin{dcases}
    	h_1^*&\textup{if } L_1(w_{1,i,j}x+b_{1,i}) = h_1,\\
    	h_2^*&\textup{if } L_1(w_{1,i,j}x+b_{1,i}) = h_2,\\
    	\end{dcases}
    	\end{aligned}\\
    	\ \forall i\in d_1, j\in d_0.
\end{equation}
Looking at \autoref{eq:prop1_bin}, we notice that in network $B^*$, if $\mathbf{b}^*_1 = \mathbf{b}_1 + \delta$ and $\mathbf{W}_1^* = \mathbf{W}_1$, then \autoref{eq:prop1_lay1} holds.
\begin{equation}\label{eq:prop1_bin}
\begin{aligned}
\ \forall i\in d_1, j\in d_0 : f^*(\underset{w_{1,i,j}^*}{\underbrace{w_{1,i,j}}}x+\underset{b_{1,i}^*}{\underbrace{b_{1,i}+\delta}}) &= \begin{aligned}
    	\begin{dcases}
    	h_1^*&\textup{if } w_{1,i,j}x+b_{1,i}+\delta < t^*,\\
    	h_2^*&\textup{otherwise},\\
    	\end{dcases}
    	\end{aligned}\\
 &= \begin{aligned}
    \begin{dcases}
    h_1^*&\textup{if }w_{1,i,j}x+b_{1,i} < t,\\
    h_2^*&\textup{otherwise},\\
    \end{dcases}
    \end{aligned}\\
&= \begin{aligned}
    \begin{dcases}
    h_1^*&\textup{if } L_1(w_{1,i,j}x+b_{1,i}) = h_1,\\
    h_2^*&\textup{if } L_1(w_{1,i,j}x+b_{1,i}) = h_2,\\
    \end{dcases}
    \end{aligned}
\end{aligned}
\end{equation}

We now have to set weights and biases of layer $L_2$ in network $B^*$, such that
$$
\mathbf{W}_{2}f(\mathbf{x})+\mathbf{b}_{2} = \mathbf{W}^*_{2}f^*(\mathbf{x})+\mathbf{b}^*_{2}\ .\\
\label{eq_bin}
$$
Doing so, we would have $B(\mathbf{x}) = B^*(\mathbf{x})$ for the single-hidden layer case. Let $\mathbf{h}$ and $\mathbf{h}^*$ be column vectors, the respective output of $f(\mathbf{x})$ and $f^*(\mathbf{x})$ for a given input $\mathbf{x}$.
$$
\begin{aligned}
\ \forall i\in d_2: \mathbf{w}_{2,i}f(\mathbf{x})+b_{2,i} &= \mathbf{w}_{2,i}\mathbf{h}+b_{2,i}\\
&= \mathbf{w}_{2,i}(a\mathbf{h}^*-k)+b_{2,i}\\
&= a\mathbf{w}_{2,i}\mathbf{h}^*-k\underset{\hat{w}_2}{\underbrace{\left(\sum_jw_{2,i,j}\right)}}+b_{2,i}\\
&= a\mathbf{w}_{2,i}f^*(\mathbf{x})-k\hat{w}_2+b_{2,i}\\
&= \mathbf{w}_{2,i}^*f^*(\mathbf{x})+b^*_{2,i}\\
\end{aligned}\\
$$
We now have that $B(\mathbf{x}) = B^*(\mathbf{x})$ holds for single-hidden layer networks if $\{\mathbf{W}_i^*,\mathbf{b}_i^*\}^{2}_{i=1}$ are set as follows:
\begin{quote}
$
\begin{aligned}
\mathbf{W}_{i}^* &= \begin{aligned}
                \begin{dcases}
                \mathbf{W}_i&\textup{if } i = 1\ ,\\
                a\mathbf{W}_i&\textup{if } i = 2\ ,\\
                \end{dcases}
                \end{aligned}\\
                &= \begin{aligned}
                \begin{dcases}
                \mathbf{W}_i&\textup{if } i = 1\ ,\\
                \frac{h_1-h_2}{h_1^*-h_2^*}\mathbf{W}_i&\textup{if } i = 2\ ,
                \end{dcases}
                \end{aligned}
\end{aligned}
$
\end{quote}
\begin{quote}
$
\begin{aligned}
\mathbf{b}_{i}^* &= \begin{aligned}
                \begin{dcases}
                \mathbf{b}_{i}+\delta&\textup{si } i = 1\ ,\\
                \mathbf{b}_{i}-k\hat{w}_2&\textup{if } i = 2\ ,\\
                \end{dcases}
                \end{aligned}\\
                &= \begin{aligned}
                \begin{dcases}
                \mathbf{b}_{i}+(t^*-t)&\textup{if } i = 1\ ,\\
                \mathbf{b}_{i}-\frac{(h_1^*-h_1)(h_1-h_2)}{h_1^*-h_2^*}\left(\sum_j\mathbf{w}_{i,j}\right)&\textup{if } i = 2\ .
                \end{dcases}
                \end{aligned}
\end{aligned}
$
\end{quote}

The generalisation to deep BANNs is made by considering that in order to obtain $B^*(\xbf) = B(\xbf)\ \forall \xbf \in \mathbb{R}^{d_0}$, we must have $L_{1:k}(\mathbf{x})=L^*_{1:k}(\mathbf{x}) \forall k \in \{1,\dots,l\},\mathbf{x}\in\mathbf{R}^{d_0}$. Until now, we have considered how to modify weights and biases $\{\mathbf{W}_i^*,\mathbf{b}_i^*\}^{2}_{i=1}$ so that $\mathbf{W}_2L_1(\mathbf{x})+\mathbf{b}=\mathbf{W}_2^*L_1(\mathbf{x})+\mathbf{b}^*$ given a same input $\mathbf{x}$. By iteratively modifying the weights and biases of the network the same way we did for a single-hidden layer BANN, we obtain
$$
\begin{aligned}
\mathbf{x} &= \mathbf{x}\\
\Rightarrow \mathbf{W}_2\underset{L_1(\mathbf{x})}{\underbrace{f(\mathbf{W}_1\mathbf{x}+\mathbf{b}_1)}}+\mathbf{b}_2 &= a\mathbf{W}_2\underset{L_1^*(\mathbf{x})}{\underbrace{f^*(\mathbf{W}_1\mathbf{x}+\mathbf{b}_1+\delta)}}+\mathbf{b}_2-k\hat{w}_2\\
\Rightarrow \mathbf{W}_3f(\mathbf{W}_2L_1(\mathbf{x})+\mathbf{b}_2)+\mathbf{b}_3 &= a\mathbf{W}_3f^*(\underset{\mathbf{W}^*_2}{\underbrace{a\mathbf{W}_2}}L_1^*(\mathbf{x})+\underset{\mathbf{b}^*_2}{\underbrace{\mathbf{b}_2-k\hat{w}_2+\delta}})+\mathbf{b}_3-k\hat{w}_3\\
\end{aligned}
$$
By that last recursion, we obtain the following weights and biases $\{\Wbf_k^*,\bbf_k^*\}_{k=1}^l$ for any network of depth $l,l\in\mathbb{N}\setminus\{0,1\}$:
\begin{quote}
$
\begin{aligned}
\mathbf{W}_{i}^* &= \begin{aligned}
                \begin{dcases}
                \mathbf{W}_i&\textup{if } i = 1\ ,\\
                a\mathbf{W}_i&\textup{if } i \in \{2,\dots,l\}\ ,\\
                \end{dcases}
                \end{aligned}\\
                &= \begin{aligned}
                \begin{dcases}
                \mathbf{W}_i&\textup{if } i = 1\ ,\\
                \frac{h_1-h_2}{h_1^*-h_2^*}\mathbf{W}_i&\textup{if } i \in \{2,\dots,l\}\ ,\\
                \end{dcases}
                \end{aligned}
\end{aligned}
$
\end{quote}
\bigbreak
\begin{quote}
$
\begin{aligned}
\mathbf{b}_{i}^* &= \begin{aligned}
                \begin{dcases}
                \mathbf{b}_{i}+\delta&\textup{si } i = 1\ ,\\
                \mathbf{b}_{i}-k\left(\sum_j\mathbf{w}_{i,j}\right)+\delta&\textup{if } i \in \{2,\cdots,l-1\}\ ,\\
                \mathbf{b}_{i}-k\left(\sum_j\mathbf{w}_{i,j}\right)&\textup{if } i = l\ ,\\
                \end{dcases}
                \end{aligned}\\
                &= \begin{aligned}
                \begin{dcases}
                \mathbf{b}_{i}+(t^*-t)&\textup{if } i = 1\ ,\\
                \mathbf{b}_{i}-\frac{(h_1^*-h_1)(h_1-h_2)}{h_1^*-h_2^*}\left(\sum_j\mathbf{w}_{i,j}\right)+(t^*-t)&\textup{if } i \in \{2,\dots,l-1\}\ ,\\
                \mathbf{b}_{i}-\frac{(h_1^*-h_1)(h_1-h_2)}{h_1^*-h_2^*}\left(\sum_j\mathbf{w}_{i,j}\right)&\textup{if } i = l\ .\\
                \end{dcases}
                \end{aligned}
\end{aligned}
$
\end{quote}
\end{proof}

\subsection{Proof of Propositions \ref{prop:bin_bnd} and \ref{prop:bin_bnd_classif}}

\paragraph{Proposition 2.1.}
\textit{Let $B$ be a BANN of depth $l$. Let $Q(B,S)$ be the mean squared error of predictor $B$ on dataset $S$ and $||\cdot||$ be the Euclidean norm. Then}
\begin{align*}
Q(B,S)
\eqdots 
\sum_{\mathclap{(\xbf, \mathbf{y})\in S}} \frac{\|B(\xbf) - \mathbf{y}\|^2}m
\geq \sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\frac{\left|\ybf_{\pbf,j}^{(l-1)}\right|}{m}\textup{Var}\left(\ybf_{\pbf,j}^{(l-1)}\right) \geq \dots \geq  \sum_{j=1}^{d_1}\sum_{\pbf\in\mathcal{L}_1}\frac{\left|\ybf_{\pbf,j}^{(1)}\right|}{m}\textup{Var}\left(\ybf_{\pbf,j}^{(1)}\right),
\end{align*}
\textit{where} $m = |S|$ \textit{and} $\ybf_{\pbf,j}^{(k)} = \{y_{i,j}\ |\  (\mathbf{x}_i,\mathbf{y}_i)\in S,\ (L_k \circ \dots \circ L_1)(\xbf_i) = \pbf\} \ \forall k\in\{1,\dots,l-1\}$.

\begin{proof}
\begin{align} \allowdisplaybreaks
Q(B,S) &= \sum_{\mathclap{(\xbf, \mathbf{y})\in S}} \frac{\|B(\xbf) - \mathbf{y}\|^2}m\\
&= \frac{1}{m}\sum_{j=1}^{d_l}\sum_{(\xbf, \mathbf{y})\in S} \left(B(\xbf) - y_j\right)^2\\
&= \frac{1}{m}\sum_{j=1}^{d_l}\sum_{(\xbf, \mathbf{y})\in S} (\mathbf{w}_{l,j}\cdot L_{1:l-1}(\mathbf{x}) + b_{l,j} - y_j)^2\\
&= \frac{1}{m}\sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:l-1}(\mathbf{x})=\mathbf{p}}}(\mathbf{w}_{l,j}\cdot\pbf + b_{l,j} - y_j)^2\\
&\geq \frac{1}{m}\sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\argmin_{\mathbf{w},b}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:l-1}(\mathbf{x})=\mathbf{p}}}(\mathbf{w}\cdot\pbf + b - y_j)^2\\
&= \sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\frac{\left|\ybf_{\pbf,j}^{(l-1)}\right|}{m}\ \frac{1}{\left|\ybf_{\pbf,j}^{(l-1)}\right|}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:l-1}(\mathbf{x})=\mathbf{p}}}\left(\overline{\ybf_{\pbf,j}^{(l-1)}} - y_j\right)^2\\
&= \sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\frac{\left|\ybf_{\pbf,j}^{(l-1)}\right|}{m}\textup{Var}\left(\ybf_{\pbf,j}^{(l-1)}\right)\ .
\end{align}

The inequality is obtained by considering different values for $\mathbf{w}$ and $b$ in each of the $|\mathcal{L}_{1:l-1}|$ regions, while they're fixed in a typical neural network. The predictor $\overline{\ybf_{\pbf,j}^{(l-1)}}$ is the constant predictor minimizing the MSE on a given task: the mean of the labels to predict. We saw in \autoref{section:BANN_hyperplanes} that hidden layers $L_2$ to $L_l$ regroup regions created by the preceding hidden layer. Using the fact that for two datasets $S_1$ and $S_2$, we have
$$
\min_{c} \left(\sum_{(\xbf, y)\in S_1} (c-y_i)^2 + \sum_{(\xbf, y)\in S_2} (c-y_i)^2\right) \geq \min_{c} \sum_{(\xbf, y)\in S_1} (c-y_i)^2 + \min_{c} \sum_{(\xbf, y)\in S_2} (c-y_i)^2\ ,
$$
we obtain that $\forall k \in \{2,\dots,l-1\}$: 
$$
\sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:k}}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:k}(\mathbf{x})=\mathbf{p}}}\left(\overline{\ybf_{\pbf,j}^{(k)}} - y_j\right)^2\geq\sum_{j=1}^{d_l}\sum_{\pbf\in\mathcal{L}_{1:k-1}}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:k-1}(\mathbf{x})=\mathbf{p}}}\left(\overline{\ybf_{\pbf,j}^{(k-1)}} - y_j\right)^2\ ,
$$
completing the proof.
\end{proof}

\autoref{prop:bin_bnd_classif} aims to bound, as it's been done in \autoref{prop:bin_bnd}, the empirical loss obtained by a binary activated neural network, but this time on binary classification tasks and with the mean 0-1 loss. The proof is very similar to the one of \autoref{prop:bin_bnd}:

\renewcommand\theprop{2.2}
\begin{prop}\label{prop:bin_bnd_classif}
Let $B$ be a BANN of depth $l$. Let $Z(B,S)$ be the mean 0-1 error of predictor $B$ on dataset $S = \{(\mathbf{x}_i, y_i)\}_{i=1}^m$, $y_i \in \{-1,1\} \forall i \in \{1,\dots,m\}$. Let $I(\cdot)$ be the indicator function, where $I(a) = 1$ if $a$ is true, 0 otherwise. Then
\begin{align*}
Z(B,S)
\eqdots 
\frac{1}{m} \sum_{\mathclap{(\xbf, y)\in S}} I(B(\mathbf{x}) \neq y)
\geq \frac{1-\sum_{\pbf\in\mathcal{L}_{1:l-1}}\frac{\left|\ybf_{\pbf}^{(l-1)}\right|}{m}\left|\overline{\ybf_{\pbf}^{(l-1)}}\right|}{2} \geq \dots \geq  \frac{1-\sum_{\pbf\in\mathcal{L}_{1}}\frac{\left|\ybf_{\pbf}^{(1)}\right|}{m}\left|\overline{\ybf_{\pbf}^{(1)}}\right|}{2},
\end{align*}
where \textit{and} $\ybf_{\pbf,j}^{(k)} = \{y_{i,j}\ |\  (\mathbf{x}_i,\mathbf{y}_i)\in S,\ (L_k \circ \dots \circ L_1)(\xbf_i) = \pbf\} \ \forall k\in\{1,\dots,l-1\}$.
\end{prop}

\begin{proof}
\begin{align} \allowdisplaybreaks
Z(B,S) &= \frac{1}{m} \sum_{\mathclap{(\xbf, y)\in S}} I(B(\mathbf{x}) \neq y)\\
&= \frac{1}{m}\sum_{(\xbf, \mathbf{y})\in S} I(\sgn(\mathbf{w}_{l}\cdot L_{1:l-1}(\mathbf{x}) + b_{l}) \neq y)\\
&= \frac{1}{m}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:l-1}(\mathbf{x})=\mathbf{p}}}I(\sgn(\mathbf{w}_{l}\cdot \mathbf{p} + b_{l}) \neq y)\\
&\geq \frac{1}{m}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\argmin_{\mathbf{w},b}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:l-1}(\mathbf{x})=\mathbf{p}}}I(\sgn(\mathbf{w}\cdot \mathbf{p} + b) \neq y)\\
&= \frac{1}{m}\sum_{\pbf\in\mathcal{L}_{1:l-1}}\sum_{\substack{(\xbf, \mathbf{y})\in S \\ L_{1:l-1}(\mathbf{x})=\mathbf{p}}}I\left(\sgn\left(\overline{\ybf_{\pbf}^{(l-1)}}\right) \neq y\right)\\
&= \frac{1}{2}\left(1-\sum_{\pbf\in\mathcal{L}_{1:l-1}}\frac{\left|\ybf_{\pbf}^{(l-1)}\right|}{m}\left|\overline{\ybf_{\pbf}^{(l-1)}}\right|\right)\ .
\end{align}

The rest of the proof follows that of \autoref{prop:bin_bnd}.
\end{proof}

We obtain similar quantity to minimize in both \autoref{prop:bin_bnd} and \autoref{prop:bin_bnd_classif}; in order to minimize a lower bound on the training loss, the hyperplanes have to be placed such that the points being grouped together have similar labels. In the classification case, it simply tells us that if points with very disimilar labels are grouped together, necessarly, classification error will be made on several points of that group. In fact, another relation can be established; we know from \autoref{eq:simple_hyperplane}, when placing a neuron, in the regression context, the quantity to minimize in order to minimize the bound from \autoref{prop:bin_bnd}; it can be rewritten as follows:
$$
\begin{aligned}
\argmin_{\mathbf{w},b}&\left(\frac{|\ybf_{-1}|}{m}\textup{Var}\left(\ybf_{-1}\right)+\frac{|\ybf_{+1}|}{m}\textup{Var}\left(\ybf_{+1}\right)\right) \\
= \argmin_{\mathbf{w},b}&\left(\sum_{\substack{(\xbf, y)\in S \\ \mathbf{w}\cdot\mathbf{x}+b<0}}\left(\overline{\ybf_{-1}} - y\right)^2+\sum_{\substack{(\xbf, y)\in S \\ \mathbf{w}\cdot\mathbf{x}+b>0}}\left(\overline{\ybf_{+1}} - y\right)^2\right) \\
= \argmin_{\mathbf{w},b}&\left(\sum_{\substack{(\xbf, y)\in S \\ \mathbf{w}\cdot\mathbf{x}+b<0}} y^2-|\ybf_{-1}|\ \overline{\ybf_{-1}}^2+\sum_{\substack{(\xbf, y)\in S \\ \mathbf{w}\cdot\mathbf{x}+b>0}} y^2-|\ybf_{+1}|\ \overline{\ybf_{+1}}^2\right) \\
= \argmin_{\mathbf{w},b}&\bigg(-|\ybf_{-1}|\ \overline{\ybf_{-1}}^2-|\ybf_{+1}|\ \overline{\ybf_{+1}}^2\bigg) \\
= \argmax_{\mathbf{w},b}&\left(\frac{|\ybf_{-1}|}{m}\ \overline{\ybf_{-1}}^2+\frac{|\ybf_{+1}|}{m}\ \overline{\ybf_{+1}}^2\right).
\end{aligned}
$$
In the classification setting, from \autoref{prop:bin_bnd_classif}, we obtain these quantities to optimize:
$$
\begin{aligned}
\argmin_{\mathbf{w},b}\ \frac{1}{2}\left(1-\frac{\left|\ybf_{-1}\right|}{m}\left|\overline{\ybf_{-1}}\right|-\frac{\left|\ybf_{+1}\right|}{m}\left|\overline{\ybf_{+1}}\right|\right) &= \argmax_{\mathbf{w},b}\ \left(\frac{\left|\ybf_{-1}\right|}{m}\left|\overline{\ybf_{-1}}\right|+\frac{\left|\ybf_{+1}\right|}{m}\left|\overline{\ybf_{+1}}\right|\right).
\end{aligned}
$$
In short, it seems like in the regression case, with MSE loss function, in order to minimize the bound from \autoref{prop:bin_bnd}, when a hyperplane is placed, one must maximize the weighted sum of the \textbf{squared mean} (or minimize the variance) of the examples falling on each side of the hyperplane, whereas in the binary classification setting, with the 0-1 loss, in order to minimize the bound from \autoref{prop:bin_bnd_classif}, when a hyperplane is placed, one must maximize the weighted sum of the \textbf{absolute mean} of the examples falling on each side of the hyperplane. These results can be generalized to a whole layer.

\subsection{Proof of \autoref{theo:lb} (i) and (ii)}

The demonstration of \autoref{theo:lb} (i) and (ii) is done in three steps: we start by showing how binary activated neural networks (BANNs) can estimate a function $g$ defined as follows: $g(x) = x^2$, $x \in [0,1]$ (\autoref{prop:building_block_1}). Then, we take that first result in order to estimate the function $g(x,y) = xy$, $|x|,|y| < m$, with $m \in \mathbb{R}^*$ fixed (\autoref{prop:building_block_2}). Finally, we use those last two results as building blocks for estimating any function from $F_{n,d}$ (\autoref{theo:lb} (i) and (ii)). The first step is original, while the last two steps are adaptations of the work from \citet{DBLP:journals/nn/Yarotsky17} to BANNs.

\begin{prop}
Let $\epsilon>0$. There exists a fully connected binary activated neural network $B$ approximating $g(x)$ = $x^2$ on $x\in\textup{[0,1]}$
such that: (i) $\forall x \in [0,1] : |g(x) - B(x)| \leq \epsilon$; (ii) $B$ has only one hidden layer; (iii) the width of that hidden layer is in O($\epsilon^{-1}$); (iv) $g(x) - B(x) = 0 \forall x \in \{0,1\}$.
\label{prop:building_block_1}
\end{prop}

\begin{proof}
The proof will be done in two parts: we show that there exists a non-continuous piece-wise linear function formed of $r$ parallel lines $g^{(r)}(x)$ possessing the characteristic (iv), such that $|g(x) - g^{(r)}(x)| \leq (2r)^{-1}\forall x \in [0,1]$; then, we show that there exists a fully connected BANN having characteristics (ii) and (iii) and being capable of constructing $g^{(r)}(x)$. By using $\epsilon = (2r)^{-1}$, the proof will be completed.

\begin{figure}
\caption{Approximation of $f(x) = x^2, x\in [0,1]$ by an ensemble of $r$ parallel lines, corresponding to the output of a binary activated neural network.}
\label{fig:approx}
\end{figure}

First, when estimating $g(x) = x^2$ on $[0,1]$ with $r$ parallel lines, in order to obtain the minimum approximation error, the lines must be such that (i) they're equidistant and (ii) that distance must be twice the maximum error. We will also need that (iii) on the bounds on the considered interval for $x$, the approximation error is null. With that configuration (see \autoref{fig:approx}), each line is, for some values of $x$, the best estimator of the function of interest.

When estimating function $g(x)$ with the line $g{(1)}(x) = 0$, the maximum obtained error is $1$; which is the same when considering the line $g^{(1)}(x) = 1$. These lines will be bounding $g(x)$: that last function is always between the two lines. When adding lines to the estimator (thus, $r > 2$), we uniformly place these lines between $g^{(1)}(x) = 0$ and $g^{(1)}(x) = 1$. The maximum error thus corresponds to $1$ divided by $2r$:~the space between each line is now $\frac{1}{r}$, thus two times the maximum error. We can then estimate $g(x) = x^2$ on $[0,1]$ with an ensemble of $r$ straight parallel lines and obtaining $|g(x) - g^{(r)}(x)| \leq \frac{1}{2r} = (2r)^{-1} \forall x \in [0,1]$.

The ensemble of lines depicted on \autoref{fig:approx} is a predictor that can be written as follows:

\begin{equation}
g^{(r)}(x) = b_0 - \sum_{i=0}^{2(r-2)} w_{2i} \mathbb{I}_{\{xw_{2i+1} - b_{i+1} > 0\}}
\label{eq:func2}
\end{equation}
\begin{equation}
g^{(r)}(x) = (b_0+0.5) - \sum_{i=0}^{2(r-2)} \frac{w_{2i}}{2} \sgn(xw_{2i+1} - b_{i+1}).
\label{eq:func3}
\end{equation}

Both Equations (\ref{eq:func2}-\ref{eq:func3}) present a particular form: it's a sum of terms being the result of the application of a binary function, each of them having as input a linear function of $x$. The output of the network presented on \autoref{fig:res} is similar to Equations (\ref{eq:func2}-\ref{eq:func3}). Indeed, it represents the output of a fully connected single-hidden layer BANN. Also, Equations (\ref{eq:func2}-\ref{eq:func3}) explicitly show how $r$ directly affects the number of binary activation functions in the equations; $r$ is thus proportional to the number of neurons in the hidden layer of the network. Also, $\epsilon$ is inversely proportional to $r$.

\begin{figure}
\caption{A shallow but wide neural network; the architecture of a BANN having an output as presented in Equations (\ref{eq:func2}
\label{fig:res}
\end{figure}

Since that configuration respects (ii) and (iii) when choosing $r = \lceil(2\epsilon)^{-1}\rceil$, the proof of \autoref{prop:building_block_1} is completed.
\end{proof}

\begin{prop}
Let $x,y \in \mathbb{R}$ and $m > 0$. There exists a fully connected BANN $\times_m' : \mathbb{R}^2 \mapsto \mathbb{R}$, such that (i) $\forall|x|, |y| \leq m$, the error $\epsilon_{\times'} = |\times'(x,y) - xy|$ is upper bounded by $\epsilon$; (ii) the depth of the network is constant; (iii) the width of the first hidden layer of the network is in O($m^2\epsilon^{-1}$); (iv) $\epsilon_{\times'} = 0$ if $x = 0 \lor y = 0$.
\label{prop:building_block_2}
\end{prop}

\begin{proof}
Let $g^{(r)}_{\delta}(x)$ be the approximation function from \autoref{prop:building_block_1}, where $\epsilon_{g^{(r)}_{\delta}} \leq \delta$. The proof will make use of the demonstration of \autoref{prop:building_block_1} and the fact that $xy = \frac{(x+y)^2}{2}-\frac{x^2}{2}-\frac{y^2}{2}$. We have that
\begin{equation}
\begin{aligned}
\times_m'(x,y) &= m^2\left(2g^{(r)}_{\delta}\left(\frac{|x+y|}{m}\right) - \frac{1}{2}g^{(r)}_{\delta}\left(\frac{|x|}{m}\right) - \frac{1}{2}g^{(r)}_{\delta}\left(\frac{|y|}{m}\right)\right)\ .\\
\end{aligned}
\label{eq:funcprod}
\end{equation}

The upper bound of $\times_m'(x,y)$ as a function of $\delta$ can easily be found:
$$
\begin{aligned}
\epsilon_{\times'(x,y)} &\leq m^2(2\delta+\frac{1}{2}\delta+\frac{1}{2}\delta)\\
&= 3m^2\delta\ .
\end{aligned}
\label{eq:errbnd}
$$

\autoref{eq:funcprod} indicates that, in the end, the network estimating the function $g(x,y) = xy$ takes the form of a weighted sum of three sub-networks. Since the global network is constituted of a constant number of sub-networks, we have that (ii) is verified. Also, \autoref{eq:funcprod} confirms that (iv) is verified, considering that \autoref{prop:building_block_1} indicates that the error, on bounds of the considered interval for the input of the network, is 0.

As for (i), it is sufficient to choose $\delta = \frac{\epsilon}{3}$ for it to be verified. Finally, since the width of each sub-network is in O($\epsilon^{-1}$), (iii) is also verified.
\end{proof}

\paragraph{\autoref{theo:lb} (i) and (ii)}
\textit{For n, d $\in$ $\mathbb{N}^*$, there exists a fully connected BANN (i) capable of estimating any function from $F_{n,d}$ with absolute error upper bounded by $\epsilon$; (ii) having a first hidden layer with width in O$\left(c\epsilon^{-(d/n+1)}\right)$, where $c$ is a constant depending on $n$ and $d$, (iii) and having a constant depth.}

\begin{proof}
Let $f$ be a function from $F_{n,d}$. The proof will be done in two main parts: first, we approximate $f$ by $f_1$, a combination of Taylor polynomials of order $n - 1$ and one-dimensional piecewise linear functions while finding the error bound. We finish the proof by estimating $f_1$ by $g_1$, a binary activated feed-forward neural network, using \autoref{prop:building_block_2}.

We will need the use of a partition of unity on $\textup{[0,1]}^d$ : $\sum\nolimits_{\mathbf{m}} \psi_\mathbf{m} \equiv 1$. Here, $\mathbf{x} \in \textup{[0,1]}^d$, $\mathbf{m} = (m_1,\cdots, m_d) \in \{0,1,\cdots,N\}^d$ and $N \in \mathbb{N}^*$. $\psi_\mathbf{m}(\mathbf{x})$ is defined as follows:
$$
\psi_\mathbf{m}(\mathbf{x}) = \prod_{k=1}^dh(3Nx_k-3m_k)\ ,
$$
and 
$$
h(x) = \begin{aligned}
    	\begin{dcases}
    	1,&|x|\leq 1\\
    	2-|x|,&1<|x|<2\\
    	0,&|x|\geq 2
    	\end{dcases}
    	\end{aligned}\ .
$$

\begin{figure}
\caption{Figure taken from \citet{DBLP:journals/nn/Yarotsky17}
\label{part}
\end{figure}

See \autoref{part} for an example of $\psi_\mathbf{m}(\mathbf{x})$. Note that $||h||_\infty = 1$ and $||\psi_\mathbf{m}||_\infty = 1\ \forall \mathbf{m}$; that is: $$\underset{\mathbf{x}}{\textup{sup}}\ h = \underset{\mathbf{x}}{\textup{sup}}\ \psi_\mathbf{m} = 1\ \forall \mathbf{m}\,.$$ Also, we have that 
\begin{equation}
\begin{aligned}
\textup{supp}\ \psi_\mathbf{m} \subset \{\mathbf{x} : |x_k - \frac{m_k}{N}| < \frac{1}{N} \forall k\}\ ,
\label{theo_func_3}
\end{aligned}
\end{equation}
since supp $\psi_\mathbf{m} = \{\mathbf{x} \in \textup{[0,1]}^d | \psi_\mathbf{m}(\mathbf{x}) \neq 0\}$ and that $\psi_\mathbf{m}(\mathbf{x}) = 0$ when at least one component $x_k$ is such that $|x_k - \frac{m_k}{N}| \geq \frac{2}{3N}$.

Now, for all $\mathbf{m}$, consider the order $n-1$ Taylor polynomial for the function $f$ at $x = \frac{\mathbf{m}}{N}$:
$$
P_\mathbf{m}(\mathbf{x}) = \sum_{\mathbf{n}:|\mathbf{n}|<n} \frac{D^\mathbf{n}f}{\mathbf{n}!} _{\mathbf{x}=\frac{\mathbf{m}}{N}}\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\ .
$$

We will use the usual conventions $\mathbf{n}!$ = $\prod_{k=1}^d n_k!$ and $(\mathbf{x}-\frac{\mathbf{m}}{N})^\mathbf{n} = \prod_{k=1}^d (x_k-\frac{m_k}{N})^{n_k}$. We get an approximation $f_1$ of $f$ using $P_\mathbf{m}$ and $\psi_\mathbf{m}$ as follows:
$$
f_1 = \sum_{\mathbf{m}\in\{0,\cdots,N\}^d} \psi_\mathbf{m}P_\mathbf{m}\ .\\
$$

Therefore, $f_1$ is a sum of $(N+1)^d$ terms. Each of them is a multiplication of a linear term ($\psi_\mathbf{m}$) and a sum of $d^n$ terms ($P_\mathbf{m}$). Let's bound the approximation error of $f_1$:

\begin{equation}
\begin{aligned}
|f(\mathbf{x})-f_1(\mathbf{x})| &= \left|\sum_\mathbf{m}\psi_\mathbf{m}(\mathbf{x})(f(\mathbf{x})-P_\mathbf{m}(\mathbf{x}))\right|\\[10pt]
&\leq \sum_{\mathbf{m}}|\psi_\mathbf{m}(f(\mathbf{x})-P_\mathbf{m}(\mathbf{x}))|\\[10pt]
&\leq \sum_{\mathbf{m}:\left|x_k-\frac{m_k}{N}\right|<\frac{1}{N}\forall k} |\psi_\mathbf{m}(f(\mathbf{x})-P_\mathbf{m}(\mathbf{x}))|\\[10pt]
&\leq \sum_{\mathbf{m}:\left|x_k-\frac{m_k}{N}\right|<\frac{1}{N}\forall k} |f(\mathbf{x})-P_\mathbf{m}(\mathbf{x})|\\[10pt]
&\leq 2^d \underset{\mathbf{m}:\left|x_k-\frac{m_k}{N}\right|<\frac{1}{N}\forall k}{\textup{max}}  |f(\mathbf{x})-P_\mathbf{m}(\mathbf{x})|\\[10pt]
&\leq \frac{2^dd^n}{n!}\left(\frac{1}{N}\right)^n \underset{\mathbf{n}:|\mathbf{n}|\leq n}{\textup{max}} \underset{\mathbf{x}\in \textup{[0,1]}^d}{\textup{ess sup}}\ |D^\mathbf{n}f(\mathbf{x})|\\[10pt]
&\leq \frac{2^dd^n}{n!}\left(\frac{1}{N}\right)^n\ .\\
\end{aligned}
\label{theo_func_4}
\end{equation}

Considering \autoref{theo_func_4}, the first inequality simply uses a basic property of the absolute value function. The second one emerges from the fact that we now use a larger space for $\mathbf{m}$ than before (as illustrated in \autoref{theo_func_3}). The third inequality uses the fact that 
$\sup_\mathbf{x} \psi_\mathbf{m} = 1$. The fourth inequality illustrates that there are at most $2^d$ terms in the summation and bounds by the maximum value of $|f(\mathbf{x})-P_\mathbf{m}(\mathbf{x})|$. It is so because for a fixed value $\mathbf{x}$, at most $2^d$ combinations of $\mathbf{m}$ exist such that $\left|x_k-\frac{m_k}{N}\right|<\frac{1}{N}\forall k$. The fifth one uses a standard bound for the Taylor remainder. The last one uses the fact that we only consider functions not in the Sobolev space $\mathcal{W}^{n,\infty}(\textup{[0,1]}^d)$ but in the unit ball $F_{n,d}$; it uses \autoref{eq:unit_ball}.\\

The final manipulation of the first part of the proof is to choose wisely the value of $N$; when choosing
$$
\begin{aligned}
\frac{\epsilon}{2} &= \frac{2^dd^n}{n!}\left(\frac{1}{N}\right)^n\\
\Rightarrow N &= \left\lceil\left(\frac{n!}{2^{d+1}d^n}\frac{\epsilon}{2}\right)^{-1/n}\right\rceil\ ,
\end{aligned}
$$
it follows that 
$$
||f-f_1||_\infty \leq \frac{\epsilon}{2}\ .
$$

Note that, because of how $F_{n,d}$ is defined (see \autoref{eq:unit_ball}), every coefficient of the polynomials $P_\mathbf{m}$ (that is, every term $\frac{D^\mathbf{n}f}{\mathbf{n}!} |_{\mathbf{x}=\frac{\mathbf{m}}{N}}$) is uniformly bounded by 1 $\forall f \in F_{n,d}$. Therefore, expanding $f_1$ leads to the following:
\begin{equation}
f_1(\mathbf{x}) = \sum_{\mathbf{m}\in\{0,1,\cdots,N\}^d}\sum_{\mathbf{n}:|\mathbf{n}|<n} a_{\mathbf{m},\mathbf{n}}\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n},\ |a_{\mathbf{m},\mathbf{n}}|\leq 1\ .
\label{theo_func_5}
\end{equation}

Since every $a_{\mathbf{m},\mathbf{n}}$ is constant, $f_1$ is a linear combination of at most (as mentioned before) $d^n(N+1)^d$ terms $\psi_\mathbf{m}(\mathbf{x})(\mathbf{x}-\frac{\mathbf{m}}{N})^\mathbf{n}$\ . Each of these terms is a product of at most $d+n-1$ piece-wise linear univariate factors: $d$ functions $h(3Nx_k-3m_k)$ which compose $\psi_\mathbf{m}$ and at most $n-1$ linear expressions $x_k-\frac{m_k}{N}$ which compose $P_\mathbf{m}(\mathbf{x})$.

We now approximate $f_1$ with a sum of $d^n(N+1)^d$ binary activated feed-forward neural networks, each of them approximating a term $\psi_\mathbf{m}(\mathbf{x})(\mathbf{x}-\frac{\mathbf{m}}{N})^\mathbf{n}$\ by applying recursively the function $\times'_{m,\delta}$ from \autoref{prop:building_block_2}, with an arbitrary parameter $m$ and maximum error rate $\delta$ to be chosen later, leading us to the following:
$$
\Tilde{f}_{\mathbf{m},\mathbf{n}} = \times'_{m,\delta}(h(3Nx_1-3m_1),\times'_{m,\delta}(h(3Nx_2-3m_2),\cdots,\times'_{m,\delta}(x_k-\frac{m_k}{N},\cdots)\cdots))\ .
\label{theo_func_6}
$$
The value of $m$ is then chosen so that when using $\times'_{m,\delta}(c,d)$, both $|c|<m$ and $|d|<m$ are respected. That is, we must find the largest number that could possibly be given as an argument to one of the functions $\times'$ that appears in $\Tilde{f}_{\mathbf{m},\mathbf{n}}$. Since $\sup_{\mathbf{x}}= \sup_{\mathbf{x}} \psi_\mathbf{m} = 1\ \forall \mathbf{m}$, the first application of a function $\times'$ will be on arguments bounded by $1$. On the second application, the first argument is bounded by $1$ and the other, which is the result of the first application, bounded by $1+\delta$, and so on. On the last application, the first argument is bounded by 1 while the other one is bounded by $1+(d+n-3)\delta$. Since $\delta < 1$, choosing $m = d+n-2$ covers the worst case possible. Also, since there are $d+n-2$ recursive calls, the depth of the network estimating $\Tilde{f}_{\mathbf{m},\mathbf{n}}$ is constant, but depends on $d$ and $n$ (depth $\propto d+n$).

Since $\Tilde{f}_{\mathbf{m},\mathbf{n}}$ contains $d+n-2$ applications of $\times'$, whose error is bounded by $\delta$, the maximum approximation error for $\left|\Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\right|$ is given by $(d+n-2)\delta$.

Now, all we need to do is define the full approximation of $f$ by $\Tilde{f}$ and estimate the approximation errorlast two
\begin{equation}
\Tilde{f} = \sum_{\mathbf{m}\in\{0,1,\cdots,N\}^d}\sum_{\mathbf{n}:|\mathbf{n}|<n} a_{\mathbf{m},\mathbf{n}}\Tilde{f}_{\mathbf{m},\mathbf{n}}\\
\label{theo_func_7}
\end{equation}

\begin{equation}
\begin{aligned}
|\Tilde{f}(\mathbf{x})-f_1(\mathbf{x})| &= \left|\sum_{\mathbf{m}\in\{0,1,\cdots,N\}^d}\sum_{\mathbf{n}:|\mathbf{n}|<n} a_{\mathbf{m},\mathbf{n}}\left(\Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\right)\right|\\
&= \left|\sum_{\mathbf{m}:\mathbf{x}\in \textup{supp} \psi_\mathbf{m}}\sum_{\mathbf{n}:|\mathbf{n}|<n} a_{\mathbf{m},\mathbf{n}}\left(\Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\right)\right|\\
&\leq \sum_{\mathbf{m}:\mathbf{x}\in \textup{supp} \psi_\mathbf{m}}\sum_{\mathbf{n}:|\mathbf{n}|<n}\left| a_{\mathbf{m},\mathbf{n}}\left(\Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\right)\right|\\
&\leq \sum_{\mathbf{m}:\mathbf{x}\in \textup{supp} \psi_\mathbf{m}}\sum_{\mathbf{n}:|\mathbf{n}|<n}\left| \left(\Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\right)\right|\\
&\leq 2^d \underset{\mathbf{m}:x\in \textup{supp}
\psi_\mathbf{m}}{\textup{max}}\sum_{\mathbf{n}:|\mathbf{n}|<n}\left| \left(\Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x}-\frac{\mathbf{m}}{N}\right)^\mathbf{n}\right)\right|\\
&\leq 2^dd^n(d+n-2)\delta\ .
\label{theo_func_8}
\end{aligned}
\end{equation}
In \autoref{theo_func_8}, the second equality illustrates that $\mathbf{x} \notin supp \psi_\mathbf{m} \Rightarrow \Tilde{f}_{\mathbf{m},\mathbf{n}}-\psi_\mathbf{m}(\mathbf{x})\left(\mathbf{x} - \frac{\mathbf{m}}{N}\right)^\mathbf{n} = 0$, coming from (iv) of \autoref{prop:building_block_2}. The first inequality comes from a simple property of the absolute value function. The second one illustrates the fact that $|a_{\mathbf{m},\mathbf{n}}| \leq 1$, as pointed out in \autoref{theo_func_5}. The third inequality consider that for any $\mathbf{x}$, there are at most $2^d \psi_\mathbf{m}(\mathbf{x})$ that are not equal to 0. The fourth one considers that every $\mathbf{n}$ such that $|\mathbf{n}|<n$ leads to the worst case (there are $d^n$ of them), where each one is equal to $(d+n)$.

Finally, if we choose
$$
\delta = \frac{\epsilon}{2^{d+1}d^n(d+n-2)}\ ,\\
$$
for every $\times'_{m,\delta}$, then $||\Tilde{f}-f_1||_\infty \leq \epsilon/2$ and (iii) is respected. By the triangular inequality, we have 
$$
\begin{aligned}
||\Tilde{f}-f||_\infty &\leq ||\Tilde{f}-f_1||_\infty + ||f_1-f||_\infty\\
&= \epsilon/2 + \epsilon/2\\
&= \epsilon\ .\\
\end{aligned}
$$

As illustrated in \autoref{theo_func_7}, $\Tilde{f}$ can be implemented by a weighted sum of $(N+1)^dd^n$ parallel binary activated feed-forward neural networks; as shown in \autoref{prop:building_block_2}, each of these network has width in $O\left(\frac{1}{\delta}\right)$. Therefore, the total network now has a width that is $(N+1)^dd^n$ times bigger. Let's not forget that we chose $\delta = \frac{\epsilon}{2^{d+1}d^n(d+n)}$; the new width of the first layer is in O$\left(\frac{N^dd^n2^{d+1}d^n(d+n)}{\epsilon}\right)$ = O$\left(\frac{(2N)^dd^{2n}(d+n)}{\epsilon}\right)$. The upper bound of the depth will remain constant (a function of $d+n$), since doing an aggregation of weighted network will only add one layer.

Since we already fixed $N$, we can expand our width upper bound that way:
$$
\begin{aligned}
N^d\frac{d^{2n}(d+n)}{\epsilon} &= \left\lceil\left(\frac{n!}{2^dd^n}\frac{\epsilon}{2}\right)^{-1/n}\right\rceil^d\frac{d^{2n}(d+n)}{\epsilon}\\
&\approx  \left(\frac{n!}{2^dd^n}\frac{\epsilon}{2}\right)^{-d/n}\frac{d^{2n}(d+n)}{\epsilon}\\
&=  \left(\frac{n!}{2^{d+1}d^n}\right)^{-d/n}d^{2n}(d+n)\epsilon^{-(d/n+1)}\\
&=  c\epsilon^{-(d/n+1)}\ ,
\end{aligned}
\label{funcd}
$$
where $c$ is a constant depending on fixed numbers; that is: $n$ and $d$. Therefore, the width is in O($c\epsilon^{-(d/n+1)}$), concluding the proof.
\end{proof}

\subsection{Proof of \autoref{theo:lb} (iii)}

The following proof is inspired by the work of \citet{DBLP:journals/nn/Yarotsky17}. We will need the following technical lemma:

\begin{lemma}
Let a $>$ 0. Then: x $\geq$ 2a log(a) $\Rightarrow$ x $\geq$ a log(x). It follows that a necessary condition for the inequality x $<$ a log(x) to hold is that x $<$ 2a log(a).\\
\label{lemma1}
\end{lemma}

\paragraph{\autoref{theo:lb} (iii)} \textit{For d, n $\in$ $\mathbb{N}^*$, there exists a binary activated network (i) which is capable of estimating any function from $F_{n,d}$ with maximum error $\epsilon$; (ii) for which the number of weights $W$ is $\Omega\left(c\epsilon^{-d/n}\textup{ln}^{-1}\left(\epsilon^{-1}\right)\right)$, where $c$ is a constant depending on $n$ and $d$.}

\begin{proof}
Let $H$ be a class of Boolean functions on [0,1]$^d$. We are interested in the case where $H$ is the family of functions obtained when applying a threshold function to a binary network with fixed architecture but variable weights. We will need both the concept of VC-dimension and the fact that, in our case:

\begin{equation}
\textup{VCdim}(H) \in O(W ln(W))\ ,
\label{vc}
\end{equation}

where $W$ is the number of weights in our network. Consider N, a positive integer to be chosen later. Now, select $S$, a set of $N^d$ points $\mathbf{x}_1,\cdots,\mathbf{x}_{N^d}$ in the hypercube $[0,1]^d$ such that the distance between any two points from $S$ is not less than $\frac{1}{N}$. Given any $\mathbf{y} = [y_1,\cdots,y_{N^d}] \in \mathbb{R}^{N^d}$, we can construct a smooth function $f$ satisfying $f(\mathbf{x}_m) = y_m \forall m$ by setting
$$
f(\mathbf{x}) = \sum_{m=1}^{N^d} y_m\phi(N(\mathbf{x}-\mathbf{x}_m))\ ,
$$

with some $C^\infty$ function $\phi : \mathbb{R}^d \rightarrow \mathbb{R}$ such that $\phi(\mathbf{0}) = 1$ and $\phi(\mathbf{x}) = 0$ if $|\mathbf{x}| > 1/2$. Now, we need to make sure that $f \in F_{n,d}$:
$$
\begin{aligned}
\underset{\mathbf{x}}{\textup{max}}|D^\mathbf{n}f(\mathbf{x})| &= \underset{\mathbf{x}}{\textup{max}}\left|D^\mathbf{n}\left(\sum_{m=1}^{N^d} y_m\phi(N(\mathbf{x}-\mathbf{x}_m))\right)\right| \\\\
&= \underset{\mathbf{x}}{\textup{max}}\left|\sum_{m=1}^{N^d} y_mD^\mathbf{n}\phi(N(\mathbf{x}-\mathbf{x}_m))\right| \\\\
&= \underset{\mathbf{x}}{\textup{max}}\left|\underset{m}{\textup{max}}\ y_mD^\mathbf{n}\phi(N(\mathbf{x}-\mathbf{x}_m))\right| \\\\
&= \underset{m}{\textup{max}}\ |y_m|\ \underset{\mathbf{x}}{\textup{max}}\left|D^\mathbf{n}\phi(N(\mathbf{x}-\mathbf{x}_m))\right| \\
&= N^{|\mathbf{n}|}\underset{m}{\textup{max}}\ |y_m|\underset{\mathbf{x}}{\textup{max}}\left|D^\mathbf{n}\phi(\mathbf{x}-\mathbf{x}_m)\right| \\
&= N^{|\mathbf{n}|}\underset{m}{\textup{max}}\ |y_m|\underset{\mathbf{x}}{\textup{max}}\left|D^\mathbf{n}\phi(\mathbf{x})\right|\ .\\
\end{aligned}
$$

The third equality comes from the fact that there are no points in $S$ that are less than $\frac{1}{N}$-far from each other; therefore, $\phi$ will yield $1$ for only one term of the summation, which corresponds to the $max$ of the function. Thus, the $x$ maximizing $y_mD^\mathbf{n}\phi(N(\mathbf{x}-\mathbf{x}_m))$ corresponds to the $m$ maximizing $y_m$.

We got to make sure that $f \in F_{n,d}$: having $\underset{m}{\textup{max}} |y_m| \leq cN^{-n}$, where $c = (\textup{max}_{\mathbf{n}:|\mathbf{n}|\leq n}\textup{max}_x|D^\mathbf{n}\phi(\mathbf{x})|)^{-1}$ is the first step. We now set $\epsilon = \frac{cN^{-n}}{3}$. We suppose that there exists a binary activated feed-forward neural network that can approximate any function from $F_{n,d}$ with maximum error $\epsilon$. We will denote by $\nu(\mathbf{x},\mathbf{w})$ the output of the network for input $\mathbf{x}$ and weights $\mathbf{w}$. Also, let $\mathbf{z} \in \{0,1\}^{N^d}$ be any set of boolean $N^d$ values. We set $y_m = z_mcN^{-n} \forall m \in \{1,2,\cdots,N^d\}$. Therefore, $f \in F_{n,d}$.

By assumption, there is then a vector of weights $\mathbf{w}$ = $\mathbf{w}_\mathbf{z}$ such that for all $m$, we have $|\nu(\mathbf{x}_m,\mathbf{w}_\mathbf{z})-y_m|\leq\epsilon$ and 
$$
\nu(\mathbf{x}_m,\mathbf{w}_\mathbf{z}) \begin{aligned}
                                        \begin{dcases}
                                    	\geq cN^{-n}-\epsilon > cN^{-n}/2,&z_m = 1\\
                                    	\leq \epsilon < cN^{-n}/2,&z_m = 0\ ,
                                    	\end{dcases}
                                    	\end{aligned}
$$

so the thresholded network $\nu(\mathbf{x}_m,\mathbf{w}_\mathbf{z})_1 = 1_{\{\nu(\mathbf{x}_m,\mathbf{w}_\mathbf{z}) > cN^{-n}/2\}}$ outputs $z_m \forall m \in \{1,2,\cdots,N^d\}$. We thus conclude that $z_m$ is shattered and that 
$$
VCdim(\nu_1) \geq N^d\ .
$$

We previously defined $\epsilon = \frac{cN^{-n}}{3}$; by replacing $N$ in the inequality, we obtain
$$
VCdim(\nu_1) \geq \left(\frac{3\epsilon}{c}\right)^{-d/n}\ .
$$

Using both \autoref{lemma1} and \autoref{vc}, since $W > 0$ and $\left(\frac{3\epsilon}{c}\right)^{-d/n} > 1$, we can bound $W$ as follows:
$$
\begin{aligned}
W ln(W) &\geq \left(\frac{3\epsilon}{c}\right)^{-d/n}\\
\Rightarrow 2W ln(W) &> \left(\frac{3\epsilon}{c}\right)^{-d/n}\\
\Rightarrow W ln\left(\left(\frac{3\epsilon}{c}\right)^{-d/n}\right)&> \left(\frac{3\epsilon}{c}\right)^{-d/n}\\
\Rightarrow W &> \left(\frac{3\epsilon}{c}\right)^{-d/n}ln^{-1}\left(\left(\frac{3\epsilon}{c}\right)^{-d/n}\right)\\
\Rightarrow W &> \left(\frac{3\epsilon}{c}\right)^{-d/n}\left(-\frac{d}{n}ln\left(\frac{3\epsilon}{c}\right)\right)^{-1}\\
\Rightarrow W &> c^*\epsilon^{-d/n}ln^{-1}\left(\epsilon^{-1}\right)
\end{aligned}
$$

which concludes our proof.
\end{proof}

\subsection{Proof of \autoref{prop:BGN_min}}

\paragraph{Proposition 3.} \textit{In} \autoref{algo:bgn}: $(c_t, d_t) = \argmin_{c, d}\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)^2$.

\begin{proof}
Let's find the values of $c$ and $d$ in $\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)^2$ for which the partial derivatives  are zeros. We denote these values $\hat{c}$ and $\hat{d}$ respectively.

$$
\begin{aligned}
\allowdisplaybreaks
\frac{\partial}{\partial c}\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)^2 &= -2\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)\\
\Rightarrow 0 &= -2\sum_i \left(y_i^{(t)}-\hat{c}\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-\hat{d}\right)\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)\\
\Rightarrow 0 &= \sum_i y_i^{(t)}\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-n\hat{c}-\sum_i \hat{d}\cdot\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)\\
\Rightarrow \hat{c} &= \frac{1}{n}\sum_i (y_i^{(t)}-\hat{d})\cdot\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)\\
\end{aligned}
$$

$$
\begin{aligned}
\frac{\partial}{\partial d}\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)^2 &= -2\sum_i \left(y_i^{(t)}-c\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-d\right)\\
\Rightarrow 0 &= -2\sum_i \left(y_i^{(t)}-\hat{c}\ \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)-\hat{d}\right)\\
\Rightarrow \hat{c} &= \frac{\sum_i y_i^{(t)}-n\hat{d}}{\sum_i \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)}\\
\end{aligned}
$$

In order to lighten the notation, we write $\sum_{\pm1}$ to denote $\sum_{i:\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)=\pm1}$.

$$
\begin{aligned}
\frac{1}{n}\sum_i (y_i^{(t)}-\hat{d})\sgn(\mathbf{w}_t\cdot\xbf_i+b_t) &= \frac{\sum_i y_i^{(t)}-n\hat{d}}{\sum_i \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)}\\
\Rightarrow \frac{(\sum_{+1}y_i^{(t)}-\sum_{-1}y_i^{(t)})-\hat{d}(\sum_{+1}1-\sum_{-1}1)}{\sum_{+1}1+\sum_{-1}1} &= \frac{(\sum_{+1}y_i^{(t)}+\sum_{-1}y_i^{(t)})-\hat{d}(\sum_{+1}1+\sum_{-1}1)}{\sum_{+1}1-\sum_{-1}1}\\
\Rightarrow \hat{d}[(\sum_{+1}1+\sum_{-1}1)^2-(\sum_{+1}1-\sum_{-1}1)^2] &= (\sum_{+1}y_i^{(t)}+\sum_{-1}y_i^{(t)})(\sum_{+1}1+\sum_{-1}1)\\
&-(\sum_{+1}y_i^{(t)}-\sum_{-1}y_i^{(t)})(\sum_{+1}1-\sum_{-1}1)\\
\Rightarrow 4\hat{d}\sum_{+1}1\sum_{-1}1 &= 2\sum_{+1}y_i^{(t)}\sum_{-1}1+2\sum_{-1}y_i^{(t)}\sum_{+1}1\\
\Rightarrow \hat{d} &= \frac{1}{2}\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}+\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)
\end{aligned}
$$

$$
\begin{aligned}
\hat{c} &= \frac{\sum_i y_i^{(t)}-n\hat{d}}{\sum_i \sgn(\mathbf{w}_t\cdot\xbf_i+b_t)}\\
\Rightarrow \hat{c} &= \frac{\sum_{+1} y_i^{(t)}+\sum_{-1} y_i^{(t)}}{\sum_{+1}1-\sum_{-1}1}-\hat{d}\frac{\sum_{+1}1+\sum_{-1}1}{\sum_{+1}1-\sum_{-1}1}\\
\Rightarrow \hat{c} &= \frac{\sum_{+1} y_i^{(t)}+\sum_{-1} y_i^{(t)}}{\sum_{+1}1-\sum_{-1}1}-\frac{1}{2}\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}+\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)\frac{\sum_{+1}1+\sum_{-1}1}{\sum_{+1}1-\sum_{-1}1}\\
\Rightarrow \hat{c}(\sum_{+1}1-\sum_{-1}1) &= \sum_{+1}y_i^{(t)}+\sum_{-1}y_i^{(t)}-\frac{1}{2}\frac{\sum_{+1}y_i^{(t)}\sum_{-1}1+\sum_{-1}y_i^{(t)}\sum_{+1}1}{\sum_{+1}1\sum_{-1}1}(\sum_{+1}1+\sum_{-1}1)\\
\Rightarrow \hat{c}(\sum_{+1}1-\sum_{-1}1) &= \frac{1}{2}\sum_{+1}y_i^{(t)}+\frac{1}{2}\sum_{-1}y_i^{(t)}-\frac{1}{2}\frac{\sum_{+1}y_i^{(t)}\sum_{-1}1}{\sum_{+1}1}-\frac{1}{2}\frac{\sum_{-1}y_i^{(t)}\sum_{+1}1}{\sum_{-1}1}\\
\Rightarrow \hat{c}(\sum_{+1}1-\sum_{-1}1) &= \frac{1}{2}\frac{\sum_{+1}y_i^{(t)}(\sum_{+1}1-\sum_{-1}1)}{\sum_{+1}1}-\frac{1}{2}\frac{\sum_{-1}y_i^{(t)}(\sum_{+1}1-\sum_{-1}1)}{\sum_{-1}1}\\
\Rightarrow \hat{c} &= \frac{1}{2}\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}-\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}\right)\\
\end{aligned}
$$
\end{proof}

\subsection{Proof of \autoref{prop:BGN_dec}}

\paragraph{Proposition 4.} \textit{In \autoref{algo:bgn}, we have }
\begin{center}
$Q(BGN_{t-1},S) - Q(BGN_{t},S) = (c_t+d_t)(c_t-d_t) > 0, \ \forall t \in \mathbb{N}\setminus\{0,1\}.$
\end{center}

\begin{proof}
Let's use the shortcut notation $= \sum_{\pm1}$ to denote $\sum_{i:\sgn(\mathbf{w}_t\cdot\xbf_i+b_t)=\pm1} $.  
$$
\begin{aligned}
d_t-c_t &= \frac{1}{2}\left(\frac{\sum_{+1} y_i}{\sum_{+1} 1}+\frac{\sum_{-1} y_i}{\sum_{-1} 1}\right) - \frac{1}{2}\left(\frac{\sum_{+1} y_i}{\sum_{+1} 1}-\frac{\sum_{-1} y_i}{\sum_{-1} 1}\right)\\
&= \frac{\sum_{-1} y_i}{\sum_{-1} 1}\ ,\\\
d_t+c_t &= \frac{1}{2}\left(\frac{\sum_{+1} y_i}{\sum_{+1} 1}+\frac{\sum_{-1} y_i}{\sum_{-1} 1}\right) + \frac{1}{2}\left(\frac{\sum_{+1} y_i}{\sum_{+1} 1}-\frac{\sum_{-1} y_i}{\sum_{-1} 1}\right)\\
&= \frac{\sum_{+1} y_i}{\sum_{+1} 1}\ .
\end{aligned}
$$

We first show that $Q(BGN_{t-1},S) - Q(BGN_{t},S) > 0, \ \forall t \in \mathbb{N}\setminus\{0,1\}$.

$$
\begin{aligned}
&\ \quad Q(BGN_{t-1},S) - Q(BGN_{t},S)\\
&= \frac{1}{m}\sum_{i=1}^{m}\left(y_i-BGN_{t-1}(\mathbf{x}_i)\right)^2 - \frac{1}{m}\sum_{i=1}^{m}\left(y_i-BGN_{t}(\mathbf{x}_i)\right)^2\\
&= \frac{1}{m}\sum_{i=1}^{m}\left(y^{(t)}_i\right)^2 - \frac{1}{m}\sum_{i=1}^{m}\left(y^{(t)}_i-c_th_t(\mathbf{x}_i)-d_t\right)^2\\
&= \frac{1}{m}\sum_{i=1}^{m}\left(y^{(t)}_i\right)^2 - \frac{1}{m}\sum_{i=1}^{m}\left(\left(y^{(t)}_i\right)^2-2y^{(t)}_i(c_th_t(\mathbf{x}_i)+d_t)+(c_th_t(\mathbf{x}_i)+d_t)^2\right)\\
&= \frac{1}{m}\sum_{i=1}^{m}\left(2y^{(t)}_i(c_th_t(\mathbf{x}_i)+d_t)-(c_th_t(\mathbf{x}_i)+d_t)^2\right)\\
&= \textcolor{red}{c_t\frac{2}{m}\sum_{i=1}^{m}y^{(t)}_ih_t(\mathbf{x}_i)+d_t\frac{2}{m}\sum_{i=1}^{m}y^{(t)}_i}\textcolor{blue}{-\frac{1}{m}\sum_{i=1}^{m}c_t^2\underset{=1\forall\mathbf{x}}{\underbrace{h_t^2(\mathbf{x}_i)}}-\frac{1}{m}\sum_{i=1}^{m}d_t^2}\textcolor{brown}{-c_td_t\frac{2}{m}\sum_{i=1}^{m}h_t(\mathbf{x}_i)}\\
&= \frac{1}{m}\Bigg[\textcolor{red}{\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}-\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)\left(\sum_{+1}y_i-\sum_{-1}y_i\right)+\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}+\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)\left(\sum_{+1}y_i+\sum_{-1}y_i\right)}\\
&\ \quad\textcolor{blue}{-\frac{1}{4}\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}-\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)^2\left(\sum_{+1}1+\sum_{-1}1\right)-\frac{1}{4}\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}+\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)^2\left(\sum_{+1}1+\sum_{-1}1\right)}\\
&\ \quad\textcolor{brown}{-\frac{1}{2}\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}-\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}+\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)\left(\sum_{+1}1-\sum_{-1}1\right)}\Bigg]\\
&= \frac{1}{m}\Bigg[\textcolor{red}{2\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\sum_{+1}1}+2\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\sum_{-1}1}}\textcolor{blue}{-\frac{1}{2}\left(\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\left(\sum_{+1}1\right)^2}+\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\left(\sum_{-1}1\right)^2}\right)\left(\sum_{+1}1+\sum_{-1}1\right)}\\
&\ \quad\textcolor{brown}{-\frac{1}{2}\left(\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\left(\sum_{+1}1\right)^2}-\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\left(\sum_{-1}1\right)^2}\right)\left(\sum_{+1}1-\sum_{-1}1\right)}\Bigg]\\
&= \frac{1}{m}\Bigg[2\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\sum_{+1}1}+2\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\sum_{-1}1}-\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\sum_{+1}1}-\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\sum_{-1}1}\Bigg]\\
&= \frac{1}{m}\Bigg[\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\sum_{+1}1}+\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\sum_{-1}1}\Bigg]>0\\
\end{aligned}
$$

Since we have that 
$$
\left(\sum_{+1}y_i^{(t)}\right)^2\left(\sum_{-1}1\right)+\left(\sum_{-1}y_i^{(t)}\right)^2\left(\sum_{+1}1\right)+\left(\sum_{+1}y_i^{(t)}\right)\left(\sum_{-1}y_i^{(t)}\right)\left(\sum_{+1}1+\sum_{-1}1\right)=0,
$$

then 
$$
\begin{aligned}
Q(BGN_{t-1},S) - Q(BGN_{t},S) &= \frac{1}{m}\Bigg[\frac{\left(\sum_{+1}y_i^{(t)}\right)^2}{\sum_{+1}1}+\frac{\left(\sum_{-1}y_i^{(t)}\right)^2}{\sum_{-1}1}\Bigg]\\
&= \frac{1}{m}\Bigg[\frac{\left(\sum_{+1}y_i^{(t)}\right)^2\left(\sum_{-1}1\right)+\left(\sum_{+1}y_i^{(t)}\right)^2\left(\sum_{-1}1\right)}{\left(\sum_{+1}1\right)\left(\sum_{-1}1\right)}\Bigg]\\
&= -\frac{1}{m}\Bigg[\frac{m\left(\sum_{+1}y_i^{(t)}\right)\left(\sum_{+1}y_i^{(t)}\right)}{\left(\sum_{+1}1\right)\left(\sum_{-1}1\right)}\Bigg]\\
&= -\left(\frac{\sum_{+1}y_i^{(t)}}{\sum_{+1}1}\right)\left(\frac{\sum_{-1}y_i^{(t)}}{\sum_{-1}1}\right)\\
&= -(d_t+c_t)(d_t-c_t)\\
&= (c_t+d_t)(c_t-d_t)\\
\end{aligned}
$$
\end{proof}

\end{document}