Quantile Inverse Optimization: Improving Stability in Inverse Linear Programming

TL;DR

This paper introduces a new inverse linear programming method focused on stability, addressing sensitivity issues in data-driven inverse LP by ensuring solutions are robust to data imperfections and applicable in real-world scenarios.

Contribution

The paper proposes a novel inverse LP approach that guarantees stability under data noise, formulates it as a mixed-integer program, and develops heuristics for efficient solutions.

Findings

Enhanced stability of inverse LP solutions under data noise

Effective application in diet and transshipment problems

Development of heuristics for large-scale problems

Abstract

Inverse linear programming (LP) has received increasing attention due to its potential to generate efficient optimization formulations that can closely replicate the behavior of a complex system. However, inversely inferred parameters and corresponding forward solutions from the existing inverse LP method can be highly sensitive to noise, errors, and uncertainty in the input data, limiting its applicability in data-driven settings. We introduce the notion of inverse and forward stability in inverse LP and propose a novel inverse LP method that determines a set of objective functions that are stable under data imperfection and generate solutions close to the relevant subset of the data. We formulate the inverse model as a mixed-integer program and elucidate its connection to bi-clique problems, which we exploit to develop efficient heuristics. We also show how this method can be used for online learning. We numerically evaluate the stability of the proposed method and demonstrate its practical use in the diet recommendation and transshipment applications.

Figures

Tables

Table 1. Table 1 : MQIO results via MIP solver, 𝐃 ¯ ¯ 𝐃 \bar{\mathbf{D}} -Alg-Exact, 𝐃 ¯ ¯ 𝐃 \bar{\mathbf{D}} -Alg-Heuristic, and 𝐃 ~ ~ 𝐃 \tilde{\mathbf{D}} -Alg-Heuristic.

		Objective Function Value					Time (s)
$(n,m,K)$	$\tau$	MIP	$\overline{𝐃}$-Exact	$\overline{𝐃}$-Heur	$\widetilde{𝐃}$-Heur		MIP	$\overline{𝐃}$-Exact	$\overline{𝐃}$-Heur	$\widetilde{𝐃}$-Heur
(15,100,35)	3.0	*14.00	14.50	11.50	11.00		5000.00	12.42	1.32	0.18
	3.5	15.00	15.00	13.00	10.50		2441.31	26.61	0.72	0.20
	4.0	15.00	15.00	12.00	9.50		28.62	5.77	0.60	0.20
(15,100,200)	3.0	*12.50	12.50	11.00	11.50		5000.00	103.03	2.77	0.90
	3.5	15.00	15.00	14.00	10.50		766.59	77.94	2.35	1.05
	4.0	15.00	15.00	15.00	13.00		426.08	2.18	2.18	0.88
(50,300,35)	3.0	*26.50	37.50	30.50	25.50		5000.00	2000.15	43.26	0.77
	3.5	*28.00	44.50	36.50	27.00		5000.00	3511.62	34.28	0.83
	4.0	*24.50	49.00	41.00	28.50		5000.00	2181.47	26.31	0.85
(50,300,200)	3.0	*0.00	33.00	27.00	20.00		5000.00	1393.53	145.27	4.23
	3.5	*0.00	40.50	29.00	26.00		5000.00	2811.66	182.90	4.75
	4.0	*0.00	47.50	33.00	27.00		5000.00	2705.12	140.94	5.33
(50,300,500)	3.0		27.00	27.00	10.50		5000.00	384.28	384.28	10.23
	3.5		35.00	29.50	13.50		5000.00	775.12	375.80	12.21
	4.0		43.50	37.00	15.50		5000.00	3027.96	186.59	12.84

Equations

$$\text{FO}(\mathbf{c}):\underset{\mathbf{x}}{\mathrm{minimize}}\{{\mathbf{c}}^{\prime }\mathbf{x}|\mathbf{Ax}\ge \mathbf{b}\},$$

$$\underset{\mathbf{c},\{{\mathbf{ϵ}}^k{\}}_{k\in \mathcal{K}},\mathbf{y}}{\mathrm{minimize}}$$

$${\mathbf{A}}^{\prime }\mathbf{y}=\mathbf{c},$$

$$\mathbf{y}\ge 0,$$

$$\mathbf{A}({\hat{\mathbf{x}}}^k-{\mathbf{ϵ}}^k)\ge \mathbf{b},\forall k\in \mathcal{K},$$

$${\mathbf{c}}^{\prime }({\hat{\mathbf{x}}}^k-{\mathbf{ϵ}}^k)={\mathbf{b}}^{\prime }\mathbf{y},\forall k\in \mathcal{K},$$

$$\Vert \mathbf{c}{\Vert }_p=1,$$

$$\underset{\tilde{\mathcal{X}}}{\min }\{d(\hat{\mathcal{X}},\tilde{\mathcal{X}})|\hat{\mathcal{C}}\cap \tilde{\mathcal{C}}=\varnothing \},$$

$$\underset{\mathbf{x}\in {\mathcal{X}}^{\ast }(\hat{\mathbf{c}})}{\max }\{d(\hat{\mathcal{X}},\mathbf{x})\},$$

$$\text{QIO}(\mathcal{K},\tau ,\theta ):\underset{\mathbf{c},\{{\mathbf{ϵ}}^k{\}}_{k\in \mathcal{K}},\mathbf{y},\mathcal{S}}{\mathrm{minimize}}$$

$${\mathbf{A}}^{\prime }\mathbf{y}=\mathbf{c},$$

$$\mathbf{y}\ge 0,$$

$$\mathbf{A}({\hat{\mathbf{x}}}^k-{\mathbf{ϵ}}^k)\ge \mathbf{b},\forall k\in \mathcal{S},$$

$${\mathbf{c}}^{\prime }({\hat{\mathbf{x}}}^k-{\mathbf{ϵ}}^k)={\mathbf{b}}^{\prime }\mathbf{y},\forall k\in \mathcal{S},$$

$$\Vert {\mathbf{ϵ}}^k{\Vert }_{\ell }\le \tau ,\forall k\in \mathcal{S},$$

$$|\mathcal{S}|\ge \theta K,$$

$$\mathcal{S}\subseteq \mathcal{K},$$

$$\Vert \mathbf{c}{\Vert }_p=1.$$

$$\underset{\mathbf{v},\mathbf{u},\{{\mathbf{ϵ}}^k{\}}_{k\in \mathcal{K}}}{\mathrm{minimize}}$$

$$b_i\le {{\mathbf{a}}^i}^{\prime }({\hat{\mathbf{x}}}^k-{\mathbf{ϵ}}^k)\le b_i+M_1(1-v_i),\forall i\in \mathcal{I},\forall k\in \mathcal{K},$$

$$\Vert {\mathbf{ϵ}}^k{\Vert }_{\ell }\le \tau +M_2(1-u_k),\forall k\in \mathcal{K},$$

$$\sum _{k=1}^K{u}_k\ge \theta K,$$

$$v_i,u_k\in \{0,1\},\forall i\in \mathcal{I},\forall k\in \mathcal{K},$$

$$\mathbf{MQIO}(\mathcal{K},\tau,\theta):\quad\underset{\mathbf{v},\mathbf{u},\{\boldsymbol{\epsilon}^{k}\}_{k\in\mathcal{K}}}{\text{maximize}}\;\Big{\{}\,\sum_{i\in\mathcal{I}}v_{i}\,\Big{|}\,\text{\eqref{eq:IO_mip.2}--\eqref{eq:IO_mip.5}}\Big{\}}.$$

$$\underset{\mathbf{v},\mathbf{ϵ}}{\mathrm{maximize}}$$

$$b_i\le {{\mathbf{a}}^i}^{\prime }({\hat{\mathbf{x}}}^k-\mathbf{ϵ})\le b_i+M(1-v_i),\forall i\in \mathcal{I},$$

$$\Vert \mathbf{ϵ}{\Vert }_{\ell }\le \tau ,$$

$$v_i\in \{0,1\},\forall i\in \mathcal{I}.$$

