* Patel and Menon [2008] N. Patel, S. Menon,  Simulation of spray–turbulence–flame interactions in a lean direct injection combustor,  Combustion and Flame 153 (2008) 228–257.
  * Balakrishnan et al. [2010] K. Balakrishnan, D. Nance, S. Menon,  Simulation of impulse effects from explosive charges containing metal particles,  Shock Waves 20 (2010) 217–239.
  * Toro [2009] E. F. Toro,  The hll and hllc riemann solvers,  in: Riemann solvers and numerical methods for fluid dynamics, Springer, 2009, pp. 315–344.
  * Schmidmayer et al. [2021] K. Schmidmayer, J. Cazé, F. Petitpas, E. Daniel, N. Favrie, Modelling interactions between waves and diffused interfaces, 2021. URL: https://hal.archives-ouvertes.fr/hal-03387818, working paper or preprint.
  * Shyue and Xiao [2014] K.-M. Shyue, F. Xiao,  An eulerian interface sharpening algorithm for compressible two-phase flow: the algebraic thinc approach,  Journal of Computational Physics 268 (2014) 326–354.
  * Gryngarten and Menon [2013] L. D. Gryngarten, S. Menon,  A generalized approach for sub-and super-critical flows using the local discontinuous galerkin method,  Computer Methods in Applied Mechanics and Engineering 253 (2013) 169–185.
  * Jain et al. [2021] S. S. Jain, M. C. Adler, J. R. West, A. Mani, P. Moin, S. K. Lele,  Assessment of diffuse-interface methods for compressible multiphase fluid flows and elastic-plastic deformation in solids,  arXiv preprint arXiv:2109.09729 (2021).
  * Poinsot and Lele [1992] T. J. Poinsot, S. K. Lele,  Boundary conditions for direct simulations of compressible viscous flows,  Journal of Computational Physics 101 (1992) 104–129.
  * Sridharan et al. [2015] P. Sridharan, T. L. Jackson, J. Zhang, S. Balachandar,  Shock interaction with one-dimensional array of particles in air,  Journal of Applied Physics 117 (2015) 075902.
  * Ghia et al. [1982] U. Ghia, K. Ghia, C. Shin,  High-Re solutions for incompressible flow using the navier-stokes equations and a multigrid method,  Journal of Computational Physics 48 (1982) 387–411.
  * White and Corfield [2006] F. M. White, I. Corfield, Viscous fluid flow, volume 3, McGraw-Hill New York, 2006.
  * Igra et al. [2002] D. Igra, T. Ogawa, K. Takayama,  A parametric study of water column deformation resulting from shock wave loading,  Atomization and Sprays 12 (2002).
  * Chen [2008] H. Chen,  Two-dimensional simulation of stripping breakup of a water droplet,  AIAA Journal 46 (2008) 1135–1143.
  * Rayleigh et al. [1879] L. Rayleigh, et al.,  On the capillary phenomena of jets,  Proc. R. Soc. London 29 (1879) 71–97.
  * Gallagher et al. [2017] T. Gallagher, M. Akiki, S. Menon, V. Sankaran, V. Sankaran,  Development of the generalized maccormack scheme and its extension to low mach number flows,  International Journal Numerical Methods in Fluids 85 (2017) 165–188.
  * Patel [2007] N. V. Patel, Simulation of hydrodynamic fragmentation from a fundamental and an engineering perspective, Ph.D. thesis, Georgia Institute of Technology, 2007.
  * Temkin and Kim [1980] S. Temkin, S. S. Kim,  Droplet motion induced by weak shock waves,  Journal of Fluid Mechanics 96 (1980) 133–157.
  * Pilch and Erdman [1987] M. Pilch, C. Erdman,  Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop,  International Journal of Multiphase Flow 13 (1987) 741–757.
  * Salvadori et al. [2022] M. Salvadori, A. Panchal, D. Ranjan, S. Menon,  Numerical study of detonation propagation in H2-Air with kerosene droplets,  in: AIAA SciTech 2022 Forum, 2022, pp. AIAA–2022–0394.
  * Kindracki [2015] J. Kindracki,  Experimental research on rotating detonation in liquid fuel–gaseous air mixtures,  Aerospace Science and Technology 43 (2015) 445–453.
  * Gogulya et al. [2004] M. Gogulya, M. Makhov, A. Y. Dolgoborodov, M. Brazhnikov, V. Arkhipov, V. Shchetinin,  Mechanical sensitivity and detonation parameters of aluminized explosives,  Combustion, Explosion, and Shock Waves 40 (2004) 445–457.
  * Kailasanath et al. [1985] K. Kailasanath, E. Oran, J. Boris, T. Young,  Determination of detonation cell size and the role of transverse waves in two-dimensional detonations,  Combustion and Flame 61 (1985) 199–209.
  * Pelanti [2021] M. Pelanti,  Arbitrary-rate relaxation techniques for the numerical modeling of compressible two-phase flows with heat and mass transfer,  arXiv preprint arXiv:2108.00556 (2021).


# MARIO Eval: Evaluate Your Math LLM with your Math LLM
–A mathematical dataset evaluation toolkit

Boning Zhang , Chengxi Li∗, Kai Fan∗
Zhejiang University, Alibaba Group
<EMAIL_ADDRESS><EMAIL_ADDRESS>
Equal contribution. This work was done when the first author’s Internship at
Alibaba.Corresponding Author

###### Abstract

Large language models (LLMs) have been explored in a variety of reasoning
tasks including solving of mathematical problems. Each math dataset typically
includes its own specially designed evaluation script, which, while suitable
for its intended use, lacks generalizability across different datasets.
Consequently, updates and adaptations to these evaluation tools tend to occur
without being systematically reported, leading to inconsistencies and
obstacles to fair comparison across studies. To bridge this gap, we introduce
a comprehensive mathematical evaluation toolkit that not only utilizes a
python computer algebra system (CAS) for its numerical accuracy, but also
integrates an optional LLM, known for its considerable natural language
processing capabilities. To validate the effectiveness of our toolkit, we
manually annotated two distinct datasets. Our experiments demonstrate that the
toolkit yields more robust evaluation results compared to prior works, even
without an LLM. Furthermore, when an LLM is incorporated, there is a notable
enhancement. The code for our method will be made available at
https://github.com/MARIO-Math-Reasoning/math_evaluation.

##  1 Introduction

