Revisiting Tropical Polynomial Division: Theory, Algorithms and
  Application to Neural Networks
Tropical geometry has recently found several applications in the analysis of
neural networks with piecewise linear activation functions. This paper presents
a new look at the problem of tropical polynomial division and its application
to the simplification of neural networks. We analyze tropical polynomials with
real coefficients, extending earlier ideas and methods developed for
polynomials with integer coefficients. We first prove the existence of a unique
quotient-remainder pair and characterize the quotient in terms of the convex
bi-conjugate of a related function. Interestingly, the quotient of tropical
polynomials with integer coefficients does not necessarily have integer
coefficients. Furthermore, we develop a relationship of tropical polynomial
division with the computation of the convex hull of unions of convex polyhedra
and use it to derive an exact algorithm for tropical polynomial division. An
approximate algorithm is also presented, based on an alternation between data
partition and linear programming. We also develop special techniques to divide
composite polynomials, described as sums or maxima of simpler ones. Finally, we
present some numerical results to illustrate the efficiency of the algorithms
proposed, using the MNIST handwritten digit and CIFAR-10 datasets.
Ioannis Kordonis, Petros Maragos
# Revisiting Tropical Polynomial Division: Theory, Algorithms and Application to Neural Networks.

###### Abstract

Tropical geometry has recently found several applications in the analysis of neural networks with piecewise linear activation functions. This paper presents a new look at the problem of tropical polynomial division and its application to the simplification of neural networks. We analyze tropical polynomials with real coefficients, extending earlier ideas and methods developed for polynomials with integer coefficients. We first prove the existence of a unique quotient-remainder pair and characterize the quotient in terms of the convex bi-conjugate of a related function. Interestingly, the quotient of tropical polynomials with integer coefficients does not necessarily have integer coefficients. Furthermore, we develop a relationship of tropical polynomial division with the computation of the convex hull of unions of convex polyhedra and use it to derive an exact algorithm for tropical polynomial division. An approximate algorithm is also presented, based on an alternation between data partition and linear programming. We also develop special techniques to divide composite polynomials, described as sums or maxima of simpler ones. Finally, we present some numerical results to illustrate the efficiency of the algorithms proposed, using the MNIST handwritten digit and CIFAR-10 datasets.

 Tropical geometry, Piecewise linear neural networks, Tropical polynomial division, Neural network compression

## I Introduction

Tropical geometry is a relatively new research field combining elements and ideas from polyhedral and algebraic geometry [1]. The underlying algebraic structure is the max-plus semiring (also known as tropical semiring), where the usual addition is replaced by maximization and the standard multiplication by addition. Tropical polynomials (also known as max-polynomials) are the polynomials in the max-plus algebra and have a central role in tropical geometry. This work deals with the division of tropical polynomials and presents some applications in neural network compression.

A link between tropical geometry and machine learning was recently developed [2, 3]. An important application is the analysis of neural networks with piecewise linear activations (e.g., ReLU). In this front, [3, 4, 5, 6, 7] use a tropical representation of neural networks to describe the complexity of a network structure, defined as the number of its linear regions, or describe their decision boundaries. [8] presents an algorithm to extract the linear regions of deep ReLU neural networks. Papers [9, 10, 11] deal with the problem of tropical polynomial division and its use in the simplification of neural networks. The analysis is, however, restricted to polynomials with integer coefficients. Article [12] uses tropical representations of neural networks and employs polytope approximation tools to simplify them. Another class of networks represented in terms of tropical algebra is morphological neural networks [13, 14, 15, 16, 17, 18]. Other applications of tropical algebra and geometry include the tropical modeling of classical algorithms in probabilistic graphical models [19, 20] and piecewise linear regression [21, 22]. Neural networks employing the closely related log-sum-exp nonlinearity were presented in [23, 24].

Another problem very much related to the current work is the factorization of tropical polynomials. This problem was first studied in [25, 26, 27, 28], for the single-variable case, and in [29, 30] for the multivariate case. The problem of tropical rational function simplification was studied in [31]. Another related problem is approximation in tropical algebra (see, e.g., [32]).

_Contribution:_ This work deals with the problem of tropical polynomial division and its applications in simplifying neural networks. The analysis generalizes previous works [9, 10, 11] to include polynomials with real coefficients. We first introduce a new notion of division and show that there is a unique quotient-remainder pair. Furthermore, we show that the quotient is equal to the convex bi-conjugate of the difference between the dividend and the divisor. Then, the quotient is characterized in terms of the closed convex hull of the union of a finite collection of convex (unbounded) polyhedra. This characterization leads to a simple exact algorithm based on polyhedral geometry. We propose an efficient approximate scheme inspired by an optimization algorithm for convex, piecewise-linear fitting [33]. The approximate algorithm alternates between data clustering and linear programming. Then, we focus on developing algorithms for dividing composite polynomials expressed as the sum of simpler ones. The division results are then applied to simplify neural networks with ReLU activation functions. The resulting neural network has fewer neurons with maxout activations and a reduced total number of parameters. Finally, we present numerical results for the MNIST handwritten digit recognition and the CIFAR-10 image classification problem. As a baseline for comparison, we use the structured L1 pruning without retraining. The proposed method is very competitive in the large compression regime, that is, in the case where there are
very few parameters remaining after compression.

The rest of the paper is organized as follows: Section II presents some preliminary algebraic and geometric material. In Section III, we develop some theoretical tools for tropical polynomial division. An exact division algorithm is presented in Section IV, and an approximate algorithm in Section V. Section VI-B deals with the division of composite tropical polynomials. Finally, VII gives some numerical examples. The proofs of the theoretical results are presented in the Appendix.

## II Preliminaries, Background and notation

Max-plus algebra is a modification of the usual algebra of the real numbers, where the usual summation is substituted by maximum and the usual multiplication is substituted by summation [34, 35, 36, 37]. Particularly, max-plus algebra is defined on the set \(\mathbb{R}_{\text{max}}=\mathbb{R}\cup\{-\infty\}\) with the binary operations "\(\vee\)" and "\(+\)". The operation "\(+\)" is the usual addition. The operation "\(\vee\)" stands for the maximum, i.e. for \(x,y\in\mathbb{R}_{\text{max}}\), it holds \(x\lor y=\text{max}\{x,y\}\). For more than two terms we use"\(\vee\)", i.e., \(\bigvee_{i\in I}x_{i}=\sup\{x_{i}:i\in I\}\). We will also use the '\(\wedge\)' notation to denote the minimum, i.e., \(x\wedge y=\text{min}\{x,y\}\). A unified view of max-plus and several related algebras, was developed in [38], in terms of weighted lattices.

A tropical polynomial is a polynomial in the max-plus algebra. Particularly, a tropical polynomial function [39, 40, 41] is defined as \(p:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\text{max}}\) with

\[p(\mathbf{x})=\bigvee_{i=1}^{m_{p}}(\mathbf{a}_{i}^{T}\mathbf{x}+b_{i}), \tag{1}\]

where \(\mathbf{a}_{i}\in\mathbb{R}^{n}\) and \(b_{i}\in\mathbb{R}_{\text{max}}\). In this paper, we will use the terms tropical polynomial and tropical polynomial function interchangeably.

For the tropical polynomial \(p(\mathbf{x}),\) the Newton polytope is defined as

\[\text{Newt}(p)=\text{ conv}\{\mathbf{a}_{1},\ldots,\mathbf{a}_{m_{p}}\}, \tag{2}\]

where "conv" stands for the convex hull, and extended Newton polytope as

\[\text{ENewt}(p)=\text{ conv}\{[\mathbf{a}_{1}^{T}\ b_{1}]^{T},\ldots,[\mathbf{a}_{m_{p }}^{T}b_{m_{p}}]^{T}\}. \tag{3}\]

The values of the polynomial function \(p\) depend on the upper hull of the extended Newton polytope [1] (see also [4]). A polyhedron is subset of \(\mathbb{R}^{n}\) which can be represented either as the set of solutions of a finite collection of linear inequalities (\(\mathcal{H}-\)representation)

\[P=\{\mathbf{x}\in\mathbb{R}^{n}:\mathbf{A}\mathbf{x}\leq\mathbf{b}\},\]

or in terms of a set of vertices \(\mathbf{v}_{1},\ldots,\mathbf{v}_{l}\in\mathbb{R}^{n}\) and a set of rays \(\mathbf{w}_{1},\ldots,\mathbf{w}_{s}\in\mathbb{R}^{n}\) as

\[P=\left\{\sum_{j=1}^{l}\lambda_{j}\mathbf{v}_{j}+\sum_{j=1}^{s}\mu_{j}\mathbf{w}_{j}: \lambda_{j}\geq 0,\sum_{j=1}^{s}\lambda_{j}=1,\mu_{j}\geq 0\right\}, \tag{4}\]

(\(\mathcal{V}-\)representation). The Minkowski-Weyl (resolution) theorem states that each polyhedron admits both representations [42]. Figure 1 illustrates the Minkowski-Weyl resolution of an unbounded polyhedron. the Furthermore, there are several well-known algorithms for representation conversion (e.g., [42], Ch. 9). Bounded polyhedra are called polytopes. It is often convenient to describe polytopes in terms of extended representations (e.g., [43]). An extension of a polytope \(P\subset\mathbb{R}^{n}\) is a polyhedron \(Q\subset\mathbb{R}^{n^{\prime}}\), with \(n^{\prime}>n\), along with a linear projection \(F:R^{n^{\prime}}\rightarrow\mathbb{R}^{n}\) such that \(F(Q)=P\).

For a function \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{\infty\}\) the epigraph is defined as

\[\text{epi}(f)=\{(\mathbf{x},z)\in\mathbb{R}^{n+1}:z\geq f(\mathbf{x})\}.\]

The convex conjugate of \(f\) is a function \(f^{\star}:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\{+\infty,-\infty\}\) given by (see [44])

\[f^{\star}(\mathbf{a})=\sup\{\mathbf{a}^{T}\mathbf{x}-f(\mathbf{x}):\mathbf{x}\in\mathbb{R}^{n}\}.\]

For a pair of polyhedra \(P_{1},P_{2}\subset\mathbb{R}^{n}\) the Minkowski sum is defined as

\[P_{1}\oplus P_{2}=\{x+y:x\in P_{1},y\in P_{2}\}.\]

The following properties hold (see [4])

\[\text{Newt}(p_{1}+p_{2}) =\text{Newt}(p_{1})\oplus\text{Newt}(p_{2}),\] \[\text{Newt}(p_{1}\lor p_{2}) =\text{conv}(\text{Newt}(p_{1})\cup\text{Newt}(p_{2})).\]

Similar relations hold for the extended Newton polytopes as well.

We then introduce a dual function of a tropical polynomial, called Extended Newton Function (ENF). For a tropical polynomial \(p\), the ENF is defined as:

\[\text{ENF}_{p}(\mathbf{a})=\sup\{b\in\mathbb{R}:[\mathbf{a}^{T}\ b]^{T}\in\text{ENewt} (p)\},\]

where \(\mathbf{a}\in\mathbb{R}^{n}\) and \(\sup(\emptyset)=-\infty\). This function describes the upper hull of the extended Newton polytope. Thus, the ENF fully characterizes the tropical polynomial function \(p\).

**Proposition 1**: _Let \(p_{1},p_{2}\) be two tropical polynomials. Then \(p_{1}(\mathbf{x})\leq p_{2}(\mathbf{x})\), for all \(\mathbf{x}\in\mathbb{R}^{n}\), if and only if \(\text{ENF}_{p_{1}}(\mathbf{a})\leq\text{ENF}_{p_{2}}(\mathbf{a})\), for all \(\mathbf{a}\in\mathbb{R}^{n}\)._

_Proof_: See Appendix A.

Fig. 1: _Illustration of the Minkowski-Weil theorem. The unbounded polyhedron \(P\) (green color in the figure) has two representations: \(\mathcal{H}-\)representation as the intersection of halfspaces \(H_{1},\ldots,H_{4}\) and \(\mathcal{V}-\)representation generated by vertices \(v_{1},v_{2},v_{3}\) and rays \(w_{1},w_{2}\)._
The significance of ENFs is that we can often express inequality \(\text{ENF}_{p_{1}}\leq\text{ENF}_{p_{2}}\), using a finite number of inequalities in terms of the tropical polynomial coefficients of \(p_{1}\), \(p_{2}\).

## III Tropical Polynomial Division Definition and Existence

We first define the tropical polynomial division.

**Definition 1**: _Let \(p,d\) be tropical polynomials. We define the quotient and the remainder of the division of \(p\) by \(d\)._

1. _A tropical polynomial_ \(q\) _is the quotient if_ \[p(\boldsymbol{x})\geq q(\boldsymbol{x})+d(\boldsymbol{x}),\qquad\text{for all } \boldsymbol{x}\in\mathbb{R}^{n},\] (5) _and_ \(q\) _is the maximum polynomial satisfying this inequality. Particularly, for any tropical polynomial_ \(\tilde{q}\)_, such that_ \(p(\boldsymbol{x})\geq\tilde{q}(\boldsymbol{x})+d(\boldsymbol{x})\)_, for all_ \(\boldsymbol{x}\in\mathbb{R}^{n}\)_, it holds_ \(\tilde{q}(\boldsymbol{x})\leq q(\boldsymbol{x})\)_, for all_ \(\boldsymbol{x}\in\mathbb{R}^{n}\)_._
2. _A tropical polynomial_ \(r\) _is the remainder of the division of_ \(p\) _by_ \(d\) _with quotient_ \(q\) _if_ \[p(\boldsymbol{x})=(q(\boldsymbol{x})+d(\boldsymbol{x}))\lor r(\boldsymbol{x}),\qquad\text{for all }\boldsymbol{x}\in\mathbb{R}^{n},\] (6) _and_ \(r\) _is the minimum tropical polynomial satisfying this equality. Particularly, for any tropical polynomial_ \(\tilde{r}\)_, such that_ \(p(\boldsymbol{x})=(q(\boldsymbol{x})+d(\boldsymbol{x}))\lor\tilde{r}( \boldsymbol{x})\)_, it holds_ \(\tilde{r}(\boldsymbol{x})\geq r(\boldsymbol{x})\)_, for all_ \(\boldsymbol{x}\in\mathbb{R}^{n}\)_._