$${\bar{\mathbf{D}}}_{ki}=\{\begin{array}{ll}1,& \text{mbox}ife=(k,i)\in \mathbb{E},\\ 0,& \text{mbox}otherwise.\end{array}$$

$$\underset{\tau,\mathbf{c},\{\boldsymbol{\epsilon}^{k}\}_{k\in\mathcal{K}},\mathbf{y},\mathcal{S}}{\text{minimize}}\big{\{}\tau\ \big{|}\ \text{\eqref{eq:IO1.2}--\eqref{eq:IO1.9}}\big{\}}.$$

$$\underset{\mathbf{x}}{\text{maximize}}\;\Big{\{}\,\lVert\mathbf{x}-\hat{\mathbf{x}}^{k}\rVert_{\infty}\,\Big{|}\,\mathbf{x}\in\mathcal{X}\Big{\}}.$$

$$\underset{\{{\mathbf{x}}^k{\}}_{k\in \mathcal{K}},\{{\mathcal{S}}^i{\}}_{i\in \hat{\mathcal{I}}}}{\mathrm{minimize}}$$

$$D^i({\mathbf{x}}^k)=\underset{\mathbf{z}\in {\mathcal{X}}_i}{\min }\{\Vert {\mathbf{x}}^k-\mathbf{z}{\Vert }_{\ell }\},\forall i\in \hat{\mathcal{I}},k\in \mathcal{K},$$

$$D^i({\mathbf{x}}^k)\ge \tau +\sigma ,\forall k\in {\mathcal{S}}^i,i\in \hat{\mathcal{I}},$$

$$|{\mathcal{S}}^i|\ge \lfloor (1-\theta )K\rfloor +1,i\in \hat{\mathcal{I}}.$$

$$\underset{\mathbf{v},\mathbf{u},\mathbf{c},\boldsymbol{\lambda},\{\boldsymbol{\epsilon}^{k}\}_{k\in\mathcal{K}}}{\text{maximize}}\;\Big{\{}\,\sum_{i\in\mathcal{I}}v_{i}\,\Big{|}\,\mathbf{D}\mathbf{c}\leq\mathbf{d},\,\mathbf{c}=\sum_{i\in\mathcal{I}}\lambda_{i}\mathbf{a}^{i},\,\mathbf{0}\leq\boldsymbol{\lambda}\leq\mathbf{v},\,\text{\eqref{eq:IO_mip.2}--\eqref{eq:IO_mip.5}}\Big{\}},$$

$$\text{Clique(Dˉ,θ)}:\underset{{\mathbf{v}}^{clq},{\mathbf{u}}^{clq}}{\mathrm{maximize}}$$

$$u_k^{clq}+v_i^{clq}\le 1,\forall (k,i)\in \mathcal{K}\times \mathcal{I}\text{ such that }{\bar{\mathbf{D}}}_{ki}=0,$$

$$\sum _{k\in \mathcal{K}}{u}_k^{clq}\ge \theta K,$$

$$v_i^{clq},u_k^{clq}\in \{0,1\},\forall i\in \mathcal{I},\forall k\in \mathcal{K}.$$

$$\underset{\mathbf{\alpha },\mathbf{ϵ}}{\mathrm{minimize}}$$

$$b_i\le {{\mathbf{a}}^i}^{\prime }({\hat{\mathbf{x}}}^k-\mathbf{ϵ})\le b_i+M{\alpha }_i,\forall i\in \mathcal{I},$$

$$\Vert \mathbf{ϵ}{\Vert }_{\ell }\le \tau ,$$

$$\mathbf{\alpha }\ge 0.$$

$${\tilde{\mathbf{D}}}_{ki}=\{\begin{array}{ll}1,& \text{mbox}if{\alpha }_i^k=0,\\ 0,& \text{mbox}otherwise,\end{array}$$

Full text

\DoubleSpacedXI\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough\TITLE

Quantile Inverse Optimization: Improving Stability in Inverse Linear Programming

\RUNAUTHOR

Shahmoradi and Lee

\RUNTITLE

Quantile Inverse Optimization

\ARTICLEAUTHORS\AUTHOR

Zahed Shahmoradi \AFFIndustrial Engineering, University of Houston, 4722 Calhoun Road, Houston, TX 77204, USA, \EMAILzshahmoradi@uh.edu \AUTHORTaewoo Lee \AFFIndustrial Engineering, University of Houston, 4722 Calhoun Road, Houston, TX 77204, USA, \EMAILtlee6@uh.edu

\ABSTRACT

Inverse linear programming (LP) has received increasing attention due to its potential to infer efficient optimization formulations that can closely replicate the behavior of a complex system. However, inversely inferred parameters and corresponding forward solutions from the existing inverse LP methods can be highly sensitive to noise, errors, and uncertainty in the input data, limiting their applicability in data-driven settings. We introduce the notion of inverse and forward stability in inverse LP and propose a novel inverse LP method that determines a set of objective functions that are stable under data imperfection and generate forward solutions close to the relevant subset of the data. We formulate the inverse model as a large-scale mixed-integer program (MIP) and elucidate its connection to biclique problems, which we exploit to develop efficient algorithms that solve much smaller MIPs instead to construct a solution to the original problem. We numerically evaluate the stability of the proposed method and demonstrate its use in the diet recommendation and transshipment applications.

\KEYWORDS

Inverse optimization, Inverse linear programming, Online learning

1 Introduction

Given a set of observed decisions as input data, inverse optimization infers parameters of a “forward” optimization problem, e.g., objective function coefficients, that make the given decisions optimal or near-optimal. By doing so, inverse optimization allows the forward problem to capture the preferences or utilities of the decision maker (DM) and reproduce the decisions accordingly. Inverse optimization has recently found numerous applications including finance (Bertsimas et al. 2012), transportation (Bertsimas et al. 2015), electricity market (Saez-Gallego and Morales 2017, Birge et al. 2017), incentive design (Aswani et al. 2019), and healthcare (Erkin et al. 2010, Lee et al. 2013).

Various inverse optimization models have been developed for different types of forward problems including network optimization (e.g., Ahuja and Orlin (2001), Heuberger (2004)), linear programs (Ahuja and Orlin 2001, Troutt et al. 2008, Chan et al. 2014, Ghate 2015), conic programs (Iyengar and Kang 2005), convex programs (Keshavarz et al. 2011), integer and mixed-integer programs (MIPs) (Schaefer 2009, Wang 2009, Lamperski and Schaefer 2015), and multicriteria optimization (Chan and Lee 2018, Naghavi et al. 2019). The underlying assumption of these studies is that the choice of the inverse model depends on which type of the forward problem is to be used for reproducing the decisions; e.g., if the user wants to derive a linear programming (LP) formulation that can replace a complex decision model, inverse LP is used to fit a linear objective function to the observed decisions.

As the goal of inverse optimization is to fit a model to given data, recent studies have naturally revealed and leveraged the connection between inverse optimization and regression. Aswani et al. (2018) highlight that the residuals used in their inverse model are similar to those in regression and propose a loss function analogous to the sum of squared errors in ordinary least squares. Bertsimas et al. (2015) propose nonparametric inverse optimization motivated by kernel methods in nonparametric regression. Chan et al. (2019) study the similarity between inverse LP and linear regression and propose a goodness-of-fit measure for inverse LP, analogous to the R-squared measure in regression.

Inverse LP has received a particular attention among different types of inverse problems due to its potential to find an LP formulation—whether the true decision-generating system is linear or not—which is easy to solve while closely replicating the system. While inverse optimization for more complex (e.g., non-convex) forward problems might lead to a better fit, the resulting forward problems can be computationally burdensome. Inverse LP has also been used for inferring LP formulations customized to different DMs or problem instances, which can replace a complex formulation in a personalized or distributed manner, e.g., different LP cancer therapy planning formulations inferred for different types of patients, thus enabling personalized treatment modeling (Boutilier et al. 2015).

However, such potential of inverse LP often does not translate well into situations where the observed decisions are subject to noise, measurement errors, and uncertainty. In particular, objective functions (or cost vectors) obtained by the existing inverse LP methods can be very sensitive to small changes in the data. For example, suppose that the DM’s true decisions turn out slightly different from the observed ones due to measurement errors, or some of them are simply outliers due to the DM’s inconsistent behavior. The cost vector inferred by the existing methods often changes substantially in response to such errors and outliers, leading to substantially different forward solutions; also, the forward solutions are often far removed from the observed decisions, making it hard to predict the performance of the inferred cost vector (see Section 2 for detailed illustration). Such instability can limit the applicability of inverse LP in data-driven settings.

In this study, we propose a novel inverse LP framework that addresses the instability issues caused by noise, errors, and uncertainty in data. We formally define stability measures in inverse optimization and propose a new inverse model that improves on the previous methods. Inspired by least quantile linear regression, our model infers a cost vector based on a quantile statistic of optimality errors associated with the observed decisions. Furthermore, we aim to find a set of cost vectors that are guaranteed to generate solutions within a specified distance from the relevant subset of the data. The presentation of our method in this paper focuses on the setting where the constraint matrix of the forward problem remains the same while observed decisions may vary. This setting can find various real world applications, e.g., diet problems where nutritional factors for each food do not change and network optimization problems where the geographical configuration (i.e., nodes and arcs) remains the same. We also discuss how this method can be extended when such assumptions do not hold.

1.1 Relevant Literature

Recent studies in inverse optimization focus on data-driven settings where a large (potentially noisy) dataset collected over a period of time or from many DMs is used for inferring the objective functions. Keshavarz et al. (2011) formulate an inverse model based on relaxed KKT conditions and impute a convex objective function that minimizes the KKT residuals with respect to the input solutions. Bertsimas et al. (2015) consider inverse variational inequality with noisy data and find model parameters that minimize the optimality gap associated with the data. Similarly, Esfahani et al. (2018) develop a distributionally robust inverse optimization model to infer an objective function from noisy data. Aswani et al. (2018) introduce the notion of risk consistency in inverse optimization and propose a model that finds an objective function that replicates the data in a statistically consistent manner. Inverse optimization has also been used for online learning where the inferred objective function is updated adaptively as new data are available over time (Bärmann et al. 2017, Dong et al. 2018). Although these inverse convex programming frameworks can be specialized to inverse LPs, they are rather focused on generic convex programs with assumptions that preclude the above-mentioned instability issues (e.g., strictly convex feasible regions) and thus are not designed to address these particular issues in inverse LP.

Although several recent inverse LP methods accommodate data that may not be optimal, these works focus on the development of closed-form solutions or efficient algorithms under often limiting assumptions (Chan et al. 2014, Babier et al. 2018, Ghobadi et al. 2018). As a result, how the obtained cost vectors actually work in the presence of data imperfection has received little attention. In particular, the previous studies exploit the polyhedral nature of the underlying forward LP and develop algorithms that find a cost vector orthogonal to one of the constraint vectors. While such algorithms are efficient, the resulting cost vector is highly sensitive to outliers or small data shifts.

Our modeling approach to address the instability issues follows in spirit the line of work in the robust linear regression literature. Robust regression aims to infer a model that is stable against outliers or data shifts by increasing stability measure of the fitted model, known as the breakdown point, which is defined as the fraction of the data that can be altered arbitrarily without making the model arbitrarily bad. Various robust regression techniques have been proposed, including least median of squares (Rousseeuw 1984), least trimmed squares (Hubert et al. 2008, Rousseeuw and Van Driessen 2006), and least quantile of squares methods (Bertsimas and Mazumder 2014). In these techniques, a parameter indicating the fraction of the data to be “trimmed” is pre-specified based on the user’s application-specific knowledge and preferences, and model-fitting is done by minimizing a certain quantile error statistic, excluding the residuals associated with data points that are deemed to be outliers.

1.2 Contributions and Organization of the Paper

Our contributions are as follows.

We introduce the notion of forward and inverse stability in inverse optimization and show that the previous inverse LP methods are often unstable in the face of data noise, uncertainties, and outliers. We then propose a new inverse LP method that improves on the previous methods in terms of inverse and forward stability measures and provides bounds for both measures. 2. 2.

We formulate the new inverse LP model as a large-size MIP and formally characterize the set of all feasible cost vectors. By establishing a new connection between the inverse model and a class of biclique problems, we develop efficient algorithms that solve many smaller MIPs (often in parallel) instead of directly solving the large MIP. We show that one of the algorithms is exact with a condition that is easy to check, while other heuristics are useful for very large instances. 3. 3.

We demonstrate two settings where the proposed inverse model is relevant: (i) the user observes decisions from an unknown system and aims to infer an LP formulation that can make new decisions similar to the observations; (ii) the user provides an LP formulation to the DMs as a signal and collects their decisions as a response, which are then used for inferring their objective functions for the LP. For the former, we demonstrate our method in the diet recommendation problem to quantify an individual’s diet preferences from noisy and inconsistent data. For the latter, we use the transshipment problem to infer costs perceived by the DMs from a sequence of noisy datasets.

The rest of the paper is organized as follows. In Section 2, we illustrate how the previous inverse LP can be sensitive to data noise and outliers and formalize the notion of stability in inverse optimization. In Section 3, we present the new inverse LP model that addresses the instability issues as well as its MIP reformulation. In Section 4, we discuss the connection between the proposed model and biclique problems and propose efficient solution approaches. In Section 5, we examine the performance of our model in terms of stability using various LP instances and demonstrate its use in the diet recommendation and transshipment problems. We conclude in Section 6. Unless otherwise stated, proofs are in the appendix.

Throughout this paper, vectors are column vectors and ${\mathbf{a}}^{\prime}$ denotes the transpose of vector $\mathbf{a}$, $\mathbf{e}$ denotes the vector of ones, cone$(\cdot)$ denotes the set of conic combinations of given vectors, and cone${}_{+}(\cdot)$ denotes the set of strict conic combinations of given vectors.

2 Preliminaries

2.1 Forward Linear Program

We consider the following forward optimization (FO) problem:

[TABLE]

where $\mathbf{c}\in\mathbb{R}^{n},\mathbf{x}\in\mathbb{R}^{n}$, $\mathbf{A}\in\mathbb{R}^{m\times n}$, and $\mathbf{b}\in\mathbb{R}^{m}$. Let $\mathcal{I}=\{1,\ldots,m\}$ index the constraints, $\mathcal{J}=\{1,\ldots,n\}$ index the variables, and $\mathbf{a}^{i}\in\mathbb{R}^{n}$ be a (column) vector corresponding to the $i$-th row of $\mathbf{A}$. Let $\mathcal{X}$ be the set of feasible solutions for the FO problem, assumed bounded, full-dimensional and free of redundant constraints, and $\mathcal{X}_{i}=\{\mathbf{x}\in\mathcal{X}\,|\,{\mathbf{a}^{i}}^{\prime}\mathbf{x}=b_{i}\}$, $i\in\mathcal{I}$. Without loss of generality, we assume that $\mathbf{a}^{i}$ for each $i\in\mathcal{I}$ is normalized a priori such that $\|\mathbf{a}^{i}\|_{p}=1$ for some $p\geq 1$.

2.2 Previous Inverse Linear Programming Method

Consider a set of $K$ data points (or observations) $\hat{\mathcal{X}}=\{\mathbf{\hat{x}}^{1},\ldots,\mathbf{\hat{x}}^{K}\}$ with the index set $\mathcal{K}=\{1,\ldots,K\}$. We make no assumption on the feasibility or optimality of the observations for the forward problem (1). Previous inverse LP methods with sub-optimal or infeasible observations aim to find a $\mathbf{c}$ vector that can generate a forward optimal solution that is closest to the observations $\{\mathbf{\hat{x}}^{k}\}_{k\in\mathcal{K}}$ under some distance metric, thus closely “reproducing” the data (Bertsimas et al. 2015, Aswani et al. 2018, Chan et al. 2019). In LP, using general $\ell$-norm as a distance metric, such an inverse model can be written as follows (Aswani et al. 2018, Babier et al. 2018):

[TABLE]

where $\ell\geq 1$, $p\geq 1$, and $\boldsymbol{\epsilon}^{k}$ denotes a perturbation vector for observation $\hat{\mathbf{x}}^{k}$. Given observations $\{\hat{\mathbf{x}}^{k}\}_{k\in\mathcal{K}}$, the above problem finds a $\mathbf{c}$ vector that makes each perturbed solution $\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k}$ satisfy dual feasibility (2b)–(2c), primal feasibility (2d), and strong duality (2e), thus rendering it optimal, while the perturbations (i.e., “optimality errors” associated with the observations) are minimized in $\ell$-norm in the objective function (2a). The normalization constraint (2f) prevents a trivial, all-zero cost vector from being feasible.

Note that existing inverse convex programming models (e.g., Aswani et al. (2018)) can be written equivalently as (2) when the underlying forward problem is exactly (1). While the above problem is non-convex, recent works propose an efficient, exact algorithm for the problem by exploiting the solution structure (Babier et al. 2018, Chan et al. 2019), namely that an optimal cost vector $\mathbf{c}^{*}$ is orthogonal to one of the hyperplanes defining polyhedron $\mathcal{X}$ (i.e., $\mathbf{c}^{*}=\mathbf{a}^{i}$ for some $i\in\mathcal{I}$).

2.3 Instability Issues in Inverse LP

When the dataset $\hat{\mathcal{X}}$ contains noise or measurement errors, the previous inverse model (2) can be unstable in terms of both the cost vector $\mathbf{c}^{*}$ it produces as well as the forward solution that $\mathbf{c}^{*}$ generates (i.e., $\mathbf{x}^{*}\in\argmin\textbf{FO}(\mathbf{c}^{*})$). The following examples illustrate such instability issues.

Example 2.1

Consider the forward LP: $\underset{\mathbf{x}}{\text{minimize}}\;\{\mathbf{c}^{\prime}\mathbf{x}\,|\,\mathbf{A}\mathbf{x}\geq\mathbf{b}\}$ where $\mathbf{a}^{1}=\begin{bmatrix}0\\ -1\end{bmatrix},\;\mathbf{a}^{2}=\begin{bmatrix}-1\\ 0\end{bmatrix},\;\mathbf{a}^{3}=\begin{bmatrix}0\\ 1\end{bmatrix},$ and $\mathbf{a}^{4}=\begin{bmatrix}1\\ 0\end{bmatrix}$ are the rows of $\mathbf{A}$ (written as column vectors) and $\mathbf{b}=[-2.5,-2.5,0,0]^{\prime}$, with four data points $\hat{\mathbf{x}}^{1}=\begin{bmatrix}2\\ 2.3\end{bmatrix}$, $\hat{\mathbf{x}}^{2}=\begin{bmatrix}2.2\\ 2.3\end{bmatrix}$, $\hat{\mathbf{x}}^{3}=\begin{bmatrix}2.2\\ 2\end{bmatrix}$, and $\hat{\mathbf{x}}^{4}=\begin{bmatrix}2\\ 2\end{bmatrix}$ (see Figure 1(a)). The previous inverse model (2) with $\lVert\cdot\rVert_{\ell}=\lVert\cdot\rVert_{\infty}$ finds $\mathbf{c}^{*}=\mathbf{a}^{1}=\begin{bmatrix}0\\ -1\end{bmatrix}$ as the unique optimal cost vector with the objective value of $1.4$ ($\boldsymbol{\epsilon}^{1}=\begin{bmatrix}0\\ -0.2\end{bmatrix}$, $\boldsymbol{\epsilon}^{2}=\begin{bmatrix}0\\ -0.2\end{bmatrix}$, $\boldsymbol{\epsilon}^{3}=\begin{bmatrix}0\\ -0.5\end{bmatrix}$, and $\boldsymbol{\epsilon}^{4}=\begin{bmatrix}0\\ -0.5\end{bmatrix}$). Suppose that $\hat{\mathbf{x}}^{2}$ and $\hat{\mathbf{x}}^{3}$ are shifted to $\begin{bmatrix}2.3\\ 2.2\end{bmatrix}$ and $\begin{bmatrix}2.3\\ 1.9\end{bmatrix}$, respectively (shown in gray). Even with these small shifts, the optimal cost vector changes from $\mathbf{c}^{*}=\begin{bmatrix}0\\ -1\end{bmatrix}$ (which now leads to the objective value of 1.6) to a substantially different vector $\mathbf{c}^{*}=\mathbf{a}^{2}=\begin{bmatrix}-1\\ 0\end{bmatrix}$ with the objective value of 1.4 ($\boldsymbol{\epsilon}^{1}=\begin{bmatrix}-0.5\\ 0\end{bmatrix}$, $\boldsymbol{\epsilon}^{2}=\begin{bmatrix}-0.2\\ 0\end{bmatrix}$, $\boldsymbol{\epsilon}^{3}=\begin{bmatrix}-0.2\\ 0\end{bmatrix}$, and $\boldsymbol{\epsilon}^{4}=\begin{bmatrix}-0.5\\ 0\end{bmatrix}$).

Example 2.2

Consider the same initial data points from Example 2.1. Suppose that an outlier $\hat{\mathbf{x}}^{5}=\begin{bmatrix}2.2\\ 0.3\end{bmatrix}$ is introduced (see Figure 1(b)). The optimal cost vector then changes from $\mathbf{c}^{*}=\mathbf{a}^{1}=\begin{bmatrix}0\\ -1\end{bmatrix}$ to a substantially different vector $\mathbf{c}^{*}=\mathbf{a}^{2}=\begin{bmatrix}-1\\ 0\end{bmatrix}$ with the objective value of $1.9$ ($\boldsymbol{\epsilon}^{1}=\begin{bmatrix}-0.5\\ 0\end{bmatrix}$, $\boldsymbol{\epsilon}^{2}=\begin{bmatrix}-0.3\\ 0\end{bmatrix}$, $\boldsymbol{\epsilon}^{3}=\begin{bmatrix}-0.3\\ 0\end{bmatrix}$, $\boldsymbol{\epsilon}^{4}=\begin{bmatrix}-0.5\\ 0\end{bmatrix}$, and $\boldsymbol{\epsilon}^{5}=\begin{bmatrix}-0.3\\ 0\end{bmatrix}$). The previous solution now has the objective value of $1.4\,+\,\rVert\begin{bmatrix}0\\ -2.2\end{bmatrix}\lVert_{\infty}=3.6$, hence no longer optimal.

Example 2.3

Consider the five data points from Example 2.2 where the optimal cost vector is $\mathbf{c}^{*}=\mathbf{a}^{2}=\begin{bmatrix}-1\\ 0\end{bmatrix}$. While a desirable forward solution is supposed to be close to the relevant subset of the data without the outlier, Figure 1(c) shows that using this cost vector for the forward LP can lead to a solution $\mathbf{x}^{*}$ that is far from the majority of the dataset.

Examples 2.1 and 2.2 show that the previous inverse model is sensitive to outliers or small data shifts, which we refer to as being inverse-unstable. On the other hand, Example 2.3 shows that solving the forward problem with a cost vector from the previous model can lead to a solution that is (unexpectedly) far from the relevant subset of the data, which we refer to as being forward-unstable.

We note that these instability issues are common in general inverse LP settings. For example, if a data-generating system produces errors 10% of the time, the inferred cost vector can still be heavily dragged toward these errors spread far from the majority of the data, regardless of the size of the dataset, leading to a cost vector that would have not been chosen if only the pristine data (the remaining 90%) was used. Also, when small shifts occur in multiple data points according to the same distribution (e.g., the DM’s behavior changing gradually over time), the cost vector can also suddenly change to a substantially different vector, leading to a substantially different forward solution. In general, these instability issues depend heavily on the spatial distribution of the data in relation to polyhedron $\mathcal{X}$, e.g., the distance between each data point and each facet of $\mathcal{X}$.

2.4 Stability Measures for Inverse LP

In this subsection, we formally define a notion of stability in inverse LP and propose measures that we use to assess the stability of an inverse LP model. Given a dataset $\hat{\mathcal{X}}$ and a certain inverse LP model, let $\hat{\mathcal{C}}$ be the set of cost vectors obtained by the model.

•** Inverse stability:** We propose to measure inverse stability of an inverse LP model (or solutions to the model) by the minimum tolerable data shift in $\hat{\mathcal{X}}$ until the model loses all of its initial solutions $\hat{\mathcal{C}}$. That is, if we let $\tilde{\mathcal{X}}$ denote the shifted data and $\tilde{\mathcal{C}}$ denote the set of cost vectors obtained by the model with $\tilde{\mathcal{X}}$, inverse stability is measured by

[TABLE]

where $d$ is some distance function; e.g., $d(\hat{\mathcal{X}},\tilde{\mathcal{X}})=\max\{\|\hat{\mathbf{x}}-\tilde{\mathbf{x}}\|_{\ell}\,|\,\hat{\mathbf{x}}\in\hat{\mathcal{X}},\tilde{\mathbf{x}}\in\tilde{\mathcal{X}}\}$ or $d(\hat{\mathcal{X}},\tilde{\mathcal{X}})=\sum_{k\in\mathcal{K}}\|\hat{\mathbf{x}}^{k}-\tilde{\mathbf{x}}^{k}\|_{\ell},\ell\geq 1,$ if $\hat{\mathcal{X}}$ and $\tilde{\mathcal{X}}$ are in one-to-one correspondence. That is, ideally, a stable inverse model should maintain some cost vectors when reasonably small changes occur in the data, as quantified by the above measure. Note that this measure is analogous to the stability measure in regression, i.e., the finite-sample breakdown point, defined as the fraction of the data that can be altered without spoiling the inferred model completely (Yohai 1987).

•** Forward stability:** Suppose we select $\hat{\mathbf{c}}\in\hat{\mathcal{C}}$ and find a set of forward optimal solutions $\mathcal{X}^{*}(\hat{\mathbf{c}})=\argmin\,\textrm{{FO}}(\hat{\mathbf{c}})$. How unstable this cost vector $\hat{\mathbf{c}}$ can be is assessed by how far a forward solution $\mathbf{x}\in\mathcal{X}^{*}(\hat{\mathbf{c}})$ can be from the given observations $\hat{\mathcal{X}}$, i.e., the worst-case distance between $\mathcal{X}^{*}(\hat{\mathbf{c}})$ and $\hat{\mathcal{X}}$:

[TABLE]

where the distance function $d$ can be defined similarly as above. That is, to guarantee forward stability in inverse LP, an inverse model should identify a cost vector that is guaranteed to produce an optimal solution that is reasonably close to the observations.

The inverse and forward stability measures often conflict with each other. For example, a cost vector whose elements are all zeros is perfectly inverse-stable, yet using this vector for the forward problem likely leads to a solution far from the given data (i.e., highly forward-unstable). Also, one might attempt to include as many cost vectors as possible in $\hat{\mathcal{C}}$ just to make it less vulnerable to data shifts (i.e., more inverse-stable), yet which of these vectors is truly forward-stable then becomes less clear. Another challenge is that neither of the above measures lends itself to a tractable inverse optimization problem. In the next section, we provide a tractable inverse LP framework that improves both inverse and forward stability as quantified by the above measures.

3 Models

In this section, we first propose a new inverse LP model that improves the inverse stability of the previous model (2). We then analyze the solution structure of the model, which we exploit to reformulate the model as an MIP and characterize the set of all feasible cost vectors. Finally, we introduce a specific objective function for the MIP which can serve as a proxy for the forward stability measure and leads to a cost vector with improved forward stability.

3.1 Inverse LP for Quantile Statistics: Improving Inverse Stability

As discussed, inverse instability is caused by the vulnerability of inferred cost vectors to data shifts or outliers. Our modeling strategy to address this is inspired by the close relationship between the least squares method in regression and inverse optimization: the previous inverse model (2) is similar to the least squares method in that it fits the LP model by minimizing the sum of optimality errors (Chan et al. 2019). In regression, one way to address data imperfection is to minimize the $\theta$ quantile statistic where $\theta\in[0,1]$ (i.e., ($\theta\times 100$)-th percentile) instead of the sum of the squared errors, which is known as the least quantile of squares method (Koenker and Hallock 2001, Bertsimas and Mazumder 2014) or least trimmed quantile regression (Rousseeuw 1984, Rousseeuw and Van Driessen 2006).

We adopt this idea and propose a generalized inverse framework where a cost vector is inferred based on the $\theta$ quantile optimality error, i.e., associated with a data point that induces the $\lceil\theta K\rceil$-th smallest optimality error with respect to the cost vector. Instead of minimizing this $\theta$ quantile error as was done in Bertsimas and Mazumder (2014), we find a set of cost vectors such that the $\lceil\theta K\rceil$-th smallest optimality error is no greater than a certain threshold $\tau$. That is, such cost vectors render a relevant subset of the data (of cardinality $\lceil\theta K\rceil$) within the threshold optimality error. This problem, which we call the quantile inverse optimization (QIO) problem, can be written as follows:

[TABLE]

Constraints (5b)–(5c) represent dual feasibility and (5d)–(5e) enforce primal feasibility and strong duality associated with the perturbed solutions $\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k}$ for a subset of the data $\mathcal{S}\subseteq\mathcal{K}$ whose cardinality is enforced to be no less than $\theta K$ by (5g). Constraint (5f) then ensures that optimality error for each chosen observation is within $\tau$. Clearly, these constraints ensure that the $\lceil\theta K\rceil$-th smallest optimality error is no greater than $\tau$. Note that this problem is written as a feasibility problem—we later introduce an objective function that leads to a specific subset of the feasible cost vectors (see Section 3.2). We call cost vectors that are feasible for (5b)–(5i) “inverse-feasible” cost vectors. Similar to quantile-based regression, the choice of $\theta$ depends on the application context or user preference; e.g., $\theta$ can be chosen based on the user’s belief about the fraction of outliers and can also be adjusted post-hoc depending on the result from the model. Or, it can be set to 50% if one is interested in the median error-based inverse problem. The above problem will be referred to as $\textbf{QIO}(\mathcal{K},\tau,\theta)$ or $\textbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$ interchangeably, depending on which is more convenient for the context.

A strong advantage of this modeling framework is that it generalizes the previous model (2)—given the same dataset $\hat{\mathcal{X}}$, there exists $\tau$ such that $\textbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$ with $\theta=100\%$ produces the same set of feasible cost vectors as (2). Furthermore, given such $\tau$, as $\theta$ decreases (i.e., the feasible region expands), the set of feasible cost vectors from $\textbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$ contains the cost vectors from the previous model. That is, in the presence of data shift or additional outliers, $\textbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$ does not lose all of its initial inverse-feasible solutions until the previous model does; i.e., it is at least as stable as the previous model as quantified by inverse stability measure (3). We formalize this idea as follows.

Proposition 3.1

Let $\hat{\mathcal{C}}(\theta)$ be the set of inverse-feasible cost vectors for $\mathbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$. Given $\hat{\mathcal{X}}$ and $\tau$, inverse stability measure (3) for $\hat{\mathcal{C}}(\theta)$ is non-decreasing as $\theta$ decreases.

We later further show that a solution to $\textbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$ provides a lower bound on the inverse stability measure. In the following, we first analyze the solution structure of $\textbf{QIO}(\hat{\mathcal{X}},\tau,\theta)$.

3.1.1 Solution Structure and MIP Reformulation.

The QIO problem is non-convex due to constraint (5e) and thus is hard to solve. We analyze the solution structure of this problem, which leads to an MIP reformulation and a lower bound on the inverse stability measure.

Proposition 3.2

If $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ is feasible, then there exists a feasible ${\mathbf{c}}$ for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ such that ${\mathbf{c}}=\mathbf{a}^{i}$ for some $i\in\mathcal{I}$.

Proposition 3.2 implies that the feasibility of the QIO model can be checked by evaluating at most $m$ constraint vectors, i.e., $\mathbf{a}^{i},i=1,\ldots,m$.

For example, revisiting Example 2.1 (Figure 1(a)), for the QIO problem with $\|\cdot\|_{\ell}=\|\cdot\|_{\infty}$, $\theta=0.5$, and $\tau=0.3$, both $\mathbf{c}=\mathbf{a}^{1}$ and $\mathbf{c}=\mathbf{a}^{2}$ are feasible: $\mathbf{c}\!=\!\mathbf{a}^{1}$ leads to $\|\boldsymbol{\epsilon}^{1}\|_{\infty}\!=\!\|\boldsymbol{\epsilon}^{2}\|_{\infty}\!=\!0.2$ and $\mathbf{c}\!=\!\mathbf{a}^{2}$ leads to $\|\boldsymbol{\epsilon}^{2}\|_{\infty}\!=\!\|\boldsymbol{\epsilon}^{3}\|_{\infty}\!=\!0.3$, both satisfying the threshold $\tau=0.3$.

Proposition 3.2 also suggests that if there is no $\mathbf{c}=\mathbf{a}^{i}$ feasible for (5) for any $i\in\mathcal{I}$, then there exists no cost vector that can make at least $\lceil\theta K\rceil$ observations within $\tau$-optimality. In this case, users could decrease $\theta$ or increase $\tau$ to make the model feasible. Or, $\tau$ can be set by adding a reasonable margin to the minimum possible value of $\tau$ that keeps the model feasible (see Appendix 7.1 for more details). With an appropriate threshold $\tau$, the set of inverse-feasible cost vectors that are orthogonal to some hyperplanes defining $\mathcal{X}$ can be found efficiently by evaluating each hyperplane $i$ (see Algorithm 1).

\OneAndAHalfSpacedXI

Note that Algorithm 1 returns the set of $\mathbf{a}^{i}$’s that are inverse-feasible, i.e., $\mathcal{A}$, which is only a subset of the set of all inverse-feasible solutions.

We extend the results in Proposition 3.2 and Algorithm 1 and derive the following result, which further identifies other inverse-feasible cost vectors.

Lemma 3.3

If there exists $\bar{\mathbf{x}}^{k}\in\cap_{i\in\bar{\mathcal{A}}}\mathcal{X}_{i}$ for some set $\bar{\mathcal{A}}\subseteq\!\mathcal{A}$ such that $\|\bar{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}\!\leq\tau$ for all $k\in{\bar{\mathcal{S}}}\subseteq\mathcal{K}$ where $|{\bar{\mathcal{S}}}|\geq\theta K$, then $\mathbf{c}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$ such that $\|\mathbf{c}\|_{p}=1$ is inverse-feasible.

