Expressive Monotonic Neural Networks
The monotonic dependence of the outputs of a neural network on some of its
inputs is a crucial inductive bias in many scenarios where domain knowledge
dictates such behavior. This is especially important for interpretability and
fairness considerations. In a broader context, scenarios in which monotonicity
is important can be found in finance, medicine, physics, and other disciplines.
It is thus desirable to build neural network architectures that implement this
inductive bias provably. In this work, we propose a weight-constrained
architecture with a single residual connection to achieve exact monotonic
dependence in any subset of the inputs. The weight constraint scheme directly
controls the Lipschitz constant of the neural network and thus provides the
additional benefit of robustness. Compared to currently existing techniques
used for monotonicity, our method is simpler in implementation and in theory
foundations, has negligible computational overhead, is guaranteed to produce
monotonic dependence, and is highly expressive. We show how the algorithm is
used to train powerful, robust, and interpretable discriminators that achieve
competitive performance compared to current state-of-the-art methods across
various benchmarks, from social applications to the classification of the
decays of subatomic particles produced at the CERN Large Hadron Collider.
Ouail Kitouni, Niklas Nolte, Michael Williams
# Expressive Monotonic Neural Networks

###### Abstract

The monotonic dependence of the outputs of a neural network on some of its inputs is a crucial inductive bias in many scenarios where domain knowledge dictates such behavior. This is especially important for interpretability and fairness considerations. In a broader context, scenarios in which monotonicity is important can be found in finance, medicine, physics, and other disciplines. It is thus desirable to build neural network architectures that implement this inductive bias provably. In this work, we propose a weight-constrained architecture1 with a single residual connection to achieve exact monotonic dependence in any subset of the inputs. The weight constraint scheme directly controls the Lipschitz constant of the neural network and thus provides the additional benefit of robustness. Compared to currently existing techniques used for monotonicity, our method is simpler in implementation and in theory foundations, has negligible computational overhead, is guaranteed to produce monotonic dependence, and is highly expressive. We show how the algorithm is used to train powerful, robust, and interpretable discriminators that achieve competitive performance compared to current state-of-the-art methods across various benchmarks, from social applications to the classification of the decays of subatomic particles produced at the CERN Large Hadron Collider.