Observe that the set of tropical polynomials \(q\) satisfying (5) is closed under the pointwise maximum. That is, assume that \(q_{1},q_{2}\) are tropical polynomials such that \(p(\boldsymbol{x})\geq d(\boldsymbol{x})+q_{1}(\boldsymbol{x})\) and \(p(\boldsymbol{x})\geq d(\boldsymbol{x})+q_{2}(\boldsymbol{x})\), then it holds \(p(\boldsymbol{x})\geq d(\boldsymbol{x})+(q_{1}(\boldsymbol{x})\lor q_{2}( \boldsymbol{x}))\). This property motivates us to determine the monomials

\[q_{M}^{\boldsymbol{a},b}(\boldsymbol{x})=\boldsymbol{a}^{T}\boldsymbol{x}+b,\]

that satisfy (5). For each monomial coefficient \(\boldsymbol{a}\), define the maximum value of \(b\) that satisfy \(p(\boldsymbol{x})\geq d(\boldsymbol{x})+q_{M}^{\boldsymbol{a},b}(\boldsymbol{ x})\) as

\[l(\boldsymbol{a})=\sup\{b:q_{M}^{\boldsymbol{a},b}(\boldsymbol{x})+d( \boldsymbol{x})\leq p(\boldsymbol{x}),\text{for all }x\in\mathbb{R}^{n}\}. \tag{7}\]

Note that \(l(\boldsymbol{a})\) can be written as

\[l(\boldsymbol{a}) =\sup\{b:b\leq p(\boldsymbol{x})-d(\boldsymbol{x})-\boldsymbol{a }^{T}\boldsymbol{x},\text{for all }x\in\mathbb{R}^{n}\}\] \[=\inf_{\boldsymbol{x}\in\mathbb{R}^{n}}\{p(\boldsymbol{x})-d( \boldsymbol{x})-\boldsymbol{a}^{T}\boldsymbol{x}\}.\]

Denoting by \(f(\boldsymbol{x})\) the difference \(p(\boldsymbol{x})-d(\boldsymbol{x})\), we have

\[l(\boldsymbol{a})=-f^{\star}(\boldsymbol{x}),\]

where \(f^{\star}\) is the convex conjugate of \(f\). A candidate for the quotient is

\[q(\boldsymbol{x})=\sup_{\boldsymbol{a}\in\mathbb{R}^{n}}\{\boldsymbol{a}^{T} \boldsymbol{x}+l(\boldsymbol{a})\}.\]

Note that \(q\) is the bi-conjugate of \(f\). Indeed

\[q(\boldsymbol{x})=\sup_{\boldsymbol{a}\in\mathbb{R}^{n}}\{\boldsymbol{a}^{T} \boldsymbol{x}-f^{\star}(\boldsymbol{a})\}=f^{\star\star}(\boldsymbol{x}). \tag{8}\]

The following proposition shows that \(q(\boldsymbol{x})\) is indeed the quotient of the division.

**Proposition 2**: _For any pair of tropical polynomials \(p,d\) there is a quotient-remainder pair. The quotient is given by (8). Furthermore, the polynomial functions of the quotient and remainder are unique._

See Appendix B.

**Remark 1**: _The use of the convex conjugate is closely related to the slope transform [45, 46, 47]. The slope transform was used to derive an analog of frequency response of max-plus and more generally morphological systems._

**Example 1** (Tropical Polynomial Division in 1D): _Let \(p(x)=\max(-2x-1,1,x+1,3x-3)\) and \(d(x)=\max(x,2x-1)\). Equivalently, the polynomial functions \(p\) and \(d\) are written as_

\[p(x)=\begin{cases}-2x-1,&\text{if }x\leq-1\\ 1,&\text{if }-1<x\leq 0\\ x+1,&\text{if }0<x\leq 2\\ 3x-3,&\text{if }x\geq 2\end{cases},\]

\[d(x)=\begin{cases}x,&\text{if }x\leq 1\\ 2x-1,&\text{if }x>1\end{cases}.\]

_Thus,_

\[f(x)=p(x)-d(x)=\begin{cases}-3x-1,&\text{if }x\leq-1\\ 1-x,&\text{if }-1<x\leq 0\\ 1,&\text{if }0<x\leq 1\\ -x+2,&\text{if }1<x\leq 2\\ x-2,&\text{if }x\geq 2\end{cases}.\]

_The plots of \(p,d\) and \(f\) are given in Figure 2. The difference \(f\) is not convex since the slopes \(-3,-1,0,-1,1\) are not increasing. The quotient \(q\) is given as the convex bi-conjugate of \(f\), that is the largest convex function which is less than or equal to \(f\), for all \(x\). Thus, the quotient is given \(q(x)=\max(-3x+1,1-x,-0.5x+1,x-2)\). The sum \(d(x)+q(x)\) is equal to \(p(x)\), for all \(x\), except \(x\in(0,2)\). Thus, \(r(x)=x+1\), that is, the form of \(p(x)\), for \(x\in(0,2)\)._

We call the division _nontrivial_ if \(q(\boldsymbol{x})>-\infty\) and _effective_ if \(r(\boldsymbol{x})\neq p(\boldsymbol{x})\), for some \(\boldsymbol{x}\in\mathbb{R}^{n}\). It is easy to see that an effective division is also nontrivial. As shown in the previous example, these notions are not equivalent.

**Proposition 3**: _Assume that \(p\) and \(d\) are tropical polynomials. If the division of \(p\) by \(d\) has quotient \(q\) and remainder \(r\). Then:_

Fig. 2: _This figure refers to Example 1. The left plot presents the divided \(p\) and the divisor \(d\). The middle plot their difference \(f\) and the quotient \(q\). The right part shows the dividend \(p\), its approximation \(d+q\) and the remainder \(r\)._
1. The division of \(r\) by \(d\) is not effective.
2. The division of \(p\) by \(q\) has quotient \(d\) and remainder \(r\).

_Proof_ See Appendix C.

_Remark 2:_ For the Euclidean division of a positive integer \(p\) by a positive integer \(d\) with quotient \(q\) and remainder \(r\), the notion of nontrivial division corresponds to \(q\neq 0\) and effective division to \(r<p\). For the division of positive integers, we observe that these notions are equivalent. Furthermore, the converse of (a) is true. Particularly, if for some \(\bar{q},\bar{r}\in\mathbb{N}\), it holds \(p=d\bar{q}+\bar{r}\) and the division of \(r\) by \(d\) is not effective, then \(\bar{q}\) is the quotient and \(\bar{r}\) the remainder.

On the other hand, for tropical polynomials, the converse of (a) is not true. Particularly, it is possible that for a pair of tropical polynomials \(p,d\) there is another pair \(q,r\) such that (6) is satisfied, the division of \(r\) by \(d\) is not effective (or even is trivial), but \(q,r\) is not a quotient-remainder pair. For example consider \(p(x)=\max(3x,-3x),d(x)=\max(2x,-2x)\) and \(q(x)=\max(0,x),r(x)=\max(0,-3x)\). It is easy to see that the division of \(r\) by \(d\) is trivial.

Consider the set

\[C=\{\mathbf{c}\in\mathbb{R}^{n}:\text{Newt}(\mathbf{c}^{T}\mathbf{x}+d(\mathbf{x}))\subset \text{Newt}(p(\mathbf{x}))\}. \tag{9}\]

This set appears also [9, 10, 11] for the case of discrete coefficients. It is not difficult to see that \(C\) is a convex polytope. The following proposition shows that the division is non-trivial if and only if \(C\neq 0\).

_Proposition 4:_ Let \(C\) be the set defined in (9). Then,

1. \(C\) is equal to the domain of function \(l(\cdot)\) (given by (7)). That is, \[C=\text{\bf dom}(l(\mathbf{a}))=\{\mathbf{a}\in\mathbb{R}^{n}:l(\mathbf{a})>-\infty\}.\]
2. The division is non-trivial if and only if \(C\neq\emptyset\)
3. It the quotient has the form \(q(\mathbf{x})=\bigvee_{i=1}^{m_{p}}(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+\hat{b}_{i})\), then \(\hat{\mathbf{a}}_{i}\in C\), for all \(i\).

_Proof_ See Appendix D.

_Corollary 1:_ Assume that the division of a tropical polynomial \(p(\mathbf{x})=\bigvee_{i=1}^{m_{p}}(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+b_{i})\) by another tropical polynomial \(d(\mathbf{x})=\bigvee_{i=1}^{m_{d}}(\tilde{\mathbf{a}}_{i}^{T}\mathbf{x}+\tilde{b}_{i})\) is nontrivial. Then

1. It holds \[\text{span}\{\tilde{\mathbf{a}}_{1},\dots,\tilde{\mathbf{a}}_{m_{d}}\}\subset\text{ span}\{\mathbf{a}_{1},\dots,\mathbf{a}_{m_{p}}\}.\]
2. The quotient \(q(\mathbf{x})=\bigvee_{i=1}^{m_{q}}(\tilde{\mathbf{a}}_{i}^{T}\mathbf{x}+\tilde{b}_{i})\) is such that \[\text{span}\{\hat{\mathbf{a}}_{1},\dots,\hat{\mathbf{a}}_{m_{q}}\}\subset\text{ span}\{\mathbf{a}_{1},\dots,\mathbf{a}_{m_{p}}\}.\]

_Proof:_ See Appendix E.

_Remark 3:_ Corollary 1 can be used to simplify the division of two tropical polynomials. Particularly, assume that

\[\text{span}\{\tilde{\mathbf{a}}_{1},\dots,\tilde{\mathbf{a}}_{m_{d}}\}\subset\text{ span}\{\mathbf{a}_{1},\dots,\mathbf{a}_{m_{p}}\}\subset\mathbb{R}^{d}.\]

If \(\mathbf{Q}\) is a matrix the columns of which represent an orthonormal basis of the subspace \(\text{span}\{\mathbf{a}_{1},\dots,\mathbf{a}_{m_{p}}\}\), then \(\mathbf{a}_{i}\) can be expressed in terms of a reduced dimension vector \(\mathbf{a}_{i}^{T}\) as \(\mathbf{a}_{i}=\mathbf{Q}\mathbf{a}_{i}^{T}\), and \(\tilde{\mathbf{a}}_{i}\) similarly in terms of \(\tilde{\mathbf{a}}_{i}^{T}\) as \(\tilde{\mathbf{a}}_{i}=\tilde{\mathbf{Q}}\tilde{\mathbf{a}}_{i}^{T}\). Then, the division of polynomials \(p\) and \(d\) can be expressed in terms of reduced dimension polynomials \(p^{r},d^{r}\) as \(p(\mathbf{x})=p^{r}(\mathbf{Q}^{T}\mathbf{x})=p^{r}(\tilde{\mathbf{x}}),d=d^{r}(\mathbf{Q}^{T}\mathbf{ x})=d^{r}(\tilde{\mathbf{x}})\), where the dimension of \(\tilde{\mathbf{x}}\) is equal to the rank of \(\mathbf{Q}\).

This observation is certainly useful when the number of terms of the dividend \(m_{p}\) is less than the space dimension \(d\). It will be also used in Section VI-B.

## IV An Exact Algorithm For Tropical Division

We now present an algorithm to compute the quotient and remainder of a division. The algorithm uses polyhedral computations. Particularly, we use (8), and its equivalent in terms of the epigraphs (see (34) in the appendix) to compute \(q\).

The input of the algorithm is a pair of tropical polynomials

\[p(\mathbf{x})=\bigvee_{i=1}^{m_{p}}(\mathbf{a}_{i}^{T}\mathbf{x}+b_{i}),\quad d(\mathbf{x})= \bigvee_{j=1}^{m_{d}}(\tilde{\mathbf{a}}_{j}^{T}\mathbf{x}+\tilde{b}_{j}).\]

The first step is to compute a partition of \(\mathbb{R}^{n}\) in polyhedra on which \(f(\mathbf{x})=p(\mathbf{x})-d(\mathbf{x})\) is linear. To do so, for each \(i,j\), consider the polyhedron on which the terms \(i,j\) attain the maximum in \(p\) and \(d\) respectively.

\[P_{i,j}=\{x\in\mathbb{R}^{n}:\mathbf{a}_{i}^{T}\mathbf{x}+b_{i}\geq\mathbf{a }_{i}^{T}\mathbf{x}+b_{i^{\prime}},\] \[\tilde{\mathbf{a}}_{j}^{T}\mathbf{x}+\tilde{b}_{j}\geq\tilde{\mathbf{a}}_{j^ {\prime}}^{T}\mathbf{x}+\tilde{b}_{j^{\prime}},i^{\prime}=1,\dots,m_{p},j^{\prime }=1,\dots,m_{d}\}.\]

Polyhedron \(P_{i,j}\) can be compactly written in terms of a matrix \(\mathbf{A}^{i,j}\) and a vector \(\mathbf{b}^{i,j}\) as