Lemma 3.3 suggests that we can check each subset $\bar{\mathcal{A}}$ of $\mathcal{A}$ to see if a conic combination of $\mathbf{a}^{i}$’s for $i\in\bar{\mathcal{A}}$ is also inverse-feasible. Using this observation, we propose the following MIP:

[TABLE]

where $v_{i}=1$ if ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k})\geq b_{i}$ holds with equality (i.e., $\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k}\in\mathcal{X}_{i}$) for all $k\in\mathcal{K}$ and $v_{i}=0$ otherwise, and $u_{k}=1$ if observation $\hat{\mathbf{x}}^{k}$ is “chosen” and $u_{k}=0$ otherwise. Parameters $M_{1}$ and $M_{2}$ denote sufficiently large constants; Appendix 7.2 shows how to find appropriate values for them.

In formulation (6), inverse-feasible $\mathbf{a}^{i}$’s are identified by letting $v_{i}=1$; the hypothesis of Lemma 3.3 is then explicitly written as constraints to further identify a conic combination of the chosen $\mathbf{a}^{i}$’s that is also inverse-feasible (the set $\mathcal{I}$ in (6b) and (6e) can be replaced by $\mathcal{A}$ without losing any feasible solution because $v_{i}$ can be 1 only for $i\in\mathcal{A}\subseteq\mathcal{I}$ anyways; we use $\mathcal{I}$ in (6) and subsequent formulations for notational clarity).

The following result establishes the equivalence between problem (6) and QIO$(\mathcal{K},\tau,\theta)$ in terms of the set of achievable inverse-feasible cost vectors given the same dataset $\hat{\mathcal{X}}$.

Theorem 3.4

A cost vector $\mathbf{c}$ is feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ if and only if there exists a feasible solution $(\bar{\mathbf{v}}$ $\neq$ $\mathbf{0},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ for model (6) such that $\mathbf{c}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i:\bar{v}_{i}=1})$ and $\lVert\mathbf{c}\rVert_{p}=1$.

Thus, an inverse-feasible cost vector for QIO$(\mathcal{K},\tau,\theta)$ can be obtained by finding a feasible solution $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ for model (6) and creating a conic combination of $\mathbf{a}^{i}$’s for $i$ such that $\bar{v}_{i}=1$. In fact, without having to normalize $\mathbf{c}$ post-hoc, for any conic combination of such $\mathbf{a}^{i}$’s, there is an inverse-feasible cost vector for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ that generates the same forward optimal solutions. That is, if we let $\mathcal{C}_{QIO}$ be the set of all inverse-feasible $\mathbf{c}$ vectors for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ and $\mathcal{C}_{MIP}=\bigcup_{\mathbf{v}\in\hat{\mathcal{V}}}\textup{cone}(\{\mathbf{a}^{i}\}_{i:{v}_{i}=1})$ where $\hat{\mathcal{V}}$ is the set of feasible non-zero $\mathbf{v}$ vectors for model (6), we have $\bigcup_{\mathbf{c}\in\mathcal{C}_{QIO}}\text{argmin}\,\mathbf{FO}(\mathbf{c})=\bigcup_{\mathbf{c}\in\mathcal{C}_{MIP}}\text{argmin}\,\mathbf{FO}(\mathbf{c})$. The following result formally characterizes the set of all inverse-feasible solutions for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ (proof is straightforward from the proof of Theorem 3.4 and thus is omitted).

Corollary 3.5

Let $\Pi(\hat{\mathcal{X}})$ denote the feasible region of model (6). Then the set of all inverse-feasible cost vectors, $\hat{\mathcal{C}}$, can be characterized by $\hat{\mathcal{V}}=\text{Proj}_{\mathbf{v}}(\Pi(\hat{\mathcal{X}}))\setminus\{\mathbf{0}\}$. That is, $\hat{\mathcal{C}}=\bigcup_{\mathbf{v}\in\hat{\mathcal{V}}}\textup{cone}(\{\mathbf{a}^{i}\}_{i:v_{i}=1}).$

3.1.2 Lower Bound on the Inverse Stability Measure.

Finally, we show that a solution to (6) provides a lower bound on inverse stability measure (3) for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$. To do so, we make the following assumptions: (i) data noise in $\hat{\mathcal{X}}$ is in the form of point-wise shift, i.e., a shifted data point is expressed as $\tilde{\mathbf{x}}^{k}=\hat{\mathbf{x}}^{k}+\boldsymbol{\delta}^{k}$ for each $k\in\mathcal{K}$; (ii) $\displaystyle d(\hat{\mathcal{X}},\tilde{\mathcal{X}})=\sum_{k\in\mathcal{K}}\|\hat{\mathbf{x}}^{k}-\tilde{\mathbf{x}}^{k}\|_{\ell}$ for some $\ell\geq 1$.

Proposition 3.6

Let $\xi^{*}$ be the inverse stability value for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$. Let $(\bar{\mathbf{v}}\neq\mathbf{0},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ be a feasible solution to (6), $\bar{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\bar{v_{i}}=1\}$, $\displaystyle d_{ik}^{*}=\min_{\mathbf{x}\in\mathcal{X}_{i}}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\},i\in\bar{\mathcal{I}},k\in\mathcal{K}$, and $\displaystyle d_{i[k]}^{*}$ denote the $k$-th largest value of $\displaystyle d_{ik}^{*},k\in\mathcal{K}$. Then we have \displaystyle\max_{i\in\bar{\mathcal{I}}}\bigg{\{}\sum_{k=1}^{\lfloor(1-\theta)K\rfloor+1}\max(0,\tau-d_{i[k]}^{*})\bigg{\}}\leq\xi^{*}.

Given $\mathcal{K}$ and $\tau$, we note that the lower bound increases as $\theta$ decreases, which reinforces the role of $\theta$ in increasing inverse stability. Our numerical results show that the lower bound is non-trivial and close to the numerically estimated inverse stability measure.

While problem (6) leads to the set of cost vectors with improved inverse stability, which cost vector among them to use for the forward problem is an important consideration to improve forward stability. In the next subsection, we introduce an objective function for problem (6) to find a $\mathbf{c}$ vector with improved forward stability.

3.2 Finding the Maximal Dimension Inverse-Feasible Set: Improving Forward Stability

Because it is hard to model forward stability measure (4) as a tractable objective function, we propose a surrogate, tractable function that can be used as an objective function for the QIO model. Recall that the forward stability measure for a cost vector $\mathbf{c}$ represents the worst-case distance between $\mathcal{X}^{*}\!(\mathbf{c})=\underset{\mathbf{x}}{\text{argmin}}\,\textbf{FO}(\mathbf{c})$ and the data. We compute forward stability only in terms of the data points that are actually used by the model and consider $\mathbf{c}$ to be forward-stable if the worst-case distance is no greater than the error threshold $\tau$ for the relevant subset of the data of cardinality $\lceil\theta K\rceil$ (i.e., excluding those deemed to be outliers). We first make the following observation.

Proposition 3.7

