Tightening the Difference-of-Convex Formulation for the Piecewise Linear Approximation in General Dimensions

1. Introduction

A continuous piecewise linear (CPWL) function is a continuous function defined on a compact set which can be partitioned into a finite set of affine domains within which the function is affine. CPWL functions are a fundamental tool in mathematics, engineering, and data science. This is due to their ability to approximate any continuous nonlinear function defined on a compact set to an arbitrary level of precision using a finite number of linear pieces (Huang 2020). CPWL functions are also an essential tool when fitting a data set for which the analytical expression of the nonlinear relationship is unknown (Rebennack and Krasko 2020). In optimization models, using CPWL functions to represent complex nonlinear relationships enables the use of mature mixed-integer linear programming (MILP) solvers, which are often significantly faster than nonlinear solvers (Vielma et al. 2010, Vielma 2015).

Not only can CPWL functions be represented in MILP models, but the optimal CPWL fitting problem itself can also be expressed as a MILP problem (Toriello and Vielma 2012). The objective function may vary depending on the application. For instance, given a fixed number of affine functions, one may wish to determine the CPWL approximation that minimizes an approximation error metric—for example, the maximum or average error (Rebennack and Krasko 2020). Alternatively, one may wish to identify the CPWL approximation with the minimal number of affine pieces (Kong and Maravelias 2020).

MILP formulations for the CPWL fitting problem have been extensively investigated for univariate functions and one-dimensional data sets. Frenzen et al. (2010) studied optimal partitioning of the range of the univariate function to be approximated and provided a quantitative result regarding how the number of segments depends on the function and the approximation error. Rebennack and Krasko (2020) developed a finitely convergent method to generate CPWL approximations for one-dimensional data sets or continuous univariate functions identifying the minimal number of breakpoints required to meet an acceptable approximation error. Similarly, Kong and Maravelias (2020) introduced a novel method of enforcing continuity at breakpoints through a set of linear constraints. The two formulations are compared in detail in Warwicker and Rebennack (2022), and it is shown that the formulation from Rebennack and Krasko (2020) is generally faster in practice. The work from Warwicker and Rebennack (2023) extends the work from Warwicker and Rebennack (2022) and uses a Benders decomposition approach to improve the solution time of the MILP problem. Ploussard (2024) introduces a tightening method that significantly shortens the MILP solution time of the CPWL fitting problem when hierarchically minimizing the number of linear pieces and approximation error. More recently, Warwicker and Rebennack (2025) presented an efficient approach for approximating continuous univariate functions with a CPWL function by embedding a linear-time algorithm into an adaptive discretization framework.

For bivariate functions and two-dimensional data sets, most existing methods constrain the affine domains of the CPWL approximation to be simplices. For example, D’Ambrosio et al. (2010) partition the CPWL function domain into rectangles that are later split into triangles along predefined diagonals. Toriello and Vielma (2012) use a similar method but select a set of splitting diagonals that minimize the approximation error. Rebennack and Kallrath (2015) use a different triangulation strategy that iteratively partitions the triangles covering the approximation domain into smaller ones until a target approximation error is reached. Duguet and Ngueveu (2022) propose a similar iterative refinement, but without imposing the affine domains to be simplices. Bärmann et al. (2023) introduce an approximation algorithm for the optimal piecewise linear interpolation of the continuous bivariate function xy over rectangular domains.

The literature addressing the optimal CPWL fitting problem for data sets of higher dimensions is significantly more limited (Misener and Floudas 2010, Rebennack and Kallrath 2015). Besides, most of the existing methods essentially extend the simplex-based partitioning strategy used for two dimensions, which result in computational challenges due to the exponential growth of the number of possible simplex-based partitions as the data dimension increases (Hughes and Anderson 1996). Kazda and Li (2021) introduced the first MILP formulation able to identify the optimal CPWL approximation of data sets of any dimension for a given maximum approximation error. They achieve this by leveraging the difference-of-convex (DC) representation of CPWL functions (Kripfganz and Schulze 1987), which allows for the affine domains partitioning the CPWL function to be implicitly defined by the difference of the pointwise maxima of two sets of affine functions. The strength of the DC representation lies in its ability to describe any CPWL function with affine domains of arbitrary shape. However, the computational cost of the MILP problem rapidly increases with the dimension of the data set, the number of data points, and the number of affine functions.