\[P_{i,j}=\{\mathbf{x}\in\mathbb{R}^{n}:\mathbf{A}^{i,j}\mathbf{x}\geq\mathbf{b}^{i,j}\}.\]

Note that some \(P_{i,j}\)'s can be empty.

The second step is to compute the polyhedra that represent the epigraph of \(f(\mathbf{x})=p(\mathbf{x})-d(\mathbf{x})\) for each \(P_{i,j}\). The epigraphs are given by:

\[E_{i,j}=\{(\mathbf{x},z)\in\mathbb{R}^{n+1}:\mathbf{A}^{i,j}\mathbf{x}\geq \mathbf{b}^{i,j},\] \[\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad z \geq(\mathbf{a}_{i}-\tilde{\mathbf{a}}_{j})^{T}\mathbf{x}+b_{i}-\tilde{b}_{j}\}.\]

For each polyhedron \(E_{i,j}\), we compute the corresponding \(\mathcal{V}-\)representation, consisting of a set of vertices \(V_{i,j}\) and a set of rays \(R_{i,j}\), satisfying (4). Note that the epigraph of \(f\) is given by \(\cup_{i,j}E_{i,j}\).

The epigraph of the quotient \(q\) is given by the closed convex hull of the epigraph of \(f\) (see (34) in the appendix). Next, we consider the union of vertices \(V=\cup_{i,j}V_{i,j}\), and rays \(R=\cup_{i,j}R_{i,j}\), and consider the polyhedron \(E\) generated by \(V\) and \(R\) (as in (4)). Due to Proposition 2, \(E\) is the epigraph of the quotient \(q(\mathbf{x})\). Then, compute the \(\mathcal{H}-\)representation of \(E\)

\[E=\{(\mathbf{x},z)\in\mathbb{R}^{n+1}:[\mathbf{A}^{E,\mathbf{x}}\ \mathbf{a}^{E,z}][\mathbf{x}^{T}\ z]^{T}\geq\mathbf{b}^{E}\}, \tag{10}\]

for appropriate matrix \(\mathbf{A}^{E,\mathbf{x}}\) and vectors \(\mathbf{a}^{E,z},\mathbf{b}^{E}\). Assume that \(E\neq\mathbb{R}^{n}\). We will prove in Proposition 5 that the components of \(\mathbf{a}^{E,z}\) are positive. Thus,

\[E= \left\{(\mathbf{x},z)\in\mathbb{R}^{n+1}:z\geq-\frac{1}{[\mathbf{a}^{E,z}] }[\mathbf{A}^{E,\mathbf{x}}]_{l}\mathbf{x}+\right.\] \[\quad\quad\quad\left.+[\mathbf{b}^{E}]_{l}/[\mathbf{a}^{E,z}]_{l},\ l=1, \dots,L\right\},\]

where \(L\) is the number of rows of \(\mathbf{A}^{E,\mathbf{x}}\). Therefore,