If a feasible solution $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ for problem (6) satisfies $\displaystyle\sum_{i\in\mathcal{I}}\bar{v}_{i}=n$, then $\bar{\mathbf{c}}\in\textup{cone}_{+}(\{\mathbf{a}^{i}\}_{i:\bar{v}_{i}=1})$ satisfies $\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\bar{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\tau,\forall k:\bar{u}_{k}=1$.

Proposition 3.7 implies that if $\mathbf{c}$ is a strict conic combination of $n$ $\mathbf{a}^{i}$’s chosen by (6) (i.e., such that $\bar{v}_{i}=1$), then $\mathbf{FO}(\mathbf{c})$ has a unique optimal solution that is within $\tau$-distance from all chosen data points (Mangasarian 1979, Tavaslıoğlu et al. 2018), and thus the forward stability measure is guaranteed to be no greater than the pre-specified threshold. Furthermore, the following result shows that as we assign more 1’s to the elements of $\mathbf{v}$ in problem (6) (i.e., adding more basis vectors to the conic combination), the resulting $\mathbf{c}$ vector improves both forward and inverse stability.

Proposition 3.8

Given $\hat{\mathcal{X}}$, let $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ and $(\tilde{\mathbf{v}},\bar{\mathbf{u}},\{\tilde{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ be two distinct feasible solutions to problem (6). Let $\bar{\mathcal{C}}=\textup{cone}(\{\mathbf{a}^{i}\}_{i:\bar{v}_{i}=1})$ and $\tilde{\mathcal{C}}=\textup{cone}(\{\mathbf{a}^{i}\}_{i:\tilde{v}_{i}=1})$. If $\bar{\mathbf{v}}\geq\tilde{\mathbf{v}}$, then

1. (i)

$\bar{\mathcal{C}}$* is at least as forward-stable as $\tilde{\mathcal{C}}$; i.e., for any $\tilde{\mathbf{c}}\in\tilde{\mathcal{C}}$ there exists $\bar{\mathbf{c}}\in\bar{\mathcal{C}}$ such that $\underset{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\underset{\mathbf{x}\in\mathcal{X}^{*}(\tilde{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\},\forall k\in\mathcal{K}$.* 2. (ii)

$\bar{\mathcal{C}}$* is at least as inverse-stable as $\tilde{\mathcal{C}}$; i.e., $\underset{\mathcal{X}\subset\mathbb{R}^{K\times n}}{\min}\{d(\hat{\mathcal{X}},\mathcal{X})\,|\,\bar{\mathcal{C}}\cap\mathcal{C}(\mathcal{X})=\emptyset\}\geq\underset{\mathcal{X}\subset\mathbb{R}^{K\times n}}{\min}\{d(\hat{\mathcal{X}},\mathcal{X})\,|\,\tilde{\mathcal{C}}\cap\mathcal{C}(\mathcal{X})=\emptyset\}$ where $\mathcal{C}(\mathcal{X})$ denotes the set of inverse-feasible cost vectors for $\mathcal{X}$.*

Propositions 3.7 and 3.8 imply that we can add the objective function of maximizing the number of non-zero $v_{i}$’s to problem (6) to find which cost vector to eventually use for the forward problem to further improve forward stability without compromising inverse stability. This leads us to the following mixed-integer quantile inverse optimization (MQIO) model:

[TABLE]

Our numerical results show that even when the optimal value of (7) is not $n$ but close to $n$, the MQIO model is more forward-stable than the previous inverse model. Computing forward stability is inherently challenging as it involves finding the maximum distance between the dataset and some face of polyhedron $\mathcal{X}$. An upper bound for the forward stability measure may be computed if some geometric information about the polyhedron is given, e.g., the minimum volume ellipsoid or other size measure for some faces of $\mathcal{X}$, yet obtaining such information itself is not easy in general. In Appendix 7.4, we show that an optimal solution to the MQIO model leads to an upper bound on the forward stability measure if such information is available.

Appendix 7.5 provides extensions to the MQIO model that can accommodate situations where $(\mathbf{A},\mathbf{b})$ varies over different data points, $(\mathbf{A},\mathbf{b})$ is subject to uncertainty, or the cost vector is constrained. The MQIO problem is a large-size MIP in general and thus is computationally challenging. In the next section, we provide efficient solution approaches by exploiting the problem structure.

4 Solution Approaches

In this section, we develop efficient solution approaches to the MQIO problem. In Section 4.1, we first show that the MQIO problem is NP-hard by establishing its equivalence to a biclique problem. We then provide a formal biclique representation of a feasible solution to the QIO (hence MQIO) problem in Section 4.2, which indicates that algorithms designed for biclique problems can be adapted to solve the MQIO problem. Finally, we propose an exact algorithm and heuristics in Section 4.3.

4.1 Problem Complexity and Connection to Biclique Problems

We show that MQIO$(\mathcal{K},\tau,\theta)$ is NP-hard by establishing its connection to biclique problems, which we later use to develop alternative algorithms for the problem. To do so, we first show that there exists an optimal solution to $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ with exactly $\lceil\theta K\rceil$ observations selected.

Lemma 4.1

If $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ is feasible, then there exists an optimal solution $(\mathbf{v}^{*},\mathbf{u}^{*},\{\boldsymbol{\epsilon}^{k*}\}_{k\in\mathcal{K}})$ such that $\sum_{k=1}^{K}u^{*}_{k}=\lceil\theta K\rceil$.

Proposition 4.2

The problem $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ is NP-hard.

To prove Proposition 4.2 we first introduce the maximum $\kappa$-subset intersection (MSI) problem. Consider a ground set $\mathcal{G}=\{g_{1},\ldots,g_{N}\}$ and a set of its subsets $\mathscr{R}=\{\mathcal{R}_{1},\ldots,\mathcal{R}_{Q}\}$, i.e., $\mathcal{R}_{q}\subseteq\mathcal{G},q=1,\ldots,Q$. Given a positive integer $\kappa$, MSI finds $\kappa$ subsets in $\mathscr{R}$ whose intersection has maximum cardinality. MSI problems are known to be NP-hard (Dawande et al. 2001, Xavier 2012).

Proof 4.3

Proof of Proposition 4.2: Consider the following feasibility problem for each $k\in\mathcal{K}$.

[TABLE]

Let $(\mathbf{v}^{k},\boldsymbol{\epsilon}^{k})$ be a feasible solution to (8) with respect to $\hat{\mathbf{x}}^{k}$ and $\mathcal{A}_{k}=\{i\in\mathcal{I}\,|\,v_{i}^{k}=1\}$ be a subset of $\mathcal{I}$ for each $k\in\mathcal{K}$, and construct the set of subsets $\mathscr{A}=\{\mathcal{A}_{1},\ldots,\mathcal{A}_{K}\}$. Finding exactly $\lceil\theta K\rceil$ subsets in $\mathscr{A}$ such that their intersection has maximum cardinality is equivalent to solving MSI with $\kappa=\lceil\theta K\rceil$. Note that for each $k\in\mathcal{K}$ there can be multiple solutions for $\mathbf{v}^{k}$ satisfying (8), each leading to a different $\mathcal{A}_{k}$. As a result, $\mathscr{A}$ may not be unique. Thus, $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ is equivalent to solving the MSI problem multiple times with different $\mathscr{A}$’s and finding the maximum cardinality. Therefore, $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ is at least as hard as MSI. \Halmos

Dawande et al. (2001)** show that a general MSI problem can be reformulated as a version of the biclique problem (called the maximum one-sided edge cardinality problem).

4.2 Biclique Representation of the Solution

Motivated by its connection to the MSI problem, we also cast the MQIO problem as a biclique problem with a bipartite graph constructed as follows. We first create a node $k$ in $V_{1}$ for each data point $\mathbf{\hat{x}}^{k}$ and a node $i$ in $V_{2}$ for each $\mathbf{a}^{i}$. With every node $k\in V_{1}$ and $i\in V_{2}$ we associate an edge $e=(k,i)\in\mathbb{E}$ if $v_{i}^{k}=1$ where $\mathbf{v}^{k}$ denotes a feasible solution to problem (8) with respect to $\mathbf{\hat{x}}^{k}$. Then we define a bi-adjacency matrix $\bar{\mathbf{D}}\in\{0,1\}^{K\times m}$ associated with graph $G=(V_{1}\cup V_{2},\mathbb{E})$ as follows.

[TABLE]

Matrix $\bar{\mathbf{D}}$ can be built by solving problem (8) for each $k\in\mathcal{K}$ to find a feasible solution $\mathbf{v}^{k}$ and letting the $k$-th row of $\bar{\mathbf{D}}$ be ${\mathbf{v}^{k}}^{\prime}$. The following result shows how an inverse-feasible solution can be found from the $\bar{\mathbf{D}}$ matrix above.

Proposition 4.4

If there exists $\bar{\mathbf{D}}$ that satisfies (9) and has an all-one submatrix whose rows and columns correspond to $\bar{\mathcal{S}}$ and $\bar{\mathcal{A}}$, respectively, where $\bar{\mathcal{S}}\subseteq\mathcal{K}$, $|\bar{\mathcal{S}}|\geq\theta K$, and $\bar{\mathcal{A}}\subseteq\mathcal{I}$, then there exists a solution $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ where $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$.

*Proposition 4.4 suggests that if we can find a matrix $\bar{\mathbf{D}}$ that has an all-one submatrix with at least $\lceil\theta K\rceil$ rows, there exists a corresponding inverse-feasible cost vector for the QIO (hence MQIO) problem. The following example further illustrates this. **

---

Example 4.5

Consider Example 2.3 again (see Figure 2(a)), and the QIO model with the following parameters: $\tau=1$, $\theta=0.8$, $\|\cdot\|_{\ell}=\|\cdot\|_{\infty}$, and $\|\cdot\|_{p}=\|\cdot\|_{\infty}$. Let $(\bar{\mathbf{v}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ be a feasible solution to problem (8) for each data point $k$: we have $\bar{\mathbf{v}}^{k}=[1,1,0,0]^{\prime}$ for $k=1,2,3,4,$ with $\bar{\boldsymbol{\epsilon}}^{1}=\begin{bmatrix}-0.5\\ -0.2\end{bmatrix}$, $\bar{\boldsymbol{\epsilon}}^{2}=\begin{bmatrix}-0.3\\ -0.2\end{bmatrix}$, $\bar{\boldsymbol{\epsilon}}^{3}=\begin{bmatrix}-0.3\\ -0.5\end{bmatrix}$, and $\bar{\boldsymbol{\epsilon}}^{4}=\begin{bmatrix}-0.5\\ -0.5\end{bmatrix}$, respectively, and $\bar{\mathbf{v}}^{5}=[0,1,1,0]^{\prime}$ with $\bar{\boldsymbol{\epsilon}}^{5}=\begin{bmatrix}-0.3\\ 0.3\end{bmatrix}$. Then we can construct a bipartite graph, shown in Figure 2(b), where each edge identifies $(i,k)$ such that $\bar{v}_{i}^{k}=1$, and the corresponding $\bar{\mathbf{D}}$ matix in Figure 2(c). Note that this $\bar{\mathbf{D}}$ has an all-one submatrix with rows and columns corresponding to $\bar{\mathcal{S}}=\{1,2,3,4\}$ and $\bar{\mathcal{A}}=\{1,2\}$, respectively (the shaded rectangle), and $|\bar{\mathcal{S}}|=4\geq\theta K$. Now consider a cost vector $\bar{\mathbf{c}}=0.5\mathbf{a}^{1}+0.5\mathbf{a}^{1}=\begin{bmatrix}-0.5\\ -0.5\end{bmatrix}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$. We can see that $\bar{\mathbf{c}}$ is inverse-feasible because $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{1},\ldots,\bar{\boldsymbol{\epsilon}}^{5}\},\bar{\mathbf{y}},\bar{\mathcal{S}})$ where $\bar{\mathbf{y}}=[0.5,0.5,0,0]^{\prime}$ is feasible for QIO($\mathcal{K},\tau,\theta$).

*Furthermore, Proposition 7 implies that the number of columns in the all-one submatrix corresponds to the objective function value of $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$, which is to be maximized. Thus, it is important to find a submatrix with as many columns as possible. Building on these observations, in the next section, we propose specific algorithms to find an all-one submatrix of $\bar{\mathbf{D}}$ with at least $\lceil\theta K\rceil$ rows and as many columns as possible.

---

4.3 Algorithms

While the matrix $\bar{\mathbf{D}}$ is constructed using $\{\mathbf{v}^{k}\}_{k\in\mathcal{K}}$ obtained by solving problem (8) for each data point $k$, because there may be multiple solutions for (8), there may be multiple $\bar{\mathbf{D}}$ matrices achievable. Thus, it is important to find a $\bar{\mathbf{D}}$ matrix that can lead to an all-one submatrix of a desirable size as described above, which eventually leads to a desirable feasible solution to $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$. In this subsection, we now present our main algorithm that finds such a $\bar{\mathbf{D}}$ matrix and corresponding inverse-feasible cost vectors.**

*Our first idea is that we solve problem (8) with a weighted objective function for each $k$ such that resulting solutions $\{\mathbf{v}^{k}\}_{k\in\mathcal{K}}$ are similar to each other, thus leading to similar rows of $\bar{\mathbf{D}}$. For the first data point, i.e., $k=1$, we solve problem (8) with respect to $\hat{\mathbf{x}}^{1}$ with the objective function $\displaystyle\sum_{i\in\mathcal{I}}v_{i}$, which will find a solution $\mathbf{v}^{1}$ with as many ones as possible. For $k\geq 2$, we solve problem (8) with respect to $\hat{\mathbf{x}}^{k}$ with the weighted objective function $\displaystyle\sum_{i\in\mathcal{I}}w_{i}^{k}v_{i}$ where $\mathbf{w}^{k}=\boldsymbol{\vartheta}^{k}/\|\boldsymbol{\vartheta}^{k}\|$ and $\displaystyle\vartheta^{k}_{i}=\sum_{\ell=1}^{(k-1)}v^{\ell}_{i}$; i.e., we assign a higher weight to column $i$ for which $v_{i}^{\ell}$ was assigned 1 more often for the previous data points $\ell=1,\ldots,k-1$. The resulting vector ${\mathbf{v}^{k}}^{\prime}$ is then inserted to the $k$-th row of $\bar{\mathbf{D}}$. Note that because the weights are updated in a sequential manner, the order of the data points can impact the resulting $\bar{\mathbf{D}}$. We assume without loss of generality that the data points are sorted in the increasing order of the sum of the distances to the rest of the dataset (i.e., $\displaystyle\sum_{k\in\mathcal{K}}\|\hat{\mathbf{x}}^{1}-\hat{\mathbf{x}}^{k}\|_{\ell}\leq\ldots\leq\sum_{k\in\mathcal{K}}\|\hat{\mathbf{x}}^{K}-\hat{\mathbf{x}}^{k}\|_{\ell}$); we find this ordering particularly effective because starting with data points that are closer to the rest of the dataset leads to more similar $\mathbf{v}^{k}$’s in the subsequent iterations.

---

*Next, since the construction of $\bar{\mathbf{D}}$ depends on the solution $\mathbf{v}^{1}$ for problem (8) with respect to the first data point $\hat{\mathbf{x}}^{1}$ and this problem may have multiple solutions, it is important to consider different solutions for this problem to place different “seeds” for $\bar{\mathbf{D}}$ (i.e., different first rows of $\bar{\mathbf{D}}$), which allows to explore different inverse solutions. To do so, once an optimal solution to problem (8) (with the added objective function) for the first data point $\hat{\mathbf{x}}^{1}$ is found, we impose a cut to exclude the solution from the feasible region and re-solve the problem repeatedly.

---

*Finally, once a certain $\bar{\mathbf{D}}$ matrix is found, its all-one submatrix (i.e., a biclique in the bipartite graph corresponding to $\bar{\mathbf{D}}$) with at least $\lceil\theta K\rceil$ rows and as many columns as possible can be found (if exists) by solving the associated maximum clique problem. While any existing algorithm for a general maximum clique problem would work, we use an MIP formulation to solve the problem exactly, which we call Clique($\bar{\mathbf{D}},\theta$); this formulation can be found in Appendix 8.1.

---

*This algorithm is formally presented in Algorithm 2. In Lines 4-9, problem (8) is solved with the weighted objective function and cuts. To generate different first rows of $\bar{\mathbf{D}}$ (i.e., different $\mathbf{v}^{1}$’s), a cut defined in Line 11 is added to the problem in Line 5. The following result shows that, under a certain condition, this algorithm finds inverse-feasible cost vectors that correspond to the optimal solution to MQIO($\mathcal{K},\tau,\theta$), hence exact. We call this algorithm $\bar{\mathbf{D}}$-Alg-Exact. **

---

Proposition 4.6

Let $(\mathbf{v}^{*},\mathbf{u}^{*},\{\boldsymbol{\epsilon}^{k*}\}_{k\in\mathcal{K}})$ be an optimal solution to $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ and $z^{*}=\sum_{i\in\mathcal{I}}v^{*}_{i}$. Let $\mathcal{V}_{z}=\{\mathbf{v}^{1}\in\{0,1\}^{m}\,|\,\exists(\mathbf{v}^{1},\boldsymbol{\epsilon}^{1})$ satisfying (8b)–(8d) with respect to $\hat{\mathbf{x}}^{1}\!,\;{\mathbf{v}^{1}}^{\prime}\mathbf{e}=z\}$. If $u^{*}_{1}=1$, then Algorithm 2 returns $z^{Alg}=z^{*}$ and $\mathbf{v}^{Alg}=\mathbf{v}^{*}$ in $r\leq\displaystyle\sum_{z=z^{*}}^{n}|\mathcal{V}_{z}|$ iterations.

*Once the algorithm is run, the condition $u_{1}^{*}=1$ can be easily checked by whether the first row of $\bar{\mathbf{D}}_{r}$ is also included in its all-one submatrix where $r$ is the iteration at which the most recent $\mathbf{v}^{Alg}$ was achieved. Our numerical results show that the condition is met in all instances for which the algorithm was run completely (not interrupted by time limit). If this condition is not met, although the algorithm may still generate a reasonably good solution, the user may reorder the data points to place a different data point in the position of $k=1$ and rerun the algorithm. ***

\OneAndAHalfSpacedXI

Algorithm 2 without cuts, i.e., finding only one $\bar{\mathbf{D}}$, would be faster, albeit not exact, which we call $\bar{\mathbf{D}}$-Alg-Heuristic. We also propose an even faster heuristic called $\tilde{\mathbf{D}}$-Alg-Heuristic in Appendix 8.2, which solves a linear relaxation of (8) and uses weighted $\ell_{1}$ minimization to make $\mathbf{1}\!-\!\mathbf{v}$ sparse where $\mathbf{v}\in[0,1]^{m}$. The proposed algorithm framework can lead to many different variants depending on how the weights are updated and how cuts are generated. For example, the algorithms can be modified to accommodate a huge dataset by applying the sequential weight updating only to the first few MIPs (say $\bar{K}$ of them), which leads to the weight parameters $\mathbf{w}^{\bar{K}}$, and solving the remaining ($K\!-\!\bar{K}$) MIPs in parallel with the same weights $\mathbf{w}^{\bar{K}}$. Also, the algorithms can be modified to infer cost vectors in an online manner based on datasets that become available through different time points and $(\mathbf{A},\mathbf{b})$ parameters changing over time; Appendix 8.3 shows more details about this algorithm.**

5 Numerical Results

We examine the performance of the algorithms proposed for the MQIO model using various-sized LP instances. We also assess forward and inverse stability of the MQIO model as well as the previous inverse LP method using the measures proposed in Section 2. Finally, we demonstrate the stability performance of the MQIO model in the diet recommendation and transshipment applications.

5.1 Performance of the Algorithms

We evaluate the performance of the proposed algorithms for randomly generated LP instances with $n\in\{15,50\}$, $m\in\{100,300\}$, $K\in\{35,200,500\}$, and $\tau\in\{3,3.5,4\}$. An instance is defined by a tuple $(n,m,K,\tau)$ and each instance was solved twice—with $\theta=0.75$ and $0.85$. For each instance, as a pre-processing step, we excluded $\mathbf{a}^{i}$’s that were not inverse-feasible a priori from the construction of $\mathbf{c}$ (see Algorithm 1).

Table 1 shows the results for the MQIO problem (7) obtained by the Gurobi solver (Gurobi Optimization 2020) (labeled “MIP”), $\bar{\mathbf{D}}$-Alg-Exact, $\bar{\mathbf{D}}$-Alg-Heuristic, and $\tilde{\mathbf{D}}$-Alg-Heuristic,* averaged over two results with $\theta=0.75$ and $0.85$. Recall that the objective function value represents the number of basis vectors ($\mathbf{a}^{i}$’s) to construct the set of inverse-feasible cost vectors. All optimization problems were solved by Gurobi 7.5 with a 4-core 3.6 GHz processor and 32 GB memory.*

*For small instances where $(n,m)=(15,100)$, the exact MIP formulation found an optimal solution within the time limit of 5,000 seconds for most cases. Note that the solution time for MIP decreases as $\tau$ goes up—we conjecture that this is because the problem becomes easier as the threshold distance $\tau$ becomes more “generous.” ** The exact algorithm achieved the same objective values as the MIP values within significantly less time. Both heuristics found objective values reasonably close to the exact values and were even faster. The computational benefit of the proposed algorithms becomes more clear for larger instances; the MIP solver failed to solve any instance with $(n,m)=(50,300)$. As expected, there is a clear trade-off between the algorithms. In all instances, $\bar{\mathbf{D}}$-Alg-Exact returned the greatest objective values but was the slowest, whereas $\tilde{\mathbf{D}}$-Alg-Heuristic was the fastest algorithm yet the solution quality was not as good as the other two algorithms. Although the results suggest that both heuristics can be used for very large instances, since $\tilde{\mathbf{D}}$-Alg-Heuristic solves LP relaxations for each $\hat{\mathbf{x}}^{k}$ instead of MIPs, we suggest using this algorithm for instances with very large $n$ and $m$.

---

5.2 Stability of MQIO

5.2.1 Forward Stability.

To assess the forward stability of our proposed MQIO model, we created an instance with $(n,m,K,\tau)=(15,100,35,3)$ and $\theta=0.75$, for which we knew the true optimal objective value was $\sum_{i\in\mathcal{I}}v^{*}_{i}=15$. To see the effect of the number of $\mathbf{a}^{i}$’s used to form an inverse solution (i.e., $\sum_{i\in\mathcal{I}}v_{i}$) on forward stability, we solved the inverse problem on this instance repeatedly, each with a constraint $\sum_{i\in\mathcal{I}}v_{i}\leq h$, $h\in\{1,2,\ldots,15\}$. Each problem was solved exactly via the MQIO formulation. Note that the case where $h=1$ can be considered similar to the previous inverse LP model (2) in that it finds a single cost vector identical to one of the $\mathbf{a}^{i}$’s. For each $h$, once the set of inverse-feasible cost vectors, $\mathcal{C}_{h}$, was found, we randomly generated 50 cost vectors $\mathbf{c}\in\mathcal{C}_{h}$ and computed $d(\hat{\mathcal{X}},\mathbf{x}^{*}(\mathbf{c}))=\underset{k\in\mathcal{S}_{h}}{\max}\,\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}^{*}(\mathbf{c})\|_{\infty}\}$, where $\mathcal{S}_{h}$ denotes the set of chosen data points in iteration $h$ and $\mathbf{x}^{*}(\mathbf{c})\in\argmin\textbf{FO}(\mathbf{c})$, as an estimate for the forward stability measure.

Figure 3(a) shows the resulting distances for each $h\in\{1,2,\ldots,15\}$: each box-plot represents 50 values of $d(\hat{\mathcal{X}},\mathbf{x}^{*}(\mathbf{c}))$ from 50 different cost vectors. We observe that when $h$ is low (i.e., few $\mathbf{a}^{i}$’s are used for generating a $\mathbf{c}$ vector) the distance $d$ is significantly greater than $\tau$, hence forward-unstable. In other words, if this cost vector is used for the forward problem, the resulting solutions may not be close to the input data. On the other hand, as $h$ increases (i.e., as more $\mathbf{a}^{i}$’s are used) the cost vector produces a solution that is closer to the data, eventually within the threshold distance $\tau$ (indicated by the horizontal dashed line). The overall decreasing trend of the distance $d$ supports the idea of maximizing $\sum_{i\in\mathcal{I}}v_{i}$ as a surrogate objective function to maximize forward stability.

5.2.2 Inverse Stability.

To evaluate the inverse stability of MQIO, we investigated the effect of $\theta$ and $\tau$ on how sensitive the inverse solution set is to changes in the dataset $\hat{\mathcal{X}}$. Again, an instance with $(n,m,K)=(15,100,35)$ was considered. First, we considered the MQIO model with different $\tau$ values increasing from 2 to 3 by 0.1 (the objective value of MQIO increases monotonically as $\tau$ increases). For each $\tau$ value, we first solved the MQIO model with the given dataset $\hat{\mathcal{X}}$ to find the initial inverse-feasible set $\hat{\mathcal{C}}$. ** We then shifted $\lfloor(1-\theta)K\rfloor+1$ data points in $\hat{\mathcal{X}}$: $\{\hat{\mathbf{x}}^{k}+\boldsymbol{\gamma}^{k}\}_{k\in\mathcal{K}}$ where $\boldsymbol{\gamma}^{k}$ is uniformly distributed in $[0,\Gamma]^{n}$ if $\hat{\mathbf{x}}^{k}$ is chosen to be shifted, and $\boldsymbol{\gamma}^{k}=\mathbf{0}$ otherwise. Given this shifted data, denoted by $\tilde{\mathcal{X}}(\Gamma)$, we re-solved the MQIO problem to find the new inverse-feasible set $\tilde{\mathcal{C}}$. ** We repeated this process while increasing $\Gamma$ until $\tilde{\mathcal{C}}$ had no common cost vector shared by the initial set $\hat{\mathcal{C}}$, i.e., $\hat{\mathcal{C}}\cap\tilde{\mathcal{C}}=\emptyset$. When such a $\Gamma$ value was reached, we computed the distance $\displaystyle d(\hat{\mathcal{X}},\tilde{\mathcal{X}}(\Gamma))=\sum_{k\in\mathcal{K}}\|\hat{\mathbf{x}}^{k}-\tilde{\mathbf{x}}^{k}(\Gamma)\|_{\infty}$ as an estimate for the inverse stability measure in Section 2.

Figure 3(b) shows the result of this experiment for all $\tau$ values and $\theta=0.75,0.85$ and $1$. For a fixed value of $\theta$, the distance $d(\hat{\mathcal{X}},\tilde{\mathcal{X}}(\Gamma))$ increases as $\tau$ increases (i.e., more $\mathbf{a}^{i}$ vectors are included in the inverse-feasible set), which means the model becomes more tolerant to changes in the dataset, hence more inverse-stable. For a fixed value of $\tau$, Figure 3(b) also shows that increasing $\theta$ decreases inverse stability, reinforcing the idea that allowing for no or few outliers can lead to an “impatient” inverse model that can easily switch its solutions affected by such data imperfection.

*Note that the estimated inverse stability values in Figure 3(b) are in fact upper bounds on the true inverse stability values. To show this, we refer the reader to formulation (12) in Appendix 7.3, whose optimal value is the true inverse stability value for the MQIO model. Since no initial cost vector in $\hat{\mathcal{C}}$ remains inverse-feasible given the shifted data $\tilde{\mathcal{X}}$, $\tilde{\mathcal{X}}$ is a feasible solution to formulation (12), with the objective value of $d(\hat{\mathcal{X}},\tilde{\mathcal{X}})$. Because (12) is a minimization problem, its optimal value (i.e., the true inverse stability value) cannot be greater than $d(\hat{\mathcal{X}},\tilde{\mathcal{X}})$. Figure 3(b) shows that the lower bounds on the inverse stability measure (based on Proposition 3.6) are non-trivial, and since the estimated values are upper bounds, the difference between the lower bounds and the true inverse stability values must be even less than shown in Figure 3(b). ***

Finally, we conducted a similar experiment to assess the inverse stability of the previous inverse model (2). We solved model (2) with $\hat{\mathcal{X}}$ to find an initial inverse solution $\hat{\mathbf{c}}$, re-solved the problem repeatedly with $\tilde{\mathcal{X}}(\Gamma)$ with increasing $\Gamma$ values until the resulting cost vector was different from $\hat{\mathbf{c}}$, and computed the distance between the original dataset and the shifted one, which is shown as a thick $\boldsymbol{\times}$ marker in Figure 3(b). The result is almost identical to the lowest inverse stability achieved by the MQIO model (i.e., with $\theta=1$ and smallest $\tau$).

5.3 Diet Recommendation

Literature on diet prediction/recommendation using historical data largely focuses on a “direct” replication of diets where the goal is to create diets that are closest to the old diets in terms of some distance function (e.g., $\ell_{2}$ distance between new and old diet vectors) (Perignon et al. 2016, Darmon et al. 2006). However, such an approach may not work if there is any change in the constraints of the underlying diet problem such as changes in nutritional requirements or available foods. In this case, learning the objective function that represents one’s preferences instead and creating new diets “indirectly” is a more robust and transferable way of learning the individual’s diet behavior (similar motivations can be found in learning driving behaviors in Abbeel and Ng (2004)**). Moreover, by doing so the model can find an optimal diet with respect to the inferred preferences, whereas replicating the diets directly can inherit undesirable diet patterns from the past.

We use the database from the National Health and Nutrition Examination Surveys (NHANES), which includes nutritional requirements and nutrition facts per serving for each food type, and build a diet problem for a subset of the foods and nutrients from the database to keep the experiment simple. We classify the foods into food “types” (see Table 10.1 in Appendix 10.1). We assume that once the number of servings for each food type is determined, more detailed decisions (e.g., specific menus) can be made by dietitians based on this. We believe this is a realistic consideration as there exist a myriad of different menus. We first solve the diet problem with some arbitrary cost vector, assumed to reflect true preferences of an individual, and find $\mathbf{x}^{*}$ that represents the number of servings for food type $i$ per day. In practice, if the data come from some unknown system, the assumed forward LP feasible region may not “fit” the data well. In this case, the user needs to tweak some parameters of the forward LP (e.g., constraint parameters or the list of variables), similar to linear regression.

Given $\mathbf{x}^{*}$ with no noise, both the MQIO and previous models can find the true cost vector and generate the same $\mathbf{x}^{*}$ as a forward solution. However, as one’s diet behavior is assumed inconsistent and noisy, to examine the impact of data noise on both models we generate multiple diets $\mathbf{x}^{*}+\boldsymbol{\gamma}$* where $\boldsymbol{\gamma}$ is a normal-distributed random noise vector $\boldsymbol{\gamma}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$ where $\mathbf{I}$ is the identity matrix of an appropriate dimension, which form the input dataset $\hat{\mathcal{X}}$.*

We apply the MQIO model as well as the previous inverse model (2) with the same dataset $\hat{\mathcal{X}}$ as input and compare the resulting recommended diets. For MQIO, we find a set of inverse-feasible cost vectors $\mathcal{C}^{*}$ and solve the forward diet problem 20 times each with a different $\mathbf{c}$ vector selected from $\mathcal{C}^{*}$ to generate a set of diet recommendations. Since the previous inverse model finds a single cost vector, we solve the diet problem with this vector only once. Figure 4(a) shows diet recommendations obtained by the MQIO model (circular markers) are much closer to the input data (diamond-shaped markers; $\sigma=0.2$ to create the noisy dataset $\hat{\mathcal{X}}$) than those from the previous inverse model are (squared markers), reinforcing the improved stability in the MQIO model. Also note that the MQIO model leads to multiple such diet recommendations based on the multiple cost vectors.

*In Figure 4(b), we further increase the noise in the input data by increasing $\sigma$ from 0.2 to 1. While the previous model generates the same recommendation, the MQIO model “adapts to” the increased variability in the observations and generates more diverse recommendations for some food types; e.g., food types 5–7. In summary, the MQIO model generates recommendations that are consistent with the individual’s past behavior and stable in the face of data noise; thus, the objective function found by the model can better predict one’s eating behavior. *

*Finally, we show that the MQIO framework is amenable to the case where data point are collected at different time points and the goal of the inverse problem is to infer cost vectors based on more recent data points. Let $\hat{\mathbf{x}}^{t}$ be a data point collected at time $t=1,\ldots,T$. We let the threshold parameter $\tau_{t}$ vary over time: the more recent $\hat{\mathbf{x}}^{t}$ is, the smaller $\tau_{t}$ we impose (i.e., the closer we want to make this data point to optimality). Figure 5 shows the result of such a modified model with time-stamped diet observations. Darker diamond markers represent more recent observations. As more recent data were collected, the modified MQIO model generates recommendations (circular markers) that are closer to the more recent ones. This result suggests that our MQIO framework can offer modeling flexibility that can lend itself to adaptive settings where the preference function can be updated over time. *

5.4 Transshipment Problem

We demonstrate the use of our inverse model for online learning in the transshipment application, where production and shipment costs are inferred adaptively as batches of decisions are observed sequentially over time. We consider a transshipment problem with one distribution node, two supply nodes, and two demand nodes (see Appendix 10.2). Suppose at each time $t$, demand $\mathbf{d}^{t}$ is revealed and the DM makes a set of decisions on production ($\mathbf{x}^{(p)}$) and transshipment ($\mathbf{x}^{(t)}$) to minimize the total cost ${\mathbf{c}^{(p)}}^{\prime}\mathbf{x}^{(p)}+{\mathbf{c}^{(t)}}^{\prime}\mathbf{x}^{(t)}$ where $\mathbf{c}^{(p)}$ and $\mathbf{c}^{(t)}$ denote production and transshipment costs, respectively. Let $\mathbf{c}=[{\mathbf{c}^{(p)}};{\mathbf{c}^{(t)}}]$ and $\mathbf{x}=[{\mathbf{x}^{(p)}};{\mathbf{x}^{(t)}}]$. Let $\hat{\mathcal{X}}_{t}$ denote the set of decisions observed at time $t$. We assume that there is a true cost vector $\mathbf{c}^{true}$ (shown in Appendix 10.2), and generate a dataset $\hat{\mathcal{X}}_{t}$ by adding noise to an optimal solution with respect to $\mathbf{c}^{true}$ and randomly generated demand $\mathbf{d}^{t}\sim U[0,1.1]$.

We use an online learning extension to $\bar{\mathbf{D}}$-Alg-Heuristic (see Appendix 8.3) to infer the set of cost vectors $\mathcal{C}_{t}$ at each time $t$ from the datasets collected through time $t$. We consider $T=600$ (i.e., 600 sequential batches of data) and $\theta=0.75$. At each iteration $t=1,\ldots,600$, we randomly select 20 cost vectors $\mathbf{c}^{t}\in\mathcal{C}_{t}$ and solve the corresponding forward problem with demand $\mathbf{d}^{t}$ to obtain a forward solution $\bar{\mathbf{x}}^{t}$. We assess the performance of the algorithm at each time $t$ in terms of the distance between our forward solution $\bar{\mathbf{x}}^{t}$ and data batch $\hat{\mathcal{X}}_{t}$, i.e., $d(\hat{\mathcal{X}}_{t},\bar{\mathbf{x}}^{t})=\underset{\mathbf{x}\in\hat{\mathcal{X}}_{t}}{\max}\,\|\mathbf{x}-\bar{\mathbf{x}}^{t}\|_{\infty}$.

*Figure 6(a) shows the average distance $\sum_{t=1}^{\bar{t}}d(\hat{\mathcal{X}}_{t},\bar{\mathbf{x}}^{t})/\bar{t}$ at iteration $\bar{t}$ over 20 trials, i.e., 20 randomly selected cost vectors. Each light-color line shows each distance per trial and the dark line shows the average. As $\bar{t}$ grows, the average distance between the dataset and the forward solution converges. In fact, the inferred set of cost vectors remains unchanged from iteration 415 onwards, and importantly, this set includes $\mathbf{c}^{true}$. That is, as more data are collected, the algorithm successfully finds a set of cost vectors that represent the preferences encoded in the true cost vector and generate forward solutions close to the given data (hence forward-stable). ** Additional results on the comparison between the objective function values achieved by our forward solution $\bar{\mathbf{x}}^{t}$ and that from the given data are available in Appendix 10.2.

---

We also test the algorithm with different values of $\theta\in\{0.75,0.85,1\}$. Figure 6(b) shows that as $\theta$ decreases, the algorithm takes more iterations to converge, which is intuitive as less data points are considered relevant at each iteration, requiring the algorithm to collect more data until it finally arrives at a certain set of cost vectors. When $\theta=1$, although the distance converges more quickly, the resulting inverse set does not include the true cost vector because the algorithm with $\theta=1$ is more susceptible to data noise. This again supports the need to explicitly account for potential outliers in inverse LP. Finally, Figure 6(b) also shows the performance of the previous inverse model (2) in this setting: at each iteration $\bar{t}$ we find a cost vector by solving (2) with all data points up to time $\bar{t}$. This model fails to find the true cost vector at all iterations. Also, the distance between its corresponding forward solution and the data is greater than that from the MQIO model.

The performance of the online learning algorithm (i.e., Algorithm 4) also depends on the variability of demand parameters. When the interval of the demand distribution increases, the true vector is found more slowly. We conjecture that this is because increased variability in demands causes the feasible region of the forward problem to vary more, rendering some constraints redundant and thus the algorithm cannot detect some of the $\mathbf{a}^{i}$’s needed to construct the true cost vector.

6 Conclusion

*In this paper we developed a new inverse LP method that can capture noise, errors, and uncertainty in the input data and infer cost vectors that are more stable than those obtained by the previous methods. We formulated the model as a large-scale MIP and developed efficient algorithms by exploiting its connection to the well-known biclique problems. Our inverse method was demonstrated in the diet recommendation and transshipment applications where past data used as input can be noisy and inaccurate. Many directions for future research exist. For example, preferences of a large group of DMs can be clustered into smaller representative groups using the quantile statistic-based inverse optimization approach. ** A goodness-of-fit measure that is amenable to the quantile inverse model can be useful for assessing the validity of the assumed forward LP, which remains our future work. Extending the inverse nonlinear programming techniques to accommodate the quantile-based framework is also an interesting future work, which can allow to assess our proposed inverse model in a broader context.

---

{APPENDICES}

7 Supplemental Materials for the Inverse Model

7.1 Finding the Threshold Parameter $\tau$

The threshold parameter $\tau$ for the QIO model (i.e., (5)) can be set by adding a reasonable margin to the minimum possible value that keeps the model feasible, which can be found by solving the following problem where $\tau$ is now a variable:

[TABLE]

The above problem (10) finds a cost vector such that the $\theta$ quantile optimality error is minimized; this is similar to the structure of the least quantile method in the regression context in Bertsimas and Mazumder (2014)**. The following result shows the solution structure of an optimal cost vector for the above problem (proof can be found in Appendix 9).

Proposition 7.1

*If problem (10) is feasible, then there exists an optimal ${\mathbf{c}}$ for (10) such that ${\mathbf{c}}=\mathbf{a}^{i}$ for some $i\in\mathcal{I}$. *

Proposition 7.1 suggests that the search for an optimal cost vector for (10) can be done by evaluating each hyperplane $i=1,\ldots,m$, i.e., solving formulation (23) for each $i$ and each $k$ (see proof of Proposition 7.1 in Appendix 9) and finding $i$ that induces the minimum $\tau^{*}$; formulation (23) is a convex optimization problem and thus is straightforward to solve.

*To further illustrate the implication of Proposition 7.1, consider the forward feasible region and the initial data points in Example 2.1 (Figure 1(a)) and problem (10) with $\theta=0.5$ and $\|\cdot\|_{\ell}=\|\cdot\|_{\infty}$. Cost vector $\mathbf{c}=\mathbf{a}^{1}$ is the optimal solution to problem (10), leading to the objective value of $\tau^{*}=0.2$ with the optimality errors $\|\boldsymbol{\epsilon}^{1}\|_{\infty}=\|\boldsymbol{\epsilon}^{2}\|_{\infty}=0.2$; the objective function values for (10) achieved by $\mathbf{c}=\mathbf{a}^{2}$, $\mathbf{c}=\mathbf{a}^{3}$, and $\mathbf{c}=\mathbf{a}^{4}$ are $0.3$, $2$, and $2$, respectively. ***

7.2 Finding Big $M$ Parameters

For each $k\in\mathcal{K}$, consider the following formulation that finds $\mathbf{x}\in\mathcal{X}$ that has maximum $\ell_{\infty}$ distance from $\hat{\mathbf{x}}^{k}$:

[TABLE]

Let $d_{k}$ denote the optimal value of (11) with respect to $\hat{\mathbf{x}}^{k}$. Given that $\mathcal{X}$ is bounded, $d_{k}$ is also bounded. To find $d_{k}$, we first reformulate (11) as $\displaystyle\displaystyle\max_{j\in\mathcal{J}}\;\underset{\mathbf{x}\in\mathcal{X}}{\max}\;\{|x_{j}-\hat{x}_{j}^{k}|\}$; for each $j\in\mathcal{J}$, if $\bar{\mathbf{x}}$ denotes an optimal solution to the inner max problem, then {\max}\{|\bar{x}_{j}-\hat{x}_{j}^{k}|\}=\max\Big{\{}\bar{x}_{j}-\hat{x}_{j}^{k},\hat{x}_{j}^{k}-\bar{x}_{j}\Big{\}}. As a result, we have \displaystyle d_{k}=\max_{j\in\mathcal{J}}\;\max\Big{\{}\max_{\mathbf{x}\in\mathcal{X}}\;\{x_{j}-\hat{x}_{j}^{k}\},\max_{\mathbf{x}\in\mathcal{X}}\;\{\hat{x}_{j}^{k}-x_{j}\}\Big{\}}, i.e., $d_{k}$ can be obtained by solving $2n$ LPs and selecting the largest objective value. Note that $M_{2}=\underset{k\in\mathcal{K}}{\max}\;\{d_{k}\}$ is a sufficiently large number valid for constraint (6c). More efficiently, one may use a constraint-specific parameter $M_{2k}=d_{k}$ for each $k\in\mathcal{K}$.

To find $M_{1}$, we use Cauchy-Schwarz and Minkowski inequalities as follows. For each $i\in\mathcal{I}$ and $k\in\mathcal{K}$, by Cauchy-Schwarz inequality, ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k})\leq\lVert\mathbf{a}^{i}\rVert_{2}(\lVert\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k}\rVert_{2})$. Using Minkowski inequality, we have $\lVert\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k}\rVert_{2}\leq\lVert\mathbf{\hat{x}}^{k}\rVert_{2}+\lVert\boldsymbol{\epsilon}^{k}\rVert_{2}$ for each $k\in\mathcal{K}$. Moreover, for each vector $\boldsymbol{\epsilon}^{k}\in\mathbb{R}^{n}$ we have $\lVert\boldsymbol{\epsilon}^{k}\rVert_{2}=\sqrt{\sum_{i=1}^{n}\lvert\epsilon_{i}^{k}\rvert^{2}}\leq\sqrt{\sum_{i=1}^{n}(\underset{i}{\max}\{\lvert\epsilon_{i}^{k}\rvert\})^{2}}=\sqrt{n\times(\underset{i}{\max}\{\lvert\epsilon_{i}^{k}\rvert\})^{2}}=\sqrt{n}\lVert\boldsymbol{\epsilon}^{k}\rVert_{\infty}$. It follows that, ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k})-b_{i}\leq\lVert\mathbf{a}^{i}\rVert_{2}(\lVert\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k}\rVert_{2})-b_{i}\leq\lVert\mathbf{a}^{i}\rVert_{2}(\lVert\mathbf{\hat{x}}^{k}\rVert_{2}+\lVert\boldsymbol{\epsilon}^{k}\rVert_{2})-b_{i}\leq\lVert\mathbf{a}^{i}\rVert_{2}(\lVert\mathbf{\hat{x}}^{k}\rVert_{2}+\sqrt{n}\lVert\boldsymbol{\epsilon}^{k}\rVert_{\infty})-b_{i}$. From the previous paragraph, we have $\lVert\boldsymbol{\epsilon}^{k}\rVert_{\infty}\leq d_{k}$ for each $k\in\mathcal{K}$. Hence, choosing M_{1}=\underset{i\in\mathcal{I}}{\max}\{\lVert\mathbf{a}^{i}\rVert_{2}\}\big{(}\underset{k\in\mathcal{K}}{\max}\{\lVert\hat{\mathbf{x}}^{k}\rVert_{2}\}+\sqrt{n}\,\underset{k\in\mathcal{K}}{\max}\{d_{k}\}\big{)}+\underset{i\in\mathcal{I}}{\max}\{-b_{i}\}) is valid for each constraint of $\eqref{eq:IO_mip.4}$. As a more efficient and constraint-specific alternative, one may choose $M_{1ik}=\lVert\mathbf{a}^{i}\rVert_{2}(\lVert\mathbf{\hat{x}}^{k}\rVert_{2}+\sqrt{n}d_{k})-b_{i}$ for constraint $\eqref{eq:IO_mip.4}$ for each $i\in\mathcal{I}$ and $k\in\mathcal{K}$.