Deep learning techniques offer an alternative way to identify CPWL approximations (Huang 2020). In fact, the function represented by a neural network (NN) using ReLU activation function is a CPWL function, and the NN loss function represents the CPWL approximation error (DeVore et al. 2021). To the best of our knowledge, the only existing methods to identify CPWL approximations of data sets in multiple dimensions without explicitly specifying the partitioning of affine domains are the ReLU NN approach and the DC MILP approach introduced by Kazda and Li (2021). Training an NN to identify a good CPWL approximation is significantly faster than solving an MILP problem, contributing to NN growing popularity in recent years. However, the MILP approach provides three key advantages over the ReLU NN approach:

• The NN approach uses a gradient-based method that identifies local optima, whereas the MILP approach is guaranteed to yield the CPWL function with the smallest approximation error.

• Contrary to the MILP approach, the NN approach does not allow for the explicit enforcement of a target maximum approximation error.

• Because of the complex relationship between the structure (width and depth) of an NN and the number of affine pieces of the CPWL functions they can model, it is difficult for NNs to identify good CPWL approximations with a manageable number of affine pieces (Chen et al. 2023). Conversely, the MILP approach allows the user to identify CPWL functions with a predefined, or minimal, number of linear pieces. This is a critical feature when these CPWL functions are subsequently embedded in MILP formulations to model nonlinear relationships, where the number of binary variables required to model CPWL functions typically increases with the number of affine pieces that compose them, making the MILP problem harder to solve. In other words, although an NN can quickly identify a good-quality CPWL approximation of a nonlinear relationship, the resulting approximation may render the MILP formulation in which it is embedded computationally intractable.

To address the high computational cost of solving the MILP problem, authors Kazda and Li (2024) introduce an iterative linear programming (LP) approach that identifies a valid CPWL approximation given a target maximum error. However, the identified CPWL solution is not necessarily optimal. A priori, identifying the global optimum of the CPWL fitting problem can only be achieved by solving an MILP problem. Therefore, it is critical to improve the MILP formulation of the CPWL fitting problem to reduce its computational cost. The method described in Ploussard (2024) achieves a tightening of the MILP problem that significantly improves the solution time, but the method only applies to one-dimensional data sets. This raises the need to develop MILP tightening methods that can be extended to data sets of any dimensions.

This paper aims to address this research gap. To achieve this, we introduce a critical class of CPWL interpolations called “well-behaved” and demonstrate that any CPWL interpolation has a well-behaved version. This important result is then leveraged to identify six tightening strategies that can be combined to significantly reduce the solution time of the MILP CPWL fitting problem.

The five main contributions of this paper are the following. First, we formalize the concept of well-behaved CPWL interpolations, a class of CPWL interpolations in which each affine piece interpolates a number of points greater than the dimension of the interpolated data set. Second, we demonstrate that any CPWL interpolation of a data set has a well-behaved version, which serves as a key theoretical foundation for our tightening method. Third, we leverage this result to introduce six tightening strategies for the MILP formulation of the CPWL fitting problem. Fourth, we analyze the theoretical impact of each tightening strategy together with the time complexity of the preprocessing step involved. Fifth, we assess the practical impact of multiple combinations of the tightening strategies. The proposed tightening strategies enable the identification of an optimal CPWL approximation in a relatively low computation time (i.e., up to 20 times faster than the existing literature).

The remainder of the article is organized as follows: Section 2 introduces the theoretical background of the CPWL-fitting problem. Section 3 describes the proposed MILP formulation. Section 4 formalizes the concept of well-behaved CPWL fitting and demonstrates the fact that any CPWL interpolation or approximation of a data set has a well-behaved version. Section 5 introduces the six tightening strategies of the MILP formulation. Section 6 introduces the case studies and assesses the performances of several combinations of tightening strategies. Section 7 presents the conclusions.