\[q(\mathbf{x})=\bigvee_{l=1}^{L}\left(-\frac{1}{[\mathbf{a}^{E,z}]_{l}}[\mathbf{A}^{E,\mathbf{x} }]
The procedure is illustrated in Figure 3.

Having computed the quotient, we formulate the sum \(\tilde{p}(\mathbf{x})=d(\mathbf{x})+q(\mathbf{x})\). For each \(P_{i,j}\) we choose a point \(\mathbf{x}_{i,j}\) in its relative interior, for example

\[\mathbf{x}_{i,j}=\frac{1}{|V_{i,j}|}\sum_{v\in V_{i,j}}v+\sum_{w\in R_{i,j}}w. \tag{12}\]

Let \(I\) be the set of indices \(i\) such that there is an index \(j\) satisfying \(p(\mathbf{x}_{i,j})>\tilde{p}(\mathbf{x}_{i,j})\). Then,

\[r(\mathbf{x})=\bigvee_{i\in I}(\mathbf{a}_{i}^{T}\mathbf{x}+b_{i}). \tag{13}\]

The computation of the quotient and the remainder is summarized in Algorithm 1.

```
0: Input: Tropical Polynomials \(p(\mathbf{x})\) and \(d(\mathbf{x})\)
1: Compute a partition of \(\mathbb{R}^{n}\) into polyhedra \(P_{i,j}\), \(i=1,\ldots,m_{p}\), \(j=1,\ldots,m_{d}\), and their interior points \(\mathbf{x}_{i,j}\) according to (12)
2: Compute the epigraphs \(E_{i,j}\)
3: For each epigraph, compute the resolution into a set of vertices \(V_{i,j}\) and a set of rays \(R_{i,j}\), using a representation conversion algorithm
4: Compute the \(\mathcal{H}-\)representation of the polyhedron generated by the union of the set of vertices \(V=\cup_{i,j}V_{i,j}\), and the union of the set of rays \(R=\cup_{i,j}R_{i,j}\), in the form (10)
5: If \(L>0\) then \(q(\mathbf{x})\) is given by (11). Otherwise, \(q(\mathbf{x})=-\infty\).
6: Compute the set of indices \(i\) such that there is a \(j\) with \(p(\mathbf{x}_{i,j})-d(\mathbf{x}_{i,j})-q(\mathbf{x}_{i,j})>0\). Compute \(r(\mathbf{x})\) according to (13)
7: Return \(q(\mathbf{x})\), \(r(\mathbf{x})\)
```

**Algorithm 1** Exact Polynomial Division

**Proposition 5**: _If \(E\neq\mathbb{R}^{n}\), then the output of the algorithm is indeed the quotient and the remainder of the division. If \(E=\mathbb{R}^{n}\), then the quotient is \(q(\mathbf{x})=-\infty\), for all \(x\in\mathbb{R}^{n}\) and the remainder is equal to the dividend, i.e., \(r=p\)._

_Proof_ See Appendix F.

_Remark 4:_ The algorithm involves first computing the resolution of a number of polyhedra and then the computation of the \(\mathcal{H}-\)representation of the convex hull of their union. The computational complexity of these operations is exponential on their input size (e.g. [48]). Therefore, we expect that the algorithm will be usable only for small examples.

_Example 2:_ Let us divide the tropical polynomial \(p(x,y)=\max(0,3x+3y,6x)\), by \(d(x,y)=\max(x,x+y,2x+y)\), using Algorithm 1. The quotient is

\[q(x,y)=\max(1.5x+1.5y,3x,0),\]

and \(d(x)+q(x)\) is given by \(\max(x,x+y,2.5x+2.5y,3.5x+2.5y,5x+y,4x)\).

The linearity regions of \(f(x,y)\) and \(q(x)\) are shown in Figure 4. Furthermore, Figure 4 illustrates the Newton polygons for \(p,d\) and \(d+q\).

## V Approximate Algorithms For Tropical Division

We now present an approximate algorithm for computing the quotient of a tropical division, assuming that the quotient has a predefined maximum number of terms \(\tilde{m}_{q}\). The algorithm mimics the procedure developed in [33].

The quotient has the form \(q(\mathbf{x})=\bigvee_{i=1}^{m_{q}}(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+\hat{b}_{i})\). From the proof of Proposition 2 we know that \(q(\mathbf{x})\) is the maximum tropical polynomial function, all the terms of which satisfy \(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+\hat{b}_{i}\leq f(\mathbf{x})\), for all \(\mathbf{x}\in\mathbb{R}^{n}\). We will utilize this property to formulate an optimization problem on a set of sample points \(\mathbf{x}_{1},\ldots,\mathbf{x}_{N}\in\mathbb{R}^{n}\). Given the set of sample points and the corresponding set of the values of function \(f(\cdot)\), i.e., \(f(\mathbf{x}_{1}),\ldots,f(\mathbf{x}_{N})\), the optimization problem is formulated as:

\[\underset{\hat{\mathbf{a}}_{i},\hat{b}_{i},i=1,\ldots,\hat{m}_{q}}{ \text{maximize}} \sum_{j=1}^{N}\bigvee_{i=1}^{\tilde{m}_{q}}(\hat{\mathbf{a}}_{i}^{T} \mathbf{x}_{j}+\hat{b}_{i})\] (14) subject to \[\hat{\mathbf{a}}_{i}^{T}\mathbf{x}_{j}+\hat{b}_{i}\leq f(\mathbf{x}_{j}), \forall i,j\] \[\hat{\mathbf{a}}_{i}\in C,\ \ i=1,\ldots,\tilde{m}_{q}\]

where \(C\) is the set defined in (9). Since \(C\) is a polytope, the last constraint in (16) can be written as a set of linear inequalities. The details are presented in subsection V-A.

The algorithm alternates between two phases. Phase 1 partitions the data. Particularly, assuming a set of solutions

Fig. 3: _Illustration of the algorithm applied to Example 1. Function \(f(x)=p(x)-d(x)\) is linear on each of the subsets \(P_{1},\ldots,P_{b}\). The epigraphs are generated as follows: \(E_{1}\) by vertex \(V_{1}\) and rays \(R_{1}\) and \(R_{2}\), epigraph \(E_{2}\) by vertices \(V_{1},V_{2}\) and ray \(R_{2}\), epigraph \(E_{3}\) by vertices \(V_{2},V_{3}\) and ray \(R_{2}\), epigraph \(E_{4}\) by vertices \(V_{3},V_{4}\) and rays \(R_{2}\), and finally epigraph \(E_{5}\) by vertex \(V_{4}\) and rays \(R_{2},R_{3}\). Their closed convex hull \(E\) is generated by vertices \(V_{1},V_{2},V_{4}\) and rays \(R_{1},R_{3}\). The quotient \(q(\mathbf{x})\) is computed by \(\mathcal{H}-\)representation of \(E\)._

Fig. 4: _The linear regions of \(f(x,y)\) and \(q(x,y)\). The two functions coincide on the bold half-lines. The figure also illustrates the Newton polygons of \(p(x,y)\), \(d(x,y)\), and \(d(x,y)+q(x,y)\)._
\((\hat{\mathbf{a}}_{i},b_{i})\), \(i=1,\ldots,\tilde{m}_{q}\), we partition the samples \(\mathbf{x}_{j}\) into sets \(I_{1},\ldots,I_{\tilde{m}_{q}}\) with \(j\in I_{i}\), if

\[\hat{\mathbf{a}}_{i}^{T}\mathbf{x}_{j}+\hat{b}_{i}\geq\hat{\mathbf{a}}_{i^{\prime}}^{T}\mathbf{x} _{j}+\hat{b}_{i^{\prime}},\ \ \text{for all }i^{\prime}. \tag{15}\]

Phase 2, solves the linear optimization problem:

\[\underset{\hat{\mathbf{a}}_{i},b_{i},i=1,\ldots,\tilde{m}_{q}}{\text{ maximize}} \sum_{i=1}^{\tilde{m}_{q}}\hat{\mathbf{a}}_{i}^{T}\mathbf{x}_{j}+\hat{b}_{i}\] (16) subject to \[\hat{\mathbf{a}}_{i}^{T}\mathbf{x}_{j}+\hat{b}_{i}\leq f(\mathbf{x}_{j}),\ \ \forall i,j\] \[\hat{\mathbf{a}}_{i}\in C,\ \ i=1,\ldots,\tilde{m}_{q}\]

We then simplify (16) in two ways. First, it may contain redundant inequalities. To remove the redundant inequalities we compute the lower convex hull of the set \(\{(\mathbf{x}_{j},f(\mathbf{x}_{j})):j=1,\ldots,N\}\) with respect to its last component, and denote by \(J\) the set of indices \(j\), for which \((\mathbf{x}_{j},f(\mathbf{x}_{j})\) belongs to the lower convex hull.

Second, problem (16) decomposes into a set of \(\tilde{m}_{q}\) linear programs

\[\underset{\hat{\mathbf{a}}_{i},b_{i}}{\text{maximize}} \mathbf{s}_{i}^{T}\hat{\mathbf{a}}_{i}+N_{i}\hat{b}_{i}\] (17) subject to \[(\mathbf{x}_{j})^{T}\ 1]\begin{bmatrix}\hat{\mathbf{a}}_{i}\\ \hat{b}_{i}\end{bmatrix}\leq f(\mathbf{x}_{j}),\ \ j\in J\ \,\] \[\hat{\mathbf{a}}_{i}\in C\]

where \(\mathbf{s}_{i}=\sum_{j\in I_{i}}\mathbf{x}_{j}\), \(N_{i}=|I_{i}|\).

The computations are summarized in Algorithm 2.

```
0: Input: Tropical Polynomials \(p(\mathbf{x})\), \(d(\mathbf{x})\), the degree of the approximate quotient \(\tilde{m}_{q}\), the number of sample points \(N\), the maximum number of iterations
1: Compute the inequalities corresponding to set \(C\) using the formulas in Appendix V-A
2: Sample the points \(\mathbf{x}_{j}\), for \(j=1,\ldots,N\), and compute the values of \(f(\mathbf{x}_{j})=p(\mathbf{x}_{j})-d(\mathbf{x}_{j})\)
3: Compute the lower convex hull of the set \(\{(\mathbf{x}_{j},f(\mathbf{x}_{j})):j=1,\ldots,N\}\), and use the samples belonging to it to formulate the first constraint in (16)
4: Initialize the partition \(I_{1},\ldots,I_{\tilde{m}_{q}}\)
5: Solve Problems (17), for all \(i\)
6: Compute the partitions \(I_{i}\) that satisfy (15)
7: If the number of iterations has not exceeded the maximum number of iterations go to Step 5
8: Return \(\hat{q}(\mathbf{x})=\sqrt{\hat{m}_{i=1}^{\tilde{m}_{q}}}\hat{\mathbf{a}}_{i}^{T}\mathbf{x }+\hat{b}_{i}\),
```

**Algorithm 2** Approximate Polynomial Division

We then examine the quality of approximation of \(p(\mathbf{x})\) by \(d(\mathbf{x})+q_{t}(\mathbf{x})\) on the sample set given by

\[e(t)=\sum_{j=1}^{N}[p(\mathbf{x}_{j})-d(\mathbf{x}_{j})-q_{t}(\mathbf{x}_{j})], \tag{18}\]

where \(q_{t}\) is the quotient computed after \(t\) steps of the algorithm. The following proposition shows that the quality of approximation improves as the number of steps increases.

**Proposition 6**: _The quantity \(e(t)\) is nonnegative and non-increasing with respect to \(t\)._

Proof:: See Appendix G.

**Remark 5**: _In practice some issues with Algorithm 1 may arise. Some of the classes \(I_{i}\) may end up empty. In this case, we may split another existing class randomly. Furthermore, we expect the algorithm to converge to a local optimum. We may use a multiple start method to avoid bad local optima._

**Remark 6**: _Let us comment about the relationship between Algorithm 2 and [33]. Inspired by [33], Algorithm 2 partitions the data in several clusters. The major difference is that in Algorithm 2 we use the sharper lower linear approximation of each class of data subject to the constraint that this approximation does not exceed any sample point of any class. Thus, instead of alternating between clustering and linear regression, we alternate between clustering and linear programming._

Another related work is [49] which studies the problem of computing a convex under-approximation for a given dataset by a piecewise linear function. The algorithm uses an iterative linear programming technique and has two phases. At phase one of each iteration, it considers a subset of the data set and minimizes the maximum approximation error for this subset, using a piecewise linear function with a number of terms equal to the iteration count. At phase two, it adds to this subset the data point that maximizes the approximation error. Then, it moves to the next iteration.

### _The Linear Constraints for Set \(C\)_

We next compute a set of linear inequalities describing set \(C\). Let the extreme points of \(\text{Newt}(p)\) be \(\mathbf{a}_{e,1},\ldots,\mathbf{a}_{e,k}\). Then, \(\hat{\mathbf{a}}\in C\) if and only if \(\hat{\mathbf{a}}+\text{Newt}(d(\mathbf{x}))\subset\text{Newt}(p(\mathbf{x}))\). That is,

\[\hat{\mathbf{a}}+\tilde{\mathbf{a}}_{j}\in\text{conv}\{\mathbf{a}_{e,1},\ldots,\mathbf{a}_{e,k}\},\]

for all \(j\). Equivalently if there are auxiliary variables \(\lambda_{j,l}\), such that

\[\hat{\mathbf{a}}+\tilde{\mathbf{a}}_{j}=\sum_{l=1}^{k}\mathbf{a}_{e,l}\lambda_{j,l},\ \ \ \lambda_{j,l}\geq 0,\ \ \ \ \sum_{l=1}^{k}\lambda_{j,l}=1\]

for all \(j=1,\ldots,m_{q}\), \(l=1,\ldots,k\).

Problem (17) becomes:

\[\underset{\hat{\mathbf{a}}_{i},b_{i},\lambda}{\text{maximize}} \mathbf{s}_{i}^{T}\hat{\mathbf{a}}_{i}+N_{i}\hat{b}_{i}\] (19) subject to \[[\mathbf{x}_{j}^{T}\ 1]\begin{bmatrix}\hat{\mathbf{a}}_{i}\\ \hat{b}_{i}\end{bmatrix}\leq f(\mathbf{x}_{j}),\ \ j\in J\] \[\hat{\mathbf{a}}_{i}+\tilde{\mathbf{a}}_{j}=\sum_{l=1}^{k}\mathbf{a}_{e,l} \lambda_{j,l},\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \lambda_{j,l}\geq 0\] \[\sum_{l=1}^{k}\lambda_{j,l}=1\]

## VI Division for Composite Polynomials

In this section, we present some analytical results and some algorithms for dividing tropical polynomials that are sums or maxima of simpler ones.
### _General Properties_

We will use \(Q(p,d)\) and \(R(p,d)\) to denote the quotient and remainder of the division of a tropical polynomial \(p\) by another tropical polynomial \(d\).

**Proposition 7**: _The following inequalities hold_

\[Q(p_{1}+p_{2},d) \geq Q(p_{1},d)+Q(p_{2},d)+d, \tag{20}\] \[Q(p_{1}\lor p_{2},d) \geq Q(p_{1},d)\lor Q(p_{2},d),\] (21) \[Q(p,d_{1}+d_{2}) \geq Q(p,d_{1})+Q(p,d_{2})-p,\] (22) \[Q(p,d_{1}\lor d_{2}) \geq Q(p,d_{1})\wedge Q(p,d_{2}). \tag{23}\]

_Furthermore, if \(R(p_{1},d)=R(p_{2},d)=-\infty\), then (20),(21),(22) hold as equalities. Finally, if \(s\) is a monomial, i.e., it has the form \(s(\mathbf{x})=\mathbf{a}^{T}\mathbf{x}+b\), then_

\[Q(p+s,d)=Q(p,d)+s=Q(p,d-s), \tag{24}\]

_Proof_ See Appendix H.

**Remark 7**: _Assume that polynomial \(d\) is fixed. Inequality (21) implies that function \(Q(\cdot,d)\) is nondecreasing. If \(d(\mathbf{x})\geq 0\), for all \(\mathbf{x}\in\mathbb{R}^{n}\), then function \(Q(\cdot,d)\) is also super-additive._

In the following we call a tropical polynomial described as a summation of other tropical polynomials as composite. Furthermore, we refer to a tropical polynomial in the form 1 as simple1.

Footnote 1: These notions refer to the representation of the tropical polynomials, and not to whether or not they can be written as a sum of other polynomials. That is, it is possible that a tropical polynomial with a simple representation (in the form (1)), can be factorized and written in a composite form.

### _Dividing a Composite Polynomial by a Simple_

In the following we propose algorithms to approximately divide a composite polynomial \(p(\mathbf{x})=\sum_{\nu=1}^{N}p^{\nu}(\mathbf{x})\) or \(p(\mathbf{x})=\bigvee_{\nu=1}^{N}p^{\nu}(\mathbf{x})\) by another polynomial \(d(\mathbf{x})\). The computation of the value of the composite polynomial \(p(\mathbf{x})\) is trivial, provided the values of the simpler polynomials \(p^{\nu}(\mathbf{x})\). Thus, the only difficulty in order to apply Algorithm 2 is the computation of the Newton polytope \(\text{Newt}(p)\).

It is convenient to consider the Newton polytopes \(\text{Newt}(p^{\nu})\) in terms of an extended description. Particularly, assume matrices \(\mathbf{A}^{\nu},\mathbf{F}^{\nu}\), and vectors \(\mathbf{\beta}^{\nu}\) of appropriate dimensions, such that

\[\text{Newt}(p^{\nu})=\{\mathbf{a}\in\mathbb{R}^{n}:\mathbf{a}=\mathbf{F}^{\nu}\mathbf{\alpha} ^{\nu},\mathbf{A}^{\nu}\mathbf{\alpha}^{\nu}\leq\mathbf{\beta}_{\nu}\}. \tag{25}\]

For \(p(\mathbf{x})=\sum_{\nu=1}^{N}p^{\nu}(\mathbf{x})\), the Newton polytope is given by

\[\text{Newt}(p)=\{\mathbf{F}^{1}\mathbf{\alpha}^{1}+\cdots+\mathbf{F}^{N}\mathbf{\alpha}^{N}: \mathbf{A}^{\nu}\mathbf{\alpha}^{\nu}\leq\mathbf{\beta}_{\nu},\nu=1,\ldots,N\},\]

which corresponds to the Minkowski sum \(\text{Newt}(p^{1})\oplus\cdots\oplus\text{Newt}(p^{N})\). For \(p(\mathbf{x})=\bigvee_{\nu=1}^{N}p^{\nu}(\mathbf{x})\), the Newton polytope is given by

\[\text{Newt}(p)= \{\mathbf{F}^{1}\mathbf{\alpha}^{1}+\cdots+\mathbf{F}^{N}\mathbf{\alpha}^{N}: \exists\lambda^{\nu}\geq 0,\nu=1,\ldots,N,\] \[\mathbf{A}^{\nu}\mathbf{\alpha}^{\nu}\leq\mathbf{\beta}_{\nu}\lambda^{\nu}, \ \sum_{\nu=1}^{N}\lambda^{\nu}=1\}.\]

### _Neural Networks As Differences of Composite Polynomials_

We then present some examples of neural networks represented as the difference of two composite tropical polynomials.

**Example 3** (Single hidden layer ReLU network): _This example follows [9]. Consider a neural network with \(n\) inputs, a single hidden layer with \(m_{p}\) neurons and ReLU activation functions, and a single linear output neuron. The output of the \(\nu\)-th neuron of the hidden layer is:_

\[z_{\nu}=\max(\mathbf{w}_{\nu}^{T}\mathbf{x}+\beta_{\nu},0)\]

_The output \(y\) can be written as_

\[y=\sum_{\nu=1}^{N_{1}}w_{\nu}^{2+}z_{\nu}-\sum_{\nu=N_{1}+1}^{N}w_{\nu}^{2-}z_ {\nu}+\beta^{2},\]

_where \(N_{1}+N_{2}=N\), and \(N_{1},N_{2}\) are the number of neurons that have positive weight and negative weights respectively. Without loss of generality, we assume that the neurons are ordered such that the weights of the first neurons to be positive and the weights of the last negative. The output of the neural network can be expressed as the difference of two tropical polynomials \(p_{1}(\mathbf{x})\), \(p_{2}(\mathbf{x})\), plus a constant term. Each of these tropical polynomials is written as the sum of simpler ones_

\[p_{1}(\mathbf{x}) =\sum_{\nu=1}^{N_{1}}\max(w_{\nu}^{2+}\mathbf{w}_{\nu}^{T}\mathbf{x}+w_{ \nu}^{2+}\beta_{\nu},0)\] \[=\sum_{\nu=1}^{N_{1}}\max((\mathbf{a}_{1}^{\nu})^{T}\mathbf{x}+b_{\nu},0)\] \[p_{2}(\mathbf{x}) =\sum_{\nu=N_{1}+1}^{N}\max(w_{\nu}^{2-}\mathbf{w}_{\nu}^{T}\mathbf{x}+w_ {\nu}^{2-}\beta_{\nu},0)\] \[=\sum_{\nu=N_{1}+1}^{N}\max((\mathbf{a}_{2}^{\nu})^{T}\mathbf{x}+b_{\nu},0)\]

_Note that \(p_{1}\), \(p_{2}\), expressed in their canonical form (1), can have a large number of terms, corresponding to linear regions of tropical polynomial functions \(p_{1}\), \(p_{2}\) (for an enumeration of the linear regions see [3, 4])._

The Newton polytope of the simple tropical polynomial \(\max((\mathbf{a}_{1}^{\nu})^{T}\mathbf{x}+b_{\nu},0)\) is the line segment \([0,\mathbf{a}_{1}^{\nu}]\) in the \(d\)-dimensional space. This polytope is equivalently described in the form (25) as \(\text{Newt}(p^{\nu})=\{\alpha^{\nu}\mathbf{a}_{1}^{\nu}:0\leq\alpha^{\nu}\leq 1\}\). Thus, the Newton polytope of \(p_{1}\) is given by

\[\text{Newt}(p_{1})=\{\alpha^{1}\mathbf{a}_{1}^{1}+\cdots+\mathbf{\alpha}^{1}\mathbf{a}_{1}^{ n_{1}}:0\leq\alpha^{\nu}\leq 1\},\]

which is a zonotope (see also [4]).

The constraint \(\mathbf{\hat{a}}_{i}\in C\) of problem (17) becomes

\[\mathbf{\hat{a}}_{i}+\mathbf{\bar{a}}_{j}=\sum_{l=1}^{n}\mathbf{a}^{l}\lambda_{j,l},\ \ 0\leq \lambda_{j,l}\leq 1.\]

Finally, note that if the number of neurons of the hidden layer \(N_{1},N_{2}\) is smaller than the input dimension \(n\), then the tropical polynomials \(p_{1},p_{2}\) can be expressed in dimensions \(N_{1},N_{2}\).

**Example 4** (Multi-class to binary classification with a ReLU network): _Assume a neural network with a single ReLU hidden
layer with \(N\) units and a softmax output layer with \(K\) units trained to classify samples to \(K\) classes. The equations are the following

\[z_{\nu}^{1}=\max((\mathbf{w}_{\nu}^{1})^{T}x+b_{\nu}^{1},0),\ \ z_{k}^{2}=(\mathbf{w}_{k}^{2})^{T}\mathbf{z}^{1},\ \ y_{k}=\frac{e^{z_{k}^{2}}}{\sum_{k^{ \prime}}e^{z_{k^{\prime}}^{2}}}.\]

We then assume, that for a given set of samples, we know that the correct class is either \(i_{1}\) or \(i_{2}\). The neural network can be simply reduced for binary classification, substituting the last layer with a single neuron with sigmoid activation with output given by

\[y=\sigma((\mathbf{w}_{i_{1}}^{2}-\mathbf{w}_{i_{2}}^{2})^{T}\mathbf{z}^{2}+\beta^{2}),\]

where \(\sigma(\cdot)\) is the sigmoid activation function,, i.e. \(\sigma(z)=1/(1-e^{-z})\). Note that this transformation corresponds to the application of Bayes' rule with prior probability of \(1/2\) for classes \(i_{1},i_{2}\), and likelihoods given by the output of the initial neural network.

The output of the new network can be represented as a difference of two tropical polynomials as in the previous example.

**Example 5** (Representing a multi-class network as a vector of tropical polynomials): _Consider a neural network with a single-hidden layer having ReLU activations and a linear output layer with \(K\) neurons, followed by softmax. Each output (before the softmax) \(z_{k}^{2}\) can be expressed as a difference of two tropical polynomials plus a constant term_

\[z_{k}^{2}=p_{k}^{+}(\mathbf{x})-p_{k}^{-}(\mathbf{x})+b_{k}.\]

_The output is given by_

\[y_{k}=\frac{e^{z_{k}^{2}}}{\sum_{k^{\prime}}e^{z_{k^{\prime}}^{2}}}.\]

_Let us add to each of \(z_{k}^{2}\)'s the sum of all negative polynomials. We get the tropical polynomials_

\[\tilde{z}_{k}^{2}=p_{k}^{+}(\mathbf{x})+\sum_{k^{\prime}\neq k}p_{k^{\prime}}^{-}( \mathbf{x})+b_{k}. \tag{26}\]

_Then, it is easy to see that_

\[y_{k}=\frac{e^{\tilde{z}_{k}^{2}}}{\sum_{k^{\prime}}e^{\tilde{z}_{k^{\prime}} ^{2}}}.\]

_Finally, denoting by \(z_{i}^{1}\), \(i=1,...M\) the output of the hidden layer before the ReLU, then \(\tilde{z}_{k}^{2}\) can be written in the form_

\[\tilde{z}_{k}^{2}=\sum_{i=1}^{M}\bar{w}_{k,i}\max(z_{i}^{1},0), \tag{27}\]

_for appropriate non-negative constants \(\bar{w}_{k,i}\)._

**Remark 8**: _In the tropical framework, the difference of two polynomials is considered as a tropical rational function. The addition of all the negative polynomials (denominators) corresponds to making the tropical fractions have the same denominator._

### _Division Algorithm for Composite Quotient_

We then describe an algorithm to divide a composite tropical polynomial \(p(\mathbf{x})=\sum_{i=1}^{n}\max(\mathbf{a}_{i}^{T}\mathbf{x}+b_{i},0)\), where \(n\) is the dimension of the underlying space, by the zero polynomial \(d(\mathbf{x})=0\) and search for a quotient in the form \(q(\mathbf{x})=\sum_{i=1}^{m_{q}}\max(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+\hat{b}_{i},0)\). Such polynomials occur in single hidden layer ReLU neural networks if the number of neurons in the hidden layer is smaller than the space dimension, using a dimensionality reduction transformation (e.g., QR).

We assume that the vectors \(\mathbf{a}_{i}\) are linearly independent. Denote by \(\mathbf{A}\) the square matrix of vectors \(\mathbf{a}_{i}\), i.e., \(\mathbf{A}=[\mathbf{a}_{1}\ \ldots\mathbf{a}_{n}]\), and by \(\mathbf{B}\) the vector of \(b_{i}\), i.e., \(\mathbf{B}=[b_{1}\ \ldots\ b_{n}]^{T}\).

The algorithm is based on the fact that the quotient \(q\) is the largest tropical polynomial function which is less than or equal to \(p\). Consider a set of sample points \(x^{1},\ldots,x^{N}\). Then the problem of maximizing \(q\) can be approximated by

\[\underset{\hat{\mathbf{a}}_{i},\hat{b}_{i},i=1,\ldots,m_{i_{q}}}{ \text{maximize}} \sum_{j=1}^{N}\sum_{i=1}^{m_{q}}\max(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}^{j }+\hat{b}_{i},0)\] (28) subject to \[q(\mathbf{x})\leq p(\mathbf{x}),\ \ \text{for all}\ x\in\mathbb{R}^{n}\]

The following proposition reformulates the constraint \(q(\mathbf{x})\leq p(\mathbf{x}),\ \ \text{for all}\ x\in\mathbb{R}^{n}\) in terms of the variables \(\hat{\mathbf{a}}_{i},\hat{b}_{i},i=1,\ldots,m_{q}\).

**Proposition 8**: _Let \(p(\mathbf{x}),q(\mathbf{x})\) as above. Then, \(q(\mathbf{x})\leq p(\mathbf{x})\), for all \(x\in\mathbb{R}^{n}\) if and only if_

\[\mathbf{0}\leq\mathbf{A}^{-1}\hat{\mathbf{a}}_{i},\ \ \mathbf{A}^{-1}\sum_{i=1}^{m_{n_{q}}}\hat{\mathbf{a}} _{i}\leq\mathbf{1},\ \ \hat{b}_{i}\leq\mathbf{B}^{T}\mathbf{A}^{-1}\hat{\mathbf{a}}_{i} \tag{29}\]

_for \(i=1,\ldots,m_{q}\), where \(\mathbf{0}\) is the vector having all its entries equal to zero and \(\mathbf{1}\) all entries equal to one._

**Phase 1**: _For the set of sample points \(\mathbf{x}^{1},\ldots,\mathbf{x}^{N}\) we compute \(y_{i,j}=\chi_{\hat{\mathbf{a}},\mathbf{x}^{j}+b_{i}\geq 0}\), where \(\chi\) is the indicator function (i.e., \(\chi_{\hat{\mathbf{a}},\mathbf{x}^{j}+b_{i}\geq 0}=1\) if \(\hat{\mathbf{a}}_{i}\mathbf{x}^{j}+b_{i}\geq 0\) and \(0\) otherwise). Based on the values of \(y_{i,j}\), we determine the local form of the objective in (28). Then, we compute \(\mathbf{c}_{1}\) and \(\mathbf{c}_{2}\) as_

\[\big{[}\sum_{j=1}^{N}y_{1,j}\mathbf{x}_{j}^{T}\ \ldots\ \sum_{j=1}^{N}y_{m_{n_{q}},j}\mathbf{x}_{j}^{T} \big{]},\ \big{[}\sum_{j=1}^{N}y_{1,j}\ \ldots\ \sum_{j=1}^{N}y_{m_{n_{q}},j}\big{]}.\]

**Phase 2**: _We solve the linear programming problem_

\[\underset{\hat{\mathbf{a}}_{i},\hat{b}_{i},i=1,\ldots,m_{i_{q}}}{ \text{maximize}} \mathbf{c}_{1}\begin{bmatrix}\hat{\mathbf{a}}_{1}\\ \vdots\\ \hat{\mathbf{a}}_{m_{n_{q}}}\end{bmatrix}+\mathbf{c}_{2}\begin{bmatrix}\hat{b}_{1}\\ \vdots\\ \hat{b}_{m_{n_{q}}}\end{bmatrix}\] (30) subject to \[\mathbf{0}\leq\mathbf{A}^{-1}\hat{\mathbf{a}}_{i},\ \ \ \hat{b}_{i}\leq\mathbf{B}^{T}\mathbf{A}^{-1}\hat{\mathbf{a}}_{i}\] \[\mathbf{A}^{-1}\sum_{i=1}^{m_{q}}\hat{\mathbf{a}}_{i}\leq\mathbf{1},\]

_and denote its solution by \(\hat{\mathbf{a}}_{i}^{\text{new}},b_{i}^{\text{new}}\). We set_

\[\hat{\mathbf{a}}_{i}\leftarrow(1-\rho)\hat{\mathbf{a}}_{i}+\hat{\mathbf{a}}_{i}^{\text{new}}\rho, \tag{31}\] \[\hat{b}_{i}\leftarrow(1-\rho)b_{i}+\rho b_{i}^{\text{new}},\]

_where \(\rho\in(0,1)\) and go to Phase 1._

**Remark 9**: _The algorithm corresponds to the application of Frank-Wolfe (conditional gradient) method (see, e.g., [50]) to (28)._
Let us note that this linear programming problem (30) is analytically solvable. Indeed, observing that the coefficients of \(\hat{b}_{i}\)'s are nonnegative, we can use the second inequality to compute the the optimal value of \(\hat{b}_{i}\)'s. Substituting back to the problem, we end up with \(n\) independent linear programs. The feasible set of each of these problems has \(\tilde{m}_{q}+1\) vertices. We obtain the optimal solution of the problem comparing the value of the linear objective on these vertices.

**Remark 10**: _Let us comment on the difference of the algorithm of this subsection with the algorithm of Subsection VI-B. The algorithm of this subsection provides a solution a composite tropical polynomial, which corresponds to a (compressed) neural network with ReLu activations, while the algorithm of Subsection VI-B provides maxout solutions._

### _Dividing a Vector of Composite Polynomials_

In this subsection, we propose a simplified algorithm for dividing a vector of tropical polynomials by the zero polynomial. The motivation comes from Example 5. Consider a set of tropical polynomials in the form (27).

We then describe a simplified version of the the linear optimization part of Algorithm 2. Consider a set of sample points \(\boldsymbol{z}^{1}(1),\ldots,\boldsymbol{z}^{1}(N)\). Observe that, with this formulation, both constraints (17) are satisfied by any \(\hat{\boldsymbol{a}}_{k}\), such that its \(i-\)th component \(\hat{a}_{i,k}\) satisfies \(0\leq\hat{a}_{i,k}\leq\bar{w}_{i,k}\) and \(\hat{b}_{k}=0\). Then, the solution of (17), for the \(l-th\) term of the division of the \(k-\)th polynomial is approximated by the optimal solution of

\[\begin{split}\underset{\hat{\boldsymbol{a}}_{k}^{l}}{\text{maximize}}& \boldsymbol{s}_{l,k}^{T}\hat{\boldsymbol{a}}_{k}^{l}\\ \text{subject to}& 0\leq\hat{a}_{i,k}^{l}\leq\bar{w}_{i,k}, \;\;i=1,\ldots,M\end{split}. \tag{32}\]

The solution of (32) is given by:

\[\hat{a}_{i,k}^{l}=\begin{cases}\bar{w}_{i,k}&\text{if}\;\;\;s_{i,l,k}>0\\ 0&\text{otherwise}\end{cases}\]

## VII Numerical Results

This section, presents some numerical examples for the tropical division algorithms. First, we present some simple examples to build some intuition for the behaviour of Algorithm 2. Then, some applications of tropical division to the compression of neural networks are presented. We focus on MNIST handwritten and CIFAR-10 datasets.

### _Numerical Examples for Algorithm 2_

As a first example, we present the application of Algorithm 2 to the tropical polynomial division of Example 2. We choose \(200\) sample points distributed according to the normal distribution with unit variance \(\mathcal{N}(0,I_{2})\). The algorithm converges in two steps and gives an approximate solution2.

Footnote 2: The code for the numerical experiments is available online at [https://github.com/jkordonis/TropicalML](https://github.com/jkordonis/TropicalML).

\[\hat{q}(x,y)= \max(1.5038x+1.4962y+0.0182,\] \[0.0003x+0.0288,3x+0.0241),\]

which is very close to the actual quotient (see Example 2). Running again the algorithm with larger sample sizes we conclude that increasing the number of sample points reduces the error.

As a second example we divide a random tropical polynomial with \(m_{p}=128\) terms in \(n=3\) dimensions with another having \(m_{d}=2\) terms. First, we run Algorithm 2, for \(\tilde{m}_{q}=5\), with multiple initial partitions. The evolution of the error for the several runs is shown in Figure 5. The results illustrate the need for a multi-start method. Indeed, for different random initial partitions the algorithm converges to different local minima. The quality of the approximation \(e\) after 10 steps of the algorithm is given in Figure 6 for a varying number of terms \(\tilde{m}_{q}\). This result shows that we may derive a good approximation of the division using a small number of terms (e.g. \(8\)).

### _MNIST Dataset_

This example considers the MNIST handwritten digits data set. It consists of \(60000\) training examples and \(10000\) test examples of \(28\times 28\) gray scale images, that represent digits \(0-9\). We start with a single hidden layer neural network with 100 hidden units with ReLU activations and \(10\) softmax output units. The network is trained in the original training data set using standard techniques (Adam optimizer on cross-entropy loss, batch size of 128). We then design a smaller network discriminating between digits \(3\) and \(5\). The input space has \(28\cdot 28=784\) dimensions.

We first use the technique described in Example 4 to reduce the output layer of the original network to a single neuron. Then, using the technique described in Example 3, we express

Fig. 5: _The total error \(e\) for executions of Algorithm 2, for different initial partitions._

Fig. 6: _The total error \(e\) for divisions having different number of terms for the approximate quotient_
the output (before the sigmoid) as the difference of two tropical polynomials. These polynomials are expressed as the summation of \(100\) terms in the form \(\max(\mathbf{\alpha}_{l}^{r}\mathbf{x},0)\), where \(l=1,2\). We then apply a reduced QR decomposition to the matrices \(\mathbf{A}_{1}=[\mathbf{\alpha}_{1}^{1}\ \ldots,\ \mathbf{\alpha}_{1}^{100}]\) and \(\mathbf{A}_{2}=[\mathbf{\alpha}_{2}^{1}\ \ldots,\ \mathbf{\alpha}_{2}^{100}]\), writing them as \(\mathbf{A}_{1}=\mathbf{Q}_{1}\mathbf{A}_{1}^{T}\) and \(\mathbf{A}_{2}=\mathbf{Q}_{2}\mathbf{A}_{2}^{T}\). The \(\mathbf{A}_{l}^{r}\) matrices have as columns \(\mathbf{\alpha}_{l}^{r,y}\), i.e., an expression of the vectors \(\mathbf{\alpha}_{l}^{r}\) in the subspace of their span. For each polynomial the input has been transformed as \(\mathbf{x}_{l}=\mathbf{Q}_{l}^{T}\mathbf{x}\) (recall Remark 3).

We divide each of these polynomials with the zero polynomial \(d=0\), applying two different variants of the Algorithm 2. In both cases, we use as sample points the first \(200\) training points. The first variant (modification of Section VI-B), the approximate quotients have the form \(q_{l}(\mathbf{x}_{l})=\bigvee_{i=1}^{m_{q}}(\hat{\mathbf{a}}_{i}^{l}\mathbf{x}_{l}+\hat{b }_{i}^{l})\), \(l=1,2\). Then, the output of the original network can be approximated as

\[y^{\text{app}}=\sigma(q_{1}(\mathbf{Q}_{1}^{T}\mathbf{x})-q_{2}(\mathbf{Q}_{2}^{T}\mathbf{x}) +\beta^{2})\]

The computation of the approximate output corresponds to a neural network an input layer, a hidden layer with two maxout units, a linear layer with a single unit (computing the difference of the two maxout units) and an output layer with a single sigmoid unit.

The second variant, uses approximates the quotients in the form \(\sum_{i=1}^{m_{q}}\max(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+\hat{b}_{i},0)\) and uses the modification of subsection VI-D. The compressed network in this case is simply a smaller ReLU network.

Table I presents the error rate on the test set for the original neural network, as well as the computed maxout and ReLU networks. We compare the reults when each maxout unit has \(3\), \(5\), and \(10\) terms respectively and each of the tropical polynomials (positive and negative) has \(3\), \(5\), and \(10\) ReLU units. It also presents the percentage of parameters that remain after the compression. As a baseline for comparison, we use structured L1 pruning without retraining (see e.g. [51]). That is, this method deletes the neurons having weights with small L1 norm. We compare reduced neural networks with the same number of parameters. We present two results, the binary comparison of digits \(3\) and \(5\) and an average of the \(45=10\cdot 9/2\) different binary comparisons.

The performance of the algorithms is shown in Figure 7. We observe that the algorithm of Section VI-B produce superior results compared to the algorithm of Section VI-D. Furthermore, both algorithms outperform L1 structured pruning, especially in the high compression regime. Let us note that the linear programming problem of Section VI-D is analytically solvable, and thus the corresponding algorithm is faster.

**Remark 11**: _It is worth noting substituting \(\mathbf{x}_{l}\) into \(q_{l}\), we get_

\[q_{1}(\mathbf{Q}_{1}^{T}\mathbf{x})=\bigvee_{i=1}^{\hat{m}_{q}}(\hat{\mathbf{a}}_{i}^{l} Q_{l}^{T}\mathbf{x}+\hat{b}_{i}^{l}).\]

_Thus we do not need to store \(Q_{1},Q_{2}\) but only the vectors \(\hat{\mathbf{a}}_{i}^{l}Q_{l}^{T}\), and the scalars \(\hat{b}_{i}^{l}\), \(\beta^{2}\), for \(l=1,2\) and \(i=1,\ldots,m_{q}\)._

**Remark 12**: _Let us note that for computing the reduced networks, we used only the first 200 samples of the training data._

In the following, we continue exploring the algorithm of Section VI-B to more complex examples.

### _CIFAR-10 Dataset-Binary_

This example concerns CIFAR-10 dataset which contains small (\(32\times 32\)) color images. The training set consists of 50000 images of 10 classes (5000 samples for each class) and the test set of 10000 images (1000 samples per class). We first train a VGG-like network [52], having three blocks consisting of two convolution layers with ReLU activation followed by a \(2\times 2\) max-pooling layer with padding, a dense layer with 1024 neurons, and an output layer with 10 neurons. The convolution layers in the first, second, and third block have 32, 64, and 128 neurons, respectively. The network was trained using standard techniques (Adam optimizer for the cross-entropy loss, batch normalization and dropout for regularization and data augmentation). Let us note that the input of the dense layer has dimension \(4\cdot 4\cdot 128=2048\).

We then design a simplified neural network discriminating between pairs of classes (e.g., "Automobile" and "Truck").

Fig. 7: _The average error rate of the compressed networks for the binary MNIST classification, compared with L1 structured pruning. Both techniques lead to better accuracies compared to L1 pruning. Furthermore, maxout networks clearly outperform both the L1 structured pruning and the ReLU compressed network, especially in the large compression regime._
To do so, we approximate the dense layer with 1024 neurons using the techniques of the previous section. Particularly, using the technique of Example 4 we describe the output of the neural network as the difference of two tropical polynomials in \(2048\) dimensions. These polynomials are then divided by the zero polynomial applying the Algorithm 2, with the modification of Section VI-B, and using as sample points the first \(500\) training points.

Table II presents the error rate on the test set when the maxout units on the dense layer have \(3\), and \(5\) terms respectively. Three tests results are presented: on pairs Automobile-Truck (A-T), Cat-Dog (C-D) and the average of the 45=\(9\cdot 10/2\) possible binary comparisons. The results obtained with the tropical division algorithm are compared with the L1 structured pruning. Table II also presents the number of parameters of the dense layer of the networks. We observe that the dense layer can be compressed to less than \(1\%\) of its original size with minimum performance loss.

**Remark 13**: _The proposed scheme is very competitive for the large compression regime, that is, in the case where there are very few parameters remaining. Furthermore, it can be performed assuming access to a small portion of the training data (500 out of 50000)._

### _CIFAR-10 Multiclass_

We then compress the multi-class neural network trained in the previous subsection for the CIFAR-10 dataset. As the input to the compression algorithm, we consider the output of the convolutional part of the neural network with dimension 2048.

We first represent the logits (the output before the softmax) as a vector of tropical polynomials using (26) in Example 5. We then simplify each of the tropical polynomials using the simplified algorithm of Section VI-E and a random sample of \(500\) training points. The results are summarized in Table III. Additionally, we present an improved prediction scheme, where the \(K\times\tilde{m}_{q}\) values \((\tilde{\mathbf{a}}_{k}^{l})^{T}\mathbf{z}^{1}\), with \(k=1,\ldots,K\) and \(l=1,\ldots,\tilde{m}_{q}\) are fed to a small single hidden layer ReLU neural network with 100 units. We refer to this small neural network as head.

Note that the additional parameters for the improved prediction are very few. Indeed for the case of \(\hat{m}_{q}=3,5,7\) we have \(4110\), \(6110\) and \(8110\) parameters respectively.

We observe again that the division algorithm works well on the large compression regime. Furthermore, in this example there is no improvement when we use more terms in the quotient polynomials. This is probably due to stacking in local optima.

We then exploit the good behaviour of the binary classification algorithm and apply it to the multi-class problem. To do so, we use a subset of the tropical polynomials obtained for binary classification. In each binary comparison, the division algorithm gives two tropical polynomials and the \(9\times 10/2=45\) binary comparisons give a total of \(90\) tropical polynomials. The value of a subset of the tropical polynomials is fed to a single hidden layer network with ReLU activations, as described in Figure 8.

Figure 9 presents the error of the compressed networks obtained using the multiple binary division algorithm of the previous paragraph. To produce these results we used a random set of \(m\) tropical polynomials with \(m=10,20,40,60\).

We observe that the proposed method works better than L1-structured pruning, when the number of remaining parameters is small.

**Remark 14**: _Note that in all the examples, the tropical division compression algorithms used only a small subset of the training data. This can be useful when simplifying a network where only parts of the training dataset are available._

**Remark 15**: _An interesting point of future research is to investigate whether producing additional sample points with the usual data augmentation techniques improves the final result._

\begin{table}
\begin{tabular}{c||c|c|c} Network & Orig. 2 classes & \(\hat{m}_{q}=5\) & \(\hat{m}_{q}=3\) \\ \hline Err. L1 A-T & 4.04\% & 25.8\% & 46.45\% \\ \hline Err. Div. X-A & 4.05\% & **4.6\%** & **4.4\%** \\ \hline Err. L1 C-D & 13.7\% & 49.3\% & 50\% \\ \hline Err. Div. C-D & 13.7\% & **15.4\%** & **16.85\%** \\ \hline Err. LiV avg. & 24.2\% & 25.4\% & 43.4\% & 43.62\% \\ \hline Err. Div. Avg. & 27.4\% 25.7\% & **3.85\%/3.6\%** & **4.06\%** & **3.37\%** \\ \hline \# Param. & \(2.1\cdot 10^{6}\) & \(2\cdot 10^{4}\) & \(1.2\cdot 10^{4}\) \\ \hline \% of Param. & 100\% & 0.95\% & 0.57\% \\ Remaining & 100\% & & 0.95\% \\ \end{tabular}
\end{table} TABLE II: Results for the CIFAR-10 dataset. Error rate of the original network, the reduced maxout networks and the reduced networks obtained using the L1 structured pruning algorithm. We present the Automobile vs Truck case, the Cat vs Dog case, and the average of the \(45=10\cdot 9/2\) different binary comparisons. The results are indexed using A-T, C-D, and Avg.

\begin{table}
\begin{tabular}{c|c|c|c|c} Network & Original & \(\hat{m}_{q}=7\) & \(\hat{m}_{q}=5\) & \(\hat{m}_{q}=3\) \\ \hline Error L1 & 14.18\% & 45\% & 57.3\% & 71.7\% \\ \hline Error Div. & 14.18\% & 41\(\pm\)1.8 \% & 40.1 \(\pm\) 1.7\% & 42.5 \(\pm\) 2\% \\ \hline Err. Div. Improved & 14.18\% & **26.6\% 0.3 \%** & **26.4 \(\pm\) 0.3\%** & **28.1 \(\pm\) 0.1\%** \\ \hline \# Param. & \(2.1\cdot 10^{6}\) & \(1.4\cdot 10^{6}\) & \(10^{5}\) & \(6.6\cdot 10^{4}\) \\ \hline \% of Param. & & & & \\ Remaining & 100\% & & & \\ \end{tabular}
\end{table} TABLE III: Results for the compression on the CIFAR-10 dataset. Average error rates for multiple executions. Comparison of the original division algorithm, the improved division algorithm and the L1 pruning. The number of parameters include also the parameters of the head.

Fig. 8: The structure of the proposed reduced network. Each pink box consists of linear neurons obtained by successive implementations of the division algorithm.
## VIII Conclusion

This work proposed a new framework for tropical polynomial division and its application to neural networks. We showed the existence of a unique quotient-remainder pair and characterized the quotient as the convex bi-conjugate of the difference of the dividend from the divisor. This characterization led to an exact algorithm based on polyhedral geometry. We also proposed an approximate algorithm based on the alternation between clustering and linear programming. We then focused on dividing composite tropical polynomials and proposed two modified algorithms, one producing simple quotients and the other composite. Finally, we applied tropical polynomial division techniques to simplify neural networks with ReLU activations. The resulting neural networks have either maxout or ReLU activations, depending on the division algorithm chosen. The results are promising and compare favorably with a simple baseline (L1 pruning).

There are several directions for future research. First, we may combine L1 regularization to the problem (19) to induce sparse solutions and improve further neural network compression. Another direction is to study compression problems for neural networks with many outputs, using a single division. A possible tool in this direction is the Cayley trick [53, 54]. Finally, the study of alternative optimization algorithms is of certain interest.

## Acknowledgment

The authors are grateful to Dr. George Retsinas for his valuable suggestions and comments.

### _Proof of Proposition 1_

The proof follows closely the proof of Proposition 3.1.6 in [39]. Observe that

\[p_{i}(x) =\max_{\mathbf{a}\in\mathbb{R}^{n}}\{[\mathbf{x}^{T}\ 1][\mathbf{a}^{T}\ \text{ ENF}_{p_{i}}(\mathbf{a})]^{T}\}\] \[=\max_{\mathbf{a}\in\mathbb{R}^{n}}\{[\mathbf{x}^{T}\ 1][\mathbf{a}^{T}\ \text{ ENF}_{p_{i}}(\mathbf{a})]^{T}\}. \tag{33}\]

Which implies that if \(\text{ENF}_{p_{1}}(\mathbf{a})\leq\text{ENF}_{p_{2}}(\mathbf{a})\) then \(p_{1}(\mathbf{x})\leq p_{2}(\mathbf{x}),\) for all \(\mathbf{x}\in\mathbb{R}^{n}\). On the other hand, assume that \(p_{1}(\mathbf{x})\leq p_{2}(\mathbf{x})\), and that for some \(\mathbf{a}_{0}\) it holds \(\text{ENF}_{p_{1}}(\mathbf{a}_{0})>\text{ENF}_{p_{2}}(\mathbf{a}_{0})\). Then, the point \([\mathbf{a}_{0}^{T}\ \text{ENF}_{p_{1}}(\mathbf{a}_{0})]\) does not belong to the convex set \(C=\{[\mathbf{a}^{T}\ z]\in\mathbb{R}^{n+1}:z\leq\text{ENF}_{p_{2}}(\mathbf{a})\}\). Thus, there is a \([c_{x}^{T}\ c_{z}]\in\mathbb{R}^{n+1}\) such that \(\mathbf{c}_{x}^{T}\mathbf{a}_{0}+c_{z}\text{ENF}_{p_{1}}(\mathbf{a}_{0})>\mathbf{c}_{x}^{T}\bm {a}+c_{z}z\), for all \([\mathbf{a}^{T}\ z]\in C\). Since \(C\) is unbounded below, \(c_{z}>0\). Furthermore, since \([\mathbf{a}^{T}\ \text{ENF}_{p_{2}}(\mathbf{a})]\in C\) it holds

\[\mathbf{a}_{0}^{T}\mathbf{x}_{0}+\text{ENF}_{p_{1}}(\mathbf{a}_{0})>\mathbf{a}^{T}\mathbf{x}_{0}+ \text{ENF}_{p_{2}}(\mathbf{a}),\]

for all \(\mathbf{a}\in\mathbb{R}^{n}\), where \(\mathbf{x}_{0}=\mathbf{c}_{x}/c_{z}\). Taking the maximum with respect to \(\mathbf{a}\) and using (33) in both sides we have:

\[p_{1}(\mathbf{x}_{0})\geq\mathbf{a}_{0}^{T}\mathbf{x}_{0}+\text{ENF}_{p_{1}} (\mathbf{a}_{0})> \max_{\mathbf{a}\in\mathbb{R}^{n}}\{\mathbf{a}^{T}\mathbf{x}_{0}+\text{ENF}_{p_{ 2}}(\mathbf{a})\}\] \[=p_{2}(\mathbf{x}_{0}),\]

which is a contradiction.

### _Proof of Proposition 2_

It is not difficult to see that if \(\tilde{q}(\mathbf{x})=\bigvee_{i=1}^{m_{q}}(\tilde{\mathbf{a}}_{i}^{T}\mathbf{x}+\tilde{b }_{i})\), and \(\tilde{q}\) satisfies (5), then \(\tilde{q}(\mathbf{x})\leq q(\mathbf{x})\), for all \(\mathbf{x}\in\mathbb{R}^{n}\). Indeed, since

\[\tilde{\mathbf{a}}_{i}^{T}\mathbf{x}+\tilde{b}_{i}\leq\tilde{q}(\mathbf{x})\leq p(\mathbf{x})- d(\mathbf{x}),\]

it holds \(\tilde{b}_{i}\leq l(\tilde{\mathbf{a}}_{i})\). Thus,

\[\tilde{\mathbf{a}}_{i}^{T}\mathbf{x}+\tilde{b}_{i}\leq\tilde{\mathbf{a}}_{i}^{T}\mathbf{x}+l( \mathbf{a}_{i})\leq q(\mathbf{x}).\]

Hence, \(\tilde{q}(\mathbf{x})\leq q(\mathbf{x})\), for all \(\mathbf{x}\in\mathbb{R}^{n}\).

To prove that \(q\) is a quotient, it remains to show that it is a tropical polynomial. The epigraph of \(q\) is the closed convex hull of the epigraph of \(f\)[44]

\[\text{epi}(q)=\overline{\text{conv}(\text{epi}(f))}. \tag{34}\]

On the other hand, \(f\) is a piecewise linear function. Thus, its epigraph is a union of a finite number of convex polyhedra. Indeed if \(A_{1},\ldots,A_{m_{f}}\) its linear regions (i.e., intersections of linear regions of \(p\) and \(d\)) then \(A_{i},\,i=1,...,m_{f}\) are convex polyhedra and

\[\text{epi}(f)=\bigcup_{i=1}^{m_{f}}\{(\mathbf{x},t)\in\mathbb{R}^{n+1}:\mathbf{x}\in A _{i},t\geq f(\mathbf{x})\}.\]

Observe that \(\text{epi}(q)\) is a convex polyhedron and \(q\) is a piecewise linear convex function. Thus, \(q\) is a tropical polynomial. From the definition it is obvious that the quotient is unique.

Let \(q\) be the quotient and assume that for some tropical polynomial \(\tilde{r}\) it holds

\[p(\mathbf{x})=(q(\mathbf{x})+d(\mathbf{x}))\vee\tilde{r}(\mathbf{x}),\qquad\text{for all }\mathbf{x}\in\mathbb{R}^{n}.\]

Assume that for some point \(\mathbf{x}_{0}\) it holds \(p(\mathbf{x}_{0})>d(\mathbf{x}_{0})+q(\mathbf{x}_{0})\) and that \(p\) has the linear form \(\mathbf{a}_{i0}^{T}\mathbf{x}+b_{i0}\) in a neighborhood of \(\mathbf{x}_{0}\). Then, \(\tilde{r}(\mathbf{x})=p(\mathbf{x})\) in this neighborhood. Furthermore, since \(\tilde{r}\) is convex, we have

\[\tilde{r}(\mathbf{x})\geq\bigvee_{i\in I}\mathbf{a}_{i}\mathbf{x}+b_{i}=r(\mathbf{x}),\]

where \(I=\{i\in\{1,\ldots,m_{p}\}:\exists\mathbf{x}\in\mathbb{R}^{n}\text{ with }p(\mathbf{x})>d(\mathbf{x})+q(\mathbf{x}),\text{ and }p(\mathbf{x})=\mathbf{a}_{i}\mathbf{x}+b_{i}\}\). It is not difficult to see that \(r\) is a remainder.

Fig. 9: _The error rate of the compressed neural networks obtained using the multiple division algorithm. The blue and yellow lines use divisions with \(\tilde{m}_{q}=3\) and \(\tilde{m}_{q}=5\), respectively. The horizontal axis represents the percentage of the remaining parameters. Finally, the green line corresponds to using random vectors instead of the quotients._
### _Proof of Proposition 3_

a) Let \(\bar{q}\) and \(\bar{r}\) be the quotient and the remainder of the division of \(r\) by \(d\). Then

\[r(\mathbf{x})=(d(\mathbf{x})+\bar{q}(\mathbf{x}))\vee\bar{r}(\mathbf{x}) \tag{35}\]

Combining with

\[p(\mathbf{x})=(d(\mathbf{x})+q(\mathbf{x}))\lor r(\mathbf{x}),\]

we get

\[p(\mathbf{x})=(d(\mathbf{x})+q(\mathbf{x}))\vee(d(\mathbf{x})+\bar{q}(\mathbf{x})) \vee\bar{r}(\mathbf{x})\] \[p(\mathbf{x})=(d(\mathbf{x})+q(\mathbf{x})\vee\bar{q}(\mathbf{x}))\vee\bar{r}( \mathbf{x}). \tag{36}\]

Equation (35) implies that \(\bar{r}(\mathbf{x})\leq r(\mathbf{x})\). This fact combined with the fact that \(q(\mathbf{x})\) is the quotient of the division of \(p(\mathbf{x})\) by \(d(\mathbf{x})\) and (36) implies that \(\bar{q}(\mathbf{x})\lor q(\mathbf{x})=q(\mathbf{x})\). Hence, since \(r(\mathbf{x})\) is the remainder of the division of \(p(\mathbf{x})\) by \(d(\mathbf{x})\) it holds \(r(\mathbf{x})\leq\bar{r}(\mathbf{x})\). Hence, \(r(\mathbf{x})=\bar{r}(\mathbf{x})\), for all \(\mathbf{x}\in\mathbb{R}^{n}\)

b) The proof is trivial.