The $M$ parameter in formulation (8), with respect to each $k\in\mathcal{K}$, can be calculated similarly to $M_{1}$ with $d_{k}$ replaced by $\tau$, i.e., M=\underset{i\in\mathcal{I}}{\max}\{\lVert\mathbf{a}^{i}\rVert_{2}\}\big{(}\lVert\hat{\mathbf{x}}^{k}\rVert_{2}+\tau\,\sqrt{n}\big{)}+\underset{i\in\mathcal{I}}{\max}\{-b_{i}\}).

7.3 Computing the Inverse Stability Measure

Recall that inverse stability measure (3) is defined as the minimum distance of data shift such that all inverse-feasible cost vectors are no longer feasible for the inverse model with respect to the shifted data. Given a dataset $\hat{\mathcal{X}},$ let $\hat{\mathcal{C}}$ be the set of inverse-feasible cost vectors for QIO($\hat{\mathcal{X}},\tau,\theta$). As Theorem 3.4 suggests, every inverse-feasible cost vector $\mathbf{c}\in\hat{\mathcal{C}}$ is a conic combination of some $\mathbf{a}^{i}$’s that are also inverse-feasible. The following result shows that, to check whether or not all inverse-feasible cost vectors in $\hat{\mathcal{C}}$ are no longer feasible under the shifted data, it is sufficient to check whether or not all inverse-feasible $\mathbf{a}^{i}$’s whose conic hull forms $\hat{\mathcal{C}}$ are no longer feasible.

Lemma 7.2

Let $\hat{\mathcal{C}}$ be the set of inverse-feasible cost vectors for QIO$(\hat{\mathcal{X}},\tau,\theta)$ and $\hat{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\mathbf{a}^{i}\in\hat{\mathcal{C}}\}$. Let $\tilde{\mathcal{X}}$ denote the shifted dataset and $\tilde{\mathcal{C}}$ the set of inverse-feasible cost vectors for QIO$(\tilde{\mathcal{X}},\tau,\theta)$. If no $\mathbf{a}^{i},i\in\hat{\mathcal{I}}$, is inverse-feasible for QIO$(\tilde{\mathcal{X}},\tau,\theta)$, then no $\mathbf{c}\in\hat{\mathcal{C}}$ is inverse-feasible for QIO$(\tilde{\mathcal{X}},\tau,\theta)$.

Lemma 7.2 suggests that to compute the inverse stability measure for $\hat{\mathcal{C}}$ given the original dataset $\hat{\mathcal{X}}$, all we need is to find a minimally shifted set of the data points, $\tilde{\mathcal{X}}$, such that no $\mathbf{a}^{i},i\in\hat{\mathcal{I}}$, is inverse-feasible for QIO($\tilde{\mathcal{X}},\tau,\theta$). Given $\hat{\mathcal{I}}$, to make $\mathbf{c}=\mathbf{a}^{i}$ infeasible for QIO($\tilde{\mathcal{X}},\tau,\theta$) for all $i\in\hat{\mathcal{I}}$, we need to ensure that there are more than $\lfloor(1-\theta)K\rfloor$ data points in $\tilde{\mathcal{X}}$ whose minimum distance from each $\mathcal{X}_{i}$, $i\in\hat{\mathcal{I}}$, is strictly greater than $\tau$. This observation leads to the following formulation.

[TABLE]

*In the above formulation, variable $\mathbf{x}^{k}$ is considered a shifted data point from $\hat{\mathbf{x}}^{k}$. In constraint (12b), $D^{i}(\mathbf{x}^{k})$ computes the minimum distance between some shifted data point $\mathbf{x}^{k}$ and $\mathcal{X}_{i}$, which is then ensured to be no less than $\tau+\sigma$ for at least $\lfloor(1-\theta)K\rfloor+1$ data points by constraints (12c) and (12d) where $\sigma$ denotes an infinitesimal positive constant. The objective function minimizes the sum of distances between the original data points and shifted ones. As computing the exact inverse stability value would have required $D^{i}(\mathbf{x}^{k})>\tau$ in place of constraint (12c), the optimal value of the above problem is only infinitesimally greater than the exact inverse stability value, and we treat this value as the true inverse stability value (the inverse stability measure is defined on the one-dimensional open set). Note that problem (12) is always feasible as one can let $\mathbf{x}^{k},k\in\mathcal{K}$, be arbitrarily far from each $\mathcal{X}_{i},i\in\hat{\mathcal{I}}$. Problem (12) is non-convex and thus it is still challenging to compute its optimal value exactly. However, we use this formulation to compute the lower bound for the true inverse stability value, which is easy to compute (see Proposition 3.6 and its proof).

---

7.4 Upper Bound for the Forward Stability Measure

In this section, we show that an optimal solution to the MQIO model leads to an upper bound on forward stability measure (4). We make the following assumption: for each $i\in\mathcal{I}$, there exists a constant $\rho_{i}$ such that $\|\mathbf{x}^{1}-\mathbf{x}^{2}\|_{\ell}\leq\rho_{i}$ for any $\mathbf{x}^{1},\mathbf{x}^{2}\in\mathcal{X}_{i}$ and $\ell\geq 1$. The following result shows that if such $\rho_{i}$ values are known a priori, an upper bound for the forward stability measure for a cost vector found by the MQIO model can be computed.

Proposition 7.3

Let $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}\}_{k\in\mathcal{K}})$ be a feasible solution to $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$, $\bar{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\bar{v}_{i}=1\}$, $\bar{\mathcal{S}}=\{k\in\mathcal{K}\,|\,\bar{u}_{k}=1\}$, and $\bar{\mathbf{c}}=\sum_{i\in\bar{\mathcal{I}}}\lambda_{i}\mathbf{a}^{i}$ where $\lambda_{i}>0$ for all $i\in\bar{\mathcal{I}}$. Then, $\displaystyle\max_{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\underset{\mathbf{x}\in{\cap}_{i\in\bar{\mathcal{I}}}\mathcal{X}_{i}}{\min}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}+\displaystyle\min_{i\in\bar{\mathcal{I}}}\;\{\rho_{i}\}$ for all $k\in\bar{\mathcal{S}}$.

Because $\displaystyle\max_{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$ is bounded above for each $k\in\bar{\mathcal{S}}$, the forward stability measure given any distance function $d$ in (4), e.g., $\displaystyle\max_{k\in\bar{\mathcal{S}}}\max_{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$, is also bounded above.**

*Given the same inverse-feasible cost vector $\bar{\mathbf{c}}$ in Proposition 7.3, if $\ell=\infty$, the exact forward stability measure (4) can be written as \displaystyle\underset{k\in\bar{\mathcal{S}}}{\max}\;\Big{\{}\displaystyle\max_{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}\;\{\|\mathbf{x}-\hat{\mathbf{x}}^{k}\|_{\infty}\}\Big{\}}. The optimal value of the inner max problem for each $k$ can be obtained by solving $2n$ LPs: $z^{(1)}_{kj}=\underset{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}{\max}\;\{x_{j}-\hat{x}^{k}_{j}\}$ and $z^{(2)}_{kj}=\underset{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}{\max}\;\{\hat{x}^{k}_{j}-x_{j}\}$ for each $j\in\mathcal{J}=\{1,\ldots,n\}$ and selecting the largest value among them, i.e., $\displaystyle\max_{j\in\mathcal{J}}\,\{z^{(1)}_{kj},z^{(2)}_{kj}\}$. Thus, we have the exact forward stability measure: \displaystyle\underset{k\in\bar{\mathcal{S}}}{\max}\;\Big{\{}\displaystyle\max_{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}\;\{\|\mathbf{x}-\hat{\mathbf{x}}^{k}\|_{\infty}\}\Big{\}}=\displaystyle\max_{k\in\bar{\mathcal{S}}}\;\max_{j\in\mathcal{J}}\,\{z^{(1)}_{kj},z^{(2)}_{kj}\}.

---

7.5 Extensions to the MQIO model

In Section 3, we assume the constraint parameters $(\mathbf{A},\mathbf{b})$ remain the same over different observations. However, the proposed formulations can also apply when $\mathbf{b}$ varies over $k$ (i.e., $\mathbf{b}^{k}$ for observation $k$) as they only exploit the structure of $\mathbf{A}$. Should $\mathbf{A}$ vary over $k$ (i.e., $\mathbf{A}^{k}$ for observation $k$), an inverse-feasible cost vector can be obtained by replacing constraint (6b) by $b_{i}^{k}\leq{\mathbf{a}^{ki}}^{\prime}(\mathbf{\hat{x}}^{k}-\boldsymbol{\epsilon}^{k})\leq b_{i}^{k}+M_{1}(1-v_{i}^{k}),\forall i\in\mathcal{I},\forall k\in\mathcal{K}$, constructing observation-specific conic hulls $\mathcal{C}_{k}$ using $\mathbf{v}^{k}$ similarly as in Section 3, and finding the intersection of the multiple conic hulls (e.g., via the following additional constraint: $\mathbf{c}=\sum_{i\in\mathcal{I}}\lambda_{i}^{k}\mathbf{a}^{ki}$ and $\mathbf{0}\leq\boldsymbol{\lambda}^{k}\leq\mathbf{v}^{k}$ for all $k\in\mathcal{S}$).

*Additionally, we note that the MQIO model is also more stable than the previous inverse LP method when the parameters ($\mathbf{A}$,$\mathbf{b}$) are subject to uncertainty (i.e., the realized parameters may turn out different from the assumed ones). Suppose that both the previous inverse LP and our MQIO models are used to infer cost vectors given some “nominal” parameters $(\bar{\mathbf{A}},\bar{\mathbf{b}})$. Since the previous model finds a cost vector that is identical to one of the $\mathbf{a}^{i}$’s, a slight change in the corresponding $\mathbf{a}^{i}$ can make the inferred cost vector no longer inverse-feasible. On the other hand, if we let $\bar{\mathcal{C}}$ be the set of inverse-feasible cost vectors from the original MQIO model, there exists $\mathbf{c}\in\bar{\mathcal{C}}$ that is still inverse-feasible for the realized MQIO problem (with $(\tilde{\mathbf{A}},\tilde{\mathbf{b}})\neq(\bar{\mathbf{A}},\bar{\mathbf{b}})$) as long as there exists $\mathbf{x}^{k}$ within $\tau$-distance from at least $\lceil\theta K\rceil$ data points such that $\mathbf{x}^{k}\in\cap_{i\in\tilde{\mathcal{I}}}\mathcal{X}_{i}$ for some set $\tilde{\mathcal{I}}\subseteq\mathcal{I}$ and $\bar{\mathcal{C}}\cap\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\tilde{\mathcal{I}}})\neq\emptyset$. In general, if $(\mathbf{A},\mathbf{b})$ is subject to uncertainty, the MQIO problem lends itself to traditional robust optimization or stochastic programming techniques as it is a linear MIP.

---