2. Background

Let $S={({\mathit{x}}_i,z_i)}_{i=1,\dots ,N}={(x_{i,1},\dots ,x_{i,d},z_i)}_{i=1,\dots ,N}$ be a set of N points of ${\mathbb{R}}^{d+1}$. We assume that the points ${({\mathit{x}}_i)}_{i=1,\dots ,N}$ are in general position in ${\mathbb{R}}^d$—that is, no subset of $d+1$ points is contained within a single hyperplane (Edelsbrunner et al. 1986). Let $P_{{\mathbb{R}}^d}:{\mathbb{R}}^{d+1}\to {\mathbb{R}}^d$ be the projection defined as $P_{{\mathbb{R}}^d}(x_1,\dots ,x_d,z)=(x_1,\dots ,x_d)$. Let $S_{{\mathbb{R}}^d}={({\mathit{x}}_i)}_{i=1,\dots ,N}$ be the projection of S by $P_{{\mathbb{R}}^d}$. Let $D=\mathit{\text{Conv}}(S_{{\mathbb{R}}^d})$ be the convex hull of $S_{{\mathbb{R}}^d}$.

Note that, contrary to the more common definition of CPWL functions from Geißler et al. (2012), we do not assume that the $D_k$ form a partition of D. The intersection of two affine domains is not required to have an empty interior. In addition, an affine domain can be empty, and two affine domains can be equal as long as the corresponding affine pieces are equal. This implies that the set of affine functions representing the CPWL function is not unique. However, apart from this, Definition 1 describes exactly the same mathematical object as the CPWL functions defined in the literature. The flexibility of our definition enables us to compare CPWL functions to one another, which will be useful in the following section.

The proof of Proposition 1 is provided in Kripfganz and Schulze (1987). The above formulation is called the DC representation of f. Note that the DC representation of a CPWL function is not unique.

The proof of Proposition 2 is given in Kripfganz and Schulze (1987) and stems from the DC representation of f.

3. MILP Formulation

We aim to identify the optimal CPWL fitting of a set of points S. More specifically, we aim to identify the CPWL ε-approximation of S that minimizes an objective Q subject to a given error tolerance $\epsilon$. For example, Q can be the maximum fitting error, the average fitting error, or the number of affine pieces.

The following MILP formulation is largely based on the formulation presented in Kazda and Li (2021), where $f({\mathit{x}}_i)$, $f^{+}({\mathit{x}}_i)$, $f^{-}({\mathit{x}}_i)$, $e_i$, ${\mathit{a}}_j^c$, and $b_j^c$ are continuous variables; ${\delta }_{i,j}^c$ are binary variables; and $M_i^c$ are big-M parameters.