With appropriate prompts or external tools, Large language models (LLMs) have
attained human parity in various tasks. Nonetheless, mathematical reasoning
remains a formidable challenge for LLMs, necessitating a systematic evaluation
to accurately compare their performance on such tasks. However, this field
suffers from poor automatic evaluation that is neither robust nor complete
enough for accurate evaluation on math problem answers, owing to their
complexity and diversity. Recent LLMs reasoning methods (ToRA [4] and MathChat
[12]) lack unified math evaluation standards even on the same dataset such as
MATH [4]. Such discrepancy between different evaluation scripts may not
accurately reflect their true reasoning capability. Consequently, we are
motivated to design a more comprehensive automatic evaluation toolkit that
encompasses various math concepts within answers. Our aim is to establish a
convenient and standardized evaluation framework to support future research in
mathematical reasoning.

Upon conducting a thorough review of existing automatic math evaluation
methods, we pinpointed several key shortcomings. Traditionally, the assessment
of mathematical answers has heavily relied on simplistic methods such as
direct string comparisons or simple rules, inadequate to address complex
situations. As illustrated in Figure 1, identical answer expression may imply
different math concepts, while different expressions may be actually
equivalent under certain conditions.

Figure 1: Most previous evaluation tools judge correctness solely based on the answer, while ours also takes into account the answer type of the question implied. Type | Code Friendly Definition
---|---
Real | A number can be used to measure a continuous one-dimensional quantity.
Complex | A number can be written in the form of $a+bi$, $a,b\in\mathcal{R}$ and $i$ is the imaginary unit.
Set | A collection of different elements which are typically mathematical objects of any kind.
Interval(s) | A interval is a set of real numbers that lie between two fixed endpoints without any gaps.
| Intervals represent the intersection or union of several interval.
Vector | A collection of ordered elements which are typically mathematical objects of any kind.
Matrix | A rectangular array of numbers, symbols, or expressions, arranged in rows and columns.
Expression | A combination of symbols that is well-formed according to rules that depend on the
| context, including at least one unknown variable, but no equal sign, _e.g._
, $x^{2}+y^{2}-1$
Function | A function from a set $X$ to a set $Y$ is an assignment of an element of $Y$ to each element
| of $X$, _e.g._ , $f(x,y)=x^{2}+y^{2}-1$.
Equation | A equation states that two quantities are the same, in the form of $A=B$.
| In our case, at least one of $A$ and $B$ should be Expression, _e.g._ ,
$x^{2}+y^{2}=1$.
Inequality | A relation which makes a non-equal comparison between two numbers or expressions,
|  _e.g._ , $x^{2}+y^{2}<1$.
Others | Do not belong to above 10 types.
Table 1: The type definitions lack complete rigor by design, with the
intention to align these types with those defined in Python or SymPy as Table
2.

In light of such limitations, we propose a novel two-stage math evaluation
toolkit that features an optional integration with LLM. The initial stage
involves identifying the type of answer based on our predefined set of
mathematical concepts. Subsequently, the second stage invokes the relevant
type-specific function to evaluate the equivalence between the expected and
predicted answers. When integrating the LLM, it has the potential to
significantly enhance the precision of both the classification of answer types
and the determination of answer equivalence. We evaluate the equivalence of
answers using our toolkit in various configurations: without LLM, with LLM,
and LLM-only. This evaluation reveals that the hybrid approach (with LLM)
effectively leverages the numerical precision offered by Python packages and
the natural language comprehension capabilities of the LLM. In summary, our
key contributions are as follows:

  1. 1.

We construct a mathematical evaluation toolkit with an optionally integratable
LLM. It covers a variety of mathematical reasoning answer types.

  2. 2.

We curated three mathematical answer evaluation datasets, including original
questions and answers, model predicted answers, and our manually annotated
answer types and equivalences.

  3. 3.

We plan to make our datasets and evaluation toolkit publicly available to
support and advance the efforts of the research community.

##  2 Main Framework

In this section, we outline the design pattern of our mathematical evaluation
toolkit. It mainly includes two modules. One is a type classification module
that is utilized to ascertain the type of the expected answer. The other is an
evaluation module that is crafted to assess the equivalence between the
predicted response and the expected answer.

###  2.1 Type Definitions