Finally, should the cost vector be constrained, say $\mathbf{D}\mathbf{c}\leq\mathbf{d}$, the MQIO formulation can be rewritten as follows:

[TABLE]

where the first three constraints are added to ensure that the resulting $\mathbf{c}$ satisfies the cost vector constraint and is still a conic combination of some $\mathbf{a}^{i}$’s.

8 Supplemental Materials for the Algorithms

8.1 The Biclique Problem Formulation

Let $\bar{\mathbf{D}}\in\{0,1\}^{K\times m}$ be the bi-adjacency matrix defined as (9). The following MIP finds an all-one submatrix of $\bar{\mathbf{D}}$ (i.e., a biclique) with at least $\lceil\theta K\rceil$ rows (i.e., data points) and as many columns ($\mathbf{a}^{i}$’s) as possible.

[TABLE]

*In (14), $u_{k}^{clq}$ and $v_{i}^{clq}$ are 1 if node $k\in V_{1}$ and $i\in V_{2}$ are selected to be in the clique, respectively, and 0 otherwise. The objective function maximizes the number of nodes selected from $V_{2}$ (i.e., columns of the all-one submatrix of $\bar{\mathbf{D}}$). Constraint (14b) ensures that at most one of the nodes can be selected from $k\in V_{1}$ and $i\in V_{2}$ if there is no edge corresponding to $(k,i)$, and constraint (14c) guarantees that at least $\lceil\theta K\rceil$ nodes from $V_{1}$ (i.e., rows of the submatrix) are selected. Once problem (14) is solved, the submatrix can be retrieved by selecting rows and columns of $\bar{\mathbf{D}}$ corresponding to $\{k\in\mathcal{K}\,|\,u_{k}^{clq}=1\}$ and $\{i\in\mathcal{I}\,|\,v_{i}^{clq}=1\}$, respectively. ***

8.2 Alternative Heuristic

The algorithms proposed in Section 4.3 can still be computationally burdensome because the construction of each $\bar{\mathbf{D}}$ requires solving the MIP problem (8) $K$ times. To address this, we propose an LP relaxation of problem (8) using the idea of weighted $\ell_{1}$ minimization (Candes et al. 2008). For each $k\in\mathcal{K}$ we solve the following LP formulation:

[TABLE]

Then we construct a matrix $\tilde{\mathbf{D}}\in\{0,1\}^{K\times m}$:

[TABLE]

*where $\boldsymbol{\alpha}^{k}$ is a feasible solution to (15) with respect to $\mathbf{\hat{x}}^{k}$. Note that the construction of $\tilde{\mathbf{D}}$ is more efficient than that of $\bar{\mathbf{D}}$. ** The following result shows the relationship between a solution to problem (15) and a solution to problem (8). **

---

Lemma 8.1

Given an observation $\hat{\mathbf{x}}^{k}$, let $(\bar{\boldsymbol{\alpha}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ be a feasible solution to problem (15) and $\bar{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\bar{\alpha}^{k}_{i}=0\}$. Then, there exists a feasible solution $(\bar{\mathbf{v}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ to problem (8) with respect to $\hat{\mathbf{x}}^{k}$ where $\bar{v}^{k}_{i}=1$ for $i\in\bar{\mathcal{I}}$ and $\bar{v}^{k}_{i}=0$ for $i\in\mathcal{I}\setminus\bar{\mathcal{I}}$.

Next, similar to Proposition 4.4, the following result shows how a feasible solution to the QIO (hence MQIO) model can be retrieved from $\tilde{\mathbf{D}}$.

Proposition 8.2

If there exists $\tilde{\mathbf{D}}$ that satisfies (16) and has an all-one submatrix whose rows and columns correspond to $\bar{\mathcal{S}}$ and $\bar{\mathcal{A}}$, respectively, where $\bar{\mathcal{S}}\subseteq\mathcal{K}$, $|\bar{\mathcal{S}}|\geq\theta K$, and $\bar{\mathcal{A}}\subseteq\mathcal{I}$, then there exists a solution $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ such that $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$.

We use the same instance provided in Example 4.5 to illustrate the result of Proposition 8.2.

Example 8.3

Let $(\bar{\boldsymbol{\alpha}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ denote the optimal solution to problem (15) associated with each data point $k=1,\ldots,5$: we have $\bar{\boldsymbol{\alpha}}^{k}=[0,\,0,\,\bar{\alpha}_{3}^{k}>0,\,\bar{\alpha}_{4}^{k}>0]^{\prime}$ for $k=1,2,3,4,$ (with $\bar{\boldsymbol{\epsilon}}^{1}=\begin{bmatrix}-0.5\\ -0.2\end{bmatrix}$, $\bar{\boldsymbol{\epsilon}}^{2}=\begin{bmatrix}-0.3\\ -0.2\end{bmatrix}$, $\bar{\boldsymbol{\epsilon}}^{3}=\begin{bmatrix}-0.3\\ -0.5\end{bmatrix}$, and $\bar{\boldsymbol{\epsilon}}^{4}=\begin{bmatrix}-0.5\\ -0.5\end{bmatrix}$), and $\bar{\boldsymbol{\alpha}}^{5}=[\bar{\alpha}_{1}^{5}>0,\,0,\,0,\,\bar{\alpha}_{4}^{5}>0]^{\prime}$ (with $\bar{\boldsymbol{\epsilon}}^{5}=\begin{bmatrix}-0.3\\ 0.3\end{bmatrix}$). We then construct a bipartite graph by defining node $i$ for $\mathbf{a}^{i}$, node $k$ for $\hat{\mathbf{x}}^{k}$, and edge $(k,i)$ where $\bar{\alpha}_{i}^{k}=0$, which leads to the same bipartite graph in Figure 2(b). As a result, the corresponding bi-adjacency matrix $\tilde{\mathbf{D}}$ is identical to the $\bar{\mathbf{D}}$ matrix in Figure 2(c). Similar to Example 4.5, we can obtain an inverse-feasible cost vector from the all-one submatrix of $\tilde{\mathbf{D}}$.

Similar to the algorithms in Section 4.3, there can be multiple $\tilde{\mathbf{D}}$’s because there may be multiple optimal solutions to problem (15). We propose a heuristic, which we call $\tilde{\mathbf{D}}$-Alg-Heuristic, via a similar weighting approach to find a matrix $\tilde{\mathbf{D}}$ and the same clique problem to find the all-one submatrix of $\tilde{\mathbf{D}}$ (i.e., Clique($\tilde{\mathbf{D}},\theta$)). The pseudo-code for this heuristic can be found in Algorithm 3.

8.3 Application to Online Learning

The assumption we made in Section 3 is that the data points are available all at once in advance from a fixed forward LP formulation. However, in many application domains, data may become available through different time points in separate batches. Moreover, some parameters of the underlying forward LP formulation may change at each time point as well. In this subsection, we show how our inverse method can be extended for learning the cost vectors adaptively over time in an online manner.

Suppose at each time $t\in\mathcal{T}=\{1,\ldots,T\}$, the DM observes an external signal $\mathbf{b}=\mathbf{b}^{t}$ as the right-hand-side of the forward problem and makes a set of decisions $\hat{\mathcal{X}}_{t}$ with the index set $\mathcal{K}_{t}$ of cardinality $K_{t}$. We assume that the decisions can be noisy. The following formulation is a modification to MQIO that finds inverse-feasible cost vectors for the entire collection of batches $\bigcup_{t\in\mathcal{T}}\hat{\mathcal{X}}_{t}$:

[TABLE]

The above formulation still aims to solve the inverse problem in the traditional setting, which means it has to wait until all data points are collected, leading to a large-size MIP. Instead, we propose an online learning approach using a modified version of $\bar{\mathbf{D}}$-Alg-Heuristic in Algorithm 4.

In this algorithm, we start with an empty $\bar{\mathbf{D}}$ matrix and update it at each time $t$ by finding $\mathbf{v}^{*}\in\underset{\mathbf{v}}{\text{argmax}}\{{\mathbf{w}^{t}}^{\prime}\mathbf{v}\,|\,\eqref{eq:IO_mip.2}$–$\eqref{eq:IO_mip.5}\}$* with $\hat{\mathcal{X}}_{t}$ as input data (where $\mathbf{w}^{t}$ is an appropriate weight vector at time $t$) and inserting ${\mathbf{v}^{*}}^{\prime}$ into the next empty row of $\bar{\mathbf{D}}$. At each time $t$, given the updated $\bar{\mathbf{D}}$ we find its largest all-one submatrix and let $\mathbf{v}^{t}$ be a binary vector with 1’s corresponding to the columns of the submatrix. Recall that each column of this matrix corresponds to each $\mathbf{a}^{i}$; thus, $\mathbf{v}^{t}$ indicates which $\mathbf{a}^{i}$’s should be used for creating the set of inverse-feasible cost vectors at time $t$, i.e., $\mathcal{C}_{t}=\textup{cone}(\{\mathbf{a}^{i}\}_{i:v_{i}^{t}=1})$. We show in the numerical results section that with a sufficiently large $T$, there exists $\mathbf{c}\in\mathcal{C}_{T}$ that is also inverse-feasible for the above MIP (17). Again, with such a large $T$ and an excessive number of data points, the MIP problem (17) would have been computationally extremely challenging, whereas the online learning approach is more efficient as it attempts to solve the problem in a distributed manner.*

9 Proofs

Proof 9.1

Proof of Proposition 3.1: Given $\hat{\mathcal{X}}$, $\tau$, $\theta_{1}$, and $\theta_{2}$, let $\hat{\mathcal{C}}(\theta_{1})$, $\hat{\mathcal{C}}(\theta_{2})$, $\tilde{\mathcal{C}}(\theta_{1})$, and $\tilde{\mathcal{C}}(\theta_{2})$ denote the set of cost vectors obtained by $\mathbf{QIO}(\hat{\mathcal{X}},\tau,\theta_{1})$, $\mathbf{QIO}(\hat{\mathcal{X}},\tau,\theta_{2})$, $\mathbf{QIO}(\tilde{\mathcal{X}},\tau,\theta_{1})$, and $\mathbf{QIO}(\tilde{\mathcal{X}},\tau,\theta_{2})$, respectively, where $\tilde{\mathcal{X}}$ denotes some shifted dataset. To show that inverse stability measure (3) for $\hat{\mathcal{C}}(\theta)$ is non-decreasing as $\theta$ decreases, we need to show that given $\theta_{1}$ and $\theta_{2}$, if $\theta_{2}\leq\theta_{1}$, then $\underset{\tilde{\mathcal{X}}\subset\mathbb{R}^{K\times n}}{\min}\{d(\hat{\mathcal{X}},\tilde{\mathcal{X}})|\ \hat{\mathcal{C}}(\theta_{2})\cap\tilde{\mathcal{C}}(\theta_{2})=\emptyset\}\geq\underset{\tilde{\mathcal{X}}\subset\mathbb{R}^{K\times n}}{\min}\{d(\hat{\mathcal{X}},\tilde{\mathcal{X}})|\ \hat{\mathcal{C}}(\theta_{1})\cap\tilde{\mathcal{C}}(\theta_{1})=\emptyset\}$.  Given $\theta_{2}\leq\theta_{1}$, we have $\hat{\mathcal{C}}(\theta_{1})\subseteq\hat{\mathcal{C}}(\theta_{2})$ and $\tilde{\mathcal{C}}(\theta_{1})\subseteq\tilde{\mathcal{C}}(\theta_{2})$. Let $\hat{\mathcal{C}}_{\text{diff}}=\hat{\mathcal{C}}(\theta_{2})\setminus\hat{\mathcal{C}}(\theta_{1})$ and $\tilde{\mathcal{C}}_{\text{diff}}=\tilde{\mathcal{C}}(\theta_{2})\setminus\tilde{\mathcal{C}}(\theta_{1})$. Then,

[TABLE]

which completes the proof. \Halmos

Proof 9.2

Proof of Proposition 3.2: Let $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ be a feasible solution to $\textrm{{QIO}}(\mathcal{K},\tau,\theta)$. Note that $\bar{y}_{i}>0$ for at least one $i\in\mathcal{I}$. Let $\bar{\mathcal{I}}=\{i\in\mathcal{I}\ \arrowvert\ \bar{y}_{i}>0\}$. From (5b), $\bar{\mathbf{c}}=\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}\mathbf{a}^{i}$. From (5e), for each $k\in\bar{\mathcal{S}}$, $\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}{\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})=\sum_{i\in\bar{\mathcal{I}}}b_{i}\bar{y}_{i}$ and thus $\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}({\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})-b_{i})=0$. Because $\bar{y}_{i}>0$, we have ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})-b_{i}=0,\forall i\in\bar{\mathcal{I}}$. Now, pick any arbitrary $\tilde{i}\in\bar{\mathcal{I}}$ and let $\tilde{\mathbf{y}}=\mathbf{e}_{\tilde{i}}$ where $\mathbf{e}_{\tilde{i}}$ denotes the $\tilde{i}$-th unit vector. Then, from (5b), we can construct a new cost vector $\tilde{\mathbf{c}}=\mathbf{a}^{\tilde{i}}$. By replacing $\bar{\mathbf{c}}$ with this $\tilde{\mathbf{c}}$ in (5e), (5e) becomes ${\mathbf{a}^{\tilde{i}}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})-b_{\tilde{i}}=0,\forall k\in\bar{\mathcal{S}}$, which still holds because $\tilde{i}\in\bar{\mathcal{I}}$. It is clear that $(\tilde{\mathbf{c}}=\mathbf{a}^{\tilde{i}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\tilde{\mathbf{y}}=\mathbf{e}_{\tilde{i}},\bar{\mathcal{S}})$ also satisfies constraints (5b)–(5i) and thus is feasible for (5). \Halmos

Proof 9.3

Proof of Lemma 3.3: Given some sets $\bar{\mathcal{A}}\subseteq\mathcal{A}$ and $\bar{\mathcal{S}}$ such that $|\bar{\mathcal{S}}|\geq\theta K$, let $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$ such that $\|\bar{\mathbf{c}}\|_{p}=1$. Recall that $\mathcal{X}_{i}=\{\mathbf{x}\in\mathcal{X}\,|\,{\mathbf{a}^{i}}^{\prime}\mathbf{x}=b_{i}\},i\in\mathcal{I}$. We want to show that if there exists $\bar{\mathbf{x}}^{k}\in\mathcal{X}_{i}$ for all $i\in\bar{\mathcal{A}}$ such that $\|\bar{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}\leq\tau$ for all $k\in\bar{\mathcal{S}}$, then there exists $\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}}$ and $\bar{\mathbf{y}}$ such that $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ is feasible for QIO$(\mathcal{K},\tau,\theta)$ (i.e., formulation (5)). Note that $\bar{\mathbf{c}}$ can be written as $\bar{\mathbf{c}}=\sum_{i\in\bar{\mathcal{A}}}\bar{\lambda}_{i}\mathbf{a}^{i}$ where $\bar{\lambda}_{i}>0$ for at least one $i\in\bar{\mathcal{A}}$. Let $\bar{\lambda}_{i}=0$ for $i\in\mathcal{I}\setminus\bar{\mathcal{A}}$. Because $\bar{\mathbf{x}}^{k}\in\mathcal{X}^{i},\forall i\in\bar{\mathcal{A}}$, for each $k\in\bar{\mathcal{S}}$ we have ${\mathbf{a}^{i}}^{\prime}\bar{\mathbf{x}}^{k}=b_{i}$ for all $i\in\bar{\mathcal{A}}$. Multiplying both sides by $\bar{\lambda}_{i}$ and summing over all $i\in\bar{\mathcal{A}}$, we have $\sum_{i\in\bar{\mathcal{A}}}\bar{\lambda}_{i}{\mathbf{a}^{i}}^{\prime}\bar{\mathbf{x}}^{k}=\sum_{i\in\bar{\mathcal{A}}}\bar{\lambda}_{i}b_{i}$, and because $\bar{\lambda}_{i}=0$ for $i\in\mathcal{I}\setminus\bar{\mathcal{A}}$ we further have $\sum_{i\in\mathcal{I}}\bar{\lambda}_{i}{\mathbf{a}^{i}}^{\prime}\bar{\mathbf{x}}^{k}=\sum_{i\in\mathcal{I}}\bar{\lambda}_{i}b_{i}$, hence $\bar{\mathbf{c}}^{\prime}\bar{\mathbf{x}}^{k}=\mathbf{b}^{\prime}\bar{\boldsymbol{\lambda}}$ for all $k\in\bar{\mathcal{S}}$. By letting $\bar{\mathbf{y}}=\bar{\boldsymbol{\lambda}}$, we have $\mathbf{A}^{\prime}\bar{\mathbf{y}}=\bar{\mathbf{c}}$ and $\bar{\mathbf{y}}\geq\mathbf{0}$, which satisfy constraints (5b) and (5c), respectively. For $k\in\bar{\mathcal{S}}$, let $\bar{\boldsymbol{\epsilon}}^{k}=\hat{\mathbf{x}}^{k}-\bar{\mathbf{x}}^{k}$; for $k\in\mathcal{K}\setminus\bar{\mathcal{S}}$, let $\bar{\boldsymbol{\epsilon}}^{k}$ be an arbitrary vector in $\mathbb{R}^{n}$. Then constraint (5d) is satisfied because $\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k}=\bar{\mathbf{x}}^{k}\in\mathcal{X},\forall k\in\bar{\mathcal{S}}$, and constraint (5e) is satisfied because $\bar{\mathbf{c}}^{\prime}(\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})=\bar{\mathbf{c}}^{\prime}\bar{\mathbf{x}}^{k}=\mathbf{b}^{\prime}\bar{\boldsymbol{\lambda}}=\mathbf{b}^{\prime}\bar{\mathbf{y}},\forall k\in\bar{\mathcal{S}}$. Constraint (5f) is satisfied because $\|\bar{\boldsymbol{\epsilon}}^{k}\|_{\ell}=\|\hat{\mathbf{x}}^{k}-\bar{\mathbf{x}}^{k}\|_{\ell}\leq\tau,\forall k\in\bar{\mathcal{S}}$. Finally, (5g), (5h), and (5i) are satisfied as it was assumed that $\bar{\mathcal{S}}\subseteq\mathcal{K}$, $|\bar{\mathcal{S}}|\geq\theta K$ and $\|\bar{\mathbf{c}}\|_{p}=1$. As a result, $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ is feasible for QIO$(\mathcal{K},\tau,\theta)$.\Halmos

Proof 9.4

Proof of Theorem 3.4: $(\Rightarrow)$ Assume $(\bar{\mathbf{c}}\neq\mathbf{0},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ is feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$, i.e., model (5). Then from (5b), $\bar{\mathbf{c}}=\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}\mathbf{a}^{i}$ where $\bar{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\bar{y}_{i}>0\}$. Construct $\bar{\mathbf{v}}\in\mathbb{R}^{|\mathcal{I}|}$ by setting $\bar{v}_{i}=1$ for all $i\in\bar{\mathcal{I}}$ and $\bar{v}_{i}=0$ for all $i\in\mathcal{I}\setminus\bar{\mathcal{I}}$. Also define $\bar{\mathbf{u}}\in\mathbb{R}^{|\mathcal{K}|}$ such that $\bar{u}_{k}=1$ for $k\in\bar{\mathcal{S}}$ and $\bar{u}_{k}=0$ for $k\in\mathcal{S}\setminus\bar{\mathcal{S}}$. Let $\boldsymbol{\epsilon}^{k}=\bar{\boldsymbol{\epsilon}}^{k}$ for $k\in\bar{\mathcal{S}}$ and $\boldsymbol{\epsilon}^{k}=\hat{\mathbf{x}}^{k}-\mathbf{x}$ for $k\in\mathcal{K}\setminus\bar{\mathcal{S}}$ where $\mathbf{x}\in\mathcal{X}$ satisfying ${\mathbf{a}^{i}}^{\prime}\mathbf{x}=b_{i},\forall i\in\bar{\mathcal{I}}$. We want to show that $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ is feasible for model (6). From (5b) and (5e), we have $\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}{\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})=\sum_{i\in\bar{\mathcal{I}}}b_{i}\bar{y}_{i}$, or equivalently $\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}({\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})-b_{i})=0$ for each $k\in\bar{\mathcal{S}}$. Since $\bar{y}_{i}>0$ for all $i\in\bar{\mathcal{I}}$, for each $k\in\bar{\mathcal{S}}$ we have ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})=b_{i}$ for $i\in\bar{\mathcal{I}}$ (i.e., $i:\bar{v}_{i}=1$) and ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})\geq b_{i}$ for $i\in\mathcal{I}\setminus\bar{\mathcal{I}}$ (i.e., $i:\bar{v}_{i}=0$). Also, for each $k\in\mathcal{K}\setminus\bar{\mathcal{S}}$, based on how $\boldsymbol{\epsilon}^{k}$ is defined above, we also have ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})=b_{i}$ for all $i\in\bar{\mathcal{I}}$ and ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})\geq b_{i}$ for all $i\in\mathcal{I}$. Therefore, $\{{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}}$ and $\bar{\mathbf{v}}$ together satisfy (6b). Because $\bar{\mathbf{u}}$ is constructed based on whether or not $\|\bar{\boldsymbol{\epsilon}}^{k}\|_{\ell}\leq\tau$, $\{{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}}$ and $\bar{\mathbf{u}}$ clearly satisfy (6c). Also, $\sum_{k=1}^{K}\bar{u}_{k}=\sum_{k\in\bar{\mathcal{S}}}\bar{u}_{k}=|\bar{\mathcal{S}}|\geq\theta K$, satisfying (6d). Thus, $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ is feasible for model (6). Finally, from (5b), $\bar{\mathbf{c}}=\mathbf{A}^{\prime}\bar{\mathbf{y}}=\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}\mathbf{a}^{i}=\sum_{i:\bar{v}_{i}=1}\bar{y}_{i}\mathbf{a}^{i}$, and therefore $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i:\bar{v}_{i}=1})$. Because $\bar{\mathbf{c}}$ is feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$, $\lVert\bar{\mathbf{c}}\rVert_{p}=1$.