Equation (1) represents the objective function to be minimized. As long as Q is a linear expression, Equations (1)–(10) formulate an MILP problem. For the rest of the paper, we will refer to this MILP problem as $\mathit{\text{MILP}}1(Q)$. Equation (2) is the DC representation of the CPWL function, where $f({\mathit{x}}_i)$, $f^{+}({\mathit{x}}_i)$, and $f^{-}({\mathit{x}}_i)$ represent the value at ${\mathit{x}}_i$ of the functions f, $f^{+}$, and $f^{-}$. Equations (3) and (4) formulate the representation of each convex function $f^c$ as the maximum of a set of affine functions $f_j^c$, where ${\mathit{a}}_j^c$ and $b_j^c$ are the linear coefficients and bias terms of each affine piece $f_j^c$. The binary variable ${\delta }_{i,j}^c$ is equal to one when the point ${\mathit{x}}_i$ belongs to the domain $D_j^c$, and the big-M parameter $M_i^c$ ensures that $f^c({\mathit{x}}_i)=f_j^c({\mathit{x}}_i)$ when ${\mathit{x}}_i\in D_j^c$. Note that the big-M parameter $M_i^c$ in Equation (4) must be large enough so that it does not constrain the term $f^c({\mathit{x}}_i)-{\mathit{a}}_j^c{}^T{\mathit{x}}_i-b_j^c$ when ${\delta }_{i,j}^c=0$. In Equation (5), we will identify a tight value for $M_i^c$. $P^{+}$ and $P^{-}$ are the numbers of affine pieces of $f^{+}$ and $f^{-}$, respectively. Equation (5) ensures that each point ${\mathit{x}}_i$ belongs to at least one affine domain of $f^{+}$ and one affine domain of $f^{-}$. Note that, contrary to the literature (e.g., Kazda and Li 2021), we use an inequality constraint instead of an equality constraint. By doing so, we allow the points ${\mathit{x}}_i$ to belong to multiple affine domains. This can occur if, for example, some affine pieces are identical or if some points are at the border between two or more affine domains. This inequality is also critical to ensure that some well-behaved CPWL solutions are feasible in $\mathit{\text{MILP}}1(Q)$, which is a property that will be discussed in the next section. Equation (6) defines the approximation error $e_i$ at each point ${\mathit{x}}_i$ as the distance between $z_i$ and the value of the CPWL function at ${\mathit{x}}_i$. Equation (7) ensures that this distance is never greater than the specified error tolerance $\epsilon$. Finally, Equations (8)–(10) define the domain of the continuous and binary variables.

Note that some optimization solvers like Gurobi (Gurobi Optimization LLC 2024a) make it possible to formulate Equation (4) without using a big-M parameter. Instead, Equation (4) can be replaced by the indicator constraint below (Gurobi Optimization LLC 2024b):

We note $\mathit{\text{MILPIC}}1(Q)$ the alternative MILP formulation using the indicator constraint, which is composed of Equations (1)–(3) and (5)–(11).

The objective can be set to minimize different error metrics. Formulating the objective with Equation (12) minimizes the average fitting error ($L_1$ norm), the discrete analogue of minimizing the area between two functions. Alternatively, using Equation (13) minimizes the maximum fitting error ($L_{\infty }$ norm), which provides a strict pointwise guarantee on the approximation quality.

Other objective functions may include the number of affine pieces of $f^{+}$, $f^{-}$, or f, as well as a combination of error metric and number of affine pieces. These alternative objective functions can be found in Electronic Companion Section EC.1.

4. Well-Behaved CPWL Interpolations

This section introduces a new class of CPWL interpolations called “well-behaved” CPWL interpolations. Loosely speaking, well-behaved CPWL interpolations are composed of affine pieces whose gradients are not excessively steep. This class of CPWL interpolations aligns with the intuitive representation of what the CPWL interpolation of a data set should look like, and they have desirable properties when modeled in $\mathit{\text{MILP}}1(Q)$. This section formalizes the definition of a well-behaved CPWL interpolation and defines what constitutes a well-behaved version of a CPWL interpolation. We then demonstrate that any CPWL interpolation has a well-behaved version. This result is critical to the tightening procedure introduced in Section 5.

Proposition 3 is a fundamental result from linear algebra (Boyd and Vandenberghe 2004). The system of linear equations ${\mathit{a}}^T{\mathit{x}}_i+b=z_i,i\in \{1,\dots ,d+1\}$, with variables $\mathit{a}\in {\mathbb{R}}^d$ and $b\in \mathbb{R}$, is exactly determined and has a unique solution due to the general position of the ${\mathit{x}}_i$. If the number of points is less than $d+1$, the system is underdetermined, and there is an infinite number of possible linear interpolations. If the number of points is greater than $d+1$, because the ${\mathit{x}}_i$ are in general position, the system becomes overdetermined, and there exists either one or no linear interpolation.

In other words, a well-behaved version of a CPWL interpolation f of S is a transformation of f where each affine piece has been adjusted, or “tilted,” to interpolate at least $d+1$ points of S, including the points it was originally interpolating (before being tilted). An illustration for $d=1$ is provided in Figure 1. In the following sections, we show that well-behaved CPWL interpolations have interesting properties that can be used to tighten the MILP formulation. To the best of the authors’ knowledge, this is the first time this class of CPWL interpolations is introduced. We prove below that for any CPWL interpolation, a well-behaved version always exists.