### _Proof of Proposition 4_

Assume \(p(\mathbf{x})\) is given by (1) and \(d(\mathbf{x})\) by

\[d(\mathbf{x})=\bigvee_{i=1}^{m_{d}}\tilde{\mathbf{a}}_{i}\mathbf{x}+\tilde{b}_{i}.\]

(a) We first show that \(C\subset\textbf{dom}(l(\mathbf{a}))\). Assume that \(\mathbf{a}\in C\). That is \(\mathbf{a}+\tilde{\mathbf{a}}_{i}\in\textbf{Newt}(p)\) for all \(i=1,\ldots,m_{d}\). Then, there are \(\lambda_{1},\ldots,\lambda_{m_{d}}\geq 0\) with \(\sum_{i=1}^{m_{p}}\lambda_{1}=1\) such that

\[\mathbf{a}+\tilde{\mathbf{a}}_{i}=\sum_{i=1}^{n_{p}}\lambda_{i}\mathbf{a}_{i}.\]

Thus,

\[(\mathbf{a}+\tilde{\mathbf{a}}_{i})^{T}\mathbf{x}\leq\bigvee_{i=1}^{n_{p}}\mathbf{a}_{i}\mathbf{x }\leq p(\mathbf{x})-\bigwedge_{i=1}^{n_{p}}b_{i}.\]