Footnote 1: [https://github.com/niklasnolte/MonotoneNorm](https://github.com/niklasnolte/MonotoneNorm)

## 1 Introduction

The need to model functions that are monotonic in a subset of their inputs is prevalent in many ML applications. Enforcing monotonic behaviour can help improve generalization capabilities (Milani Fard et al., 2016; You et al., 2017) and assist with interpretation of the decision-making process of the neural network (Nguyen Martinez, 2019). Real world scenarios include various applications with fairness, interpretability, and security aspects. Examples can be found in the natural sciences and in many social applications. Monotonic dependence of a model output on a certain feature in the input can be informative of how an algorithm works--and in some cases is essential for real-word usage. For instance, a good recommender engine will favor the product with a high number of reviews over another with fewer but otherwise identical reviews (_ceteris paribus_). The same applies for systems that assess health risk, evaluate the likelihood of recidivism, rank applicants, filter inappropriate content, _etc_.

In addition, robustness to small perturbations in the input is a desirable property for models deployed in real world applications. In particular, when they are used to inform decisions that directly affect human actors--or where the consequences of making an unexpected and unwanted decision could be extremely costly. The continued existence of adversarial methods is a good example for the possibility of malicious attacks on current algorithms (Akhtar et al., 2021). A natural way of ensuring the robustness of a model is to constrain its Lipschitz constant. To this end, we recently developed an architecture whose Lipschitz constant is constrained by design using layer-wise normalization which allows the architecture to be more expressive than the current state-of-the-art with stable and fast training (Kitouni et al., 2021). Our algorithm has been adopted to classify the decays of subatomic
particles produced at the CERN Large Hadron Collider in the real-time data-processing system of the LHCb experiment, which was our original motivation for developing this novel architecture.

In this paper, we present expressive monotonic Lipschitz networks. This new class of architectures employs the Lipschitz bounded networks from Kitouni et al. (2021) along with residual connections to implement monotonic dependence in any subset of the inputs by construction. It also provides exact robustness guarantees while keeping the constraints minimal such that it remains a universal approximator of Lipschitz continuous monotonic functions. We show how the algorithm is used to train powerful, robust, and interpretable discriminators that achieve competitive performance compared to current state-of-the-art methods across various benchmarks, from social applications to its original target application: the classification of the decays of subatomic particles produced at the CERN Large Hadron Collider.

## 2 Related Work

Prior work in the field of monotonic models can be split into two major categories.

* **Built-in and constrained monotonic architectures**: Examples of this category include Deep Lattice Networks (You et al., 2017) and networks in which all weights are constrained to have the same sign (Sill, 1998). The major drawbacks of most implementations of constrained architectures are a lack of expressiveness or poor performance due to superfluous complexity.
* **Heuristic and regularized architectures (with or without certification)**: Examples of such methods include Sill & Abu-Mostafa (1996) and Gupta et al., which penalizes point-wise negative gradients on the training sample. This method works on arbitrary architectures and retains much expressive power but offers no guarantees as to the monotonicity of the trained model. Another similar method is Liu et al. (2020), which relies on Mixed Integer Linear Programming to certify the monotonicity of piece-wise linear architectures. The method uses a heuristic regularization to penalize the non-monotonicity of the model on points sampled uniformly in the domain during training. The procedure is repeated with increasing regularization strength until the model passes the certification. This iteration can be expensive and while this method is more flexible than the constrained architectures (valid for MLPs with piece-wise linear activations), the computational overhead of the certification process can be prohibitively expensive. Similarly, Sivaraman et al. (2020) propose guaranteed monotonicity for standard ReLU networks by letting a Satisfiability Modulo Theories (SMT) solver find counterexamples to the monotonicity definition and adjust the prediction in the inference process such that monotonicity is guaranteed. However, this approach requires queries to the SMT solver during inference time for each monotonic feature, and the computation time scales harshly with the number of monotonic features and the model size (see Figure 3 and 4 in Sivaraman et al. (2020)).

Our architecture falls into the first category. However, we overcome both main drawbacks: lack of expressiveness and impractical complexity. Other related works appear in the context of monotonic functions for normalizing flows, where monotonicity is a key ingredient to enforce invertibility (De Cao et al., 2020; Huang et al., 2018; Behrmann et al., 2019; Wehenkel & Louppe, 2019).

## 3 Methods

The goal is to develop a neural network architecture representing a vector-valued function

\[f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{n},\quad d,n\in\mathbb{N}, \tag{1}\]

that is provably monotonic in any subset of its inputs. We first define a few ingredients.

**Definition 3.1** (Monotonicity).: Let \(\mathbf{x}\in\mathbb{R}^{d},\mathbf{x}_{\mathbb{S}}\equiv\mathbbm{1}_{\mathbb{S}} \odot\mathbf{x}\), and the Hadamard product of \(\mathbf{x}\) with the indicator vector \(\mathbbm{1}_{\mathbb{S}}(i)=1\) if \(i\in\mathbb{S}\) and \(0\) otherwise for a subset \(\mathbb{S}\subseteq\{1,\cdots,d\}\).

We say that outputs \(\mathbb{Q}\subseteq\{1,\cdots,n\}\) of \(f\) are monotonically increasing in features \(\mathbb{S}\) if

\[f(\mathbf{x}_{\mathbb{S}}^{\prime}+\mathbf{x}_{\mathbb{S}})_{i}\leq f(\mathbf{x}_{ \mathbb{S}}+\mathbf{x}_{\mathbb{S}})_{i}\quad\forall i\in\mathbb{Q}\text{ and }\forall\mathbf{x}_{\mathbb{S}}^{\prime}\leq\mathbf{x}_{\mathbb{S}}, \tag{2}\]

where \(\mathbb{\bar{S}}\) denotes the complement of \(\mathbb{S}\) and the inequality on the right uses the product (or component-wise) order.
**Definition 3.2** (Lip\({}^{p}\) function).: \(g:\mathbb{R}^{d}\to\mathbb{R}^{n}\) is Lip\({}^{p}\) if it is Lipschitz continuous with respect to the \(L^{p}\) norm in every output dimension, _i.e._,

\[||g(\mathbf{x})-g(\mathbf{y})||_{\infty}\leq\lambda\|\mathbf{x}-\mathbf{y}\|_{p}\quad\forall\bm {x},\mathbf{y}\in\mathbb{R}^{n}\,. \tag{3}\]

### Lipschitz Monotonic Networks (LMN)

We will henceforth and without loss of generality only consider scalar-valued functions (\(n=1\)). We start with a model \(g(\mathbf{x})\) that is Lip\({}^{1}\) with Lipschitz constant \(\lambda\). Note that the choice of \(p=1\) is crucial for decoupling the magnitudes of the directional derivatives in the monotonic features. More details on this can be found below and in Figure 1. The \(1\)-norm has the convenient side effect that we can tune the robustness requirement for each input individually.

With a model \(g(\mathbf{x})\) we can define an architecture with built-in monotonicity by adding a term that has directional derivative \(\lambda\) for each coordinate in \(\mathbb{S}\):

\[f(\mathbf{x})=g(\mathbf{x})+\lambda(\mathbf{1}_{\mathbb{S}}\cdot\mathbf{x})=g(\mathbf{x})+\lambda \sum_{i\in\mathbb{S}}x_{i}. \tag{4}\]

This residual connection \(\lambda(\mathbf{1}_{\mathbb{S}}\cdot\mathbf{x})\) enforces monotonicity in the input subset \(\mathbf{x}_{\mathbb{S}}\):

\[\frac{\partial g}{\partial x_{i}}\in[-\lambda,\lambda],\;\;\forall i \in\mathbb{N}_{1:n} \tag{5}\] \[\Rightarrow\frac{\partial f}{\partial x_{i}}=\frac{\partial g}{ \partial x_{i}}+\lambda\geq 0\;\;\forall\,\mathbf{x}\in\mathbb{R}^{n},i\in \mathbb{S}\,. \tag{6}\]

The importance of the norm choiceThe construction presented here does not work with \(p\neq 1\) constraints because dependencies between the partial derivatives may be introduced, see Figure 1. The \(p=1\)-norm is the only norm that bounds the gradient within the green square and, crucially, allows the directional derivatives to be as large as \(2\lambda\) independently. When shifting the constraints by introducing the linear term, the green square allows for all possible gradient configurations, given that we can choose \(\lambda\) freely. As a counter example, the red circle, corresponding to \(p=2\) constraints, prohibits important areas in the configuration space.

To be able to represent all monotonic Lip\({}^{1}\) functions with \(2\lambda\) Lipschitz constant, the construction of \(g(\mathbf{x})\) needs to be a universal approximator of Lip\({}^{1}\) functions. In the next section, we will discuss possible architectures for this task.

### Lip\({}^{p=1}\) approximators

Our goal is to construct a universal approximator of Lip\({}^{1}\) functions, _i.e._, we would like the hypothesis class to have two properties:

Figure 1: \(p\)-norm constrained gradients showing (red) \(p=2\) and (green) \(p=1\). The gradient of a function \(g(\mathbf{x})\) that is Lip\({}^{p=2}\) resides within the dashed red line. For a Lip\({}^{p=1}\) function, the boundary is the green dashed line. Note that \(\mathbf{x}\) is taken to be a row vector. The residual connection (in blue) effectively shifts the possible gradients to strictly positive values and thus enforces monotonicity. Note how the red solid circle does not include all possible gradient configurations. For instance, it does not allow for very small gradients in both inputs, whereas the green square includes all configurations, up to an element-wise maximum of \(2\lambda\).
1. It always satisfies Eq. 3, _i.e._, be \(\mathrm{Lip}^{1}\).
2. It is able to fit all possible \(\mathrm{Lip}^{1}\) functions. In particular, the bound in Eq. 3 needs to be attainable \(\forall\,\mathbf{x},\mathbf{y}\).

Lip\({}^{1}\) constrained modelsTo satisfy the first requirement, fully connected networks can be Lipschitz bounded by constraining the matrix norm of all weight matrices (Kitouni et al., 2021; Gouk et al., 2020; Miyato et al., 2018). We recursively define the layer \(l\) of the fully connected network of depth \(D\) with activation \(\sigma\) as

\[\mathbf{z}^{l}=\sigma(\mathbf{z}^{l-1})\mathbf{W}^{l}+\mathbf{b}^{l}, \tag{7}\]

where \(\mathbf{z}^{0}=\mathbf{x}\) is the input and \(f(\mathbf{x})=z^{D}\) is the output of the neural network. It follows that \(g(\mathbf{x})\) satisfies Eq. 3 if

\[\prod_{i=1}^{D}\|\mathbf{W}^{i}\|_{1}\leq\lambda, \tag{8}\]

and \(\sigma\) has a Lipschitz constant less than or equal to 1. There are multiple ways to enforce Eq. 8. Two existing possibilities that involve scaling by the operator norm of the weight matrix (Gouk et al., 2020) are:

\[\mathbf{W}^{i}\rightarrow\mathbf{W}^{\prime i}=\lambda^{1/D}\frac{\mathbf{W}^{i}}{\text{ max}(1,\|\mathbf{W}^{i}\|_{1})}\qquad\text{or}\qquad W^{i}\to W^{\prime i}= \frac{\mathbf{W}^{i}}{\text{max}(1,\lambda^{-1/D}\cdot\|\mathbf{W}^{i}\|_{1})}\,. \tag{9}\]

In our studies, the latter variant seems to train slightly better. However, in some cases it might be useful to use the former to avoid the scale imbalance between the neural network's output and the residual connection used to induce monotonicity. We note that in order to satisfy Eq. 8, it is not necessary to divide the entire matrix by its 1-norm. It is sufficient to ensure that the absolute sum over each column is constrained:

\[\mathbf{W}^{i}\rightarrow\mathbf{W}^{\prime i}=\mathbf{W}^{i}\text{diag}\left(\frac{1}{ \text{max}\left(1,\lambda^{-1/D}\cdot\sum_{j}|W^{i}_{jk}|\right)}\right)\,. \tag{10}\]

This novel normalization scheme tends to give even better training results in practice, because the constraint is applied in each column individually. This reduces correlations of constraints, in particular, if a column saturates the bound on the norm, the other columns are not impacted. While Eq. 10 may not be suitable as a general-purpose scheme, _e.g._ it would not work in convolutional networks, its performance in training in our analysis motivates its use in fully connected architectures and further study of this approach in future work.

In addition, the constraints in Eq. 9 and Eq. 10 can be applied in different ways. For example, one could normalize the weights directly before each call such that the induced gradients are propagated through the network like in Miyato et al. (2018). While one could come up with toy examples for which propagating the gradients in this way hurts training, it appears that this approach is what usually is implemented for spectral norm in PyTorch and TensorFlow (Miyato et al., 2018). Alternatively, the constraint could be applied by projecting any infeasible parameter values back into the set of feasible matrices after each gradient update as in Algorithm 2 of Gouk et al. (2020).

Constraining according to Eq. 8 is not the only way to enforce \(\mathrm{Lip}^{1}\). Anil et al. (2019) provide an alternative normalization scheme:

\[\|\mathbf{W}^{1}\|_{1,\infty}\cdot\prod_{i=2}^{m}\|\mathbf{W}^{i}\|_{\infty}\leq\lambda \tag{11}\]

Similarly to how the 1-norm of a matrix is a column-wise maximum, the \(\infty\)-norm of a matrix is determined by the maximum 1-norm of all rows and \(\|W\|_{1,\infty}\) simply equals the maximum absolute value of an element in the matrix. Therefore, normalization schemes similar to Eq. 10, can be employed to enforce the constraints in Eq. 11 by replacing the column-wise normalization with a row- or element-wise normalization where appropriate.
Preserving expressive powerGuaranteeing that the model is Lipschitz bounded is not sufficient, it must also able to saturate the bound to be able to model all possible \(\text{Lip}^{1}\) functions. Some Lipschitz network architectures, _e.g._ Miyato et al. (2018), tend to over constrain the model such that it cannot fit all \(\text{Lip}^{1}\) functions due to _gradient attenuation_. For many problems this is a rather theoretical issue. However, it becomes a practical problem for the monotonic architecture since it often works on the edges of its constraints, for instance when partial derivatives close to zero are required, see Figure 1. As a simple example, the authors of Huster et al. (2018) showed that ReLU networks are unable to fit the function \(f(x)=|x|\) if the layers are norm-constrained with \(\lambda=1\). The reason lies in fact that ReLU, and most other commonly used activations, do not have unit gradient with respect to the inputs over their entire domain. While monotonic element-wise activations like ReLU cannot have unit gradient almost everywhere without being exactly linear, the authors of Anil et al. (2019) explore activations that introduce non-linearities by reordering elements of the input vector. They propose **GroupSort** as an alternative to point-wise activations, and it is defined as follows:

\[\sigma_{G}(\mathbf{x}) =\text{sort}_{1:G}(\mathbf{x}_{1:G})+\text{sort}_{G+1:2G}(\mathbf{x}_{G +1:2G})+\ldots\] \[=\sum_{i=0}^{n/G-1}\text{sort}_{iG+1:(i+1)G}(\mathbf{x}_{iG+1:(i+1)G}), \tag{12}\]

where \(\mathbf{x}\in\mathbb{R}^{n}\), \(\mathbf{x}_{i:j}=\mathbf{1}_{i:j}\odot\mathbf{x}\), and \(\text{sort}_{i:j}\) orders the elements of a vector from indices \(i\) to \(j\) and leaves the other elements in place. This activation sorts an input vector in chunks (groups) of a fixed size \(G\). The GroupSort operation has a gradient of unity with respect to every input, giving architectures constrained with Eq. 8 greatly increased expressive power. In fact, Anil et al. (2019) prove that GroupSort networks with the normalization scheme in Eq. 11 are universal approximators of \(\text{Lip}^{1}\) functions. Therefore, these networks fulfill the two requirements outlined in the beginning of this section.

For universal approximation to be possible, the activation function used needs to be gradient norm preserving (GNP), _i.e._, have gradient 1 almost everywhere. Householder activations are another instance of GNP activations of which GroupSort-2 is a special case (Singla et al., 2021). The Householder activation is defined as follows:

\[\sigma(\mathbf{z})=\begin{cases}\mathbf{z}&\mathbf{zv}>0\\ \mathbf{z}(\mathbf{1}-2\mathbf{v}\mathbf{v}^{T})&\mathbf{zv}\leq 0\end{cases} \tag{13}\]

Here, \(\mathbf{z}\) is the preactivation row vector and \(\mathbf{v}\) is any column unit vector. Householder Lipschitz Networks naturally inherit the universal approximation property.

In summary, we have constructed a neural network architecture \(f(\mathbf{x})\) via Eq. 4 that can provably approximate all monotonic Lipschitz bounded functions. The Lipschitz constant of the model can be increased arbitrarily by controlling the parameter \(\lambda\) in our construction.

## 4 Experiments

"Beware of bugs in the above code, I have only proved it correct, not tried it" (Knuth). In the spirit of Donald Knuth, in this section we test our algorithm on many different domains to show that it works well in practice and gives competitive results, as should be expected from a universal approximator.

### Toy Example

Figure 2 shows a toy example where both a monotonic and an unconstrained network are trained to regress on a noisy one-dimensional dataset. The true underlying model used here is monotonic, though an added heteroskedastic Gaussian noise term can obscure this in any realization. As can be seen in Figure 2, no matter how the data are distributed at the edge of the support, the monotonic Lipschitz network is always non-decreasing outside of the support as guaranteed by our architecture. Such out-of-distribution guarantees can be extremely valuable in cases where domain knowledge dictates monotonic behavior is either required or desirable.

### Real-Time Decision-Making at 40 MHz at the LHC

Because many physical systems are modeled with well-known theoretical frameworks that dictate the properties of the system, monotonicity can be a crucial inductive bias in the physical sciences.
For instance, modeling enthalpy, a thermodynamic quantity measuring the total heat content of a system, in a simulator requires a monotonic function of temperature for fixed pressure (as is known from basic physical principles).

In this section, we describe a real-world physics application which requires monotonicity in certain features--and robustness in all of them. The algorithm described here has, in fact, been implemented by a high-energy particle physics experiment at the European Center for Nuclear Research (CERN), and is actively being used to collect data at the Large Hadron Collider (LHC) in 2022, where high-energy proton-proton collisions occur at 40 MHz.

The sensor arrays of the LHC experiments produce data at a rate of over 100 TB/s. Drastic data-reduction is performed by custom-built read-out electronics; however, the annual data volumes are still \(O(100)\) exabytes, which cannot be put into permanent storage. Therefore, each LHC experiment processes its data in real time, deciding which proton-proton collision events should be kept and which should be discarded permanently; this is referred to as _triggering_ in particle physics. To be suitable for use in trigger systems, classification algorithms must be robust against the impact of experimental instabilities that occur during data taking--and deficiencies in simulated training samples. Our training samples cannot possibly account for the unknown new physics that we hope to learn by performing the experiments!

A ubiquitous inductive bias at the LHC is that outlier collision events are more interesting, since we are looking for physics that has never been observed before. However, uninteresting outliers are frequently caused by experimental imperfections, many of which are included and labeled as background in training. Conversely, it is not possible to include the set of all possible interesting outliers _a priori_ in the training. A solution to this problem is to implement _outliers are better_ directly using our expressive monotonic Lipschitz architecture from Section 3.

Our architecture was originally developed for the task of classifying the decays of heavy-flavor particles produced at the LHC. These are bound states containing a beauty or charm quark that travel an observable distance \(\mathcal{O}(1\,\mathrm{cm})\) before decaying due to their (relatively) long lifetimes. This example uses a dataset of simulated proton-proton (\(pp\)) collisions in the LHCb detector. Charged particles recorded by LHCb are combined pairwise into decay-vertex (DV) candidates. The task concerns discriminating DV candidates corresponding to heavy-flavor decays from all other sources. Heavy-flavor DVs typically have substantial separation from the \(pp\) collision point, due to the relatively long heavy-flavor particle lifetimes, and large transverse momenta, \(p_{\mathrm{T}}\), of the component particles, due to the large heavy-flavor particle masses. The main sources of background DVs, described in Kitouni et al. (2021), mostly have small displacement and small \(p_{\mathrm{T}}\), though unfortunately they can also have extremely large values of both displacement and momentum.

Figure 2: Our monotonic architecture (green) and an unconstrained network (red) trained on two realizations (purple data points) of a one dimensional dataset. The shaded regions are where training data were absent. Each model is trained using 10 random initialization seeds. The dark lines are averages over the seeds, which are each shown as light lines. The unconstrained models exhibit overfitting of the noise, non-monotonic behavior, and highly undesirable and unpredictable results when extrapolating beyond the region occupied by the training data. Conversely, the monotonic Lipschitz models are always monotonic, even in scenarios where the noise is strongly suggestive of non-monotonic behavior. In addition, the Lipschitz constraint produces much smoother models.
Figure 3 shows a simplified version of this problem using only the two most-powerful inputs. Our inductive bias requires a monotonic increasing response in both features (detailed discussion motivating this bias can be found in Kitouni et al. (2021)). We see that an unconstrained neural network rejects DVs with increasing larger displacements (lower right corner), and that this leads to a decrease of the signal efficiency (true positive rate) for large lifetimes. The unconstrained model violates our inductive bias. Figures 3 and 4 show that a monotonic BDT (Auguste et al., 2020) approach works here. However, the jagged decision boundary can cause problems in subsequent analysis of the data. Figure 3 also shows that our novel approach from Section 3 successfully produces a smooth and monotonic response, and Figure 4 shows that this provides the monotonic lifetime dependence we desire in the efficiency.

In addition, we note that the added benefit of guaranteed Lipschitz robustness is a major advantage for many real world applications. Specifically for particle physicists, this kind of robustness directly translates to important guarantees when considering experimental instabilities.

Due to the simplicity and practicality of our method, the LHCb experiment is now using the proposed architecture for real-time data selection at a data rate of about \(40\,\)Tbit/s.

Figure 4: From Kitouni et al. (2021): True positive rate (efficiency) of each model shown in Figure 3 versus the proper lifetime of the decaying heavy-quark particle selected. The monotonic models produce a nearly uniform efficiency above a few picoseconds at the expense of a few percent lifetime-integrated efficiency. Such a trade off is desirable as explained in the text.

Figure 3: From Kitouni et al. (2021): Simplified version of the heavy-quark selection problem using only two inputs, which permits displaying the response everywhere in the feature space; shown here as a heat map with more signal-like (background-like) regions colored blue (red). The dark solid line shows the decision boundary (upper right regions are selected). Shown are (left) a standard fully connected neural network, (middle) a monotonic BDT, and (right) our architecture. The quantities shown on the horizontal and vertical axes are related to how long the particle lived before decaying and how massive the particle was, respectively.
### Public datasets with monotonic dependence

In this section, we follow as closely as possible the experiments done in Liu et al. (2020), and some experiments done in Sivaraman et al. (2020) to be able to directly compare to state-of-the-art monotonic architectures. Liu et al. (2020) studied monotonic architectures on four different datasets: COMPAS (Larson and Kirchner, 2016), BlogFeedback (Buza, 2014), LoanDefaulter (Kaggle, 2015), and ChestXRay (Wang et al., 2017). From Sivaraman et al. (2020) we compare against one regression and one classification task: AutoMPG (Dua and Graff, 2017) and HeartDisease (Gennari et al., 1989). Results are shown in Table 1.
an augmented version (i.e. with an additional monotonic feature added artificially) of CIFAR100 in Appendix A.

## 5 Limitations

We are working on improving the architecture as follows: First, common initialization techniques are not optimal for weight-normed networks (Arpit et al., 2019). Simple modifications to the weight initialization might aid convergence, especially for large Lipschitz parameters. Secondly, we are currently constrained to activation functions that have a gradient norm of \(1\) over their entire domain, such as **GroupSort**, to ensure universal approximation, see Anil et al. (2019). We will explore other options in the future. Lastly, there is not yet a proof for universal approximation for the architecture described in Eq. 8. However, it appears from empirical investigation that the networks do approximate universally, as we have yet to find a function that could not be approximated well enough with a deep enough network. We do not consider this a major drawback, as the construction in Eq. 11 does approximate universally, see Anil et al. (2019). Note that none of these limitations have any visible impact on the performance of the experiments in Section 4.

## 6 Conclusion and Future Work

We presented an architecture that provably approximates Lipschitz continuous and partially monotonic functions. Monotonic dependence is enforced via an end-to-end residual connection to a minimally Lip\({}^{1}\) constrained fully connected neural network. This method is simple to implement, has negligible computational overhead, and gives stronger guarantees than regularized models. Our architecture achieves competitive results with respect to current state-of-the-art monotonic architectures, even when using a tiny number of parameters, and has the additional benefit of guaranteed robustness due to its known Lipschitz constant. For future directions of this line of research, we plan to tackle the problems outlined in the limitation section, especially improving initialization of weight-normed networks.

\begin{table}
\begin{tabular}{c|c c} \multicolumn{3}{c}{**COMPAS**} \\ \hline Method & Parameters & \(\left\uparrow\right.\) Test Acc \\ \hline Certified & 23112 & \((68.8\pm 0.2)\%\) \\
**LMN** & **37** & \((\mathbf{69.3\pm 0.1})\%\) \\ \multicolumn{3}{c}{**LoanDefaulter**} \\ \hline Method & Parameters & \(\left\uparrow\right.\) Test Acc \\ \hline Certified & 8502 & \((65.2\pm 0.1)\%\) \\
**LMN** & **753** & \((\mathbf{65.44\pm 0.03})\%\) \\
**LMN mini** & **69** & \((\mathbf{65.28\pm 0.01})\%\) \\ \multicolumn{3}{c}{**Heart Disease**} \\ \hline Method & \(\left\uparrow\right.\) Test Acc \\ \hline COMET & \((86\pm 3)\%\) \\
**LMN** & \((89.6\pm 1.9)\%\) \\ \end{tabular} 
\begin{tabular}{c|c c} \multicolumn{3}{c}{**BlogFeedback**} \\ \hline Method & Parameters & \(\left\downarrow\right.\) RMSE \\ \hline Certified & 8492 & \(.158\pm.001\) \\
**LMN** & **2225** & \(\mathbf{.160\pm.001}\) \\
**LMN mini** & **177** & \(\mathbf{.155\pm.001}\) \\ \multicolumn{3}{c}{**ChestXRay**} \\ \hline Method & Parameters & \(\left\uparrow\right.\) Test Acc \\ \hline Certified & 12792 & \((62.3\pm 0.2)\%\) \\ Certified E-E & 12792 & \((66.3\pm 1.0)\%\) \\
**LMN** & **1043** & \((\mathbf{67.6\pm 0.6})\%\) \\
**LMN E-E** & **1043** & \((\mathbf{70.0\pm 1.4})\%\) \\ \multicolumn{3}{c}{**Auto MPG**} \\ \hline Method & \(\left\downarrow\right.\) MSE \\ \hline COMET & \((8.81\pm 1.81)\%\) \\
**LMN** & \((\mathbf{7.58\pm 1.2})\%\) \\ \end{tabular} 
\begin{tabular}{c|c c} \multicolumn{3}{c}{**Auto MPG**} \\ \hline Method & \(\left\downarrow\right.\) MSE \\ \hline COMET & \((8.81\pm 1.81)\%\) \\
**LMN** & \((\mathbf{7.58\pm 1.2})\%\) \\ \end{tabular}
\end{table}
Table 1: We compare our method (in bold) against state-of-the-art monotonic models (we only show the best) on a variety of benchmarks. The performance numbers for other techniques were taken from Liu et al. (2020) and Sivaraman et al. (2020). In the ChestXRay experiment, we train one model with frozen ResNet18 weights (second to last) and another with end-to-end training (last). While our models can generally get quite small, we can achieve even smaller models when only taking a subset of all the features. These models are denoted with “mini”.
## 7 Reproducibility Statement

All experiments with public datasets are reproducible with the code provided at [https://github.com/niklasnolte/monotonic_tests](https://github.com/niklasnolte/monotonic_tests). This code uses the package available in [https://github.com/niklasnolte/MonotoneNorm](https://github.com/niklasnolte/MonotoneNorm), which is meant to be a standalone pytorch implementation of Lipschitz Monotonic Networks. The experiments in Section 4.2 were made with data that is not publicly available. The code to reproduce those experiments can be found under [https://github.com/niklasnolte/HLT_2Track](https://github.com/niklasnolte/HLT_2Track) and the data will be made available in later years at the discretion of the LHCb collaboration.

#### Acknowledgments

This work was supported by NSF grant PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, [http://iaifi.org/](http://iaifi.org/)).