$(\Leftarrow)$* Assume $(\bar{\mathbf{v}}\neq\mathbf{0},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ is feasible for model (6). Let $\bar{\mathcal{I}}=\{i\in\mathcal{I}\ |\ \bar{v}_{i}=1\}$ and $\bar{\mathcal{S}}=\{k\in\mathcal{K}\ |\ \bar{u}_{k}=1\}$, and consider $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{I}}})$ such that $\|\bar{\mathbf{c}}\|_{p}=1$. We can write $\bar{\mathbf{c}}=\sum_{i\in\bar{\mathcal{I}}}\bar{\lambda}_{i}\mathbf{a}^{i}$ where $\bar{\lambda}_{i}\geq 0$ $\forall i\in\bar{\mathcal{I}}$. If we let $\bar{y}_{i}=\bar{\lambda}_{i}$ for all $i\in\bar{\mathcal{I}}$ and $\bar{y}_{i}=0$ for all $i\in\mathcal{I}\setminus\bar{\mathcal{I}}$ then $\bar{\mathbf{c}}=\sum_{i\in\bar{\mathcal{I}}}\bar{y}_{i}\mathbf{a}^{i}=\mathbf{A}^{\prime}\bar{\mathbf{y}}$ satisfies (5b); because $\bar{\mathbf{y}}\geq\mathbf{0}$, (5c) is also satisfied. Because ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})\geq b_{i}$ holds for all $i\in\mathcal{I},k\in\mathcal{K}$ by (6b), (5d) is satisfied for all $k\in\mathcal{K}$. Also because model (6) is feasible, ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})=b_{i}$ for all $i\in\bar{\mathcal{I}}$ and all $k$ for which $\bar{u}_{k}=1$, i.e., $k\in\bar{\mathcal{S}}$; multiplying both sides by $\bar{y}_{i}$ and summing over all $i\in\mathcal{I}$, we have $\sum_{i\in\mathcal{I}}\bar{y}_{i}{\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})=\sum_{i\in\mathcal{I}}\bar{y}_{i}b_{i}$ or equivalently $\bar{\mathbf{c}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})=\mathbf{b}^{\prime}\bar{\mathbf{y}}$, satisfying (5e). Clearly, $\{\bar{\boldsymbol{\epsilon}}_{k}\}_{k\in\mathcal{K}}$ and $\bar{\mathcal{S}}$ satisfy constraints (5f)–(5h). Therefore, $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ is feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$. *\Halmos

Proof 9.5

Proof of Proposition 3.6: For this proof, we first refer the reader to Appendix 7.3 where we show a formulation to compute the true inverse stability measure; this formulation is non-convex but we use it for proving Proposition 3.6. Note that the true inverse stability measure $\xi^{*}$ is the optimal value of problem (12) in Appendix 7.3.

We first present a variant of (12) for each $i\in\bar{\mathcal{I}}$ separately (i.e., $|\bar{\mathcal{I}}|$ of them) and show that the optimal objective value of (12), i.e., $\xi^{*}$, is bounded below by the largest optimal objective value among these $|\bar{\mathcal{I}}|$ smaller problems. Consider the following problem for each $i\in\bar{\mathcal{I}}$:

[TABLE]

Let $\zeta_{i}^{*}$ be the optimal value of the above problem for each $i\in\bar{\mathcal{I}}$. Note that unlike in formulation (12) where constraints of $D^{i}(\mathbf{x}^{k})\geq\tau+\sigma$ in (12c) must be satisfied for all $i\in\hat{\mathcal{I}}$, problem (19) requires constraint (19c) be satisfied for only one $i\in\bar{\mathcal{I}}\subseteq\hat{\mathcal{I}}$. Hence, we have $\zeta^{*}_{i}\leq\xi^{*}$ for all $i\in\bar{\mathcal{I}}$, or equivalently, $\displaystyle\max_{i\in\bar{\mathcal{I}}}\{\zeta^{*}_{i}\}\leq\xi^{*}$. Problem (19) is still non-convex and thus finding $\zeta_{i}^{*}$ is still a challenge. We next show that its optimal objective value $\zeta_{i}^{*}$ is equal to the sum of the optimal objective values of even smaller problems that can be efficiently solved. Given $i\in\bar{\mathcal{I}}$ and $k\in\mathcal{K}$, consider the following problem:

[TABLE]

*Let $\tilde{\zeta}_{ik}$ denote the optimal value of the above problem. For a given $i$, let $\tilde{\zeta}_{i[k]}$ be the $k$-th smallest value of $\tilde{\zeta}_{ik},k\in\mathcal{K}$. Then we have $\displaystyle\zeta_{i}^{*}=\sum_{k=1}^{\lfloor(1-\theta)K\rfloor+1}\tilde{\zeta}_{i[k]}$. Thus, \displaystyle\max_{i\in\bar{\mathcal{I}}}\bigg{\{}{\sum_{k=1}^{\lfloor(1-\theta)K\rfloor+1}\tilde{\zeta}_{i[k]}}\bigg{\}}=\max_{i\in\bar{\mathcal{I}}}\{\zeta_{i}^{*}\}\leq\xi^{*}.

---

Next, we show how we can find the optimal value for formulation (20), i.e., $\tilde{\zeta}_{ik}$. Let $d_{ik}^{*}=\underset{\mathbf{x}\in\mathcal{X}_{i}}{\min}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$ and $\bar{\tau}=\tau+\sigma$. Note that when $d_{ik}^{*}\geq\bar{\tau}$, we have $\tilde{\zeta}_{ik}=0$ because $\mathbf{x}^{k}=\hat{\mathbf{x}}^{k}$ satisfies both (20b) and (20c). Now we want to show that if $d_{ik}^{*}<\bar{\tau}$, $\tilde{\zeta}_{ik}=\bar{\tau}-d_{ik}^{*}$. To do so, we first show that $\tilde{\zeta}_{ik}\geq\bar{\tau}-d^{*}_{ik}$. Let $\bar{\mathbf{x}}\in\underset{\mathbf{x}\in\mathcal{X}_{i}}{\text{argmin}}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$. Suppose to the contrary that there exists an optimal solution $\tilde{\mathbf{x}}^{k}$ that attains $\tilde{\zeta}_{ik}<\bar{\tau}-d^{*}_{ik}$. By the definition of $d^{*}_{ik}$ we can rewrite this inequality as $\|\tilde{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}<\bar{\tau}-\|\hat{\mathbf{x}}^{k}-\bar{\mathbf{x}}\|_{\ell}$ or equivalently $\|\tilde{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}+\|\hat{\mathbf{x}}^{k}-\bar{\mathbf{x}}\|_{\ell}<\bar{\tau}$. Furthermore, by triangular inequality the following is true: $\|\tilde{\mathbf{x}}^{k}-\bar{\mathbf{x}}\|_{\ell}\leq\|\tilde{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}+\|\hat{\mathbf{x}}^{k}-\bar{\mathbf{x}}\|_{\ell}<\bar{\tau}$. However, since $\bar{\mathbf{x}}\in\mathcal{X}_{i}$, we have $D^{i}(\tilde{\mathbf{x}}^{k})=\underset{\mathbf{z}\in\mathcal{X}_{i}}{\min}\{\|\tilde{\mathbf{x}}^{k}-\mathbf{z}\|_{\ell}\}\leq\|\tilde{\mathbf{x}}^{k}-\bar{\mathbf{x}}\|_{\ell}<\bar{\tau}$, which violates (20c) and thus contradicts the feasibility of $\tilde{\mathbf{x}}^{k}$. Thus, it must be that $\tilde{\zeta}_{ik}\geq\bar{\tau}-d^{*}_{ik}$. All we need now is to show that there exits a feasible solution $\tilde{\mathbf{x}}^{k}$ to problem (20) that attains the objective function value of $\bar{\tau}-d^{*}_{ik}$. Note that if a given $\tilde{\mathbf{x}}^{k}$ is feasible for problem (20), we must have $D^{i}(\tilde{\mathbf{x}}^{k})=\underset{\mathbf{z}\in\mathcal{X}_{i}}{\min}\left\{\|\tilde{\mathbf{x}}^{k}-\mathbf{z}\|_{\ell}\right\}$ bounded below by $\bar{\tau}$ (due to (20c)). Consider $\tilde{\mathbf{x}}^{k}=\frac{\bar{\tau}}{d^{*}_{ik}}\hat{\mathbf{x}}^{k}-\frac{\bar{\tau}}{d^{*}_{ik}}\bar{\mathbf{x}}+\bar{\mathbf{x}}$. Clearly, $\mathbf{x}^{k}=\tilde{\mathbf{x}}^{k}$ and $\mathbf{z}=\bar{\mathbf{x}}$ together satisfy all constraints in (20) with equality because we have $\bar{\mathbf{x}}\in\mathcal{X}_{i}$ and

[TABLE]

*Moreover, $\mathbf{x}^{k}=\tilde{\mathbf{x}}^{k}$ leads to the objective value $\tilde{\zeta}_{ik}=\|\hat{\mathbf{x}}^{k}-\tilde{\mathbf{x}}^{k}\|_{\ell}=\|\hat{\mathbf{x}}^{k}+\frac{\bar{\tau}}{d^{*}_{ik}}\hat{\mathbf{x}}^{k}-\frac{\bar{\tau}}{d^{*}_{ik}}\bar{\mathbf{x}}+\bar{\mathbf{x}}\|_{\ell}=\frac{\bar{\tau}-d^{*}_{ik}}{d^{*}_{ik}}\|\hat{\mathbf{x}}^{k}-\bar{\mathbf{x}}\|_{\ell}=\bar{\tau}-d^{*}_{ik}$ and thus is an optimal solution to (20).

---

As a result, \displaystyle\max_{i\in\bar{\mathcal{I}}}\bigg{\{}\sum_{k=1}^{(\lfloor(1-\theta)K\rfloor+1)}\!\!\!\!\!\max(0,\bar{\tau}-d^{*}_{i[k]})\bigg{\}}=\max_{i\in\bar{\mathcal{I}}}\bigg{\{}\sum_{k=1}^{(\lfloor(1-\theta)K\rfloor+1)}\!\!\!\!\!\tilde{\zeta}_{i[k]}\bigg{\}}=\max_{i\in\bar{\mathcal{I}}}\{\zeta_{i}^{*}\}\leq\xi^{*}. Clearly, because $\tau<\bar{\tau}$ we have \displaystyle\max_{i\in\bar{\mathcal{I}}}\bigg{\{}\sum_{k=1}^{(\lfloor(1-\theta)K\rfloor+1)}\!\!\!\!\!\max(0,\tau-d^{*}_{i[k]})\bigg{\}}\leq\xi^{*}, as desired. \Halmos**

Proof 9.6

Proof of Proposition 3.7: Given a feasible solution $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ for (6), let $\bar{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\bar{v}_{i}=1\}$ and $\bar{\mathcal{S}}=\{k\in\mathcal{K}\,|\,\bar{u}_{k}=1\}$. Due to constraints (6b) and (6c), for each $k\in\bar{\mathcal{S}}$ we have $\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k}\in\{\mathbf{x}\in\mathcal{X}\,|\,{\mathbf{a}^{i}}^{\prime}\mathbf{x}^{k}=b_{i},\forall i\in\bar{\mathcal{I}}\}=\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}$ and $\|\bar{\boldsymbol{\epsilon}}^{k}\|_{\ell}\leq\tau$. Now consider $\bar{\mathbf{c}}\in\textup{cone}_{+}(\{\mathbf{a}^{i}\}_{i:\bar{v}_{i}=1})$. This $\bar{\mathbf{c}}$ vector leads to $\mathcal{X}^{*}(\bar{\mathbf{c}})=\arg\min\;\textbf{FO}(\bar{\mathbf{c}})=\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}$. Recall that $\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k}\in\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}$; this means that $\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k}\in\mathcal{X}^{*}(\bar{\mathbf{c}})$ for all $k\in\bar{\mathcal{S}}$. Moreover, because $\mathcal{X}$ is an $n$-dimensional polytope and $|\bar{\mathcal{I}}|=n$, $\mathcal{X}^{*}(\bar{\mathbf{c}})$ is a singleton, i.e., $\mathcal{X}^{*}(\bar{\mathbf{c}})=\{\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k}\},\forall k\in\bar{\mathcal{S}}$ and $\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k}=\hat{\mathbf{x}}^{k^{\prime}}-\bar{\boldsymbol{\epsilon}}^{k^{\prime}}$ for all $k,k^{\prime}\in\mathcal{K}$. To complete the proof, for all $k\in\bar{\mathcal{S}}$ we have $\displaystyle\max_{\mathbf{x}\in\mathcal{X}^{*}(\bar{\mathbf{c}})}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}=\|\hat{\mathbf{x}}^{k}-(\hat{\mathbf{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})\|_{\ell}=\|\bar{\boldsymbol{\epsilon}}^{k}\|_{\ell}\leq\tau$.\Halmos

Proof 9.7

Proof of Proposition 3.8: (i) Let $\bar{\mathcal{H}}=\{\mathbf{a}^{i}:\bar{v}_{i}=1\}$ and $\tilde{\mathcal{H}}=\{\mathbf{a}^{i}:\tilde{v}_{i}=1\}$. Because $\bar{\mathbf{v}}\geq\tilde{\mathbf{v}}$, we have $\tilde{\mathcal{H}}\subseteq\bar{\mathcal{H}}$ and thus $\tilde{\mathcal{C}}\subseteq\bar{\mathcal{C}}$. Given the dataset $\hat{\mathcal{X}}$, we want to show that for any $\tilde{\mathbf{c}}\in\tilde{\mathcal{C}}$ there exists $\bar{\mathbf{c}}\in\bar{\mathcal{C}}$ such that $\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\bar{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\tilde{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$, $\forall k\in\mathcal{K}$. Suppose to the contrary that there exists $\mathbf{c}^{*}\in\tilde{\mathcal{C}}$ such that $\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\mathbf{c}^{*})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}<\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\mathbf{c})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\},\forall\mathbf{c}\in\bar{\mathcal{C}}$ for some $k\in\mathcal{K}$. This means $\mathbf{c}^{*}\in\tilde{\mathcal{C}}\setminus\bar{\mathcal{C}}$ and thus $\tilde{\mathcal{C}}\setminus\bar{\mathcal{C}}\neq\emptyset$, which contradicts $\tilde{\mathcal{C}}\subseteq\bar{\mathcal{C}}$. Thus, for any $\tilde{\mathbf{c}}\in\tilde{\mathcal{C}}$ there exists $\bar{\mathbf{c}}\in\bar{\mathcal{C}}$ such that $\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\bar{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\underset{\mathbf{x}\in\mathcal{X}^{*}\!(\tilde{\mathbf{c}})}{\max}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$, $\forall k\in\mathcal{K}$, as desired.

(ii) Because $\bar{\mathbf{v}}\geq\tilde{\mathbf{v}}$, we have $\tilde{\mathcal{H}}\subseteq\bar{\mathcal{H}}$ and $\tilde{\mathcal{C}}\subseteq\bar{\mathcal{C}}$. Because $\tilde{\mathcal{C}}\subseteq\bar{\mathcal{C}}$,

[TABLE]

Proof 9.8

Proof of Lemma 4.1: Assume that there exists an optimal solution $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$ for $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ such that $\sum_{k=1}^{K}\bar{u}_{k}>\lceil\theta K\rceil$. Without loss of generality let $\sum_{k=1}^{K}\bar{u}_{k}=\lceil\theta K\rceil+r$ where $r>0$. Pick any $r$ observations for which $\bar{u}_{k}=1$ and change their values to [math], and denote this new vector by $\tilde{\mathbf{u}}$. Then we can construct a new solution $(\bar{\mathbf{v}},\tilde{\mathbf{u}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}})$, which is still feasible for $\mathbf{MQIO}(\mathcal{K},\tau,\theta)$ because $\|\bar{\boldsymbol{\epsilon}}^{k}\|_{\ell}\leq\tau+M_{2}(1-\bar{u}_{k})\leq\tau+M_{2}(1-\tilde{u}_{k}),\forall k\in\mathcal{K}$ and $\sum_{k=1}^{K}\tilde{u}_{k}=\lceil\theta K\rceil\geq\theta K$, and generates the same objective value, hence also an optimal solution. \Halmos

Proof 9.9

*Proof of Proposition 4.4:

Consider a matrix $\bar{\mathbf{D}}$ constructed as in (9). By definition, each row $k\in\mathcal{K}$ of $\bar{\mathbf{D}}$ is a binary vector ${\mathbf{v}^{k}}^{\prime}$ and there exists $\bar{\boldsymbol{\epsilon}}^{k}$ such that $(\mathbf{v}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ is feasible for (8). Assume that there exists an all-one submatrix in $\bar{\mathbf{D}}$ whose rows and columns correspond to $\bar{\mathcal{S}}$ and $\bar{\mathcal{A}}$, respectively, such that $\bar{\mathcal{S}}\subseteq\mathcal{K}$, $|\bar{\mathcal{S}}|\geq\theta K$, and $\bar{\mathcal{A}}\subseteq\mathcal{I}$. Define $\bar{\mathbf{v}}$ such that $\bar{v}_{i}=1$ for all $i\in\bar{\mathcal{A}}$ and $\bar{v}_{i}=0$ otherwise. Then we have $\bar{\mathbf{v}}\leq\mathbf{v}^{k},\forall k\in\bar{\mathcal{S}}$, and thus $(\bar{\mathbf{v}},\bar{\boldsymbol{\epsilon}}^{k})$ is also feasible for (8) for all $k\in\bar{\mathcal{S}}$. Define $\bar{\mathbf{u}}$ such that $\bar{u}_{k}=1$ for all $k\in\bar{\mathcal{S}}$ and $\bar{u}_{k}=0$ otherwise; note that $\sum_{k=1}^{K}\bar{u}_{k}\geq\theta K$. Let $\boldsymbol{\epsilon}^{k}=\bar{\boldsymbol{\epsilon}}^{k}$ for $k\in\bar{\mathcal{S}}$ and $\boldsymbol{\epsilon}^{k}=\hat{\mathbf{x}}^{k}-\mathbf{x}^{k}$ for $k\in\mathcal{S}\setminus\bar{\mathcal{S}}$ where $\mathbf{x}^{k}$ is arbitrarily chosen from $\{\mathbf{x}\in\mathcal{X}\ |\ {\mathbf{a}^{i}}^{\prime}\mathbf{x}=b_{i},\forall i\in\bar{\mathcal{A}}\}$. We will first show that the solution $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\boldsymbol{\epsilon}^{k}\}_{k\in\mathcal{K}})$ constructed as above is feasible for problem (6). For all $i\in\bar{\mathcal{A}}=\{i\,|\,\bar{v}_{i}=1\}$, we have ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})={\mathbf{a}^{i}}^{\prime}\mathbf{x}^{k}=b_{i}$ for each $k\in\mathcal{K}$. Also, for all $i\in\mathcal{I}\setminus\bar{\mathcal{A}}$ (i.e., $\{i\,|\,\bar{v}_{i}=0\}$), we have ${\mathbf{a}^{i}}^{\prime}(\hat{\mathbf{x}}^{k}-\boldsymbol{\epsilon}^{k})={\mathbf{a}^{i}}^{\prime}\mathbf{x}^{k}\geq b_{i}$ for each $k\in\mathcal{K}$. Thus, $\bar{\mathbf{v}}$ and $\{\boldsymbol{\epsilon}^{k}\}_{k\in\mathcal{K}}$ satisfy (6b). Constraint (6c) is satisfied because $\|\boldsymbol{\epsilon}^{k}\|_{\ell}=\|\bar{\boldsymbol{\epsilon}}^{k}\|_{\ell}\leq\tau$ for all $k\in\bar{\mathcal{S}}=\{k\,|\,\bar{u}_{k}=1\}$. Constraint (6d) is satisfied as we already know $\sum_{k=1}^{K}\bar{u}_{k}\geq\theta K$. As a result, $(\bar{\mathbf{v}},\bar{\mathbf{u}},\{\boldsymbol{\epsilon}^{k}\}_{k\in\mathcal{K}})$ is feasible for (6). Then, by Theorem 3.4, there must also exist a solution $(\bar{\mathbf{c}},\{{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ where $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i:\bar{v}_{i}=1})=\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$, as desired. \Halmos

Proof 9.10

Proof of Proposition 4.6: Let $P(r,k)$ denote the optimization problem solved in Line 5 of Algorithm 2 with respect to $\hat{\mathbf{x}}^{k}$, and $(\mathbf{v}^{(r,k)},\boldsymbol{\epsilon}^{(r,k)})\in\underset{(\mathbf{v},\boldsymbol{\epsilon})}{\text{argmax}}\ P(r,k)$. We first show that there exists $\bar{r}\leq\displaystyle\sum_{z=z^{*}}^{n}|\mathcal{V}_{z}|$ such that $\mathbf{v}^{(\bar{r},1)}=\mathbf{v}^{*}$ for problem $P(\bar{r},1)$ (i.e., with respect to the first data point $\hat{\mathbf{x}}^{1}$). For each iteration $r\geq 2$, each cut only prevents the solution found in the previous iteration $r-1$, i.e., $\mathbf{v}^{(r-1,1)}$, from being feasible for constraints (8b)–(8d). This means that $P(r,1)$ can be solved $\displaystyle\sum_{z=1}^{n}|\mathcal{V}_{z}|$ times with iteratively added cuts to collect all solutions satisfying (8b)–(8d) with respect to $\hat{\mathbf{x}}^{1}$. Note that because $(\mathbf{v}^{*},\boldsymbol{\epsilon}^{k*})$ satisfies (6b) for all $k\in\mathcal{K}$ and $u^{*}_{1}=1$, $(\mathbf{v}^{*},\boldsymbol{\epsilon}^{1*})$ satisfies constraints (8b)–(8d) with respect to $\hat{\mathbf{x}}^{1}$. Thus, it is guaranteed that for some $\bar{r}\leq\displaystyle\sum_{z=1}^{n}|\mathcal{V}_{z}|$ we have $\mathbf{v}^{(\bar{r},1)}=\mathbf{v}^{*}$. Furthermore, because the value function of $P(r,1)$ is non-increasing in $r$, for any $r>\displaystyle\sum_{z=z^{*}}^{n}|\mathcal{V}_{z}|$ the optimal value of $P(r,1)$ is less than $z^{*}$ and thus $\mathbf{v}^{*}$ is not achievable. Therefore, $\bar{r}$ must be no greater than $\displaystyle\sum_{z=z^{*}}^{n}|\mathcal{V}_{z}|$.

*Next, we need to show that when $\mathbf{v}^{(\bar{r},1)}=\mathbf{v}^{*}$ at iteration $\bar{r}$, the algorithm achieves $\mathbf{v}^{Alg}=\mathbf{v}^{*}$ and $z^{Alg}=z^{*}$. To do so, we first show that when $\mathbf{v}^{(\bar{r},1)}=\mathbf{v}^{*}$, Lines 4–9 of the algorithm lead to $\bar{\mathbf{D}}_{\bar{r}}$ with an all-one submatrix with $z^{*}$ columns and at least $\lceil\theta K\rceil$ rows. Let $\mathcal{I}^{*}=\{i\in\mathcal{I}\,|\,v^{*}_{i}=1\}$ and $\mathcal{S}^{*}=\{k\in\mathcal{K}\,|\,u^{*}_{k}=1\}$. For some given $r$ and $k$, $P(r,k)$ has an optimal solution with $v_{i}^{(r,k)}=0$ if $w^{k}_{i}=0$; also, based on how the weight parameters $\mathbf{w}^{k}$ are updated (see Section 4.3), if $v_{i}^{(r,k)}=0$ for $k=1,\ldots,k^{\prime}-1$, we have $w^{k^{\prime}}_{i}=0$. Generalizing this observation, if $\mathbf{v}^{(\bar{r},1)}=\mathbf{v}^{*}$ then for all $k=2,\ldots,K$ we have $w^{k}_{i}=0$ for all $i\in\mathcal{I}\setminus\mathcal{I}^{*}$ and $w^{k}_{i}>0$ for all $i\in\mathcal{I}^{*}$, and thus the maximum of the objective function of $P(\bar{r},k)$ for each $k$, i.e., ${\mathbf{w}^{k}}^{\prime}\mathbf{v}$, can be obtained when $\mathbf{v}=\mathbf{v}^{*}$. Because $(\mathbf{v}^{*},\boldsymbol{\epsilon}^{k*})$ satisfies all constraints of problem (8) as well as all cuts in $P(\bar{r},k)$ generated up to iteration $\bar{r}$ for all $k\in\mathcal{S}^{*}$, we have $(\mathbf{v}^{*},\boldsymbol{\epsilon}^{k*})\in\underset{\mathbf{v},\boldsymbol{\epsilon}}{\text{argmax}}\ P(\bar{r},k)$. This means if ${\mathbf{v}^{*}}^{\prime}$ is assigned to the first row of $\bar{\mathbf{D}}_{\bar{r}}$, it is guaranteed that ${\mathbf{v}^{*}}^{\prime}$ will also be assigned to rows $k\in\mathcal{S}^{*}$ of $\bar{\mathbf{D}}_{\bar{r}}$. Clearly, the submatrix of this $\bar{\mathbf{D}}_{\bar{r}}$ with columns and rows corresponding to $i\in\mathcal{I}^{*}$ and $k\in\mathcal{S}^{*}$, respectively, is all-one and we have $|\mathcal{I}^{*}|=z^{*}$ and $|\mathcal{S}^{*}|\geq\theta K$. The clique problem in Line 10 of the algorithm (also see formulation (14)) then takes this $\bar{\mathbf{D}}_{\bar{r}}$ as input to return its optimal solution $\bar{\mathbf{v}}^{\bar{r}}=\mathbf{v}^{*}$ with ${\bar{\mathbf{v}}^{r^{\prime}}}\mathbf{e}=z^{*}$. Since $z^{*}$ is the maximum achievable value (it is the optimal value of the exact MQIO problem), it satisfies the condition in Line 12 of Algorithm 2 and thus we have $\mathbf{v}^{Alg}=\bar{\mathbf{v}}^{\bar{r}}=\mathbf{v}^{*}$ and $z^{Alg}=z^{*}$. *

*Finally, we show that the stopping criterion is valid, i.e., the while loop condition is not violated before $\mathbf{v}^{Alg}=\mathbf{v}^{*}$ is observed. Consider a certain iteration $\bar{r}\leq\displaystyle\sum_{z=z^{*}}^{n}|\mathcal{V}_{z}|$ and suppose we have $\mathbf{v}^{Alg}\neq\mathbf{v}^{*}$ and $z^{Alg}<z^{*}$ at the end of iteration $\bar{r}$ (i.e., Lines 13–14 of Algorithm 2). For the while loop condition to be violated in the next iteration, i.e., iteration $\bar{r}+1$, we must have $|\mathcal{G}_{\bar{r}+1}|<z^{Alg}<z^{*},$ which means ${\mathbf{e}}^{\prime}\mathbf{v}^{(\bar{r},1)}<z^{*}$. However, this cannot happen because we showed earlier in this proof that for all $r\leq\displaystyle\sum_{z=z^{*}}^{n}|\mathcal{V}_{z}|$ we have ${\mathbf{e}}^{\prime}\mathbf{v}^{(r,1)}\geq z^{*}$.  \Halmos

Proof 9.11

Proof of Proposition 7.1: Proposition 3.2 implies that given a solution $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ feasible for (10), $(\tilde{\mathbf{c}}=\mathbf{a}^{\tilde{i}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\tilde{\mathbf{y}}=\mathbf{e}_{\tilde{i}},\bar{\mathcal{S}})$ is also feasible for any $\tilde{i}\in\bar{\mathcal{I}}=\{i\in\mathcal{I}\ \arrowvert\ \bar{y}_{i}>0\}$. For an arbitrary $\tilde{i}\in\bar{\mathcal{I}}$, consider the following formulation defined for each $k\in\mathcal{K}$:

[TABLE]

Let $\tilde{\boldsymbol{\epsilon}}^{k}$ be the optimal solution to (23) for each $k\in\mathcal{K}$ (here we suppress the index $\tilde{i}$ for brevity). Note that a solution constructed as $(\tilde{\mathbf{c}}=\mathbf{a}^{\tilde{i}},\{\tilde{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\tilde{\mathbf{y}}=\mathbf{e}_{\tilde{i}},\bar{\mathcal{S}})$ is also feasible for (10) because $\lVert\tilde{\boldsymbol{\epsilon}}^{k}\rVert_{\ell}\leq\lVert\bar{\boldsymbol{\epsilon}}^{k}\rVert_{\ell}\leq\tau$ for each $k\in\bar{\mathcal{S}}$. The objective of (10), i.e., minimizing $\tau$, is equivalent to minimizing $\underset{k\in\bar{\mathcal{S}}}{\max}\{\lVert\boldsymbol{\epsilon}^{k}\rVert_{\ell}\}$. Because $||\tilde{\boldsymbol{\epsilon}}^{k}||_{\ell}\leq||\bar{\boldsymbol{\epsilon}}^{k}||_{\ell}\leq\tau$ for each $k\in\bar{\mathcal{S}}$, we have $\underset{k\in\bar{\mathcal{S}}}{\max}\{\lVert\tilde{\boldsymbol{\epsilon}}^{k}\rVert_{\ell}\}\leq\underset{k\in\bar{\mathcal{S}}}{\max}\{\lVert\bar{\boldsymbol{\epsilon}}^{k}\rVert_{\ell}\}$. That is, for any given feasible solution $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ for (10) where $\bar{\mathbf{c}}$ is not identical to $\mathbf{a}^{i}$ for any $i\in\mathcal{I}$, we can always construct another feasible solution $(\tilde{\mathbf{c}}=\mathbf{a}^{\tilde{i}},\{\tilde{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\tilde{\mathbf{y}}=\mathbf{e}_{\tilde{i}},\bar{\mathcal{S}})$ without increasing the objective value of (10), which completes the proof. \Halmos

Proof 9.12

Proof of Lemma 7.2: Let $\hat{\mathcal{I}}=\{i\in\mathcal{I}\,|\,\mathbf{a}^{i}\in\hat{\mathcal{C}}\}$. Given shifted data $\tilde{\mathcal{X}}$, assume $\mathbf{c}=\mathbf{a}^{i}$ is not inverse-feasible for QIO($\tilde{\mathcal{X}},\tau,\theta$) for any $i\in\hat{\mathcal{I}}$. We want to show that there is no $\mathbf{c}\in\hat{\mathcal{C}}$ that is inverse-feasible for QIO($\tilde{\mathcal{X}},\tau,\theta$). Suppose to the contrary that there is $\bar{\mathbf{c}}\in\hat{\mathcal{C}}$ feasible for QIO($\tilde{\mathcal{X}},\tau,\theta$) and let $\bar{\mathcal{I}}$ be such that $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{I}}})$. Then, from Proposition 3.2, there must be some $\bar{i}\in\bar{\mathcal{I}}$ such that $\mathbf{c}=\mathbf{a}^{\bar{i}}$ is also feasible for QIO($\tilde{\mathcal{X}},\tau,\theta$). Furthermore, since $\bar{\mathbf{c}}\in\hat{\mathcal{C}}$, this means that $\mathbf{c}=\mathbf{a}^{\bar{i}}$ must be also feasible for QIO($\hat{\mathcal{X}},\tau,\theta$) (i.e., it must be that $\bar{i}\in\hat{\mathcal{I}}$), which is a contradiction.\Halmos

Proof 9.13

Proof of Proposition 7.3: Since $\bar{\mathbf{c}}$ is a strict conic combination of the selected $\mathbf{a}^{i}$ vectors, we have $\mathcal{X}^{*}(\bar{\mathbf{c}})=\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}$. For each $k\in\bar{\mathcal{S}}$, let $\tilde{\mathbf{x}}^{k}\in\underset{\mathbf{x}\in\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}}{\text{argmin}}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}$. Then, for all $k\in\bar{\mathcal{S}}$ and $i\in\bar{\mathcal{I}}$, we have $\displaystyle\max_{\mathbf{x}\in\mathcal{X}_{i}}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\|\tilde{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}+\displaystyle\max_{\mathbf{x}\in\mathcal{X}_{i}}\{\|\tilde{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\|\tilde{\mathbf{x}}^{k}-\hat{\mathbf{x}}^{k}\|_{\ell}+\rho_{i}$. Hence, for each $k\in\bar{\mathcal{S}}$,

[TABLE]

where the first inequality holds because we have $\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}\subseteq\mathcal{X}_{i}$ for all $i\in\bar{\mathcal{I}}$, which leads to $\displaystyle\max_{\mathbf{x}\in\underset{i\in\bar{\mathcal{I}}}{\cap}\mathcal{X}_{i}}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\}\leq\max_{\mathbf{x}\in\mathcal{X}_{i}}\{\|\hat{\mathbf{x}}^{k}-\mathbf{x}\|_{\ell}\},\forall i\in\bar{\mathcal{I}}$.\Halmos

Proof 9.14

Proof of Lemma 8.1: Because $\bar{\alpha}^{k}_{i}=0$ for all $i\in\bar{\mathcal{I}}$ and $(\bar{\boldsymbol{\alpha}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ is feasible for (15), we have ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})=b_{i}$ for all $i\in\bar{\mathcal{I}}$, ${\mathbf{a}^{i}}^{\prime}(\mathbf{\hat{x}}^{k}-\bar{\boldsymbol{\epsilon}}^{k})\geq b_{i}$ for all $i\in\mathcal{I}\setminus\bar{\mathcal{I}}$, and $\rVert\bar{\boldsymbol{\epsilon}}^{k}\lVert_{\ell}\leq\tau$. Let $\bar{v}_{i}^{k}=1-\bar{\alpha}_{i}^{k}$ for all $i\in\bar{\mathcal{I}}$ and $\bar{v}_{i}^{k}=0$ otherwise. Clearly, $(\bar{\mathbf{v}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ satisfies all constraints of (8).\Halmos

Proof 9.15

Proof of Proposition 8.2: Assume that $\tilde{\mathbf{D}}$ satisfies (16) and has an all-one submatrix whose rows and columns correspond to $\bar{\mathcal{S}}$ and $\bar{\mathcal{A}}$, respectively, where $\bar{\mathcal{S}}\in\mathcal{K},|\bar{\mathcal{S}}|\geq\theta K$, and $\bar{\mathcal{A}}\in\mathcal{I}$. By definition, for each $k\in\mathcal{K}$ there exists a feasible solution $(\bar{\boldsymbol{\alpha}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ for (15) with respect to $\hat{\mathbf{x}}^{k}$ where $\bar{\alpha}_{i}^{k}=0$ for all $i\in\bar{\mathcal{A}}$ and $k\in\bar{\mathcal{S}}$, and $\bar{\alpha}_{i}^{k}>0$ otherwise. This means that, by Lemma 8.1, for each $k\in\mathcal{K}$ there exists a corresponding feasible solution $(\bar{\mathbf{v}}^{k},\bar{\boldsymbol{\epsilon}}^{k})$ for (8) with respect to $\hat{\mathbf{x}}^{k}$ where $\bar{v}_{i}^{k}=1$ for all $i\in\bar{\mathcal{A}}$ and $k\in\bar{\mathcal{S}}$, and $\bar{v}_{i}^{k}=0$ otherwise. Construct a matrix $\bar{\mathbf{D}}\in\{0,1\}^{K\times m}$ such that $\bar{\mathbf{D}}_{ki}=1$ if $\bar{v}_{i}^{k}=1$ and $\bar{\mathbf{D}}_{ki}=0$ otherwise. Note that $\bar{\mathbf{D}}$ satisfies (9) and has an all-one submatrix whose rows and columns correspond to $\bar{\mathcal{S}}$ and $\bar{\mathcal{A}}$, respectively, where $\bar{\mathcal{S}}\subseteq\mathcal{K}$, $|\bar{\mathcal{S}}|\geq\theta K$, and $\bar{\mathcal{A}}\subseteq\mathcal{I}$. Thus, by Proposition 4.4, there exists a solution $(\bar{\mathbf{c}},\{\bar{\boldsymbol{\epsilon}}^{k}\}_{k\in\mathcal{K}},\bar{\mathbf{y}},\bar{\mathcal{S}})$ feasible for $\mathbf{QIO}(\mathcal{K},\tau,\theta)$ where $\bar{\mathbf{c}}\in\textup{cone}(\{\mathbf{a}^{i}\}_{i\in\bar{\mathcal{A}}})$.  \Halmos

10 Supplemental Materials for the Numerical Results

10.1 The Diet Problem: Data

Food items and nutrient data per serving.* *

Food Type Lower Upper

1 2 3 4 5 6 7 8 9 limit limit

Energy (KCAL) 91.53 68.94 23.51 65.49 110.88 83.28 80.50 63.20 52.16 1800.00 2500.00

Total_Fat (g) 4.95 0.71 1.80 3.48 6.84 4.41 5.80 0.94 0.18 44.00 78.00

Carbohydrate (g) 6.89 12.16 0.25 0.00 5.44 4.68 0.56 11.42 13.59 220.00 330.00

Protein (g) 4.90 3.68 1.59 7.99 6.80 5.93 6.27 2.40 0.41 56.00 NA

Fiber (g) 0.00 0.06 0.00 0.00 0.28 0.29 0.00 1.19 1.81 20.00 30.00

Vitamin C (mg) 0.01 1.76 0.00 0.00 0.17 0.16 0.00 0.02 11.19 90.00 2000.00

Vitamin B6 (mg) 0.06 0.03 0.01 0.09 0.11 0.06 0.06 0.03 0.10 1.30 100.00

Vitamin B12 (mcg) 0.67 0.39 0.09 0.65 0.11 0.63 0.56 0.00 0.00 2.40 NA

Calcium (mg) 172.09 125.72 46.24 2.21 5.90 15.03 29.00 27.21 6.14 1000.00 2500.00

Iron (mg) 0.05 0.08 0.04 0.75 0.35 0.35 0.73 0.79 0.13 8.00 45.00

Copper (mg) 0.02 0.04 0.01 0.03 0.03 0.03 0.04 0.05 0.06 0.90 10.00

Sodium (mg) 61.02 48.24 65.08 72.32 211.05 128.27 223.50 125.62 1.42 1500.00 2300.00

Vitamin A (mcg) 42.89 22.24 12.78 0.00 1.33 9.53 81.00 0.01 13.04 900.00 3000.00

Max serving 4 8 4 5 5 5 4 4 8

**

10.2 The Transshipment Problem: Problem Description and Additional Numerical Results

Figure 7 and the following formulation show the transshipment problem we consider.

[TABLE]

where $\mathcal{N}_{s}$, $\mathcal{N}_{d}$, and $\mathcal{N}_{t}$ denote the set of supply, demand, and transshipment nodes, respectively, and $\mathcal{N}$ denotes the set of all nodes. Variables $\mathbf{x}^{(p)}$ and $\mathbf{x}^{(t)}$ denote the production level at each supply node and transshipment flow on each arc, respectively. Parameters $\boldsymbol{\pi}^{(t)}$ and $\boldsymbol{\pi}^{(p)}$ represent the arc capacity and production capacity, respectively.

In addition to the numerical results presented in Section 5.4, we also compare the objective function values achieved by our forward solution $\bar{\mathbf{x}}^{t}$ (recall that $\mathbf{x}=[{\mathbf{x}^{(p)}};{\mathbf{x}^{(t)}}]$ and $\mathbf{c}=[{\mathbf{c}^{(p)}};{\mathbf{c}^{(t)}}]$) and that from the given data. We adopt the metrics from the previous study on online inverse LP (Bärmann et al. 2017); this study assumes a single noiseless (i.e., optimal) decision at each time point and shows that both $1/T\sum_{t=1}^{T}{\mathbf{c}^{true}}^{\prime}(\bar{\mathbf{x}}^{t}-\hat{\mathbf{x}}^{t})$ and $1/T\sum_{t=1}^{T}{\mathbf{c}^{t}}^{\prime}(\bar{\mathbf{x}}^{t}-\hat{\mathbf{x}}^{t})$ converge to zero as $T\rightarrow\infty$. In our case, because we consider batches of data at each time we use the following modified metrics: 1/T\big{[}\sum_{t=1}^{T}\sum_{k=1}^{K_{t}}{\mathbf{c}^{true}}^{\prime}(\bar{\mathbf{x}}^{t}-\hat{\mathbf{x}}^{k})/K_{t}\big{]} and 1/T\big{[}\sum_{t=1}^{T}\sum_{k=1}^{K_{t}}{\mathbf{c}^{t}}^{\prime}(\bar{\mathbf{x}}^{t}-\hat{\mathbf{x}}^{k})/K_{t}\big{]}. Note that $\mathbf{c}^{t}$ is a cost vector randomly selected from the inverse set at time $t$. Figure 8 shows these metrics achieved by our algorithm with noisy datasets at each iteration. Both metrics averaged over the data points also converge; however, neither of them reaches 0 because of the noise in the data.