Proposition 4 is a well-established fact from LP theory, as described in Bertsimas and Tsitsiklis (1997). The $\min (n,m)$ constraints “${\mathit{A}}_i{\mathit{x}}^{*}=b_i$” are also known as the active constraints of the solution ${\mathit{x}}^{*}$. The point ${\mathit{x}}^{*}$ lies on the boundary of U, and, if $m\ge n$, ${\mathit{x}}^{*}$ is a vertex of U.

The proof of Corollary 1 directly follows from Theorem 1 and Definition 7.

An important application of this theorem is that we can significantly reduce the set of CPWL functions to be considered in $\mathit{\text{MILP}}1(Q)$ without impacting the feasibility of the problem. Specifically, if a CPWL $\epsilon$-approximation f of S exists, a well-behaved version g of f also exists. The two CPWL solutions f and g have the same number of affine pieces and are equivalent with respect to S—that is, $f({\mathit{x}}_i)=g({\mathit{x}}_i),\forall i\in \{1,\dots ,N\}$. Therefore, we can eliminate CPWL solutions that are not well-behaved from the feasible region of $\mathit{\text{MILP}}1(Q)$ because they are redundant solutions and may have excessively steep gradients due to their underdetermined affine pieces. In the next section, we will see how this theorem allows us to derive valid tight bounds and constraints for the MILP problem.

5. Tightening the MILP Formulation

This section introduces six strategies to tighten and enhance the formulation of $\mathit{\text{MILP}}1(Q)$. First, we demonstrate that we can fix one of the affine pieces of $f^{-}$ without affecting the set of CPWL solutions. Next, we impose an ordering on the affine pieces of $f^{+}$ and $f^{-}$ that does not affect the set of CPWL solutions. We can also impose that each affine piece of $f^{+}$ and $f^{-}$ contains at least $d+1$ points, which serves as a valid tightening of the set of feasible well-behaved CPWL solutions. Alternatively, we may require that each affine piece of f contains at least $d+1$ points, which requires the use of additional variables. Finally, we identify tight values for the big-M parameters and establish tight bounds for the linear variables.

Additional strategies that leverage the convexity of the functions $f^{+}$ and $f^{-}$ and of their affine domains are included in the electronic companion, Section EC.3. However, contrary to the six strategies presented here, their impact on the feasible region is unclear, and they do not seem to have a significant impact on the solution time.

Herein, it is assumed that $N\ge d+1$.

5.1. Fixing One Affine Piece

Lemma 3 follows directly from the translation invariance of the max function—that is, $\max (a-c,b-c)=\max (a,b)-c,\hspace{1em}\forall a,b,c\in \mathbb{R}$.

Note that Equation (14) affects not only the set of possible values for ${\mathit{a}}_1^{-},b_1^{-}$, but also for all ${\mathit{a}}_j^c,b_j^c$. Also note that the fixed values of zero for ${\mathit{a}}_1^{-},b_1^{-}$ are chosen arbitrarily; any value could be imposed and serve as a valid tightening constraint for the problem. Imposing Equation (14) effectively reduces the dimensionality of the search space by $(d+1)$. Formally, Equation (14) does not introduce any additional linear constraints; instead, it introduces tight bounds on $(d+1)$ linear variables.

5.2. Sorting the Affine Pieces

Let ${\mathit{a}}_j^c=\left(\begin{array}{c}{a}_{j,1}^c\\ ⋮\\ a_{j,d}^c\end{array}\right)\in {\mathbb{R}}^d$ denote the linear coefficients of the affine piece $f_j^c$.