Furthermore,

\[\tilde{\mathbf{a}}_{i}\mathbf{x}+\tilde{b}_{i}-\tilde{b}_{i}\geq\tilde{\mathbf{a}}_{i} \mathbf{x}+\tilde{b}_{i}-\bigvee_{j=1}^{m_{d}}\tilde{b}_{j}.\]

Combining the last two equations we have

\[p(\mathbf{x})-d(\mathbf{x})-\mathbf{a}^{T}\mathbf{x}\geq\bigwedge_{i=1}^{n_{p}}b_{i}-\sum_{j= 1}^{m_{d}}\tilde{b}_{j}>-\infty.\]

Hence, \(\mathbf{a}\in\textbf{dom}(l)\).

Conversely assume that \(\mathbf{a}\in\textbf{dom}(l)\) but \(\mathbf{a}\not\in C\). Then, there is a \(\tilde{\mathbf{a}}_{i}\) such that \(\mathbf{a}+\mathbf{a}_{i}\not\in\textbf{Newt}(p)\). Furthermore, since (p) is a convex set, there is a vector \(\mathbf{v}\) such that \((\mathbf{a}+\mathbf{a}_{i})^{T}\mathbf{v}>\mathbf{a}_{i}^{T}\mathbf{v}\) for all \(i=1,\ldots,n_{p}\). Hence,