Drawing the inspiration from the way mathematical theories are formulated by
using a set of axioms, we establish our types ranging from fundamental
concepts (_e.g._ , Real) to complex composite structures (_e.g._ , Matrix).
The type definitions are delineated in Table 1. It is important to note that
these definitions are intentionally not fully rigorous; our primary aim is to
ensure consistency with the types inherent within the standard Python
language, as well as with the Python computer algebra system (CAS) package,
SymPy. For instance, the type ’Set’ in our framework specifically corresponds
to the set class in Python, represented with curly braces (“{…}"), which is
distinct from the concept of a mathematical set that also contains intervals.
Table 2 provides a comprehensive comparison between the types defined in our
framework and the built-in types in Python and SymPy. For Real and Complex, we
define the customized types, because SymPy considers real number as complex as
well.

Type | Package | Aligned type
---|---|---
Set | Python | set
Vector | Python | list
Matrix | SymPy | Matrix
Interval(s) | SymPy | Interval
Expression | SymPy | Expr
Function | SymPy | Function
Equation | SymPy | Equality
Inequality | SymPy | Relational
Table 2: Our basic design attempts to align defined types in Python language.
Algorithm 1 Design of our math evaluation

question $\mathbf{q}$, answer $\mathbf{a}$, prediction $\mathbf{p}$

1:$\mathbf{a}_{\text{type}}=\texttt{rule\\_classifier}(\mathbf{a},\mathbf{p})$

2:try

3:
$\texttt{ans}=\texttt{is\\_equiv}(\mathbf{a},\mathbf{p},\mathbf{a}_{\text{type}})$
$\triangleright$ False or Error may be re-evaluated later

4: try

5:except

6: $\texttt{ans}=\texttt{False}$

7: exceptnot ans and LLM is not None

8:$\mathbf{a}_{\text{type}}=\text{LLM}_{\text{type\\_prompt}}(\mathbf{q},\mathbf{a})$

9:try

10:
$\texttt{ans}=\texttt{is\\_equiv}(\mathbf{a},\mathbf{p},\mathbf{a}_{\text{type}})$
$\triangleright$ Only Error will be re-evaluated

11: try

12:except

13:
$\texttt{ans}=\text{LLM}_{\text{equiv\\_prompt}}(\mathbf{q},\mathbf{a},\mathbf{p})$

14: exceptans

###  2.2 Design Pattern

Based on the definitions, we initially create a rule-based type classifier
that depends solely on the expected and the predicted answers. For example, an
answer string that starts with "\beginmatrix" is classified as a Matrix type.
Subsequently, for the "Real" and "Expression" types, we develop two
fundamental equivalence functions. For all other types, their respective
equivalence functions can recursively invoke these two basic functions
according to the internal data structures.

However, as Figure 1 illustrates, relying solely on the answer string to
determine its type is not reliable, since the same string can represent vastly
different types. Inspired by the remarkable natural language understanding
capabilities of LLMs, we propose incorporating LLMs into the math evaluation
process to eliminate the confusions highlighted in the introduction. With
proper prompt, LLMs can analyze both the question and the answer to discern
the intended answer type. Furthermore, they can even directly assess the
equivalence between the answer and the prediction. To counteract the issue of
hallucination, LLMs will be employed exclusively for cases where the rule-
based method fails. The overall design is depicted in Algorithm 1.

###  2.3 Datasets and Setups

We mainly conducted experiments on two datasets, the MATH testset [5] and
GaoKao2023-Math-En (GK2023) [6], which contain 5,000 and 385 high school level
math problems, respectively. We performed supervised fine-tuning on the
DeepSeek-Math [13] base model on MATH trainset and infer on the two testsets
to extract the model-predicted answers. In addition, we downloaded the math
problem solving LLM ToRA-70B [4]. We then directly applied it to the
GaoKao2023 dataset for inference, denoting this specific test as GK2023-ToRA.
Subsequently, we manually verified the correctness of the predicted answers
for all three datasets.

For toolkit comparison, we utilized the official repositories from MATH, ToRA,
and DeepSeek-Math. We assess our evaluation toolkit from three perspectives:
equivalence accuracy, type classification accuracy, and solution accuracy.

equiv_acc
$\displaystyle=\frac{\\#\\{\text{human\\_eval}=\text{toolkit\\_eval}\\}}{\\#\text{dataset}}$
type_acc
$\displaystyle=\frac{\\#\\{\text{human\\_type}=\text{llm\\_type}\\}}{\\#\text{Dataset}}$
sol_acc
$\displaystyle=\frac{\\#\\{\text{toolkit\\_eval}=\texttt{True}\\}}{\\#\text{Dataset}}$

###  2.4 Main Results

The equivalence accuracy results are shown in Table 3. There is a significant
difference in the performance of the three open-source mathematical evaluation
toolkits when assessing the same pairs of expected and predicted answers.
Remarkably, gpt-3.5-turbo [8] appears to surpass these toolkits in terms of
making correctness judgments. In contrast, our toolkit in the configuration of
the basic design attains superior accuracy compared to gpt-3.5-turbo, and its
performance is further enhanced when integrated with it.

On the two testsets, our basic design can achieve about 97% accuracy, and the
incorporation of LLM technology provides an additional improvement of about
1%. While this improvement may appear minor, significant accuracy enhancements
are observed when integrating prior toolkits with our method outlined in
Algorithm 1. These results demonstrate that LLMs can significantly amplify the
effectiveness of existing tools. On GaoKao2023-ToRA, most of the inferred
results are incorrect, _i.e._ , the expected and predicted answers are quite
difference. Despite this, all toolkits exhibited satisfactory performance,
indicating that integrating the LLM may not be necessary. However, we can
still observe that our basic design still achieves better than ToRA toolkit,
which was specifically tailored for its output.

Toolkit | Equivalence Accuracy
---|---
MATH | GK2023 | GK2023-ToRA
LLM only | 95.03% | 93.77% | -
MATH | 92.51% | 91.39% | 92.21%
+Algorithm 1 | 95.15% | 93.18% | -
ToRA | 86.08% | 87.24% | 95.06%
+Algorithm 1 | 96.69% | 94.07% | -
DeepSeek-Math | 86.46% | 88.43% | 92.47%
+Algorithm 1 | 96.63% | 94.07% | -
basic design | 97.23% | 96.74% | 98.96%
\+ LLM type | 97.88% | 97.33% | -
\+ LLM equiv | 98.59% | 97.03% | -
Table 3: Equivalence accuracy. LLM: gpt-3.5-turbo. Numbers in red depend on LLM. Model | Accuracy on MATH testset
---|---
Type | Equivalence
gpt-3.5-turbo | 93.74% | 98.59%
Qwen-Max | 93.30% | 97.65%
Table 4: Type accuracy with different LLMs

###  2.5 Ablation Studies

Type Accuracy In the first study as illustrated in Table 4, we specially
evaluate two commercial LLMs: gpt3.5-turbo and Qwen-Max [1]. We choose them
for two primary reasons. First, their cost-effective API pricing enables the
handling of 5K problems. Because for a comprehensive assessment of the LLMs in
type classification, we utilized the entire dataset with the LLMs instead of
the procedure defined in Algorithm 1, activating LLM type prediction as
needed. Second, LLMs with larger size typically exhibit enhanced abilities in
natural language understanding and instruction-following, both of which are
crucial for accurate type determination. Meanwhile, the accuracy of
equivalence reported in Table 4 reflects the performance adhering strictly to
the procedure outlined in Algorithm 1. In general, we can conclude that better
LLMs may bring larger improvement over our basic design.

Solution Accuracy In the second study, depicted in Figure 2, we can observe
that previous evaluation toolkits exhibit a significant discrepancy (_e.g._ ,
over 10% gap in the MATH dataset) from the true accuracy, thereby making it
less reliable to access the actual mathematical reasoning capabilities of math
LLMs. Someone may argue that the metric of solution accuracy may not be
entirely accurate, as erroneous judgments over the correct and incorrect
answers can potentially cancel each other out. However, we can still observe a
strong correlation between equivalence accuracy and solution accuracy.
Nonetheless, a strong correlation exists between equivalence accuracy and
solution accuracy. For instance, in the MATH dataset and with our toolkits,
higher equivalence accuracy ($>90\%$) correlates with more robust solution
accuracy evaluations (less difference to the true accuracy).

$0.4$$0.45$$0.5$$0.55$$0.6$$0.65$Solution AccTrue AccuracyMATHToRADeepSeek-
MathOurs Figure 2: Solution accuracy with different toolkits. Left: MATH.
Right: GK2023

Case Studies We carried out three case studies using our manually annotated
test sets in the Appendix 6,7,8.

##  3 Related Works

In order to explore the reasoning ability of LLMs and their efficiency in
understanding special symbols, researchers utilize annotated data fine-tuning,
notably via ‘chain-of-thought’ data approaches [2, 10] and integration of
external knowledge sources [11]. Automated proof assistants exemplified by
GPT-f [9] further underscores this exploration. English math word problem
(MWP) corpus, such as ASDiv [7], GSM8K [3] and MATH [5] are proposed to
evaluate the math reasoning capabilities of LLMs in the past few years.
However, as the most widely used dataset, reasoning experiments on MATH lack
unified automatic evaluation methods. Our work releases a normalized math
evaluation work which reaches accuracy above $95\%$ on MATH, further reflects
the true capability of LLMs reasoning work more faithfully.

##  4 Conclusion

We presented a comprehensive mathematical evaluation toolkit that is able to
optionally incorporate LLMs. This toolkit is designed to take both advantages
of the CAS and LLM, and the experiments demonstrate our toolkit can promote
the accuracy and consistency of math evaluation and facilitating fair cross-
comparison of different studies.

## References

  * Bai et al. [2023] Jinze Bai, Shuai Bai, Yunfei Chu, Zeyu Cui, Kai Dang, Xiaodong Deng, Yang Fan, Wenbin Ge, Yu Han, Fei Huang, Binyuan Hui, Luo Ji, Mei Li, Junyang Lin, Runji Lin, Dayiheng Liu, Gao Liu, Chengqiang Lu, Keming Lu, Jianxin Ma, Rui Men, Xingzhang Ren, Xuancheng Ren, Chuanqi Tan, Sinan Tan, Jianhong Tu, Peng Wang, Shijie Wang, Wei Wang, Shengguang Wu, Benfeng Xu, Jin Xu, An Yang, Hao Yang, Jian Yang, Shusheng Yang, Yang Yao, Bowen Yu, Hongyi Yuan, Zheng Yuan, Jianwei Zhang, Xingxuan Zhang, Yichang Zhang, Zhenru Zhang, Chang Zhou, Jingren Zhou, Xiaohuan Zhou, and Tianhang Zhu.  Qwen technical report, 2023.
  * Chung et al. [2022] Hyung Won Chung, Le Hou, Shayne Longpre, Barret Zoph, Yi Tay, William Fedus, Yunxuan Li, Xuezhi Wang, Mostafa Dehghani, Siddhartha Brahma, Albert Webson, Shixiang Shane Gu, Zhuyun Dai, Mirac Suzgun, Xinyun Chen, Aakanksha Chowdhery, Alex Castro-Ros, Marie Pellat, Kevin Robinson, Dasha Valter, Sharan Narang, Gaurav Mishra, Adams Yu, Vincent Zhao, Yanping Huang, Andrew Dai, Hongkun Yu, Slav Petrov, Ed H. Chi, Jeff Dean, Jacob Devlin, Adam Roberts, Denny Zhou, Quoc V. Le, and Jason Wei.  Scaling instruction-finetuned language models, 2022.
  * Cobbe et al. [2021] Karl Cobbe, Vineet Kosaraju, Mohammad Bavarian, Mark Chen, Heewoo Jun, Lukasz Kaiser, Matthias Plappert, Jerry Tworek, Jacob Hilton, Reiichiro Nakano, Christopher Hesse, and John Schulman.  Training verifiers to solve math word problems, 2021.
  * Gou et al. [2024] Zhibin Gou, Zhihong Shao, Yeyun Gong, yelong shen, Yujiu Yang, Minlie Huang, Nan Duan, and Weizhu Chen.  ToRA: A tool-integrated reasoning agent for mathematical problem solving.  In _The Twelfth International Conference on Learning Representations_ , 2024.  URL https://openreview.net/forum?id=Ep0TtjVoap.
  * Hendrycks et al. [2021] Dan Hendrycks, Collin Burns, Saurav Kadavath, Akul Arora, Steven Basart, Eric Tang, Dawn Song, and Jacob Steinhardt.  Measuring mathematical problem solving with the math dataset, 2021.
  * Liao et al. [2024] Minpeng Liao, Wei Luo, Chengxi Li, Jing Wu, and Kai Fan.  Mario: Math reasoning with code interpreter output – a reproducible pipeline, 2024.
  * Miao et al. [2020] Shen-Yun Miao, Chao-Chun Liang, and Keh-Yih Su.  A diverse corpus for evaluating and developing english math word problem solvers.  In _Annual Meeting of the Association for Computational Linguistics_ , 2020.  URL https://api.semanticscholar.org/CorpusID:220047831.
  * OpenAI [2023] R OpenAI.  Gpt-4 technical report.  _arXiv_ , pages 2303–08774, 2023.
  * Polu and Sutskever [2020] Stanislas Polu and Ilya Sutskever.  Generative language modeling for automated theorem proving.  _ArXiv_ , abs/2009.03393, 2020.  URL https://api.semanticscholar.org/CorpusID:221535103.
  * Szegedy [2020] Christian Szegedy.  A promising path towards autoformalization and general artificial intelligence.  In _International Conference on Intelligent Computer Mathematics_ , 2020.  URL https://api.semanticscholar.org/CorpusID:220729524.
  * Thoppilan et al. [2022] Romal Thoppilan, Daniel De Freitas, Jamie Hall, Noam Shazeer, Apoorv Kulshreshtha, Heng-Tze Cheng, Alicia Jin, Taylor Bos, Leslie Baker, Yu Du, YaGuang Li, Hongrae Lee, Huaixiu Steven Zheng, Amin Ghafouri, Marcelo Menegali, Yanping Huang, Maxim Krikun, Dmitry Lepikhin, James Qin, Dehao Chen, Yuanzhong Xu, Zhifeng Chen, Adam Roberts, Maarten Bosma, Vincent Zhao, Yanqi Zhou, Chung-Ching Chang, Igor Krivokon, Will Rusch, Marc Pickett, Pranesh Srinivasan, Laichee Man, Kathleen Meier-Hellstern, Meredith Ringel Morris, Tulsee Doshi, Renelito Delos Santos, Toju Duke, Johnny Soraker, Ben Zevenbergen, Vinodkumar Prabhakaran, Mark Diaz, Ben Hutchinson, Kristen Olson, Alejandra Molina, Erin Hoffman-John, Josh Lee, Lora Aroyo, Ravi Rajakumar, Alena Butryna, Matthew Lamm, Viktoriya Kuzmina, Joe Fenton, Aaron Cohen, Rachel Bernstein, Ray Kurzweil, Blaise Aguera-Arcas, Claire Cui, Marian Croak, Ed Chi, and Quoc Le.  Lamda: Language models for dialog applications, 2022.
  * Wu et al. [2023] Yiran Wu, Feiran Jia, Shaokun Zhang, Han-Tai Li, Erkang Zhu, Yue Wang, Yin Tat Lee, Richard Peng, Qingyun Wu, and Chi Wang.  An empirical study on challenging math problem solving with gpt-4.  _ArXiv_ , abs/2306.01337, 2023.  URL https://api.semanticscholar.org/CorpusID:259063798.
  * Zhihong Shao [2024] Qihao Zhu Runxin Xu Junxiao Song Mingchuan Zhang Y.K. Li Y. Wu Daya Guo Zhihong Shao, Peiyi Wang.  Deepseekmath: Pushing the limits of mathematical reasoning in open language models, 2024.  URL https://arxiv.org/abs/2402.03300.

##  5 Appendix

##  6 Case Study on our MATH

In this case, there the same real number is presented in different form in the
answer and predicted point. Thus the evaluation toolkit must compare the
numerical equivalence between the number of coordinates to correctly assess
them.

    Example 1:
    Reason: Different Form
    Quesiton: Suppose $P$ is the point $(5,3)$ and $Q$ is the point $(-3,6)$. What
    is the midpoint of $\\overline{PQ}$?

    Ground Truth: \\left(1,\\frac{9}{2}\\right)
    prediction:   $(1.0, 4.5)$
    Human: True
    Tora: False
    MATH: False
    Ours: True

In this case, same expressions are not presented in the exact same form. The
evaluation toolkit may need to simplify the polynomial to the same level to
evaluate the them.

    Example 2:
    Reason: Expression
    Quesiton: Factor $36-4x^2$ completely.

    Ground Truth: 4(3-x)(3+x)
    prediction:   $-4(x - 3)(x + 3)$

    Human: True
    Tora: False
    MATH: False
    Ours: True

In the following case, the answer has a followed unit string, which should be
cleaned before evaluation.

    Example 3:
    Reason: Unit
    Quesiton: Diana can either invest $20,\\!000$ dollars for $4$ years with a simple
    interest rate of $6\\%$ or an interest rate of $7\\%$ which compounds quarterly.
    How many more dollars, rounded to the nearest dollar, would she get with the better
    interest rate than with the worse one?

    Ground Truth: 1599 \\text{ dollars}
    prediction:   $1,599.00$

    Human: True
    Tora: False
    MATH: False
    Ours: True

In this case, the answer is a vector, and the comparison should be applied in
element-wise manner.

    Example 4:
    Reason: Elementwise comparision
    Question: Convert the point $( 1, -1, -6 )$ in rectangular coordinates to cylindrical
    coordinates.  Enter your answer in the form $(r,\\theta,z),$ where $r > 0$ and
    $0 \\le \\theta < 2 \\pi.$

    Ground Truth: \\left( \\sqrt{2}, \\frac{7 \\pi}{4}, -6 \\right)
    Prediction:    (1.4142135623730951,5.497787143782138,-6)
    Human: True
    Tora: False
    MATH: False
    Ours: True

##  7 Case Study on our GAOKAO

In this case, the answer and the prediction are both intervals in different
form. They should be converted to a interval range to be properly compared.

    Example 1:
    Reason: Interval
    Quesiton: Given sets $M=\\{x|x+2\\geq 0\\},N=\\{x|x-1<0\\}$, find $M \\cap N$.

    Ground Truth: x \\in [-2, 1)
    prediction:   $[-2, 1)$
    Human: True
    Tora: False
    MATH: False
    Ours: True

For this case, the ground truth are expressed as a set, and it means the order
of the prediction does not matter.

    Example 2:
    Reason: Set
    Quesiton: If the universal set is $U=\\{1,2,3,4,5\\}$, and $M=\\{1,4\\},N=\\{2,5\\}$,
    find $N \\cup \\overline{M}$.

    Ground Truth: 2, 3, 5
    prediction:   5, 3, 2

    Human: True
    Tora: False
    MATH: False
    Ours: True

For this case, line functions are not presented in the exact same form. The
evaluation toolkit may need to simplify the polynomial on the right side to
the same level, then evaluate the function equation.

    Example 3:
    Reason: Expression
    Quesiton: Find the tangent line to the function $y=\\frac{e^{x}}{x+1}$ at point
    $\\left(1,\\frac{e}{2}\\right)$.

    Ground Truth: y=\\frac{e}{4}x+\\frac{e}{4}
    prediction:   $f(x) = \\frac{e(x + 1)}{4}$

    Human: True
    Tora: False
    MATH: False
    Ours: True

##  8 Case Study on MATH of ToRA

In this instance, a set of values should be answered to the question. So the
order of numbers separated by commas does not matter.

    Example 1:
    Reason: Set
    Quesiton: A line segment of length $5$ has one endpoint at $(1, 2)$ and
    the other endpoint at $(4, b)$. Find all possible values of $b$, separated by commas.

    Ground Truth: 6,-2
    prediction:   -2,6
    Human: True
    Tora: False
    MATH: False
    Ours: True


# Expressive Monotonic Neural Networks

Niklas Nolte , Ouail Kitouni††footnotemark:  , Mike Williams
The NSF AI Institute for Artificial Intelligence and Fundamental Interactions
Massachusetts Institute of Technology
Cambridge, MA 02139, USA
{nnolte<EMAIL_ADDRESS>
Equal contribution

###### 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 architecture111https://github.com/niklasnolte/MonotoneNorm 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.

##  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 & Martínez, 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-monotonicty 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

$\displaystyle f:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{n},\quad
d,n\in{\mathbb{N}},$ (1)

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

######  Definition 3.1 (Monotonicity).

Let ${\bm{x}}\in{\mathbb{R}}^{d}$,
${\bm{x}}_{\mathbb{S}}\equiv\mathbf{1}_{\mathbb{S}}\odot{\bm{x}}$, and the
Hadamard product of ${\bm{x}}$ with the indicator vector
$\mathbf{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

$\displaystyle
f({\bm{x}}_{\mathbb{S}}^{\prime}+{\bm{x}}_{\bar{{\mathbb{S}}}})_{i}\leq
f({\bm{x}}_{\mathbb{S}}+{\bm{x}}_{\bar{{\mathbb{S}}}})_{i}\quad\forall
i\in{\mathbb{Q}}\;\text{and}\;\forall{\bm{x}}_{\mathbb{S}}^{\prime}\leq{\bm{x}}_{\mathbb{S}},$
(2)

where $\bar{{\mathbb{S}}}$ denotes the complement of ${\mathbb{S}}$ and the
inequality on the right uses the product (or component-wise) order.

######  Definition 3.2 (Lipp function).

$g:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{n}$ is Lipp if it is Lipschitz
continuous with respect to the $L^{p}$ norm in every output dimension, _i.e._
,

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

###  3.1 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({\bm{x}})$ that is Lip1 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({\bm{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({\bm{x}})=g({\bm{x}})+\lambda(\mathbf{1}_{\mathbb{S}}\cdot{\bm{x}})=g({\bm{x}})+\lambda\sum_{i\in{\mathbb{S}}}{x}_{i}.$
(4)

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

$\displaystyle\frac{\partial g}{\partial{x}_{i}}$
$\displaystyle\in[-\lambda,\lambda],\ \ \forall\,i\in{\mathbb{N}}_{1:n}$ (5)
$\displaystyle\Rightarrow\frac{\partial f}{\partial{x}_{i}}$
$\displaystyle=\frac{\partial g}{\partial{x}_{i}}+\lambda\geq 0\ \
\forall\,{\bm{x}}\in{\mathbb{R}}^{n},i\in{\mathbb{S}}\,.$ (6)

##### The importance of the norm choice

The 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.

Figure 1: $p$-norm constrained gradients showing (red) $p=2$ and (green)
$p=1$. The gradient of a function $g({\bm{x}})$ that is Lipp=2 resides within
the dashed red line. For a Lipp=1 function, the boundary is the green dashed
line. Note that ${\bm{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$.

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

###  3.2 Lipp=1 approximators

Our goal is to construct a universal approximator of Lip1 functions, i.e., we
would like the hypothesis class to have two properties:

  1. 1.

It always satisfies Eq. 3, i.e., be Lip1.

  2. 2.

It is able to fit all possible Lip1 functions. In particular, the bound in Eq.
3 needs to be attainable $\forall\,{\bm{x}},{\bm{y}}$.

##### Lip1 constrained models

To 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

$\displaystyle{\bm{z}}^{l}=\sigma({\bm{z}}^{l-1}){\bm{W}}^{l}+{\bm{b}}^{l},$
(7)

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

$\displaystyle\prod_{i=1}^{D}\|{\bm{W}}^{i}\|_{1}\leq\lambda,$ (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:

$\displaystyle{\bm{W}}^{i}\rightarrow{\bm{W}}^{\prime
i}=\lambda^{1/D}\frac{{\bm{W}}^{i}}{\text{max}(1,\|{\bm{W}}^{i}\|_{1})}\qquad{\rm
or}\qquad W^{i}\rightarrow W^{\prime
i}=\frac{{\bm{W}}^{i}}{\text{max}(1,\lambda^{-1/D}\cdot\|{\bm{W}}^{i}\|_{1})}\,.$
(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:

$\displaystyle{\bm{W}}^{i}\rightarrow{\bm{W}}^{\prime
i}={\bm{W}}^{i}\text{diag}\left(\frac{1}{\text{max}\left(1,\lambda^{-1/D}\cdot\sum_{j}|{W}^{i}_{jk}|\right)}\right)\,.$
(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 Lip1. Anil et
al. (2019) provide an alternative normalization scheme:

$\displaystyle\|{\bm{W}}^{1}\|_{1,\infty}\cdot\prod_{i=2}^{m}\|{\bm{W}}^{i}\|_{\infty}\leq\lambda$
(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 power

Guaranteeing that the model is Lipschitz bounded is not sufficient, it must
also able to saturate the bound to be able to model all possible Lip1
functions. Some Lipschitz network architectures, e.g. Miyato et al. (2018),
tend to over constrain the model such that it cannot fit all Lip1 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:

$\displaystyle\sigma_{G}({\bm{x}})$
$\displaystyle=\text{sort}_{1:G}({\bm{x}}_{1:G})+\text{sort}_{G+1:2G}({\bm{x}}_{G+1:2G})+\ldots$
$\displaystyle=\sum_{i=0}^{n/G-1}\text{sort}_{iG+1:(i+1)G}({\bm{x}}_{iG+1:(i+1)G}),$
(12)

where ${\bm{x}}\in{\mathbb{R}}^{n}$,
${\bm{x}}_{i:j}=\mathbf{1}_{i:j}\odot{\bm{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 Lip1
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:

$\displaystyle\sigma({\bm{z}})=\begin{cases}{\bm{z}}&{\bm{z}}{\bm{v}}>0\\\
{\bm{z}}(\mathbf{I}-2{\bm{v}}{\bm{v}}^{T})&{\bm{z}}{\bm{v}}\leq 0\\\
\end{cases}$ (13)

Here, ${\bm{z}}$ is the preactivation row vector and ${\bm{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({\bm{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.

###  4.1 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.

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.

###  4.2 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\,{\rm 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_{\rm 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_{\rm T}$, though unfortunately they can also have
extremely large values of both displacement and momentum.

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 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 horiztonal and vertical axes are
related to how long the particle lived before decaying and how massive the
particle was, respectively. 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.

###  4.3 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 & 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 &
Graff, 2017) and HeartDisease (Gennari et al., 1989). Results are shown in
Table 1.

##### COMPAS

(Correctional Offender Management Profiling for Alternative Sanctions) refers
to a commercial algorithm used by judges and police officers to determine the
likelihood of reoffense. Larson & Kirchner (2016) discusses that the algorithm
is racially biased and provides a dataset from a two-year study of the real-
world performance of COMPAS. The task here is to predict the reoffense
probability within the next two years. The dataset has 13 features, 4 of which
have a monotonic inductive bias, and contains a total of 6172 data points.

##### BlogFeedBack

This dataset contains 54270 data points with 276 dimensions describing blog
posts. The task is to predict the number of comments following the post
publication within 24 hours. 8 of the features have a monotonic inductive
bias. Just like Liu et al. (2020), we also only consider the 90% of the data
points with the smallest targets so as to not let the RMSE be dominated by
outliers.

##### LoanDefaulter

The version of this dataset available on Kaggle was updated on a yearly basis
up to 2015. Kaggle (2015) contains a link that is, we believe, a superset of
the data used in Liu et al. (2020). Luckily, the authors have shared with us
the exact version of the dataset they used in their studies for an appropriate
comparison. The data is organized in 28 features and the task is to determine
loan defaulters. The classification score should be monotonic in 5 features:
non-decreasing in number of public record bankruptcies and Debt-to-Income
ratio, non-increasing in credit score, length of employment and annual income.

##### ChestXRay

This dataset contains tabular data and images of patients with diseases that
are visible in a chest x-ray. The task is to predict whether or not the
patient has such a disease. Just like Liu et al. (2020), we send the image
through an ImageNet-pretrained ResNet18 (He et al., 2016). The penultimate
layer output concatenated with tabular data acts as input to the monotonic
architecture. Two of the four tabular features are monotonic. In the bottom
right table in 1, there are two entries for our architecture. The E-E entry
refers to end-to-end training with ResNet18, whereas the other experiment
fixes the ResNet weights.

##### AutoMPG

(Dua & Graff, 2017) This is a dataset containing 398 examples of cars,
described by 7 numerical features and the model name. The target, MPG, is
monotonically decreasing with 3 of the features. The name is not used as a
feature.

##### HeartDisease

(Gennari et al., 1989) is a dataset of patients, described by 13 features. The
task is to determine whether or not the patient has heart disease.

As can be seen in Table 1, our Lipschitz monotonic networks perform
competitively or better than the state-of-the-art on all benchmarks we tried.

COMPAS | BlogFeedback
---|---
| Method | Parameters |  $\upuparrows$ Test Acc
---|---|---
Certified | 23112 | $(68.8\pm 0.2)\%$
LMN | 37 | $\mathbf{(69.3\pm 0.1)\%}$
| Method | Parameters |  $\downdownarrows$ RMSE
---|---|---
Certified | 8492 | $.158\pm.001$
LMN | 2225 | $\mathbf{.160\pm.001}$
LMN mini | 177 | $\mathbf{.155\pm.001}$
LoanDefaulter | ChestXRay
| Method | Parameters |  $\upuparrows$ Test Acc
---|---|---
Certified | 8502 | $(65.2\pm 0.1)\%$
LMN | 753 | $\mathbf{(65.44\pm 0.03)\%}$
LMN mini | 69 | $\mathbf{(65.28\pm 0.01)\%}$
| Method | Parameters |  $\upuparrows$ Test Acc
---|---|---
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)\%}$
Heart Disease | Auto MPG
| Method |  $\upuparrows$ Test Acc
---|---
COMET | $(86\pm 3)\%$
LMN | $\mathbf{(89.6\pm 1.9)\%}$
| Method |  $\downdownarrows$ MSE
---|---
COMET | $(8.81\pm 1.81)\%$
LMN | $\mathbf{(7.58\pm 1.2)\%}$
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”.

It is also immediately apparent that our architecture is highly expressive. We
manage to train tiny networks with few parameters while still achieving
competitive performance. Given that some of these datasets have a significant
number of features compared to our chosen network width, most parameters are
in the weights of the first layer. We manage to build and train even smaller
networks with better generalization performance when taking only a few
important features. These networks are denoted with mini in Table 1. Because
all of the presented architectures are small in size, we show practical finite
sample expressiveness for harder tasks and larger networks by achieving 100%
training accuracy on MINST, CIFAR-10, and CIFAR-100 with real and random
labels as well as 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 Lip1 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.

##  7 Reproducibility Statement

All experiments with public datasets are reproducible with the code provided
at https://github.com/niklasnolte/monotonic_tests. This code uses the package
available in 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 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/).

## References

  * Akhtar et al. (2021) Naveed Akhtar, Ajmal Mian, Navid Kardan, and Mubarak Shah.  Advances in adversarial attacks and defenses in computer vision: A survey.  _IEEE Access_ , 9:155161–155196, 2021.
  * Anil et al. (2019) Cem Anil, James Lucas, and Roger Grosse.  Sorting out Lipschitz function approximation.  In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 291–301. PMLR, 09–15 Jun 2019.  URL http://proceedings.mlr.press/v97/anil19a.html.
  * Arpit et al. (2019) Devansh Arpit, Victor Campos, and Yoshua Bengio.  How to initialize your network? robust initialization for weightnorm & resnets, 2019.
  * Auguste et al. (2020) Charles Auguste, Sean Malory, and Ivan Smirnov.  A better method to enforce monotonic constraints in regression and classification trees, 2020.
  * Behrmann et al. (2019) Jens Behrmann, Will Grathwohl, Ricky T. Q. Chen, David Duvenaud, and Jörn-Henrik Jacobsen.  Invertible residual networks, 2019.
  * Buza (2014) Krisztian Buza.  Feedback prediction for blogs.  In _Data analysis, machine learning and knowledge discovery_ , pp. 145–152. Springer, 2014.
  * De Cao et al. (2020) Nicola De Cao, Wilker Aziz, and Ivan Titov.  Block neural autoregressive flow.  In _Uncertainty in artificial intelligence_ , pp. 1263–1273. PMLR, 2020.
  * Dua & Graff (2017) Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Gennari et al. (1989) John H. Gennari, Pat Langley, and Doug Fisher.  Models of incremental concept formation.  _Artificial Intelligence_ , 40(1):11–61, 1989\.  ISSN 0004-3702.  doi: https://doi.org/10.1016/0004-3702(89)90046-5.  URL https://www.sciencedirect.com/science/article/pii/0004370289900465.
  * Gouk et al. (2020) Henry Gouk, Eibe Frank, Bernhard Pfahringer, and Michael J. Cree.  Regularisation of neural networks by enforcing lipschitz continuity, 2020\.
  * (11) Akhil Gupta, Naman Shukla, Lavanya Marla, Arinbjörn Kolbeinsson, and Kartik Yellepeddi.  How to incorporate monotonicity in deep networks while preserving flexibility?  URL https://arxiv.org/abs/1909.10662.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Huang et al. (2018) Chin-Wei Huang, David Krueger, Alexandre Lacoste, and Aaron Courville.  Neural autoregressive flows.  In _International Conference on Machine Learning_ , pp. 2078–2087. PMLR, 2018.
  * Huster et al. (2018) Todd Huster, Cho-Yu Jason Chiang, and Ritu Chadha.  Limitations of the lipschitz constant as a defense against adversarial examples, 2018.
  * Huster et al. (2019) Todd Huster, Cho-Yu Jason Chiang, and Ritu Chadha.  Limitations of the lipschitz constant as a defense against adversarial examples.  In _ECML PKDD 2018 Workshops: Nemesis 2018, UrbReas 2018, SoGood 2018, IWAISe 2018, and Green Data Mining 2018, Dublin, Ireland, September 10-14, 2018, Proceedings 18_ , pp. 16–29. Springer, 2019.
  * Kaggle (2015) Kaggle.  Lending club loan data.  2015\.  URL https://www.kaggle.com/datasets/wordsforthewise/lending-club.
  * Kitouni et al. (2021) Ouail Kitouni, Niklas Nolte, and Mike Williams.  Robust and provably monotonic networks.  _arXiv preprint arXiv:2112.00038_ , 2021.
  * (18) Donald Knuth.  Donald Knuth’s homepage.  https://www-cs-faculty.stanford.edu/~knuth/faq.html.  Accessed: 2022-09-28.
  * Larson & Kirchner (2016) SM Angwin J Larson and Lauren Kirchner.  There’s software used across the country to predict future criminals. and it’s biased against blacks, 2016.
  * Liu et al. (2020) Xingchao Liu, Xing Han, Na Zhang, and Qiang Liu.  Certified monotonic neural networks, 2020.  URL https://arxiv.org/abs/2011.10219.
  * Milani Fard et al. (2016) Mahdi Milani Fard, Kevin Canini, Andrew Cotter, Jan Pfeifer, and Maya Gupta.  Fast and flexible monotonic functions with ensembles of lattices.  In D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 29. Curran Associates, Inc., 2016.  URL https://proceedings.neurips.cc/paper/2016/file/c913303f392ffc643f7240b180602652-Paper.pdf.
  * Miyato et al. (2018) Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida.  Spectral Normalization for Generative Adversarial Networks.  _arXiv e-prints_ , art. arXiv:1802.05957, February 2018.
  * Nguyen & Martínez (2019) A. Nguyen and M. Martínez.  Mononet: Towards interpretable models by learning monotonic features.  9 2019.  doi: 10.48550/arxiv.1909.13611.  URL https://arxiv.org/abs/1909.13611v1.
  * Sill (1998) Joseph Sill.  Monotonic networks.  In M. Jordan, M. Kearns, and S. Solla (eds.), _Advances in Neural Information Processing Systems_ , volume 10. MIT Press, 1998.  URL https://proceedings.neurips.cc/paper/1997/file/83adc9225e4deb67d7ce42d58fe5157c-Paper.pdf.
  * Sill & Abu-Mostafa (1996) Joseph Sill and Yaser Abu-Mostafa.  Monotonicity hints.  In M.C. Mozer, M. Jordan, and T. Petsche (eds.), _Advances in Neural Information Processing Systems_ , volume 9. MIT Press, 1996.  URL https://proceedings.neurips.cc/paper/1996/file/c850371fda6892fbfd1c5a5b457e5777-Paper.pdf.
  * Singla et al. (2021) Sahil Singla, Surbhi Singla, and Soheil Feizi.  Improved deterministic l2 robustness on cifar-10 and cifar-100.  In _International Conference on Learning Representations_ , 2021.
  * Sivaraman et al. (2020) Aishwarya Sivaraman, Golnoosh Farnadi, Todd Millstein, and Guy Van den Broeck.  Counterexample-guided learning of monotonic neural networks.  _Advances in Neural Information Processing Systems_ , 33:11936–11948, 2020.
  * Wang et al. (2017) Xiaosong Wang, Yifan Peng, Le Lu, Zhiyong Lu, M Bagheri, and R Summers.  Hospital-scale chest x-ray database and benchmarks on weakly-supervised classification and localization of common thorax diseases.  In _IEEE CVPR_ , volume 7, 2017.  URL https://www.kaggle.com/datasets/nih-chest-xrays/sample.
  * Wehenkel & Louppe (2019) Antoine Wehenkel and Gilles Louppe.  Unconstrained monotonic neural networks.  _Advances in neural information processing systems_ , 32, 2019.
  * You et al. (2017) Seungil You, David Ding, Kevin Canini, Jan Pfeifer, and Maya Gupta.  Deep lattice networks and partial monotonic functions, 2017.
  * Zhang et al. (2021) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning (still) requires rethinking generalization.  _Communications of the ACM_ , 64(3):107–115, 2021\.

##  Appendix A Expressive power of the architecture

Robust architectures like Lipschitz constrained networks are often believed to
be much less expressive than their unconstrained counterparts Huster et al.
(2019). Here we show that our architecture is capable of (over)fitting complex
decision boundaries even on random labels in a setup simular to Zhang et al.
(2021).

We show the finite sample expressiveness of the architecture in
https://github.com/okitouni/Lipschitz-network-bench by fitting MNIST, CIFAR10,
CIFAR100 with normal and random labels to 100% training accuracy. We also
train on CIFAR100 with an additional “goodness” feature $x\in[0,1]$ to
showcase the monotonicity aspect of the architecture. This dataset is referred
to as CIFAR101 below. The synthetic monotonicity problem is currently
implemented such that samples with values above a critical threshold in the
goodness feature $x>x_{\rm crit}$ are labeled 0. An alternative implementation
is to take label 0 with probability $x$ and keep the original label (or assign
a random one) with probability $1-x$. Table 2 summarizes the setup used for
training. We use Adam with default hyper-parameters im all experiments.

Task | Width | Depth | LR | EPOCHS | Batchsize | Loss
---|---|---|---|---|---|---
MNIST | $1024$ | $3$ | $10^{-5}$ | $10^{5}$ | ALL | CE($\tau=256$)
CIFAR10 | $1024$ | $3$ | $10^{-5}$ | $10^{5}$ | ALL | CE($\tau=256$)
CIFAR100/101 | $1024$ | $3$ | $10^{-5}$ | $10^{5}$ | ALL | CE($\tau=256$)
Table 2: Training MNIST and CIFAR10/100 to 100% training accuracy with
Lipschitz networks.