It is important to note that, after sorting the affine pieces of $f^{+}$ and $f^{-}$, the procedure outlined in the proof of Theorem 2 can be applied to impose Equation (14) with no impact on the CPWL solution. Thus, Equations (14) and (15) are compatible and can be applied simultaneously with no impact on the CPWL solution. The number of possible arrangements for $P^c$ affine pieces is $P^c!$, which implies that implementing Equation (15) reduces the feasible region by a factor $P^{+}!P^{-}!$ of its original extent. Additionally, applying Equation (15) introduces $P^{+}+P^{-}-2$ constraints to the MILP problem.

5.3. Imposing $d+1$ Points per Affine Domain of $f^c$

Note that the inequality from Equation (5), which allows each point to belong to more than one affine domain of $f^{+}$ or $f^{-}$, is critical for utilizing Equation (16). If we were instead imposing a single affine domain per point, as in the formulation from Kazda and Li (2021), Equation (16) would become too restrictive. In fact, the combined requirement of having only one affine domain per point and at least $(d+1)$ points per affine domain would imply that $N\ge \max (P^{+},P^{-})(d+1)$, making the MILP problem infeasible if that condition was not satisfied. Using Equation (16) eliminates some non-well-behaved CPWL solutions from the feasible region. However, all well-behaved solutions are preserved. Owing to Theorem 1, it is proven that this does not impact the quality of the feasible CPWL solutions, as any CPWL $\epsilon$-approximation of S has a well-behaved version. Implementing Equation (16) adds $P^{+}+P^{-}$ constraints to the MILP problem.

5.4. Imposing $d+1$ Points per Affine Domain of f

Alternatively, we can impose a set of constraints that further reduces the feasible region to the set of well-behaved CPWL solutions. This approach requires the introduction of two additional sets of variables, ${\beta }_{i,j,k}$ and ${\gamma }_{j,k}$. In the following equations, the indices $i\in \{1,\dots ,N\}$, $j\in \{1,\dots ,P^{+}\}$, $k\in \{1,\dots ,P^{-}\}$, unless otherwise specified.

Equations (17)–(23) involve the use of $P^{+}{P}^{-}(N+1)$ additional variables, $P^{+}{P}^{-}(N+1)$ variable bounds, and $P^{+}{P}^{-}(4N+1)$ additional constraints. These constraints restrict the feasible space to be $\mathit{\text{CPW}}{L}^{*}(S,\epsilon ,P^{+},P^{-})$. Note that the indicator variable ${\gamma }_{j,k}$ is necessary in Equation (21) because the intersection of an affine domain of $f^{+}$ and $f^{-}$ may be empty.

5.5. Tightening the Big-M Parameters

In order to identify tight big-M values, we first need to define several sets. Let A denote the set of all affine functions $D\to \mathbb{R}$. Let ${[S]}^{d+1}$ denote the set of all subsets of S composed of $d+1$ points. Let $s={({\mathit{x}}_{i_k},z_{i_k})}_{k=1,\dots ,d+1}\in {[S]}^{d+1}$. We define the set $B(s,\epsilon )=\{{({\mathit{x}}_{i_k},z_{i_k}+e_k)}_{k=1,\dots ,d+1}:{(e_k)}_{k=1,\dots ,d+1}\in {\{-\epsilon ,+\epsilon \}}^{d+1}\}$. We define the set $A_{\epsilon ,s}=\{g\in A:g$ ε-approximates $s\}$ and $A_{\epsilon ,s}^{*}=\{g\in A:g$ interpolates $s^{\prime }\in B(s,\epsilon )\}$. We also note $A_{\epsilon }(S)={\cup }_{s\in {[S]}^{d+1}}{A}_{\epsilon ,s}$ and $A_{\epsilon }^{*}(S)={\cup }_{s\in {[S]}^{d+1}}{A}_{\epsilon ,s}^{*}$. For convenience, let $\neg c=\{\begin{array}{ll}“-”\hfill & \text{if} c=“+”\hfill \\ “+”\hfill & \text{if} c=“-”\hfill \end{array}$. Finally, given a CPWL solution $f=f^{+}-f^{-}$, we define $f_{j,k}=f_j^{+}-f_k^{-}$.