\[l(\mathbf{a})=\inf_{\mathbf{x}\in\mathbb{R}^{n}}\{p(\mathbf{x})-d(\mathbf{x})- \mathbf{a}^{T}\mathbf{x}\}\] \[\leq\lim_{\lambda\rightarrow-\infty}p(\lambda\mathbf{v})-d(\lambda \mathbf{v})-\mathbf{a}^{T}(\lambda\mathbf{v})=-\infty\]

Thus, \(\mathbf{a}\not\in\textbf{dom}(l)\), which is a contradiction.

(b) Since

\[q(\mathbf{x})=\sup_{\mathbf{a}\in\mathbb{R}^{n}}\{\mathbf{a}^{T}\mathbf{x}-f^{*}(\mathbf{a})\}\]

It holds

\[q(\mathbf{x})=-\infty\Leftrightarrow l(\mathbf{a})=-\infty,\text{ for all }\mathbf{a} \Leftrightarrow C=\textbf{dom}(l)=\emptyset\]

(c) If for some \(i\) it holds \(\hat{\mathbf{a}}_{i}\not\in C\), then \(l(\hat{\mathbf{a}}_{i})=-\infty\). That is, for all \(b\) there is an \(\mathbf{x}\in\mathbb{R}^{n}\) with \(\hat{\mathbf{a}}_{i}^{T}\mathbf{x}+b>f(\mathbf{x})\). This is particularly true for \(b=\tilde{b}_{i}\), and thus \(q(\mathbf{x})\) is not a quotient.

### _Proof of Corollary 1_

The proof of the first part is immediate from Proposition 4. To contradict, assume that the inclusion in (ii) is not true. Then, there exists a vector \(\mathbf{x}\in\text{span}\{\mathbf{a}_{1},\ldots,\mathbf{a}_{m_{p}}\}^{\perp}\) but \(\mathbf{x}\not\in\text{span}\{\hat{\mathbf{a}}_{1},\ldots,\hat{\mathbf{a}}_{m_{q}}\}^{\perp}\). Thus, there is an index \(i_{0}\) such that \(\hat{\mathbf{a}}_{i_{0}}^{T}\mathbf{x}\neq 0\). Without loss of generality assume that \(\hat{\mathbf{a}}_{i_{0}}^{T}\mathbf{x}>0\) (if it is negative, use \(-\mathbf{x}\) in the place of \(\mathbf{x}\)). Using (i) we have \(\lim_{\lambda\rightarrow\infty}(q(\lambda\mathbf{x})+d(\lambda\mathbf{x}))=\infty\), but \(\lim_{\lambda\rightarrow\infty}p(\lambda\mathbf{x})<\infty\). Hence, there is a \(\lambda\) such that \(p(\lambda\mathbf{x})<q(\lambda\mathbf{x})+d(\lambda\mathbf{x})\), which contradicts the fact that \(q\) is the quotient.

### _Proof of Proposition 5_

First we need to show that all the components of \(\mathbf{a}^{E,z}\) are positive. Due to the convexity of \(E\), either

\[\inf\{z:(\mathbf{x},z)\in E\}=-\infty,\text{ for all }\mathbf{x}\in\mathbb{R}^{n},\]

or

\[\inf\{z:(\mathbf{x},z)\in E\}>-\infty,\text{ for all }\mathbf{x}\in\mathbb{R}^{n}.\]

If \(E=\mathbb{R}^{n+1}\), then there are no non-trivial constraints, i.e., \(L=0\). In the following assume that \(L>0\), that is there are some nontrivial constraints, and \(E\neq\mathbb{R}^{n+1}\). From the definition of \(E\), for all \(\mathbf{x}\in\mathbb{R}^{n}\), it holds \(\inf\{z:(\mathbf{x},z)\in E\}<\infty\), and \(\sup\{z:(\mathbf{x},z)\in E\}=\infty\).

