{'fraction_non_alphanumeric': 0.04162494870742717, 'fraction_numerical': 0.012441526466967582, 'mean_word_length': 4.433514670471774, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inverse Mixed Integer Optimization: Polyhedral Insights and Trust Region Methods

Inverse optimization, determining parameters of an optimization problem that render a given solution optimal, has received increasing attention in recent years. While significant inverse optimization literature exists for convex optimization problems, there have been few advances for discrete problems, despite the ubiquity of applications that fundamentally rely on discrete decision-making. In this paper, we present a new set of theoretical insights and algorithms for the general class of inverse mixed integer linear optimization problems. Specifically, a general characterization of optimality conditions is established and leveraged to design new cutting plane solution algorithms. Through an extensive set of computational experiments, we show that our methods provide substantial improvements over existing methods in solving the largest and most difficult instances to date.


Introduction
Inverse optimization -inferring unobserved parameters of a (forward ) optimization problem that render a given (forward-feasible) solution optimal -has received increasing attention in recent years. Most of the literature has focused on inverse optimization for convex forward optimization problems, allowing the bi-level inverse problem to be reformulated as a singlelevel convex problem using the Karush-Kuhn-Tucker conditions. Examples include linear (e.g., Liu 1996, Ahuja andOrlin 2001), multiobjective linear (e.g. Naghavi et al. 2019), conic (e.g. Iyengar and Kang 2005) and general convex optimization problems (e.g. Xu 2010, Zhang et al. 2010). These models have been applied across a wide range of application domains (see Section 1.1.2) and have also been extended methodologically in several modern directions including estimating model parameters from multiple observed solutions while considering concepts such as statistical consistency (Aswani et al. 2018) or distributional robustness (Esfahani et al. 2018).
In contrast, there have been very few advances in inverse optimization for discrete forward problems, despite the ubiquity of methodological research in discrete optimization and real-world applications that have discrete decisions. Part of the challenge may stem from the fact that the optimality conditions for discrete problems like mixed integer linear optimization (MILO) generally do not lead to computationally tractable solution algorithms like they do for convex problems.
For example, Schaefer (2009) and Lamperski and Schaefer (2015) demonstrate how superadditive duality can be used to develop exact reformulations of bi-level inverse integer and mixed-integer optimization problems, respectively. However, the resulting single-level formulation is an exponentially large linear program that is intractable beyond problems with few variables.
The other main idea in the literature for a general-purpose inverse MILO solution method is a cutting plane algorithm. First proposed in Wang (2009), the idea is to decompose the inverse MILO problem into a master problem and a subproblem. The former is a relaxation of the inverse MILO problem that provides candidates for parameters to be inferred, whereas the latter identifies extreme points of the convex hull of the forward-feasible region that generate cuts to send back to the master problem. While it was shown that small problems could be solved efficiently, this approach does not scale well to larger problems. The computational cost of computing each extreme point is large in general, since it involves solving an instance of the forward MILO problem. The only advance to this method since its development was a proposed heuristic that parallelizes the computation of extreme points (Duan and Wang 2011).
In this paper, we develop a new cutting plane framework for solving inverse MILO problems.
We begin by providing a more general characterization of the optimality conditions of the inverse MILO problem relative to previous literature. We demonstrate that extreme points are sufficient but not necessary for characterizing inverse optimality, and that more parsimonious representations of the problem exist. Using this insight, we develop a new cutting plane method where the cut generation subroutine efficiently identifies interior forward-feasible points using trust regions. These cuts come at significantly lower computational cost compared to cuts generated from extreme points. In addition, we observe that inverse MILO problems can be solved with many fewer cuts compared to the classical cutting plane algorithm.
Our specific contributions are: 1. We present a novel characterization of optimality conditions for inverse MILO problems by introducing the concept of generator sets (Section 2). Generator sets characterize the complete family of reformulations of the MILO problem that preserve the feasible region of the inverse optimization problem. Our main result provides necessary and sufficient conditions for any set to be a generator set. We also show that there exist "small" generator sets, which are attractive from a computational perspective.
2. Motivated by our theoretical insights, we propose a new family of cutting plane methods that restrict the feasible region of the MILO problem using trust regions (Section 3). We also propose computational enhancements that further speed up the algorithm (Section 4). Our algorithms generalize the classical cutting plane approach to solving inverse MILO problems.
3. We propose an extension of the inverse MILO problem to a setting with multiple input data points and illustrate how our cutting plane algorithm can be easily extended to solve this problem (Section 5).
4. Through a comprehensive numerical study over instances drawn from the MIPLIB 2017 benchmark library (Gleixner et al. 2019), we demonstrate that our new cutting plane algorithm significantly outperforms the state-of-the-art algorithm by solving the largest and most difficult problem instances in the literature to date (Sections 6 and 7). We observe that our new cut generation methods can both substantially speed up cut generation and reduce the number of cuts required.

Related literature
Below, we review relevant inverse optimization models and solution methods, as well as the main application areas in which these models can be found.
1.1.1. Models and methods. The classical inverse optimization problem, which is to infer a set of model parameters that render a given feasible solution optimal, has been studied over a wide range of problem settings. The parameters that are inferred can be the cost vector or the constraints. While there are a few studies that focus on estimating constraint parameters (Güler and Hamacher 2010, Birge et al. 2017, Chan and Kaw 2020, the vast majority of papers focus on estimating the cost vector. Most early works focused on deriving algorithms to solve inverse models for particular forward problems (e.g., shortest path, minimum spanning tree), while more recent works consider general-purpose methods for broader classes of forward problems (e.g., mixed integer optimization). Heuberger (2004) provides an overview of the early works whereas the literature discussed in the introduction are examples of more recent works. In this paper, we focus on estimating the cost vector for general mixed integer forward problems, particularly in the form of the classical inverse model. Inverse optimization problems with "noisy" data have also been studied (e.g. Troutt et al. 2006, Keshavarz et al. 2011, Chan et al. 2014. These problems are characterized by the fact that there does not exist a set of (non-trivial) parameters that render a given solution, or a set of solutions, optimal. For these problems, the cost vector is estimated by minimizing various notions of suboptimality (e.g. Bertsimas et al. 2015, Aswani et al. 2018, Esfahani et al. 2018, Babier et al. 2021). Although we focus predominantly on the classical inverse model in the paper, we show that the insights and methods developed can be easily extended to a model for multiple "noisy" data points that resembles the models found for inverse convex optimization.
Finally, we borrow the term "trust region" from the general optimization community, a term which broadly describes the restriction of a particular search space. The use of trust region concepts have appeared in a range of different domains such as nonlinear optimization and, more relevant to our paper, decomposition algorithms. Within the latter domain, trust regions have been applied to the master problems of decomposition models for stochastic programming, and are shown to help generate stronger cuts in this context (Linderoth and Wright 2003, Zverovich et al. 2012, Rahmaniani et al. 2017. In contrast, we apply trust regions to the cut generation problem in our decomposition framework, which leads to improvements in both strength of cuts and cut generation time over a wide array of problem structures.
1.1.2. Applications. Inverse optimization has been studied in a wide variety of applications to estimate latent parameters and infer subjective preferences using observed decision data. Examples of such applications can be found in energy markets (Ruiz et al. 2013, Saez-Gallego et al. 2016, Birge et al. 2017, healthcare (Erkin et al. 2010, Chan et al. 2014, Babier et al. 2020, finance (Bertsimas et al. 2012, Utz et al. 2014, Yu et al. 2020, and transportation (Chow and Recker 2012, Chow et al. 2014, Xu et al. 2018, Zhang et al. 2018, where inverse optimization is used to provide insights into market structure, treatment design, risk aversion and route-choice preferences, respectively. This literature makes use of both classical inverse models and inverse models for "noisy" decision data. Inverse optimization is also prevalent in the bi-level optimization literature, especially in regards to pricing and incentive design problems. Recent examples of such applications include the design of tolls in traffic and commodity transportation networks (Marcotte et al. 2009, Brotcorne et al. 2011, Esfandeh et al. 2016, Kuiteing et al. 2017, price schedules and carbon taxes in energy systems (Zhou et al. 2011, Afşar et al. 2021, and profit-sharing mechanisms in carrier alliances (Agarwal andErgun 2010, Houghtalen et al. 2011). Here, the classical inverse optimization model appears as a "subproblem" within the solution process used to solve these pricing problems. Specifically, a "master problem" generates a decision "target" (e.g., a sustainable routing decision) while an inverse model computes a set of objective perturbations (e.g., tolls on roads) for which the target becomes an optimal decision for an underlying decision-maker (e.g., a transportation company). The master problem can also propose new targets if the inverse subproblem is infeasible (Afşar et al. 2021).
Despite the ubiquity of decision-making models in these application areas that fundamentally rely on making discrete choices, existing literature has focused primarily on inverse optimization for continuous forward problems. The few exceptions to date are in energy planning (Zhou et al. 2011), transportation (Chow and Recker 2012), and sustainability (Turner and Chan 2013). The first two applications were formulated as general inverse MILO problems and applied the classical cutting plane algorithm (Wang 2009). The sustainability application, an inverse knapsack problem, was reformulated using superadditive duality and solved by restricting focus to linear functions.

Notation
Throughout the paper, vectors and matrices are written in bold, while sets are defined using calligraphic letters. Subscripts denote specific elements of a vector, whereas superscripts denote different vectors. We use (·) ⊤ to denote the transpose operator. For any set S, conv(S) is its convex hull, ext(S) is its set of extreme points, and E(S) = ext(conv(S)).

Structure of Inverse MILO Problems
In this section, we study the structure of inverse MILO problems. We provide a novel characterization of inverse-feasibility over general non-convex forward optimization problems. We do this through the definition of a generator set, which characterizes the complete family of reformulations of the forward optimization problem that preserve the feasible region of the inverse optimization problem. In turn, we show how certain feasible region definitions are amenable to more efficient solution methods for the inverse optimization problem. Proofs of our results are provided in the Appendix.

Problem definition
Our inverse optimization model is based on the following MILO problem, known as the forward problem: (1) Let F(c, X ) be the optimal solution set of FP(c, X ). Elements of X and F(c, X ) are called forwardfeasible and forward-optimal, respectively. For a given forward-feasible solutionx ∈ X , let C(x, X ) define the inverse-feasible region, which is the set of all cost vectors c ∈ R n that render the solution x optimal, i.e., C(x, X ) := {c ∈ R n |x ∈ F(c, X )}. (2) The inverse optimization problem is defined as the following optimization model, which seeks to find a cost vector c that minimizes an objective function g(c), is inverse-feasible, and satisfies any context-specific constraints c ∈ P. A common choice of g(c) is The inverse optimization model with this objective seeks to find feasible cost vectors c that minimize the deviation to a given reference cost vector c 0 , as measured by the Manhattan distance.
Practically, this objective function is relevant in many applications (see Section 1.1 for details).
First, in estimation problems, c 0 can represent an initial estimate of the cost vector, and given an observed decisionx, the cost estimates can be refined (minimally) to be compatible withx. For example, c 0 may represent an estimate of network link costs used to describe route preferences (e.g., where initial estimates are based on road distances), and an observed routing decisionx can be used to refine these estimates (Burton and Toint 1992, Ahuja and Orlin 2001, Chen et al. 2021. Second, in bi-level optimization problems, c 0 represents the known cost vector of a decision-maker, and a set of objective perturbations to c 0 can be designed and prescribed such that a decision targetx becomes an optimal choice for the decision-maker. For example, the elements of c 0 may represent the cost of producing energy using different resources, in which case the vector c − c 0 would describe the value of carbon taxes or subsidies that must be imposed on each resource for a low-emission production planx to become optimal for an energy producer (Zhou et al. 2011, Rathore andJakhar 2021). In these settings, there may be context-specific constraints on c, such as non-negativity (P = R n + ), integrality (P = Z n ), variable bounds (P = [l, u] n ), or a transformation into a multi-objective problem We make this modeling choice because the model is well-motivated and because it is the formulation that is studied in all previous works on inverse mixed integer optimization (Wang 2009, Schaefer 2009, Duan and Wang 2011, Lamperski and Schaefer 2015, Bulut and Ralphs 2015. Nonetheless, we emphasize that this modeling choice comes without loss of generality. The theoretical insights and methodological contributions presented in this paper focus exclusively on difficulties posed by the inverse-feasibility constraint (equation (2)), and are thus directly applicable for any inverse MILO problem in the form of model (3). We also assume without loss of generality that the inverse optimization problem is feasible. If it is not, then the methods proposed in this paper can prove its infeasibility.
In the remainder of this section (Sections 2.2 and 2.3), we focus on the core difficulty of solving the inverse MILO problem: characterizing the inverse-feasible region C(x, X ). In contrast to convex optimization problems, where this condition can be transformed into convex constraints using the Karush-Kuhn-Tucker conditions, determining whether a givenx is optimal with respect to some c, let alone the c that optimizes the objective function of the inverse problem, requires solving the forward MILO problem, which is NP-hard in general. Furthermore, attempting to use inverse linear optimization models to approximate inverse MILO models by ignoring the integrality constraints of the forward-feasible region can result in arbitrarily bad approximations. This is further discussed and illustrated in Section EC.1 of the Electronic Companion.

Preliminary remarks
We first note three basic properties of the inverse-feasible region, then describe it using a set of linear constraints.
Remark 1. The inverse-feasible region is nonempty, since 0 ∈ C(x, X ). Remark 3. While the forward-feasible region X does not have to be bounded, the inverse-feasible region C(x, X ) will only contain objective vectors for which the forward problem is bounded. In particular, if c rec is a nonzero vector in the recession cone of conv(X ), then we have c rec / ∈ C. The reason is thatx is assumed to be forward-optimal with a finite objective value (under an unobserved objective function) and inf{(c rec ) ⊤ x | x ∈ X } = −∞. Thusx / ∈ arg min{(c rec ) ⊤ x | x ∈ X }. In other words, no such c rec can makex forward-optimal.
A natural way to explicitly formulate model (4) is to re-write c ∈ C(x, X ) using a potentially infinite set of linear constraints suggesting that a cutting plane method can be used to solve IO(c 0 ,x, X ). Indeed, Wang (2009) proposed a simple cutting plane algorithm to solve the inverse MILO problem in a finite (but exponential, in general) number of iterations by noting that it is sufficient to replace X in equation (5) with the extreme points of its convex hull.
Remark 4. (Wang 2009) The cutting plane algorithm proposed in Wang (2009) starts without constraints (6b), and iteratively adds cuts of this form through the computation of new extreme points x j . The efficiency of this approach depends heavily on the difficulty of generating cuts, which amounts to solving the forward MILO problem, as well as the number of extreme points of the convex hull of X (Bulut and Ralphs 2015). With the exception of parallelizing the search for extreme points (Duan and Wang 2011), the algorithm proposed in Wang (2009) remains the state-of-the-art in solving inverse MILO problems.

Optimality conditions
The main insight from Wang (2009) is that replacing constraints (5) with (6b) preserves the set of optimal solutions to the inverse optimization problem, i.e., C(x, E(X )) = C(x, X ). In this subsection, we characterize the full family of sets G such that C(x, G) = C(x, X ).
Definition 1. Given X , any set G ⊆ R n for which is a generator set for X . If G ⊆ X , then G is a forward-feasible generator set.
Examples of generator sets include conv(X ) and E(X ), with the latter being a forward-feasible generator set.
Next, we demonstrate there may exist a large family of such sets for X . First, note that constraint set (5) defines a polyhedral cone. Second, each half-space c ⊤ (x − x) ≤ 0 induced by a feasible point x ∈ X is not uniquely determined byx − x, but rather by any multiple of the vectorx − x. In other words, the half-space induced by x is the same as the one induced byx + λ(x −x) for any λ > 0.
The same cut can thus be induced by infinitely many other points of R n .
These observations are illustrated in Figure 1. The extreme points {x 1 , ..., x 7 } in Figure 1a induce half spaces for c ∈ R n and the intersection of these half-spaces is the cone shown in Figure 1b.
However, the same cone can be formed using the (blue) diamond points, all of which are interior forward-feasible points, plus x 1 , in Figure 1a.
The last observation above is formalized in the following lemma. Let y(λ,x, x) :=x + λ(x −x), λ > 0, define a point along the ray fromx to x.
Lemma 1 shows that the set of extreme points E(X ) is not necessary to characterize C(x, X ) and can instead be replaced with a set of interior and/or exterior points. This observation can be generalized with a necessary and sufficient condition for any collection of points to be a generator set. Let Y(x, X ) denote the polyhedral cone pointed atx generated by {y(1,x, x)} x∈X .
x 7 x 6 x 5 x 4 x 3 x 2 x 1 x 2 x 4 x 5 x 6 x 3 c 1 While Lemma 1 shows that generator sets can be constructed without using any extreme points, Theorem 1 goes further and shows that far fewer points than the number of extreme points may be sufficient to construct a generator set. As long as the set of extreme rays of Y(x, G) are the same as those of Y(x, X ), then G is a generator set -all other points in G that do not form an extreme ray are redundant in defining the inverse-feasible region. Referring back to Figure 1a, this result implies that any two points y(λ 1 ,x, x 1 ) and y(λ 7 ,x, x 7 ) with λ 1 , λ 7 > 0 are sufficient to define the generator set in that example. In Figure 1b, the inverse-feasible region is exactly the intersection of the half-spaces determined byx − x 1 andx − x 7 .
Note that for any forward-feasible point x ∈ X , the ray {y(λ,x, x)} λ≥0 is in Y(x, X ). Thus, . This is shown by the following corollary.
Analyzing forward-feasible generator sets can be particularly informative for the design of cutting plane algorithms. For example, inverse MILO problems can be solved by the cutting plane algorithm presented in Wang (2009) because E(X ) is a forward-feasible generator set. Below, we present an equivalent characterization of forward-feasible generator sets that does not rely on the notion of rays. This result provides useful intuition to aid in designing a new cutting plane algorithm. Let B(ǫ,x) denote a closed ball of radius ǫ > 0 aroundx.

Theorem 2. A set G ⊆ R n is a forward-feasible generator set if and only if G ⊆ X and there
exists an ǫ > 0 such that Theorem 2 states that a set G ⊆ X is a forward-feasible generator set if and only if the convex hull of the set of points in G and {x} contains an epsilon ball aroundx intersected with the convex hull of X . Going back to Figure 1a, G = {x 1 , x 7 } is a forward-feasible generator set since the convex hull of {x 1 , x 7 ,x} contains the part of an epsilon ball aroundx inside conv(X ), which is shaded in red.
The two key takeaways from Theorem 2 are that a forward-feasible generator set can: (i) consist of many fewer than |E(X )| points, and (ii) be generated from non-extreme points that lie much closer tox.

Cutting Plane Algorithm with Trust Regions
The theoretical insights from the previous section indicate that a new approach to generating cuts may yield a more effective solution algorithm. Instead of generating cuts by identifying extreme points of the convex hull of the forward-feasible region X , we generate cuts by identifying interior points of X , which we accomplish using trust regions. More specifically, we intersect X with a trust region, and generate cuts at interior points of X obtained as extreme points of the convex hull of this restriction of X .
Let T (x, p) be a trust region around the pointx of size p ≥ 1, defined as We simplify notation to T when the parameter p is not the focus of discussion.
Guided by Theorem 2, we define trust regions to be centered atx. The trust region, when imposed onto X , represents an attempt to capture a forward-feasible generator set. The 1-norm in equation (9) is chosen because the trust region has a polynomial number of extreme points. When intersected with X , the resulting restricted forward-feasible region T ∩ X itself is more likely to have a smaller number of extreme points compared to X . This can reduce the number of forward-feasible points that need to be identified in order to build a generator set. Furthermore, the extreme points of this smaller, restricted region may also be computationally easier to identify. These considerations will be discussed in detail below.

Model decomposition and cutting plane framework
The inverse MILO problem is decomposed into a master problem, describing a relaxation of the inverse model (4), and a cut generation subroutine that is used to iteratively tighten the master problem. The framework of our cutting plane algorithm is given in Algorithm 1. A global information set I is used to pass information between successive iterations of the subroutine, which specifically includes the trust region and outer-loop index for our purposes. It is initialized with a trust region of size p 0 centered atx, and the iteration index i = 0.

Algorithm 1 A general cutting plane algorithm for inverse MILO
Input: An inverse MILO problem instance (c 0 ,x, X ), initial trust region size p 0 Output: An inverse-optimal solution c * 1: Solve MP(x,X i ), letc i be its optimal solution 7: Run SUBROUTINE(c i ,x, X , I i ), letx i and I i+1 be its output The master problem is defined below, whereX ⊆ X denotes a finite set of forward-feasible points.
Constraints (10b) ensure thatx is optimal with respect to c overX . As mentioned in Remark 3, the inverse-feasible region C(x, X ) only contains cost vectors for which the forward optimization problem is bounded. However, due to its relaxed nature, the master problem can propose cost vectors that lead to unbounded forward problems in the subroutine. This can be prevented by constraints (10c) and (10d), as proposed in Wang (2009), which are not needed if the forwardfeasible region is bounded. Letc denote an optimal solution to the master problem MP(x,X ).

Cut generation subroutine
Givenc, the cut generation subroutine either finds a feasible pointx ∈ X satisfyingc ⊤ (x −x) > 0 (i.e.,c / ∈ C(x,X ∪ {x})), which generates a cut, or returnsx, verifying thatc ∈ C(x, X ). Note that the verification in the latter case requires solving the forward problem FP(c, X ), which means the cut generation subroutine is at least as hard.
The cut generation subroutine is presented in Subroutine 1. It consists of four key components: the cut generation problem ( Step 3) and three functions to modify the trust region (Steps 2, 5 and 8), which we describe next.

Subroutine 1 Cut generation subroutine
Input: A candidate objectivec, forward-feasible pointx, forward-feasible region X , information Output: A forward-feasible pointx and updated information set I out 1: k ← k + 1 10:

Return to
Step 2 11: else 12: Returnx, I in 3.2.1. Cut generation problem (Step 3). Given a candidate objectivec and a trust region T , the cut generation problem is the forward problem solved over T ∩ X (i.e., FP(c, T ∩ X )), which we refer to as a subregion of X . We note that the choice of 1-norm in the definition of employed trust regions ensures that the cut generation problem remains linear.
If the optimal solution of the cut generation problemx satisfies the conditionc ⊤ (x −x) > 0, thenx is returned; we call such anx a violated forward-feasible point. If no cut is found, then there are two possibilities depending on whether a trust region was used. If no trust region was used, i.e., T k = R n , thenx is returned and Algorithm 1 terminates withc. Otherwise,c renders x optimal over a subregion of X , but not necessarily over X . Thus, larger subregions need to be considered before it can be verified that no additional violated cuts exist, which leads to a trust region update.

UPDATE function (
Step 8). When a subregion T (x, p) ∩ X has been exhausted of violated cuts, we increase the size p of the trust region to δp for some δ > 1: For example, setting δ = 2 doubles the size of the trust region each time UPDATE is called. There is a natural trade-off to be considered when choosing δ. If δ is too small, then the trust regions grow slowly and may lead to the generation of subregions with no violated cuts. This will result in time-consuming, "empty" calls to the cut generation problem, i.e., calls in which the cut generation problem is solved but is unable to generate any violated forward-feasible point. On the other hand, if δ is too large, the computational advantages of solving the cut generation problem over a smaller subregion with fewer extreme points may quickly be lost.

REMOVE function (
Step 2). The REMOVE function periodically removes the existing trust region (by setting T = R n ) at regular intervals in both outer (i.e., master) iterations i and inner (i.e., subroutine) iterations k. In particular, for fixed index values i * > 0 and k * > 0, In other words, the cut generation problem is solved over X every i * iterations in the outer loop and when the inner loop reaches the iteration limit of k * . The REMOVE function guarantees finite convergence of Algorithm 1, irrespective of how the UPDATE function is defined, because an extreme point of X is returned at least every i * iterations in a finite number of (at most k * ) attempts to verify that a candidate objectivec i is inverse-feasible when no more violated cuts exist.
Note that the classical cutting plane algorithm proposed in Wang (2009) is a special case of our trust region-based approach where i * = k * = 1.

SAVE function (
Step 5). The SAVE function is called when a forward-feasible point yielding a violated cut is found by the cut generation problem. If such a point is found during an iteration where the trust region is R n (i.e., when REMOVE is called), the SAVE function saves the previous trust region. On the other hand, if the REMOVE function did not remove the trust region during the current iteration, then the current trust region is saved. This ensures that there always exists a "non-trivial" trust region (T (x, p) < ∞) saved.
Mathematically, the SAVE function is defined as: 3.2.5. Discussion and example. Trust regions can reduce the difficulty of solving the cut generation problem by reducing the size of the forward-feasible region. However, they can also lead to the generation of stronger cuts.
We use Figure 2a to illustrate this latter point. A trust region of size one (p = 1) is imposed on a two-dimensional integer forward-feasible region, creating a subregion containing two extreme points (indicated using blue diamonds), excludingx. Those two extreme points of the subregion have the same cut generation potential as the five extreme points of X labeled in gray. More specifically, the set of cuts generated by the former set of extreme points can dominate those generated by the latter. Once this trust region is exhausted, it must be updated. In Figure 2a, increasing the size of the trust region by a single unit would coverx, the only additional point needed to build a generator set. The trust region approach guarantees that at most three cuts need to be generated, or equivalently three forward-feasible points need to be identified, to solve the inverse MILO problem. On the other hand, the cutting plane algorithm without trust regions, which computes extreme points of X , may require generating many more cuts.
The REMOVE function periodically removes trust regions and provides opportunities for the algorithm to compute points outside of the existing trust region. We consider another example in Figure 2b with a trust region of T (x, 2) over a different two-dimensional integer forward-feasible region. Assuming the subregion X ∩ T (x, 2) has been exhausted of cuts, removal of this trust region would lead to the computation ofx with certainty, as it is the only extreme point of conv(X ) that can still generate a valid cut. This would terminate the algorithm at the next step. On the other hand, if the size of the trust region was increased incrementally each time without removing the trust region, many unnecessary interior points would be computed before a sufficiently large trust region capturesx. Periodic removal of trust regions thus helps to reduce the potential of being overly-conservative, and when used in conjunction with trust regions create an effective cut generation subroutine.

Algorithm Enhancements
In this section, we propose two enhancements to Subroutine 1 that can be implemented independently or together. These two enhancements are introduced to improve cut generation speed, particularly during iterations when the subregion is large.

Early-stop heuristic for the cut generation problem
As defined in Subroutine 1, the cut generation problem is solved to optimality in each iteration.
However, this can be time consuming, particularly during iterations where the trust region is large or removed entirely. Here, we introduce a solution time threshold τ . If the cut generation problem is not solved by τ , but at least one violated forward-feasible point has been found, the feasible point of maximum violation is returned. If no violated forward-feasible point has been found by τ , the cut generation problem will return the first such point found after τ . This enhancement can reduce the cut generation time and can be effective for two reasons: (i) Solving the cut generation problem may not necessarily provide the best cuts anyways, as discussed in Section 2, and (ii) the solver may have already found an optimal violated forward-feasible point in the cut generation problem but has not computed a bound strong enough to prove its optimality. When a standard MILO solver, such as CPLEX and Gurobi, is used to solve the forward problem in Subroutine 1, the early-stop heuristic can be implemented using a callback function. Returnx, I in Subroutine 1 and Subroutine 2 is that when the trust region size (p) is increased, Subroutine 2 first computes a series of low dimensional trust regions of the increased size, exhausts the corresponding subregions of violated forward-feasible points, and then considers the full dimensional trust region.

Stochastic dimensionality reduction of the trust region
More specifically, Subroutine 2 considers the full dimensional trust region T (x, p) only after h * low dimensional trust regions of size p have consecutively failed to generate a cut. The index h tracks the number of consecutive "empty" low dimensional trust regions (step 10), and is reset to zero when a cut is successfully generated within a trust region (step 7).
The trust region updating function is modified accordingly to reflect this process, and is slightly more involved than the original UPDATE function in equation (11). To generate a low dimensional trust region, we randomly sample a set S ⊂ N . Let s denote a value less than |N |, and let RAND(s) denote a function that generates a random subset of N of size s. The value s, which represents the cardinality of the set S, can be either a fixed, predetermined value or a function of the trust region size, i.e., s(p). For now, we assume the latter, and will discuss further in the following paragraph.
The new stochastic (indicated by S-) update function becomes: where δ, the trust region growth rate, is the same as in the original UPDATE function. When h < h * , a new low dimensional trust region is generated. When h = h * , the full dimensional trust region is generated instead. Finally, when the full dimensional trust region fails to identify any violated forward-feasible point, i.e., h = h * + 1, the S-UPDATE function generates a new low dimensional trust region of an increased size.
We conclude by discussing the choice of the function s(p) in equation (14), which defines the cardinality of the set S to be sampled given that the trust region is of size p. While s(p) can be defined simply to return a fixed value (less than |N |) for any input p, we define s(p) to help offset the reduction in computational efficiency associated with larger trust regions. In particular, we define s(p) to be a decreasing function of the trust region size p, i.e., we increasingly reduce the dimensionality of the trust region as its size grows. Our specific choice of the function s(p) is: In equation (15), the parameter κ ∈ (0, 1) denotes the dimensionality reduction rate: when the size of the trust region increases by 1 unit, the dimension of the trust region decreases by κ.
The parameter q ∈ (0, 1) specifies an upper bound on the dimensionality reduction, so that the dimension of the trust region will be no less than qn.

Extension to Multi-Point Inverse MILO
In this section, we present an extension of our base inverse MILO model (4). While the literature surrounding inverse MILO problems has almost exclusively focused on model (4), the literature on inverse convex optimization problems has expanded in the direction of data-driven parameter estimation, where multiple solutions can be used as input. In light of this, we propose a multipoint inverse MILO model, which takes multiple input solutions from potentially different forwardfeasible regions to generate a cost vector. We show how our trust region-based cutting plane algorithms can be directly extended to solve this model.

Problem description
Let D be the index set of D data points,x 1 , . . . ,x D , which are feasible for their respective feasible regions X 1 , . . . , X D . The following model is a natural extension of model (4) to the multi-point case: The data-driven inverse optimization literature has also considered "noisy" data cases for continuous forward problems, where it is not possible to compute a single vector c that is inverse-feasible for each and every decision input. In the light of this, we also propose the following extension of model (16), where we introduce an inverse-feasible cost vectorc d for each data pointx d ∈ X d : The regularization term λ represents a trade-off between making the cost vector c close to the prior c 0 and to the set of inverse-feasible cost vectors {c d } d∈D .

Extended multi-point cutting plane algorithm
We now show how a natural extension of Algorithm 1 can be used to solve models (16) and (17).
The key observation is that each feasible region X d in (16b) and (17b)  The master problem defines a relaxation of model (17), whereX d ⊆ X d is a finite set for each d: With multiple data points, there is a cut generation problem associated with each d ∈ D. Thus, there is a new issue that must be considered, specifically that of a breadth-first versus depth-first search: whether the master problem is called after cuts for each d are generated versus after the first cut for a single d. Let V * ∈ {1, . . . , D} be a parameter that denotes the minimum number of violated forward-feasible points, if any, that must be found until the master problem is called. The larger V * is, the more the algorithm behaves "breadth-first". In the extreme case where V * = D, the master problem is called only when the master problem candidate solution is considered for every d. This may reduce the number of master problem calls, but at the expense of potentially solving many unnecessary cut generation problems that cannot identify violated cuts. For example, this may happen when a generator set has been identified for several forward-feasible regions already but not all of them. In the other extreme where V * = 1, the master problem is called as soon as a violated forward-feasible point is found for any d ∈ D. In the algorithm, count v (count e ) denotes the number of forward-feasible regions for which a violated point has (has not) been found. The algorithm terminates when count e = D, i.e., there does not exist a single violated forward-feasible point in any forward-feasible region.

Algorithm 2 Cutting plane algorithm for multi-point inverse MILO
Input: An inverse MILO problem instance ({x d , X d } d∈D ), initial trust region size p 0 Output: An inverse-optimal solution (c * , {c * d } d∈D )

Experimental Design
In this section, we discuss the setup of our computational experiments, which include the generation of test instances, the different subroutine variants and enhancements considered, and the parameterization of each subroutine. Our computational experiments are focused on inverse MILO problems of the form presented in model (4) with P = R n . We made this modeling choice to focus the presentation and discussion on the key insights gained from applying our new cutting plane algorithms over various problem structures and forward-feasible regions. This model is also studied in all previous computational studies of inverse MILO (Wang 2009, Schaefer 2009, Duan and Wang 2011. To generate comprehensive insights into the performance of our algorithms, we apply them to a diverse set of problem instances, much larger than any set considered in the literature to date.

Test bank generation
To generate a comprehensive set of inverse MILO problems over many different problem structures, we draw from the MIPLIB 2017 benchmark library (Gleixner et al. 2019), which includes 240 MIPLIB problems. We use a subset of these problems to generate a test bank of inverse MILO problems as follows. We first consider all MIPLIB problems with fewer than 12,000 variables and constraints. This reduces our bank to 125 total MIPLIB problems. For each problem, we attempt to generate a feasible solution by solving the problem with a randomly generated cost vector. If the problem is unbounded or cannot be solved within 10 minutes using Gurobi 8.1.0, we try again with a different randomly generated cost vector. If this procedure cannot generate three feasible solutions within ten attempts, the MIPLIB problem is dropped from consideration. Otherwise, three feasible points are generated for each MIPLIB problem. This approach results in 73 MIPLIB problems being included in our test bank, each with three feasible points, totalling 219 inverse MILO problem instances. Each instance is constructed by using one of the feasible solutions asx, the corresponding original cost vector of the MIPLIB problem (i.e., not the random one used to generatex) as c 0 and the corresponding set of MIPLIB problem constraints as X . The 3 instances constructed per MIPLIB problem are labeled using the MIPLIB name and a suffix of t1, t2, and t3 (see Table EC.1).
Of the 73 included MIPLIB problems, 9 have only binary variables, 42 have binary and continuous variables, 14 have binary and integer variables, and 8 have binary, integer and continuous variables.
The number of variables in these problems ranges from 34 to 11,700 with an average of 3,185. The number of constraints ranges from 4 to 10,900 with an average of 2,315. These problems represent applications such as scheduling, production, shipment, assignment, set covering and bin packing.

Subroutine variants
To solve these inverse MILO instances, we implement and compare five different cutting plane algorithms, described below: • CP. This is the classical cutting plane approach proposed in Wang (2009), which serves as the benchmark algorithm.
• CP-ES. This is a simple extension of the CP algorithm in which we embed the early-stop heuristic presented in Section 4.1. We choose a time threshold of 5 seconds, at which time the most violated forward-feasible solution is returned (τ = 5). If no violated points have been found by 5 seconds, the subroutine returns the first one found after the 5-second threshold has been exceeded.
• CPTR. This is the cutting plane algorithm with trust regions, defined in Subroutine 1. The subroutine initializes with a trust region of size one (p 0 = 1 in equation (9)). The trust region doubles in size each time the trust region becomes redundant (δ = 2 in equation (11)). The trust region is removed every 10th cut that is generated and when no violated forward-feasible points can be found within the trust region (i * = 10 and k * = 2 in equation (12)).
• CPTR-ES. This is an extension of the CPTR algorithm in which we embed the early-stop heuristic with a threshold of 5 seconds (τ = 5), similar to the CP-ES algorithm.
• CPTR-ES-DR. This considers the CPTR algorithm with both the early-stop heuristic and the stochastic dimensionality reduction enhancement. We implement a dimensionality reduction rate of 3% (κ = 0.03 in equation (15)). The dimension of the trust region will never be lower than 80% of n (q = 0.8). Finally, we revert to the full dimensional trust region when ten consecutive low-dimensional trust regions are found to have no violated forward-feasible points (h * = 10).
All experiments were coded in Python 3.7 and optimization problems are solved using Gurobi 8.1.0 parameterized with a single thread. The experiments were conducted on a Intel Core i7-4790 processor at 3.60 GHz on a Windows 8.1 Pro. A time constraint of one hour was set for all instances.

Numerical Results
In this section, we present the computational results using the algorithms outlined in Section 6. We begin with a high-level overview of the results, summarized by Figures 3 and 4. We then provide an in-depth study of the advantages and potential trade-offs that exist when using trust regions and the additional enhancements (Section 7.1). In particular, we compare cut generation speed and cut strength across the various algorithms. We also conduct a sensitivity analysis, examining how the results change when several baseline parameters are modified (Section 7.2). We conclude with a brief summary of the main takeaways from our numerical results (Section 7.3). From the performance profiles in Figure 3, we observe that CPTR-ES and CPTR-ES-DR solve significantly more instances than CPTR or CP-ES, which in turn solve many more than CP, the baseline algorithm. The two best algorithms are also significantly faster than the others. For instance, they can each solve the same number of instances as CP (66 instances) and CP-ES (107 instances) in only 2.5% and 12% of the required time, respectively. The main takeaway from these high-level results is that either trust regions or the early-stop heuristic can lead to dramatic improvements in the solution time, compared to CP. When trust regions and early-stopping are used together, the improvements are even more significant.  Before we proceed to a more detailed discussion of these results, we acknowledge that 83 of the instances were not solved by any of the algorithms, including our best-performing ones. Ultimately, the need to solve many MILO problems is an inherent limitation of all these cutting plane algorithms. The MILO problems are used to generate cuts and validate optimality of a candidate cost vector. As these are MIPLIB problems, many of the MILO problems that need to be solved are time consuming relative to our chosen time limit. While trust regions can make certain MILO problems easier to solve, the algorithms eventually need to solve the full MILO problem to verify optimality. And we observe that many of the unsolved instances are precisely the ones in which MILO solution times increase dramatically when larger trust regions are considered.

Analysis and discussion of results
In this section, we analyze how the addition of trust regions and other enhancements lead to the The average per-instance cut generation times within different trust region sizes and with enhancements. Trust region sizes correspond to p in equation (9) and (13). Each data point is the average cut generation time within a trust region size for a particular instance. p = 64+ denotes all finite trust region sizes of p ≥ 64, whereas No TR denotes no trust region (i.e., p = ∞).
7.1.1. Early-stop heuristic. Cut generation times within CP can be very high, as evident from Figure 5. The median of the average per-instance cut generation times exceeds 200 seconds, and CP fails to generate any cut within the time limit in over 25% of the instances. When the early-stop heuristic is applied, the median is reduced to 7 seconds, close to the early-stop τ value of 5 seconds. CP-ES is thus capable of generating many more cuts, as described next.
The iteration count for instances solved by CP-ES is shown in Figure 6. First, the ability to generate cuts faster allows CP-ES to generate many more cuts, which enables it to solve many more instances. Second, we observe that compared to the instances solved by CP, the early-stop heuristic does not significantly increase the iteration count. Part of this reason is that in a number of instances, CP and CP-ES are equivalent because the early-stop feature is never used; many of the instances that can be solved by CP are naturally ones in which cuts can be computed quickly.
Nonetheless, the results in general suggest that cuts generated by the early-stop heuristic with τ = 5 are not significantly weaker. Note that lowering the value of τ will not necessarily result in faster cut generation. For example, there exist many instances in which the average cut generation time is higher than τ (as observed in Figure 5), i.e., it takes longer than τ seconds on average to find any violated forward-feasible point.
Cut generation times can also vary significantly depending on the given cost vector. In particular, many instances with an average cut generation time less than τ may still have a large number of iterations in which cut generation times are higher than τ .
7.1.2. Trust regions. Cut generation times can be substantially reduced using trust regions, as observed in Figure 5. For example, average cut generation time within trust regions of size p ≤ 8 can be orders of magnitudes lower than without trust regions. Like the early-stop heuristic, this reduction in cut generation time allows more cuts to be generated when instances are unsolved by CP, and the improved cut generation speed leads to many more instances being solved, as seen in Figure 7a. Of course, as trust region sizes increase, we observe that the cut generation times also tend to increase. For the instances solved by CP, CPTR solves those same instances in significantly fewer iterations. This observation highlights that the cuts computed by CPTR are stronger. For example, the instances that require over 250 cuts from CP require around an order of magnitude fewer cuts with CPTR. Comparing CPTR-ES to CP-ES, the reduction in iterations is similar, as seen in Figure   7b. In general, the more cuts CP and CP-ES require to solve an instance, the greater the impact of trust regions on reducing iteration count.
The reduction in cut generation time and the increase in strength of cuts make trust regions very effective. Nonetheless, the magnitude of improvement in these two factors varies across instances.
For example, reduction in cut generation time depends on how many cuts are computed within trust regions of smaller sizes. Secondly, a trust region of a particular size may lead to much faster cut generation in one instance than when the same trust region is applied to a different instance. To understand this phenomenon, it is worth noting that trust regions are, by definition, an additional set of constraints on the forward optimization problem. In some MIPLIB instances, the addition of trust regions of any size can "simplify" the forward-feasible region and lead to faster cut generation, whereas in other instances large trust regions may make the problem even more difficult to solve.
In Table 1 we show the results of adding trust regions to CP-ES over a subset of instances (selected from   generally perform better. The parameters of the dimensionality reduction enhancement can be tuned to balance this trade-off, which we discuss in the next subsection.

Sensitivity analysis
In the previous section, we analyzed results for a fixed set of baseline parameters to illustrate that performance gains can be achieved over a wide variety of problem structures without tuning. In this section, we consider several modifications of the baseline parameter values and show that even minor algorithm tuning can lead to dramatic gains. We consider modifications to most parameter values, as shown in Table 2, and report the change in performance over the baseline values. For simplicity, we consider each parameter modification independently, holding all other parameters fixed to their baseline values. We choose a convenience sample of the " t1" instances solved by both the CPTR-ES and CPTR-ES-DR algorithms using baseline parameters, resulting in a total of 42 instances.
The percentage change in total iteration count and solution time for CPTR-ES and CPTR-ES-DR, relative to the baseline parameter values, are shown in Figure 9. The results highlight that,  Table 2 List of considered algorithm parameters for the CPTR algorithm and the two enhancements.
values. However, we also show that it is possible to achieve significantly better performance if we could identify the best parameter values for each instance independently (markers denoted "best").  Comparing Figures 9a and 9b, it appears that solution time is more sensitive to the parameters for trust regions and the early-stop heuristic than the parameters for the dimensionality reduction enhancement. This result is intuitive since trust regions and the early-stop heuristic yield much larger improvements when added to any algorithm (see Figure 4).
Based on the previous results, we derive a few general insights into parameter value choices for the trust regions and early stop heuristic: • Trust region growth rate (δ). A significant gain in total solution time is observed when the trust region growth rate is high, i.e., δ = 4. When the size of trust regions are increased too rapidly, the computational gains from solving smaller cut generation problems are generally reduced.
• Frequency of trust region removal (i * ). The total solution time is higher when i * = 5 than when i * = 10 or i * = 20. This observation suggests that removing trust regions too frequently may reduce the computational gains of faster cut generation, especially when trust regions are small.
• Early-stop value (τ ). Cuts may be weaker when τ is decreased. For example, τ = 1 is the only parameter modification that leads to higher total iteration count. This is consistent with our theory, which suggests that the gain in cut strength when using trust regions comes from the computation of the extreme points of the corresponding subregion. When τ is low, an interior point of the subregion may be returned instead of an extreme point. Nonetheless, we observe through the progression of τ = 20, τ = 5 (baseline) and τ = 1 that the improvement in cut generation speed generally outweighs the increase in iteration count from a lower τ value.

Takeaways from numerical results
The three main takeaways from our numerical results are: i) trust regions improve the strength of cuts and reduce cut generation time, particularly when trust regions are small, ii) the early-stop heuristic further reduces cut generation time, whereas the dimensionality reduction can help but comes with a potential cost of reducing the strength of cuts, iii) default parameter settings seem to work well, but additional tuning of algorithm parameters can lead to significant improvements in solution time.

Conclusion
In this paper, we develop a novel class of cutting plane algorithms for solving inverse MILO problems. Our approach of using trust regions to speed up computation is simple but effective, and is well-supported by our insights on the optimality conditions of inverse MILO problems.
Through extensive computational experiments, we demonstrate that our cut generation techniques and enhancements are highly effective across a large range of inverse MILO problems.

Appendix. Proofs
Several of the proofs in this section use the result that the inverse-feasible region of the sets X , conv(X ), and E(X ) are the same, i.e., C(x, X ) = C(x, conv(X )) = C(x, E(X )). This result was previously established in Wang (2009).
Proof of Theorem 1 We prove by contradiction. (⇐) We assume Y(x, X ) = Y(x, G), but G is not a generator set. This implies that there exists some c such that either (i) c ∈ C(x, X ) but c / ∈ C(x, G) or (ii) By definition, this g cannot be in X , i.e., g / ∈ X , and thus g / ∈ conv(X ).
Furthermore, there cannot exist any x ∈ X , λ > 0 such thatx + λ(x −x) = g. This implies that g / ∈ Y(x, X ), which is a contradiction. The same argument can be applied to show that (ii) leads to a contradiction.
(⇒) We assume that G is in fact a generator set, but Y(x, X ) = Y(x, G). This implies that there exists a point y * such that either (i*) y * ∈ Y(x, X ) but y * / ∈ Y(x, G) or (ii*) y * ∈ Y(x, G) but y * / ∈ Y(x, X ). If (i*) is true, then there exists a seperating hyperplane between Y(x, G) and y * . The normal vector of this hyperplane that is pointing in the direction of Y(x, X ) defines a cost vector c such that all points c ⊤ (x − y) ≤ 0 ∀y ∈ Y(x, G), which implies that c ⊤ (x − g) ≤ 0 ∀g ∈ G. On the other hand, y * lies on the opposite side of the hyperplane, which implies that c(x − y * ) > 0. Finally, because y * ∈ Y(x, X ) and thus must be constructed by a convex combination of the rays {y(λ,x, x)} λ≥0 x ∈ X , there must exist at least one x * ∈ X such that c(x − x * ) > 0. This contradicts the initial assumption that C(x, X ) = C(x, G). The same argument can be applied to show that (ii*) leads to a contradiction.
Proof of Corollary 1 (⇒) Note that by the definition of Y(x, X ), X ⊆ Y(x, X ). If G ⊆ X is a generator set, then Y(x, G) = Y(x, X ) by Theorem 1. Thus, X ⊆ Y(x, G). These statements together yield G ⊆ X ⊆ Y(x, G). (⇐) Suppose X ⊆ Y(x, G) for some G ⊆ X . This implies that for any x ∈ X , x ∈ Y(x, G) and {y(λ,x, x)} λ≥0 ∈ Y(x, G). Therefore, Y(x, X ) ⊆ Y(x, G). Since G ⊆ X , the reverse must also be true, i.e., Y(x, G) ⊆ Y(x, X ). Thus, Y(x, G) = Y(x, X ). Applying Theorem 1, G ⊆ X must be a generator set.
Proof of Theorem 2 (⇒) We first prove the forward direction, i.e., that if G ⊆ X is a generator set, then B ǫ (x) ∩ conv(X ) must be a subset of conv(G ∪ {x}). We prove by contradiction. Suppose that there exists a point y in B ǫ (x) ∩ conv(X ) such that y / ∈ conv(G ∪ {x}). This implies that there exists a separating hyperplane between y and conv(G ∪ {x}). The pointx must lie in only of side of this hyperplane. Taking the normal vector of this hyperplane to be c, we must have that either c ⊤ (x − g) ≤ 0 ∀g ∈ G and c ⊤ (x − y) > 0, or that c ⊤ (x − g) > 0 ∀g ∈ G and c ⊤ (x − y) ≤ 0. Given that y ∈ conv(X ), this contracts the initial assumption that C(x, X ) = C(x, conv(X )) = C(x, G).
We will now prove that both sides of this equation can be simplied to obtain Y(x, X ) ⊆ Y(x, G), at which point we can apply Corollary 1 to show that G ⊆ X is a generator set. The major steps in this proof rely on the definition of the set Y(x, .), and we refer to Theorem 1 for its definition.
These steps lead to the result that Y(x, X ) ⊆ Y(x, G). Since X ⊆ Y(x, X ), and we assumed G ⊆ X , G ⊆ X ⊆ Y(x, G) and by Corollary (1), G must be a forward-feasible generator set. Let X LP denote a relaxed forward-feasible region where all integrality constraints present in X = {x | Ax ≥ b, x ∈ Z n−q × R q } are ignored. Similarly, let c LP ∈ arg min c∈P { c − c 0 : c ∈ C(x, X LP )} denote the solution to an arbitrary inverse linear optimization model where X LP is considered instead of X . Since X ⊆ X LP , the inverse-feasible region C(x, X LP ) must be a subset of C(x, X ), i.e., C(x, X LP ) ⊆ C(x, X ).
The main drawback of using inverse linear optimization over X LP is that the set C(x, X LP ) can be very different from the true inverse-feasible region C(x, X ). This means that c LP can be very far from the (or in fact any) optimal solution c * . Furthermore, there can be a large number of different regions X LP that define the same mixed integer set X when the integrality constraints On the other hand, the linear constraints used to describe this region (x 1 , x 2 ≥ 0, and the two hyperplanes indicated by orange lines) are different, resulting in different X LP and very different inverse-feasible regions for the solutionx = (1, 1). In particular, the inverse-feasible regions of (x, X LP ) become increasingly smaller from Figures (1a) to (1d), which would result in solutions c LP that become increasingly worse relative to c * for any choice of c 0 . This example, with C(x, X LP ) = {0} in Figure (1d), also shows that there is no bound on the distance between c LP and c * (or similarly, between c LP − c 0 1 and c − c 0 1 ). For these reasons, inverse linear optimization models may be poorly suited for solving inverse MILO models.

EC.2. Example of Algorithm Performance over a Single Instance
Below, we provide an example of the performance of the five algorithms over the instance assign1-5-8 t1. In Figure EC.2, the total solution time is plotted against the iteration count, and the slope of the line at any iteration reflects the cut generation speed at that particular iteration. The instance cannot be solved by CP, which has particularly high cut generation times, computing only 3 cuts is not solved by the CP algorithm, which manages only to compute 3 cuts within the time limit.
CP-ES solves the instance in 2061 seconds using 250 iterations (not shown).