Note that it is not possible to have a constraint in (10) in the form \([\mathbf{A}^{E,\mathbf{x}}]_{\mathbf{x}}\geq[b^{E}]_{l}\), that is nontrivial, i.e., \([\mathbf{A}^{E,\mathbf{x}}]_{l}\neq\mathbf{0}^{T}\). Indeed, if such a constraint existed, then for an \(\mathbf{x}\) not satisfying this constraint we would have \(\inf\{z:(\mathbf{x},z)\in E\}=\infty\), which is a contradiction. On the other hand if there were a constraint with \([\mathbf{A}^{E,\mathbf{x}}]_{\mathbf{x}}+[\mathbf{a}^{E,z}]_{l}z\geq[\mathbf{b}^{E}]_{l}\), with \([\mathbf{a}^{E,z}]_{l}<0\) then, for all \((\mathbf{x},z)\in E\) we would have \(z\leq(-[\mathbf{A}^{E,\mathbf{x}}]_{l}\mathbf{x}+[b^{E}]_{l})/[\mathbf{a}^{E,z}]_{l}\) and thus \(\sup\{z:(\mathbf{x},z)\in E\}<\infty\).

It is then easy to see that the function \(q\) given by (11) has as epigraph the set \(E\). Therefore, it is the quotient.

To determine the remainder it is sufficient to find all the terms \(i\) of \(p(\mathbf{x})\) for which \(p(\mathbf{x})>\tilde{p}(\mathbf{x})\), for an \(\mathbf{x}\in\mathbb{R}^{n}\) that the term \(i\) attains the maximum in \(p(\mathbf{x})\). The function \(p(\mathbf{x})-\tilde{p}(\mathbf{x})\) is linear on \(P_{i,j}\) for all \(i,j\). Therefore, since \(p(\mathbf{x})-\tilde{p}(\mathbf{x})\geq 0\), and \(\mathbf{x}_{i,j}\) is in the relative interior of \(P_{i,j}\) we have either \(p(\mathbf{x})-\tilde{p}(\mathbf{x})=0\), for all \(\mathbf{x}\in P_{i,j}\), or \(p(\mathbf{x}_{i,j})>\tilde{p}(\mathbf{x}_{i,j})\). Hence, the procedure described in Algorithm 1 produces the remainder.

### _Proof of Proposition 6_

Observe that

\[p(\mathbf{x}_{j})-d(\mathbf{x}_{j}
where the last inequality is a consequence of the first constraint of (16). Furthermore,

\[e(t)=\sum_{j=1}^{N}f(\mathbf{x}_{j})-\sum_{j=1}^{N}\bigvee_{i=1}^{\tilde{m}_{q}}(\hat {\mathbf{a}}_{i,t}^{T}\mathbf{x}_{j}+\hat{b}_{i,t}), \tag{37}\]

where \(\hat{\mathbf{a}}_{i,t}\) is the solution of (17) at the iteration \(t\). Each term of the last sum can be written as:

\[\bigvee_{i=1}^{\tilde{m}_{q}}(\hat{\mathbf{a}}_{i,t}^{T}\mathbf{x}_{j}+\hat{b}_{i,t}) =(\hat{\mathbf{a}}_{i^{\prime}_{j,t},t}^{T}\mathbf{x}_{j}+\hat{b}_{i^{\prime}_{j,t},t} ),\]

where \(i^{\prime}_{j,t}\) is the partition that \(j\) belongs after step \(t\), i.e., \(j\in I_{i^{\prime}_{j,t},t}\) and \(I_{i,t}\) is the partition obtained at Step 6 in iteration \(t\). Furthermore due to the optimization in (16),

\[\sum_{j=1}^{N}\hat{\mathbf{a}}_{i^{\prime}_{j,t},t}^{T}\mathbf{x}_{j}+\hat{b}_{i^{ \prime}_{j,t},t}\leq\sum_{j=1}^{N}\hat{\mathbf{a}}_{i^{\prime}_{j,t},t+1}^{T}\mathbf{ x}_{j}+\hat{b}_{i^{\prime}_{j,t},t+1}.\]

From Step 6,

\[\sum_{j=1}^{N}\hat{\mathbf{a}}_{i^{\prime}_{j,t},t+1}^{T}\mathbf{x}_{j}+\hat{b}_{i^{ \prime}_{j,t},t+1}\leq\sum_{j=1}^{N}\hat{\mathbf{a}}_{i^{\prime}_{j,t+1},t+1}^{T} \mathbf{x}_{j}+\hat{b}_{i^{\prime}_{j,t+1},t+1}.\]

Therefore, \(\sum_{j=1}^{N}\bigvee_{i=1}^{\tilde{m}_{q}}(\hat{\mathbf{a}}_{i,t}^{T}\mathbf{x}_{j}+ \hat{b}_{i,t})\) is non-increasing, and thus, \(e(t)\) is non-decreasing.

### _Proof of Proposition 7_

We first prove (20). The quotients \(Q(p_{1},d)\) and \(Q(p_{2},d)\) satisfy

\[p_{1}(\mathbf{x}) \geq d(\mathbf{x})+Q(p_{1},d)(\mathbf{x}), \tag{38}\] \[p_{2}(\mathbf{x}) \geq d(\mathbf{x})+Q(p_{2},d)(\mathbf{x}). \tag{39}\]

Thus,

\[p_{1}(\mathbf{x})+p_{2}(\mathbf{x})\geq d(\mathbf{x})+[Q(p_{1},d)+Q(p_{2},d)(\mathbf{x})+d( \mathbf{x})].\]

But \(Q(p_{1}+p_{2},d)\) is the maximum tropical polynomial satisfying this inequality. Thus it is greater than or equal to the expression in the bracket, i.e.,

\[Q(p_{1}+p_{2},d)(\mathbf{x})\geq Q(p_{1},d)(\mathbf{x})+Q(p_{2},d)(\mathbf{x})+d(\mathbf{x}).\]

To prove (21), take the maximum of (38) and (39). Using the distributive property we get

\[p_{1}(\mathbf{x})\lor p_{2}(\mathbf{x})\geq d(\mathbf{x})+[Q(p_{1},d)(\mathbf{x})\lor Q(p_{2},d)(\mathbf{x})].\]

Thus,

\[Q(p_{1}\lor p_{2},d)\geq Q(p_{1},d)(\mathbf{x})\lor Q(p_{2},d)(\mathbf{x}).\]

We then prove (22). The quotients \(Q(p,d_{1})\) and \(Q(p,d_{2})\) satisfy

\[p(\mathbf{x}) \geq d_{1}(\mathbf{x})+Q(p,d_{1})(\mathbf{x}), \tag{40}\] \[p(\mathbf{x}) \geq d_{2}(\mathbf{x})+Q(p,d_{2})(\mathbf{x}). \tag{41}\]

Thus,

\[p(\mathbf{x})\geq d_{1}(\mathbf{x})+d_{2}(\mathbf{x})+[Q(p,d_{1})(\mathbf{x})+Q(p,d_{2})(\mathbf{ x})-p(\mathbf{x})].\]

Hence,

\[Q(p,d_{1}+d_{2})(\mathbf{x})\geq Q(p,d_{1})(\mathbf{x})+Q(p_{2},d)(\mathbf{x})-p(\mathbf{x}).\]

To prove (23), take the maximum of (40) and (41)

\[p(\mathbf{x}) \geq[d_{1}(\mathbf{x})+Q(p,d_{1})(\mathbf{x})]\vee[d_{2}(\mathbf{x})+Q(p,d_{ 2})(\mathbf{x})]\] \[\geq[d_{1}(\mathbf{x})\lor d_{2}(\mathbf{x})]+[Q(p,d_{1})(\mathbf{x})\wedge Q (p,d_{2})(\mathbf{x})]\]

Therefore,

\[Q(p,d_{1}\lor d_{2})(\mathbf{x})\geq Q(p,d_{1})(\mathbf{x})\wedge Q(p,d_{2})(\mathbf{x})\]

### _Proof of Proposition 8_

The Newton polytopes of \(p\) and \(q\) are given by the zonotopes

\[\text{Newt}(p)=\{\mathbf{A}\mu:\mu\in[0,1]^{n}\}, \tag{42}\] \[\text{Newt}(q)=\{\hat{\mathbf{A}}\mathbf{\lambda}:\lambda\in[0,1]^{m_{q}}\}, \tag{43}\]

where \(\hat{\mathbf{A}}=[\hat{\mathbf{a}}_{1}\ \ldots\ \hat{\mathbf{a}}_{m\tau}]\). We start proving the following lemma.

**Lemma 1**: _Let \(\text{Newt}(p)\), \(\text{Newt}(q)\) as the above. Then \(\text{Newt}(q)\subset\text{Newt}(p)\) if and only if_

\[\mathbf{0}\leq\mathbf{A}^{-1}\hat{\mathbf{a}}_{i},\ \ \mathbf{A}^{-1}\sum_{i=1}^{m_{q}}\hat{\mathbf{a}}_{i} \leq\mathbf{1}, \tag{44}\]

_for \(i=1,\ldots,m_{q}\),_

_Proof of the lemma: \(\text{Newt}(q)\subset\text{Newt}(p)\) if and only if for any \(\mathbf{\lambda}\in[0,1]^{m_{q}}\) there exists a \(\mathbf{\mu}\in[0,1]^{n}\) such that \(\mathbf{A}\mu=\hat{\mathbf{A}}\mathbf{\lambda}\). Thus, \(\text{Newt}(q)\subset\text{Newt}(p)\) is equivalent to_

\[\mathbf{0}\leq\mathbf{A}^{-1}\hat{\mathbf{A}}\mathbf{\lambda}\leq\mathbf{1}, \tag{45}\]

_for all \(\mathbf{\lambda}\in[0,1]^{m_{q}}\)._

Direct part: Assume \(\text{Newt}(q)\subset\text{Newt}(p)\). Substituting \(\mathbf{\lambda}=\mathbf{e}_{i}\), i.e., the \(i-\)th unit vector, in (45), we get the first inequality of (44). Now with \(\mathbf{\lambda}=\mathbf{1}\) we get the second inequality of (44).

Converse part: Assume that (44) holds true. Multiplying the first inequality of (44) by \(\lambda_{i}\geq 0\) and adding, we get the first inequality of (45). Furthermore, for \(0\leq\lambda_{i}\leq 1\), we have \(\sum_{i=1}^{m_{q}}\lambda_{i}\mathbf{A}^{-1}\hat{\mathbf{a}}_{i}\leq\mathbf{A}^{-1}\sum_{i= 1}^{m_{q}}\hat{\mathbf{a}}_{i}\leq\mathbf{1}\), which implies the right part of inequality of (45). \(\Box\)

From Lemma 1, \(q(\mathbf{x})\leq p(\mathbf{x})\) for all \(\mathbf{x}\) is equivalent to \(\text{ENF}_{q}(\mathbf{a})\leq\text{ENF}_{p}(\mathbf{a})\), for all \(\mathbf{a}\). We first compute \(\text{ENF}_{p}(\mathbf{a})\). It holds

\[\text{ENF}_{p}(\mathbf{a})=\mathbf{B}^{T}\mathbf{A}^{-1}\mathbf{a}, \tag{46}\]

for \(\mathbf{a}\in\text{Newt}(p)\) and \(-\infty\) otherwise. Indeed, the extended Newton function of \(p\) is given by \(\text{ENetv}_{p}(\mathbf{a})=\max\{\sum_{i=1}^{n}\mu_{i}b_{i}:\sum_{i=1}^{n}\mu_{i} \mathbf{a}_{i}=\mathbf{a},\mu_{i}\in[0,1]\}\). However, \(\sum_{i=1}^{n}\mu_{i}\mathbf{a}_{i}=\mathbf{a}\) has a solution \(\mathbf{\mu}=[\mu_{1}\ \ldots\mu_{n}]^{T}\in[0,1]^{n}\), if and only if \(a\in\text{Newt}(p)\). Furthermore, for \(a\in\text{Newt}(p)\) the solution is unique and given by \(\mathbf{\mu}=\mathbf{A}^{-1}\mathbf{a}\). Thus \(\text{ENF}_{p}(\mathbf{a})\) is given by (46).

The first two conditions of (29) are equivalent to the inclusion \(\text{Newt}(q)\subset\text{Newt}(p)\). To prove the direct part of the proposition observe that \(\hat{b}_{i}\leq\
for all \(\hat{\mathbf{a}}\not\in\text{Newt}(q)\). It remains to show that \(\text{ENF}_{q}(\hat{\mathbf{a}})\leq\text{ENF}_{p}(\hat{\mathbf{a}})\), for all \(\hat{\mathbf{a}}\in\text{Newt}(q)\). For such an \(\hat{\mathbf{a}}\) it holds

\[\text{ENF}_{q}(\hat{\mathbf{a}}) =\max\{\sum_{i=1}^{m_{q}}\lambda_{i}\hat{b}_{i}:\sum_{i=1}^{m_{q} }\lambda_{i}\hat{\mathbf{a}}_{i}=\hat{\mathbf{a}},\lambda_{i}\in[0,1]\}\] \[\leq\max\{\sum_{i=1}^{m_{q}}\lambda_{i}\mathbf{B}^{T}\mathbf{A}^{-1}\hat{ \mathbf{a}}_{i}:\sum_{i=1}^{m_{q}}\lambda_{i}\hat{\mathbf{a}}_{i}=\hat{\mathbf{a}},\lambda _{i}\in[0,1]\}\] \[=\mathbf{B}^{T}\mathbf{A}^{-1}\hat{\mathbf{a}}=\text{ENF}_{p}(\hat{\mathbf{a}})\]
