Trade-off preservation in inverse multi-objective convex optimization

TL;DR

This paper introduces a novel inverse optimization approach for multi-objective convex problems that preserves decision maker's trade-offs, even with non-Pareto optimal inputs, and demonstrates its effectiveness on clinical prostate cancer data.

Contribution

It develops a new inverse optimization method that accounts for non-Pareto solutions and preserves trade-offs, extending existing models with a linear approximation and efficient algorithm.

Findings

Successfully applied to prostate cancer radiation therapy data

Balances trade-off preservation with computational efficiency

Encompasses many existing inverse models in the literature

Abstract

We present a new inverse optimization methodology for multi-objective convex optimization that accommodates an input solution that may not be Pareto optimal and determines a weight vector that produces a Pareto optimal solution that approximates the input solution and preserves the decision maker's intention encoded in it. We introduce a notion of trade-off preservation, which we use as a measure of similarity for approximating the input solution, and show its connection with minimizing an optimality gap. Our inverse model maintains the complexity of the traditional inverse convex models. We propose a linear approximation to the model and a successive linear programming algorithm that balance between trade-off preservation and computational efficiency, and show that our model encompasses many of the existing models from the literature. We demonstrate the proposed method using clinical data from prostate cancer radiation therapy.

Tables

Table 1. Table 1: Comparing the model of  Keshavarz et al. ( 2011 ) with different weights fixed to one for patient #1.

	KES1			KES2			KES3			KES4			KES5
OAR	${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$		${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$		${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$		${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$		${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$
Blad	0.973	0.653		0.003	0.837		0.000	0.782		0.000	0.811		0.001	0.823
Rect	0.004	1.000		0.865	0.823		0.000	0.949		0.000	0.963		0.002	0.972
LFem	0.000	1.143		0.005	0.232		0.992	0.000		0.004	0.001		0.000	0.023
RFem	0.000	0.671		0.004	0.177		0.004	0.002		0.992	0.000		0.000	0.050
Ring	0.023	0.929		0.123	0.614		0.004	0.412		0.004	0.411		0.997	0.404
KES# refers to KES with #-th weight fixed to one. Weights were renormalized so that they add up to one.

Table 2. Table 2: Comparison of the results from IOP r , LIOP r , and SLP algorithm for three patients.

		IOPr					LIOPr					SLP					FOP
Pat	OAR	${ϵ}^{\ast }$	${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$	Time		${ϵ}^{\ast }$	${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$	Time		${ϵ}^{\ast }$	${𝜶}^{\ast }$	$\frac{f_k\left({𝐱}^{\ast }\right)}{f_k\left(\widehat{𝐱}\right)}$	Time		Time
					(s)					(s)					(s)		(s)
1	Blad	0.812	0.014	0.812	1582		0.769	0.014	0.813	167		0.813	0.014	0.815	662 $[$30$]$		371
	Rect		0.935	0.812				0.936	0.812				0.935	0.812
	LFem		0.002	0.812				0.001	0.905				0.002	0.746
	RFem		0.000	0.809				0.001	0.695				0.000	0.932
	Ring		0.049	0.812				0.048	0.816				0.049	0.813
2	Blad	0.803	0.932	0.803	1384		0.748	0.937	0.802	82		0.803	0.933	0.803	1000 $[$56$]$		258
	Rect		0.002	0.803				0.008	0.754				0.002	0.802
	LFem		0.012	0.803				0.012	0.771				0.011	0.846
	RFem		0.018	0.803				0.015	0.965				0.018	0.808
	Ring		0.036	0.803				0.028	0.825				0.036	0.803
3	Blad	0.712	0.471	0.712	764		0.617	0.456	0.713	33		0.712	0.474	0.712	435 $[$30$]$		168
	Rect		0.412	0.712				0.438	0.709				0.409	0.712
	LFem		0.008	0.712				0.012	0.525				0.007	0.765
	RFem		0.033	0.712				0.027	0.889				0.033	0.689
	Ring		0.076	0.712				0.067	0.746				0.077	0.709
[$\cdot$] denotes the number of iterations.

Table 3. Table 3: Comparison of results averaged across 24 patients

	IOPr		SLP		LIOPr		KES
Var$\left(\frac{𝐟\left({𝐱}^{\ast }\right)}{𝐟\left(\widehat{𝐱}\right)}\right)$	0.004		0.009		0.076		0.938
$\left|{ϵ}^{\text{IOP}}-{ϵ}^{\ast }\right|$	0		0.001		0.050		N/A
${\Vert {𝜶}^{\text{IOP}}-{𝜶}^{\ast }\Vert }_2$	0		0.007		0.077		0.281
Solution time (s)	1,807		570		84		23
KES fixes the weight determined highest by the IOPr model to one.
${𝐱}^{\text{IOP}}$, ${ϵ}^{\text{IOP}}$ and ${𝜶}^{\text{IOP}}$ denote the optimal $𝐱$, $ϵ$ and $𝜶$ obtained by
the IOPr model, respectively.

Equations

$$\text{FOP}(\mathbf{\alpha }):\underset{\mathbf{x}}{\mathrm{minimize}}$$

$$g_l(\mathbf{x})\le 0,l=1,\dots ,L,$$

$$\mathbf{Ax}=\mathbf{b},$$

$$\underset{\mathbf{\alpha },\mathbf{\sigma },\mathbf{\pi }}{\mathrm{minimize}}$$

$$\sum _{k=1}^K{\alpha }_k{\nabla }_{\mathbf{x}}{f}_k(\hat{\mathbf{x}})+\sum _{l=1}^L{\sigma }_l{\nabla }_{\mathbf{x}}{g}_l(\hat{\mathbf{x}})-{\mathbf{A}}^{\prime }\mathbf{\pi }=0,$$

$${\sigma }_l{g}_l(\hat{\mathbf{x}})=0,l=1,\dots ,L,$$

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

$$\underset{\mathbf{x}}{\mathrm{minimize}}$$

$$(x_1-2{)}^2+(x_2-2{)}^2-1\le 0,$$

$$(4{\alpha }_1+{\alpha }_2){\hat{x}}_1+({\hat{x}}_1-2)\sigma =0,$$

$$({\alpha }_1+4{\alpha }_2){\hat{x}}_2+({\hat{x}}_2-2)\sigma =0,$$

$$(({\hat{x}}_1-2{)}^2+({\hat{x}}_2-2{)}^2-1)\sigma =0,$$

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

$$(4-\sqrt{2})(4{\alpha }_1+{\alpha }_2)-\sqrt{2}\sigma =0,$$

$$(4-\sqrt{2})({\alpha }_1+4{\alpha }_2)-\sqrt{2}\sigma =0,$$

$$\mathbf{\alpha }\ge 0,\sigma \ge 0,$$

$$1.7(4{\alpha }_1+{\alpha }_2)-0.3\sigma =0,$$

$$1.3({\alpha }_1+4{\alpha }_2)-0.7\sigma =0,$$

$$\mathbf{\alpha }\ge 0,\sigma =0,$$

$$(4{\alpha }_1+{\alpha }_2)-\sigma =0,$$

$$({\alpha }_1+4{\alpha }_2)-\sigma =0,$$

$$\mathbf{\alpha }\ge 0,\sigma =0,$$

$$\text{KES}(\hat{\mathbf{x}}):\underset{\mathbf{\alpha },\mathbf{\sigma },\mathbf{\pi },\mathbf{\delta },\mathbf{\gamma },\mathbf{\rho }}{\mathrm{minimize}}$$

$$\sum _{k=1}^K{\alpha }_k{\nabla }_{\mathbf{x}}{f}_k(\hat{\mathbf{x}})+\sum _{l=1}^L{\sigma }_l{\nabla }_{\mathbf{x}}{g}_l(\hat{\mathbf{x}})-{\mathbf{A}}^{\prime }\mathbf{\pi }=\mathbf{\delta },$$

$${\sigma }_l{g}_l(\hat{\mathbf{x}})={\gamma }_l,l=1,\dots ,L,$$

$${\pi }_j({\mathbf{a}}_j^{\prime }\hat{\mathbf{x}}-b_j)={\rho }_j,j=1,\dots ,m,$$

$${\alpha }_1=1,$$

$$\mathbf{\alpha }\ge 0,\mathbf{\sigma }\ge 0,$$

$$\underset{\mathbf{\alpha },\sigma ,\mathbf{\delta },\gamma }{\mathrm{minimize}}$$

$$1.7(4{\alpha }_1+{\alpha }_2)-0.3\sigma ={\delta }_1,$$

$$1.3({\alpha }_1+4{\alpha }_2)-0.7\sigma ={\delta }_2,$$

$$-0.42\sigma =\gamma ,$$

$${\alpha }_1=1,$$

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

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

$${\mathbf{A}}^{\prime }\mathbf{p}\le {\mathbf{C}}^{\prime }\mathbf{\alpha },$$

$${\mathbf{\alpha }}^{\prime }\mathbf{C}\hat{\mathbf{x}}=ϵ{\mathbf{b}}^{\prime }\mathbf{p},$$

$$\Vert \mathbf{\alpha }{\Vert }_1=1,$$

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

$$u_{k_{1}k_{2}}\Big{(}f_{k_{1}}(\mathbf{x}^{*})-f_{k_{1}}(\hat{\mathbf{x}})\Big{)}=\Big{(}f_{k_{2}}(\mathbf{x}^{*})-f_{k_{2}}(\hat{\mathbf{x}})\Big{)},\quad\forall k_{1},k_{2}\in\{1,\ldots,K\}.$$

$$\frac{f_{k_1}({\mathbf{x}}^{\ast })}{f_{k_1}(\hat{\mathbf{x}})}=\frac{f_{k_2}({\mathbf{x}}^{\ast })}{f_{k_2}(\hat{\mathbf{x}})},\forall k_1,k_2\in \{1,\dots ,K\}.$$

$$f_{k_1}({\mathbf{x}}^{\ast })-f_{k_1}(\hat{\mathbf{x}})=f_{k_2}({\mathbf{x}}^{\ast })-f_{k_2}(\hat{\mathbf{x}}),\forall k_1,k_2\in \{1,\dots ,K\}.$$

$$\text{IOP}(\hat{\mathbf{x}}):\underset{ϵ,\mathbf{x}}{\mathrm{minimize}}$$

$${\mu }_kϵ\ge f_k(\mathbf{x})-f_k(\hat{\mathbf{x}}),k=1,\dots ,K,$$

$$g_l(\mathbf{x})\le 0,l=1,\dots ,L,$$

$$\mathbf{Ax}=\mathbf{b}.$$

Full text

1 Introduction

Given a feasible solution to an (forward) optimization problem, the inverse optimization problem aims to determine parameter values – typically objective function parameters – that make the given solution optimal. Classical inverse approaches leverage duality to derive tractable inverse problems that retain the complexity of the underlying forward problem (Ahuja and Orlin 2001, Iyengar and Kang 2005). However, these inverse models return a trivial solution (e.g., a coefficient vector of all zeros for inverse linear optimization) if the given feasible solution is not a candidate to be optimal for the forward problem.

In general, there is no guarantee that a given solution is exactly optimal for the assumed forward problem. An implemented solution may have been adjusted after optimization, it may have been derived heuristically and is only approximately optimal, or it may simply be a noisy observation of a pristine solution. Another possibility is that the assumed forward model is itself a simplification of a complex system for which an observed solution is near optimal. When inverse problems involve such noise and uncertainty around the model and data, it is important to determine an objective function that replicates as closely as possible the intentions of the decision maker who generated the given solution. Recent studies have generalized the classical inverse models to overcome the issue of given solutions not being exactly optimal. For example, Troutt et al. (2006), Keshavarz et al. (2011), Chow and Recker (2012), Chan et al. (2014), Bertsimas et al. (2015), Aswani et al. (2017), and Chan et al. (2017) developed approximate inverse optimization models that impute model parameters that make the observed solutions minimally suboptimal. In this paper, we bring together the ideas of Keshavarz et al. (2011) and Chan et al. (2014) and develop a new inverse optimization model for multi-objective convex optimization where the given solutions are not Pareto optimal.

In multi-objective optimization, decision making is typically based in the objective space, i.e., the space where the vector of objective values resides. Deciding between different solutions in this space involves examining the trade-off in objective values between different points on the Pareto frontier. With a weighted objective function, as is common in convex multi-objective problems, these points are generated by solving the forward problem with different weight vectors, which explicitly quantify the trade-off in the objectives deemed acceptable by the decision maker. Conversely, without access to the weights, and only observing a solution on the Pareto frontier, it is possible using classical inverse optimization methods to reverse engineer the weight vector that generated the solution and therefore determine the decision maker’s intention with respect to trade-offs. However, if the observed solution is not on the Pareto frontier, approximate inverse optimization models like the ones mentioned above can be applied. Given a non-Pareto solution, such a model will return a weight vector that generates a Pareto optimal solution. We propose that an appropriate inverse optimization model should return a weight vector and corresponding Pareto optimal solution that differs from the input (non-Pareto) solution in a way that is consistent with the trade-offs the original decision maker had in mind. Mathematically, this means the Pareto optimal solution should have adjusted objective values that are “consistent” across all the objectives – a concept we refer to as trade-off preservation. In this paper, we will provide a formal definition of trade-off preservation, use it as a measure of similarity between different solutions, and develop an inverse model that preserves the trade-off encoded in the given solution.

Our work generalizes the approach of Chan et al. (2014) by considering convex multi-objective optimization problems. Also, our notion of trade-off preservation is general enough to represent various ways to characterize trade-offs across multiple objectives. For example, our model can be specialized to the duality gap minimization approaches in Chan et al. (2014) and the single-objective approximate inverse convex model in Keshavarz et al. (2011). We show that existing inverse optimization models designed for single-objective optimization which could potentially be used for multi-objective optimization, e.g., Keshavarz et al. (2011), may not take into consideration trade-offs across multiple objectives encoded in the given solution. As in Keshavarz et al. (2011) and Chan et al. (2014), we assume that a set of objectives is pre-specified. Our contributions are:

• (1)

We generalize previous inverse optimization approaches and develop a new inverse convex multi-objective optimization model that is itself a convex problem, accommodates any input solution, and determines a nonzero weight vector that preserves the trade-off encoded in the input solution. We introduce a notion of trade-off preservation that is generally applicable to multi-objective optimization and prove that some special cases of trade-off preservation are equivalent to the concept of duality gap minimization in inverse optimization.

• (2)

We propose an efficient linear approximation of the proposed inverse problem as well as a successive linear programming algorithm that bridges the exact and approximate methods. Based on the linear approximation, we propose a general inverse convex multi-objective optimization framework that encompasses many of the inverse models from the literature.

• (3)

We demonstrate the application of our inverse optimization model to a clinical treatment planning problem in prostate cancer radiation therapy. Using a clinical dataset, we show that weights that preserve the trade-off encoded by the given objective values produce treatments that maintain the clinical quality of the original treatments across all relevant metrics. We show that inverse models that are not trade-off-preserving may lead to treatments that deviate substantially from the original treatments and violate clinical acceptability criteria.

2 Background

We first define a canonical multi-objective convex optimization problem as the forward problem. Then, we briefly review the inverse optimization models from Iyengar and Kang (2005), Keshavarz et al. (2011), and Chan et al. (2014).

2.1 Forward optimization problem

Let $f_{k}:\mathbb{R}^{n}\rightarrow\mathbb{R},k=1,\ldots,K$ and $g_{l}:\mathbb{R}^{n}\rightarrow\mathbb{R},l=1,\ldots,L$ be convex functions. Let $\mathbf{x}\in\mathbb{R}^{n},\mathbf{A}\in\mathbb{R}^{m\times n},$ and $\mathbf{b}\in\mathbb{R}^{m}$. We define the forward optimization problem (FOP) as

[TABLE]

where $\alpha_{k}$ is the weight for the $k$-th objective function. Let $\mathbf{X}$ be the feasible region of (1). We assume $\boldsymbol{\alpha}\in\mathbb{R}_{+}^{K}\backslash\{\mathbf{0}\}$, $f_{k}(\mathbf{x})>0,k=1,\ldots,K$ for $\mathbf{x}\in\mathbf{X}$, and $\mathbf{A}$ has full rank. We also assume that Slater’s condition holds (Boyd and Vandenberghe 2004). We define $\Omega(\boldsymbol{\alpha})$ to be the set of optimal solutions to FOP($\boldsymbol{\alpha}$) and $\Omega:=\bigcup_{\boldsymbol{\alpha}\in\mathbb{R}_{+}^{K}\backslash\{\mathbf{0}\}}\Omega(\boldsymbol{\alpha})$. A solution $\mathbf{x}\in\mathbf{X}$ is (weakly) Pareto optimal if there is no other $\mathbf{y}\in\mathbf{X}$ such that $f_{k}(\mathbf{y})<f_{k}(\mathbf{x}),$ for all $k=1,\ldots,K$. It is known that for a convex multi-objective optimization problem, the set $\Omega$ consists of all Pareto optimal solutions (Ehrgott 2005). For any $S\subseteq\mathbf{X}$, we write $\mathbf{f}(S)=\{(f_{1}(\mathbf{x}),\ldots,f_{K}(\mathbf{x}))\,|\,\mathbf{x}\in S\}$. We denote $\mathbf{f}(\mathbf{X})$ as the feasible region in the objective space and the set $\mathbf{f}(\Omega)$ as the Pareto set.

2.2 Inverse conic optimization by Iyengar and

Kang (2005)

We begin by illustrating the approach of Iyengar and Kang (2005) using our forward problem (1). Given $K$ prespecified objectives and a solution $\hat{\mathbf{x}}\in\mathbf{X}$, assumed to be a regular point (Bazaraa et al. (2006); pp.204–207), a weight vector that makes $\hat{\mathbf{x}}$ optimal can be found by solving the following problem:

[TABLE]

Constraints in problem (2) correspond to the KKT conditions for the forward problem (1) with Lagrange multipliers $\boldsymbol{\sigma}$ and $\boldsymbol{\pi}$. If $\boldsymbol{\alpha}^{*}$ is an optimal solution that arises from solving (2), then $\hat{\mathbf{x}}\in\Omega(\boldsymbol{\alpha}^{*})$. Iyengar and Kang (2005) used an objective function $||\boldsymbol{\alpha}-\hat{\boldsymbol{\alpha}}||$ where $\hat{\boldsymbol{\alpha}}$ is a given weight vector. However, we omit it to allow for an objective of minimizing “residuals” to be introduced. Note that an arbitrary $\hat{\mathbf{x}}$ need not be in $\Omega$, in which case formulation (2) returns $\boldsymbol{\alpha}^{*}=\mathbf{0}$ as the unique solution (see Example 2.1).

Example 2.1

Consider the following bi-objective convex optimization problem:

[TABLE]

where $f_{1}(x_{1},x_{2})=4x_{1}^{2}+x_{2}^{2}$ and $f_{2}(x_{1},x_{2})=x_{1}^{2}+4x_{2}^{2}$. Constraints of the corresponding inverse problem given an input solution $\hat{\mathbf{x}}$ are:

[TABLE]

Consider the following three points: $\hat{\mathbf{x}}^{a}=((4-\sqrt{2})/2,(4-\sqrt{2})/2)\in\Omega$, $\hat{\mathbf{x}}^{b}=(1.7,1.3)\in\mathbf{X}\setminus\Omega$, and $\hat{\mathbf{x}}^{c}=(1,1)\not\in\mathbf{X}$ (see Figure 1). For $\hat{\mathbf{x}}^{a}$, constraints (4) become:

[TABLE]

which hold when $\alpha_{1}=\alpha_{2}\geq 0$. For $\hat{\mathbf{x}}^{b}$, constraints (4) become:

[TABLE]

which are satisfied only when $\alpha_{1}=\alpha_{2}=0$. For $\hat{\mathbf{x}}^{c}$, constraints (4) become:

[TABLE]

which again are satisfied only when $\alpha_{1}=\alpha_{2}=0$. In all cases, $\boldsymbol{\alpha}=\mathbf{0}$ is a feasible inverse solution. When $\hat{\mathbf{x}}\not\in\Omega$ (e.g., $\hat{\mathbf{x}}^{b}$ and $\hat{\mathbf{x}}^{c}$), $\boldsymbol{\alpha}=\mathbf{0}$ is the only feasible solution since $\hat{\mathbf{x}}$ cannot be optimal for FOP($\boldsymbol{\alpha}$) for any nonzero $\boldsymbol{\alpha}$.

2.3 Convex objective function imputation by Keshavarz

et al. (2011)

Keshavarz et al. (2011) explicitly considered noise or modeling errors in imputing an objective function from a convex model so that the observed data from an unknown, complex system are as consistent as possible with the proposed model. Given $K$ prespecified objectives denoted as a vector-valued function $\mathbf{f}\in\mathbb{R}^{K}$ and a solution $\hat{\mathbf{x}}$, they relax the KKT conditions and allow residuals, which are then minimized by solving the following problem:

[TABLE]

where $\mathbf{a}_{j}^{\prime}$ denotes the $j$-th row of $\mathbf{A}$ and $\phi$ is a function of the residuals such that $\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})=0$ if and only if $\boldsymbol{\delta}=\mathbf{0}$, $\boldsymbol{\gamma}=\mathbf{0}$, and $\boldsymbol{\rho}=\mathbf{0}$ (e.g., $\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})=||\boldsymbol{\delta}||_{2}^{2}+||\boldsymbol{\gamma}||_{2}^{2}+||\boldsymbol{\rho}||_{2}^{2}$). Constraint (8e) ensures that $\boldsymbol{\alpha}=\mathbf{0}$ is not a feasible solution and serves to implicitly normalize the resulting weight vector. The original version of formulation (8) in Keshavarz et al. (2011) accommodated multiple input data points each with their own residuals, but we simply illustrate their method with a single data point $\hat{\mathbf{x}}$. The extension to multiple input data points is straightforward. The original model also excluded the $\boldsymbol{\rho}$ residual in (8d) as the focus was on $\hat{\mathbf{x}}\in\mathbf{X}$, but we include it here to emphasize the applicability of the model to points $\hat{\mathbf{x}}\not\in\mathbf{X}$. As noted in their paper, the choice of $\phi$ impacts the solution. However, how the weights are normalized can also have a large impact, as we demonstrate in the example below and our computational results.

Example 2.2

We revisit Example 2.1 and apply formulation (8) to $\hat{\mathbf{x}}^{b}$:

[TABLE]

The optimal weight vector is $(\alpha^{*}_{1},\alpha^{*}_{2})=(1,0)$. The corresponding FOP($\boldsymbol{\alpha}^{*}$) generates the optimal solution $\mathbf{x}^{*}=(1.067,1.641)$. Comparing $\mathbf{f}(\hat{\mathbf{x}}^{b})=(13.250,9.650)$ to $\mathbf{f}(\mathbf{x}^{*})=(7.244,11.910)$ we see one objective value decreases while the other increases. Furthermore, if the normalization constraint is changed to $\alpha_{2}=1$, then the optimal weight vector changes to $\boldsymbol{\alpha}^{*}=(0,1)$ and the new objective vector $\mathbf{f}(\mathbf{x}^{*})=(11.910,7.244)$, which is completely opposite to the previous case.

2.4 Generalized inverse multi-objective linear optimization by Chan

et al. (2014)

Chan et al. (2014) considered the multi-objective linear optimization problem $\min\{\boldsymbol{\alpha}^{\prime}\mathbf{C}\mathbf{x}\;|\;\mathbf{A}\mathbf{x}=\mathbf{b},\mathbf{x}\geq\mathbf{0}\}$ as the forward problem, where $\mathbf{C}$ is a matrix composed of linear objectives $\mathbf{c}^{\prime}_{k}$ (one per row $k$). In the corresponding inverse problem, the residual to be minimized was the duality gap:

[TABLE]

where $\mathbf{p}$ denotes the dual vector corresponding to the constraint $\mathbf{A}\mathbf{x}=\mathbf{b}$ of the forward program. It was shown that this $\epsilon$-strong duality approach produces a weight vector $\boldsymbol{\alpha}^{*}$ and a corresponding forward solution $\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$ such that $\mathbf{c}_{k}^{\prime}\mathbf{x}^{*}/\mathbf{c}_{k}^{\prime}\hat{\mathbf{x}}$ is equal for each objective $k$ for which $\alpha_{k}>0$. In the next section, we propose a new inverse model that generalizes Chan et al. (2014) to general convex multi-objective optimization problems by explicitly using the notion of trade-off preservation.

3 Models

In this section, we formalize the definition of trade-off preservation and develop an inverse convex multi-objective optimization model that determines a weight vector that preserves the trade-off encoded in a given solution. We then propose approximation approaches to speed up computation. Finally, we draw connections between our inverse model and other models from the literature.

3.1 Trade-off preservation

When a given solution is not on the Pareto frontier, we propose that an inverse multi-objective optimization model should find a weight vector that generates a new solution that is Pareto optimal and preserves the trade-off encoded in the given solution. With access only to the given solution, we can only infer the trade-off from its objective values. Thus, we characterize how a trade-off is preserved by examining the direction and length of the vector that connects the Pareto optimal solution and the given solution in objective space.

We first define general trade-off preservation in terms of the component-wise relation between two solutions in objective space.

Definition 3.1** (Generalized trade-off preservation)**

Given scalars $u_{k_{1}k_{2}}$ for every pair of objectives $k_{1}$ and $k_{2}$, a solution $\mathbf{x}^{*}$ preserves the trade-off encoded in $\hat{\mathbf{x}}$ if the following equalities hold:

[TABLE]

A weight vector $\boldsymbol{\alpha}^{*}$ preserves the trade-off encoded in $\hat{\mathbf{x}}$ if $\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$.

The quantity $u_{k_{1}k_{2}}$ is a problem-specific scaling factor that captures the relative difference between objectives $k_{1}$ and $k_{2}$. For example, if one unit of measurement for objective 1 is equivalent to 10 units for objective 2, then we set $u_{12}=10$ and $u_{21}=1/10$. Then, by equating the adjusted amount of perturbation across all the objectives, the perturbed solution $\mathbf{x}^{*}$ preserves the initial trade-offs encoded in the given solution $\hat{\mathbf{x}}$. Naturally, the scaling factor should have the following properties:

1. (i)

$u_{kk}=1,\quad\forall k\in\{1,\ldots,K\},$ 2. (ii)

$u_{k_{1}k_{2}}u_{k_{2}k_{1}}=1,\quad\forall k_{1},k_{2},\in\{1,\ldots,K\}$, and 3. (iii)

$u_{k_{1}k_{2}}u_{k_{2}k_{3}}=u_{k_{1}k_{3}},\quad\forall k_{1},k_{2},k_{3}\in\{1,\ldots,K\}$.

Next, we present two specific cases of trade-off preservation that are intuitive and have a direct connection with the concept of a duality gap. The first case considers two solutions whose objective values are component-wise proportional to be trade-off preserving, which is of practical use when the objectives are measured in different units whose values lie in ranges with substantial variation. We refer to this type of trade-off preservation as relative trade-off preservation.

Definition 3.2** (Relative trade-off preservation)**

A solution $\mathbf{x}^{*}$ preserves the relative trade-off encoded in $\hat{\mathbf{x}}$ if the following equalities hold for every pair of objectives $k_{1}$ and $k_{2}$:

[TABLE]

A weight vector $\boldsymbol{\alpha}^{*}$ preserves the relative trade-off encoded in $\hat{\mathbf{x}}$ if $\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$.

Definition 3.2 says that the objective values of $\mathbf{x}^{*}$ are adjusted component-wise by the same relative amount from the objective values of $\hat{\mathbf{x}}$. Geometrically, this definition means that $\mathbf{f}(\mathbf{x}^{*})$ lies on the line joining $\mathbf{f}(\hat{\mathbf{x}})$ and the origin. Since our focus is to find a Pareto optimal solution that is trade-off preserving, we are interested in identifying where that line intersects the Pareto set, if at all. This concept is also intimately connected with duality. In particular, in the next subsection we will show that the weight vector $\boldsymbol{\alpha}^{*}$ that satisfies the relative trade-off preservation leads to the minimum relative duality gap with respect to a given solution $\hat{\mathbf{x}}\in\mathbf{X}$. The relative trade-off preservation is in fact a special case of the general definition of trade-off preservation (i.e., Definition 3.1) if $u_{k_{1}k_{2}}=f_{k_{2}}(\hat{\mathbf{x}})/f_{k_{1}}(\hat{\mathbf{x}}),\forall k_{1},k_{2}\in\{1,\ldots,K\}$. Note that these specific scaling factors also satisfy the properties described above.

Similarly, we define an absolute version of trade-off preservation, where the given solution’s objective values are perturbed by the same absolute amount. Such an approach may be preferred when the objective values are measured in the same units over the same range.

Definition 3.3** (Absolute trade-off preservation)**

A solution $\mathbf{x}^{*}$ preserves the absolute trade-off encoded in $\hat{\mathbf{x}}$ if the following equality holds for every pair of objectives $k_{1}$ and $k_{2}$:

[TABLE]

A weight vector $\boldsymbol{\alpha}^{*}$ preserves the absolute trade-off encoded in $\hat{\mathbf{x}}$ if $\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$.

Absolute trade-off preservation is a special case of general trade-off preservation in Definition 3.1 where $u_{k_{1}k_{2}}=1,\forall k_{1},k_{2}\in\{1,\ldots,K\}$. Similar to the relative case, we will show in the next subsection that absolute trade-off preservation has a direct relationship with minimizing the absolute duality gap in an inverse optimization model.

3.2 Inverse optimization models with trade-off preservation

In this subsection, we propose trade-off-preserving inverse multi-objective convex optimization models. We start with the model that takes into account the general definition of trade-off preservation and then present the models for the relative and absolute specializations.

3.2.1 Inverse optimization model with general trade-off preservation

We formulate a convex optimization model that determines a weight vector $\boldsymbol{\alpha}^{*}$ and corresponding optimal solution $\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$ that satisfy equation (11). To do so, we first define $\mu_{k}:=u_{\tilde{k}k},k=1,\ldots,K,$ where $\tilde{k}\in\{1,\ldots,K\}$ is an arbitrarily chosen reference objective. Given a set of objectives and an input solution $\hat{\mathbf{x}}$, the following formulation, which we call the inverse optimization problem (IOP), provides a necessary and sufficient condition for determining whether there exists a $\mathbf{x}^{*}\in\Omega$ that satisfies general trade-off preservation:

[TABLE]

Theorem 3.1

Let $(\epsilon^{*},\mathbf{x}^{*})$ be an optimal solution to problem (14) and $\boldsymbol{\alpha}^{*}$ be a vector of optimal Lagrange multipliers associated with constraints (14b). Then,

• (a)

There exists a solution in $\Omega$ that satisfies general trade-off preservation for $\hat{\mathbf{x}}$ if and only if there exists $(\epsilon^{*},\mathbf{x}^{*})$ satisfying the first $K$ constraints with equality.

• (b)

$\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$.

• Proof:

(a) ($\Leftarrow$) Let $(\epsilon^{*},\mathbf{x}^{*})$ be an optimal solution to problem (14) with the first $K$ constraints being tight. Then

[TABLE]

which implies that $\mathbf{x}^{*}$ satisfies trade-off preservation. All that remains is to prove $\mathbf{x}^{*}\in\Omega$. Constraints (14b) and (14c) imply that $\mathbf{x}^{*}\in\mathbf{X}$. Let $\boldsymbol{\alpha}\geq\mathbf{0}$, $\boldsymbol{\sigma}\geq\mathbf{0}$, and $\boldsymbol{\pi}$ be the Lagrange multipliers associated with the first, second, and third sets of constraints of problem (14), respectively. Consider the Lagrangian associated with problem (14):

[TABLE]

A solution ($\epsilon^{*},\mathbf{x}^{*}$) is optimal for problem (14) if and only if there exists $(\boldsymbol{\alpha},\boldsymbol{\sigma},\boldsymbol{\pi})\in\mathbb{R}_{+}^{K}\times\mathbb{R}_{+}^{L}\times\mathbb{R}^{m}$ that satisfies the following system of equations:

[TABLE]

Equations (17b) and (17d) form the KKT conditions for the original forward optimization problem (1). From (17a), $\boldsymbol{\alpha}$ is a nonzero vector. Thus, if a solution $\mathbf{x}^{*}$ is optimal for problem (14), $\mathbf{x}^{*}\in\Omega$.

($\Rightarrow$) Let $\mathbf{x}^{*}\in\Omega$ satisfy the definition of general trade-off preservation, i.e., for any two objectives $k_{1},k_{2}\in\{1,\ldots,K\}$, u_{k_{1}k_{2}}\Big{(}f_{k_{1}}(\mathbf{x}^{*})-f_{k_{1}}(\hat{\mathbf{x}})\Big{)}=\Big{(}f_{k_{2}}(\mathbf{x}^{*})-f_{k_{2}}(\hat{\mathbf{x}})\Big{)}. The KKT conditions for the forward problem imply that there exists $(\boldsymbol{\alpha}^{*},\boldsymbol{\sigma}^{*},\boldsymbol{\pi}^{*})$ that satisfies (17b) and (17d) with at least one $k$ for which $\alpha_{k}^{*}>0$. What remains is to show that (17a) and (17c) are also satisfied. Since u_{k_{1}k_{2}}\Big{(}f_{k_{1}}(\mathbf{x}^{*})-f_{k_{1}}(\hat{\mathbf{x}})\Big{)}=\Big{(}f_{k_{2}}(\mathbf{x}^{*})-f_{k_{2}}(\hat{\mathbf{x}})\Big{)},k_{1},k_{2}=1,\ldots,K, the following is true for some objective $\tilde{k}$:

[TABLE]

which implies that (17c) is satisfied. Finally, equation (17a) can be satisfied through a re-scaling of $\boldsymbol{\alpha}^{*}$, which is possible since $\boldsymbol{\alpha}^{*}$ is not identically zero. Thus, $(\epsilon^{*},\mathbf{x}^{*})$ is optimal for (14).

(b) Since $\mathbf{x}^{*}$ and $\boldsymbol{\alpha}^{*}$ satisfy (17b) and (17d), which are the KKT conditions for the forward problem (1), $\mathbf{x}^{*}\in\Omega(\boldsymbol{\alpha}^{*})$. $\square$

Theorem 3.1 suggests that solving IOP$(\hat{\mathbf{x}})$ can simultaneously identify, if they exist, a trade-off preserving solution $\mathbf{x}^{*}\in\Omega$ and a corresponding optimal weight vector $\boldsymbol{\alpha}^{*}$. Theorem 3.1 also implies that if there does not exist an optimal solution to problem (14) that satisfies the first $K$ constraints with equality, then there is no trade-off preserving solution in $\Omega$.

Note that the choice of the reference objective $\tilde{k}$ in the definition of $\mu_{k}$ simply scales the optimal objective value for problem (14), without changing the feasible region of (14). Because the functions $f_{k}(\cdot)$ and $g_{l}(\cdot)$ are convex functions, formulation (14) is a convex optimization problem.

3.2.2 Inverse optimization model with relative trade-off preservation

The following inverse model determines a weight vector $\boldsymbol{\alpha}^{*}$ and $\mathbf{x}^{*}\in\text{FOP}(\boldsymbol{\alpha}^{*})$ that satisfy relative trade-off preservation given a solution $\hat{\mathbf{x}}$:

[TABLE]

Theorem 3.1 holds for problem (19) because (19) is a special case of the general inverse model (14) where $\mu_{k}=u_{\tilde{k}k}=f_{k}(\hat{\mathbf{x}})/f_{\tilde{k}}(\hat{\mathbf{x}}),\forall k\in\{1,\ldots,K\}$. The derivation is straightforward and omitted. As mentioned previously, the definition of relative trade-off preservation provides a simple and intuitive connection with minimizing relative duality gap.

Proposition 3.1

Given $\hat{\mathbf{x}}\in\mathbf{X}$, let $(\epsilon^{*},\mathbf{x}^{*})$ be an optimal solution to $\text{\emph{IOP}}_{r}(\hat{\mathbf{x}})$. Then $1/\epsilon^{*}$ is the minimum relative duality gap with respect to $\hat{\mathbf{x}}$.

• Proof:

Consider the Lagrangian associated with $\text{IOP}_{r}(\hat{\mathbf{x}})$ (i.e., (19)):

[TABLE]

and the corresponding dual problem: $\underset{\boldsymbol{\alpha},\boldsymbol{\delta},\boldsymbol{\pi}}{\max}\;\underset{\epsilon,\mathbf{x}}{\min}\;L(\boldsymbol{\alpha},\boldsymbol{\sigma},\boldsymbol{\pi},\epsilon,\mathbf{x})$. Let $(\epsilon^{*},\mathbf{x}^{*})$ be an optimal solution to (19). Then we have $\sigma_{l}g_{l}(\mathbf{x}^{*})=0,\forall l=1,\ldots,L$ and $\mathbf{A}\mathbf{x}^{*}=\mathbf{b}$, and need $\boldsymbol{\alpha}$ to satisfy $\sum_{k=1}^{K}\alpha_{k}f_{k}(\hat{\mathbf{x}})=1$ (i.e., $\nabla_{\epsilon}L=0$), which cancel out all the terms but $\sum_{k=1}^{K}\alpha_{k}f_{k}(\mathbf{x})$ in (20). Thus,

[TABLE]

and therefore

[TABLE]

Next, consider a problem that determines the minimum relative duality gap with respect to $\hat{\mathbf{x}}$:

[TABLE]

For any feasible solution $(\tilde{\boldsymbol{\alpha}},\tilde{\mathbf{x}})$ for (23), we can find a scaled vector $\boldsymbol{\alpha}^{*}=\frac{1}{\tilde{\boldsymbol{\alpha}}^{\prime}\mathbf{f}(\hat{\mathbf{x}})}\tilde{\boldsymbol{\alpha}}$ such that ${\boldsymbol{\alpha}^{*}}^{\prime}\mathbf{f}(\hat{\mathbf{x}})=1$, $\tilde{\mathbf{x}}\in\Omega(\boldsymbol{\alpha}^{*})$, and thus $(\boldsymbol{\alpha}^{*},\tilde{\mathbf{x}})$ is feasible for (22b) with the same objective value. Therefore, the optimal objective value for (23) coincides with $1/\boldsymbol{\epsilon}^{*}$, as desired. $\square$

Note that the complementary slackness conditions associated with problem (19) provide a sufficient condition for when the relative trade-off preservation is satisfied, which is given in the next result. The proof is straightforward and omitted.

Corollary 3.1

Let $\mathbf{x}^{*}$ be an optimal solution to IOP${}_{r}(\hat{\mathbf{x}})$ (i.e., problem (19)) and $\boldsymbol{\alpha}^{*}$ be a vector of optimal Lagrange multipliers associated with constraints (19b). If $\alpha_{k}^{*}>0$ for all $k=1,\ldots,K$, then $\mathbf{x}^{*}$ satisfies relative trade-off preservation for $\hat{\mathbf{x}}$.

The following example illustrates that if there is at least one $k$ such that $\alpha_{k}^{*}=0$, then there may or may not be a solution in $\Omega$ that preserves the relative trade-off.

Example 3.1

Consider Example 2.1 again and two new solutions $\hat{\mathbf{x}}^{d}=(1.725,1.121)$ and $\hat{\mathbf{x}}^{e}=(1.789,1.096)$. Figure 2 shows the feasible region and the Pareto set in objective space (i.e., $\mathbf{f}(\mathbf{X})$ and $\mathbf{f}(\Omega)$, respectively), along with the two points $\mathbf{f}(\hat{\mathbf{x}}^{d})=(13.160,8.004)$ and $\mathbf{f}(\hat{\mathbf{x}}^{e})=(14.000,8.004)$. Solving IOPr($\hat{\mathbf{x}}^{d}$) and IOPr($\hat{\mathbf{x}}^{e}$) both return the same optimal weight vector $\boldsymbol{\alpha}^{*}=(0,1)$ and the same forward optimal solution $\mathbf{x}^{*}=\bar{\mathbf{x}}^{d}\in\Omega(\boldsymbol{\alpha}^{*})$. It is clear from Figure 2 that $\bar{\mathbf{x}}^{d}$ preserves the trade-off preference of $\hat{\mathbf{x}}^{d}$ but not $\hat{\mathbf{x}}^{e}$. Because $\mathbf{f}(\bar{\mathbf{x}}^{e})\not\in\mathbf{f}(\Omega)$, there is no non-zero weight vector that makes $\bar{\mathbf{x}}^{e}$ optimal to the forward problem. The case of $\hat{\mathbf{x}}^{d}$ also illustrates the possibility of degeneracy, i.e., $\epsilon^{*}=f_{1}(\mathbf{x}^{*})/f_{1}(\hat{\mathbf{x}}^{d})$ while $\alpha_{1}=0$.

Example 3.1 illustrates the intuitive fact that when a given $\hat{\mathbf{x}}$ is sufficiently “inferior” with respect to a particular objective $k$, the optimal weight vector determined from inverse optimization will return $\alpha^{*}_{k}=0$. In this case the relative trade-off is not perfectly preserved. However, the objectives with non-zero weights will still satisfy relative trade-off preservation.

3.2.3 Inverse optimization model with absolute trade-off preservation

Similar to problem (19) for relative trade-off preservation, a necessary and sufficient condition for determining whether there exists a solution $\mathbf{x}^{*}\in\Omega$ that satisfies absolute trade-off preservation for $\hat{\mathbf{x}}$ can be established by the following convex optimization model:

[TABLE]

Similar to the relative model, Theorem 3.1 holds for (24) because it is a special case of the general model (14) where $u_{k_{1}k_{2}}=1$ for all $k_{1},k_{2}\in\{1,\ldots,K\}$. Similar to Proposition 3.1, Proposition 3.2 suggests that a weight vector that produces a solution $\mathbf{x}^{*}$ that satisfies absolute trade-off preservation for $\hat{\mathbf{x}}$ minimizes the absolute optimality gap with respect to $\hat{\mathbf{x}}\in\mathbf{X}$. The proof of this result involves minor modifications to the proof of Proposition 3.1 and is omitted.

Proposition 3.2

Given $\hat{\mathbf{x}}\in\mathbf{X}$, let $(\epsilon^{*},\mathbf{x}^{*})$ be an optimal solution to $\text{\emph{IOP}}_{a}(\hat{\mathbf{x}})$. Then $-\epsilon^{*}$ is the minimum absolute duality gap with respect to $\hat{\mathbf{x}}$. That is, $-\epsilon^{*}=\min\{\boldsymbol{\alpha}^{\prime}\mathbf{f}(\hat{\mathbf{x}})-\boldsymbol{\alpha}^{\prime}\mathbf{f}(\mathbf{x})\,|\,\mathbf{x}\in\Omega(\boldsymbol{\alpha}),\boldsymbol{\alpha}\geq\mathbf{0}\}\geq 0$.

Finally, the following result provides a sufficient condition for when absolute trade-off preservation is satisfied. Again, the proof is straightforward and omitted.

Corollary 3.2

Let $\mathbf{x}^{*}$ be an optimal solution to problem (24) and $\boldsymbol{\alpha}^{*}$ be a vector of optimal Lagrange multipliers associated with constraints (24b). If $\alpha_{k}^{*}>0$ for all $k=1,\ldots,K$, then $\mathbf{x}^{*}$ satisfies absolute trade-off preservation for $\hat{\mathbf{x}}$.

3.3 Linearized inverse optimization models with trade-off preservation

In this section we formulate a linear approximation to IOP($\hat{\mathbf{x}}$) (i.e., formulation (14)). The reasons for presenting the linearization are two-fold. First, the linearized model helps elucidate connections with other inverse models in the literature, as we will discuss in the next section. Second, as we will demonstrate in the computational results, the linear approximation approach is computationally efficient to implement and may generate a solution that is close to the optimal solution from the exact inverse model, either as a stand-alone model or as part of a specialized solution methodology for the exact model.

The linearized inverse optimization problem (LIOP) for the general model (14) can be formulated by linearizing the first two sets of constraints of (14) around some point $\tilde{\mathbf{x}}$:

[TABLE]

Since $f_{k}$ and $g_{l}$ are convex, it is straightforward to see that formulation (25) forms an outer approximation – and thus provides a lower bound – to formulation (14). If $\tilde{\mathbf{x}}=\mathbf{x}^{*}$, where $\mathbf{x}^{*}$ is an optimal solution to problem (14), then an optimal Lagrange multiplier for constraints (25b) is also optimal for (14). Unfortunately, $\mathbf{x}^{*}$ is not known a priori. One option is simply to solve (25) with $\tilde{\mathbf{x}}=\hat{\mathbf{x}}$ (i.e., solve LIOP$(\hat{\mathbf{x}},\hat{\mathbf{x}})$, which for brevity we refer to as LIOP$(\hat{\mathbf{x}})$). Note that linearization may render problem (25) unbounded. In this case, a trust region $\mathbf{x}\in[\hat{\mathbf{x}}-\kappa\mathbf{e},\hat{\mathbf{x}}+\kappa\mathbf{e}]$ for some $\kappa$ can be added, where $\mathbf{e}$ is the vector of ones (Bazaraa et al. 2006).

Building on this idea, since (25) is a first-order approximation of (14), we can employ the well-known successive linear programming (SLP) algorithm (Zhang et al. 1985, Bazaraa et al. 2006) to solve formulation (14). The algorithm repeatedly solves the linear problem (25) using a trust region approach, generating an optimal solution $\mathbf{x}^{*}_{i}$ in iteration $i$ that is used as the input vector $\tilde{\mathbf{x}}$ in iteration $i+1$. In Section 4, we implement the SLP algorithm based on Zhang et al. (1985) – we refer the reader to that paper for details of the algorithm and a proof of convergence.

3.4 Relationships between inverse optimization models

Focusing on the general model, if $\hat{\mathbf{x}}\in\Omega$, then an optimal $\boldsymbol{\alpha}^{*}$ derived from one of IOP($\hat{\mathbf{x}}$), LIOP($\hat{\mathbf{x}}$), or KES($\hat{\mathbf{x}}$) is optimal for all of them. If $f_{k},\,k=1,\ldots,K$ and $g_{l},\,l=1,\ldots,L$ are linear, then IOP($\hat{\mathbf{x}}$) and LIOP($\hat{\mathbf{x}}$) are equivalent; particularly IOP($\hat{\mathbf{x}}$) with the relative and absolute trade-off preservation specializes to models of Chan et al. (2014) . The linearization described in the previous section provides a bridge between our inverse formulation (14) and the model of Keshavarz et al. (2011). The next result draws an equivalence between LIOP($\hat{\mathbf{x}}$) and KES($\hat{\mathbf{x}}$), thus illustrating how the model of Keshavarz et al. (2011) can be modified to take into account trade-off preservation.

Proposition 3.3

An optimal weight vector solution to KES$(\hat{\mathbf{x}})$ is identical to an optimal weight vector solution to the dual of LIOP$(\hat{\mathbf{x}})$, if $\boldsymbol{\delta}=\mathbf{0},\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})=-\boldsymbol{\gamma}^{\prime}\mathbf{e}+\boldsymbol{\rho}^{\prime}\mathbf{e}$, and $\boldsymbol{\alpha}$ is normalized by the constraint $\sum_{k=1}^{K}\mu_{k}\alpha_{k}=1$ instead of $\alpha_{1}=1$.

• Proof:

Consider the dual of LIOP$(\hat{\mathbf{x}})$:

[TABLE]

We take the inner product of constraint (26c) with $\hat{\mathbf{x}}$ to get

[TABLE]

By adding the left-hand side of equation (27) to the objective function of problem (26), we obtain the following problem:

[TABLE]

KES($\hat{\mathbf{x}}$) with $\boldsymbol{\delta}=\mathbf{0},\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})=-\boldsymbol{\gamma}^{\prime}\mathbf{e}+\boldsymbol{\rho}^{\prime}\mathbf{e}$, and $\boldsymbol{\alpha}$ normalized by the constraint $\sum_{k=1}^{K}\mu_{k}\alpha_{k}=1$ is

[TABLE]

which is equivalent to problem (28). $\square$

Note that if $\hat{\mathbf{x}}\in\mathbf{X}$, the form of $\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})$ in the statement of Proposition 3.3 can be replaced with $\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})=||\boldsymbol{\gamma}||_{1}+||\boldsymbol{\rho}||_{1}$. Proposition 3.3 states that by specifying $\boldsymbol{\delta}$, the function $\phi(\boldsymbol{\delta},\boldsymbol{\gamma},\boldsymbol{\rho})$, and the normalization constraint in KES$(\hat{\mathbf{x}})$, a decision maker can preserve the trade-off encoded in $\hat{\mathbf{x}}$ in the inverse optimization process. The following result shows how KES($\hat{\mathbf{x}}$), which was initially proposed for single-objective problems, can be interpreted in the presence of multiple objectives.

Corollary 3.3

The linear approximation of the general inverse model (25) with $\mu_{1}=1$ and $\mu_{k}=0$ for $k=2,\ldots,K$ is equivalent to KES($\hat{\mathbf{x}}$).

Corollary 3.3 indicates that KES($\hat{\mathbf{x}}$) is an approximate version of the general inverse model (14) where $\mu_{1}=1$ and $\mu_{k}=0$ for $k=2,\ldots,K$. That is, applying KES($\hat{\mathbf{x}}$) to multi-objective optimization problems puts its entire emphasis on the adjustment of the first objective and in general does not take into account the trade-off across multiple objectives. This observation is in fact reinforced by our computational results in the next section.

4 Computational results

In this section, we use data from prostate cancer radiation therapy and perform a retrospective study to demonstrate several inverse models presented in this paper. We use historical prostate cancer treatments as input into the inverse models to derive objective function weights that would approximately reproduce the treatments (via forward optimization). We compare weights generated by the different inverse models and examine how the clinical criteria are affected by the weights. While the experiments constitute a retrospective analysis with historical treatment data, we also comment on clinical implications of our results.

4.1 Background and motivation

Radiation therapy is one of the primary methods to treat many cancers including prostate cancer (Foroudi et al. 2003). In prostate cancer treatment, the primary target is the prostate, while the healthy organs to be avoided include the rectum, bladder, and femoral heads (tops of the femur bones). Beams of radiation are delivered to the patient targeting a tumor at their intersection. The patient’s anatomy is discretized into small volume elements called voxels. A radiation beam can be modeled as a set of small beamlets, whose intensities are optimized. A radiation therapy treatment planning problem is typically a multi-objective optimization problem and the traditional approach to solving these problems clinically is through a trial-and-error approach to weight selection.

In response to the time consuming, trial-and-error approach to determining the weight values, knowledge-based treatment planning is an increasingly popular approach. Such an approach uses a database of historical treatments to derive a treatment plan for a de novo patient (e.g., Wu et al. (2011), Moore et al. (2011), Chanyavanich et al. (2011), Lee et al. (2013), and Boutilier et al. (2015)). The idea is that planning effort may be reduced if we can leverage the vast knowledge and experience accumulated through creating clinically accepted plans for similar patients in the past. Naturally, the historical treatments implicitly encode information on accepted trade-offs between conflicting clinical criteria (i.e., objectives). Thus, any knowledge-based effort for determining the weight values should consider maintaining the clinical metrics (i.e., objective values) achieved by the historical treatments.

In the computational results below, we apply different inverse optimization models proposed in this paper so as to impute weight values for a convex quadratic treatment planning problem and demonstrate to what extent the clinical metrics are preserved by each of the models. Specifically, for each historical treatment, which may not be mathematically exactly optimal for the underlying treatment planning problem, the inverse models determine sets of weight values that can reproduce the plan as closely as possible. We use relative trade-off preservation in this application because the clinical metrics, while all measured in terms of dose, may take on values in largely varying ranges (e.g., max dose values in the 70’s and penalties for an organ exceeding a dose threshold in the single digits). We then show that weight values that approximate the given treatment under the consideration of trade-off preservation indeed generate a new treatment plan that is considered most similar clinically.

4.2 Forward optimization problem

Let $\mathbf{x}\in\mathbb{R}^{n}_{+}$ be a vector of beamlet intensities and $\mathbf{D}_{k}\in\mathbb{R}^{m_{k}\times n}_{+}$ be a matrix that quantifies the dose deposited to each voxel in structure $k$ from unit intensity of each beamlet, where $m_{k}$ denotes the number of voxels in structure $k$ and $n$ denotes the number of beamlets. Let $\mathcal{T}$ index the tumor. The vector $\mathbf{D}_{k}\mathbf{x}\in\mathbb{R}^{m_{k}}_{+}$ quantifies the dose to every voxel in structure $k$. The radiation therapy treatment planning problem is a forward optimization problem. We present a simple forward model below. With each healthy organ $k$, we associate an objective function $f_{k}(\mathbf{x})=||(\mathbf{D}_{k}\mathbf{x}-\theta_{k}\mathbf{e})_{+}||_{2}^{2}$, where $\theta_{k}$ denotes a dose threshold for structure $k$, and $(\cdot)_{+}$ denotes a vector with the operator $\max\{0,\cdot\}$ applied to each component.

[TABLE]

The objective penalizes delivering dose above a certain threshold for structure $k$. Formulation (30) can be written as a convex program by rewriting $f_{k}(\mathbf{x})$ as $f_{k}(\mathbf{z})=||\mathbf{z}||_{2}^{2}$, where $\mathbf{z}$ is an auxiliary decision vector satisfying $\mathbf{z}\geq\mathbf{D}_{k}\mathbf{x}-\theta_{k}\mathbf{e}$ and $\mathbf{z}\geq\mathbf{0}$. We use hard constraints (30b) to bound the lower and upper doses to organ $k$ and the tumor $\mathcal{T}$. Constraint (30c) is a stylized approach to discouraging individual beamlets from significantly exceeding the mean beamlet intensity. In the experiments below, five healthy organs-at-risk (OAR) comprise the objective function: the bladder, rectum, left and right femoral heads, and a ring of healthy tissue around the tumor (used to encourage conformity of the dose around the target). We choose $\beta=2$, $\mathbf{l}_{k}=\mathbf{0}$ for all healthy organs, $\mathbf{l}_{\mathcal{T}}=78\mathbf{e}$, $\mathbf{u}_{k}=81.9\mathbf{e}$ for all healthy organs and the tumor, $\theta_{k}=30$ for the left and right femoral heads, and $\theta_{k}=50$ for the other healthy organs, based on a protocol at Princess Margaret Cancer Centre in Toronto, Canada.

4.3 Impact of different inverse optimization models on trade-off preservation

We used 24 historical treatments delivered at Princess Margaret Cancer Centre for our input $\hat{\mathbf{x}}$ vectors. Note that a clinical treatment plan is typically designed via the repeated solution of a nonconvex optimization problem with a larger number of objectives. Thus, the historical treatment plans are not likely to be Pareto optimal for problem (30). The goal of inverse optimization here is to find a weight vector such that problem (30) generates a treatment plan that is similar to the historical one, in terms of the achieved clinical metrics (i.e., objective values).

On average across the 24 patients, $n=409$ and $m_{k}=8,157$, resulting in $6,313$ variables and $14,770$ constraints (after including auxiliary variables) of the forward problem (30). For each of the treatments $\hat{\mathbf{x}}$, we derived weights using KES($\hat{\mathbf{x}}$), IOPr($\hat{\mathbf{x}}$), LIOPr($\hat{\mathbf{x}}$), and the SLP algorithm. For KES($\hat{\mathbf{x}}$), we solved five different instances, corresponding to five different ways of normalizing the weight vector (i.e., $\alpha_{k}=1,k=1,\ldots,5$). For the SLP algorithm, parameters and termination criteria were chosen based on Zhang et al. (1985) (e.g., terminating when the $l_{2}$ norm difference between two consecutive iterates is less than 0.001) and $\hat{\mathbf{x}}$ was used as an initial solution. All the models were solved using CPLEX 12.3 on a computer with a 3.07 GHz 12-core CPU and 32 GB of RAM.

Table 1 summarizes the results of applying the five variations of the KES($\hat{\mathbf{x}}$) model to patient #1. It can be seen that the imputed weight values are heavily dependent on which objective was used for normalization in KES($\hat{\mathbf{x}}$). The weights in Table 1 suggest that the KES model gives the vast majority of the weight to the objective whose weight is used in the normalization constraint. The results were similar for the remaining patients. For each instance, $f_{k}(\mathbf{x}^{*})/f_{k}(\hat{\mathbf{x}})$ denotes the ratio of the objective value achieved by $\boldsymbol{\alpha}^{*}$ to the objective value associated with $\hat{\mathbf{x}}$ for objective $k$. In Table 1, we see that there is no well-defined pattern in the ratios. The amount of improvement varies from objective to objective and from model to model. Furthermore, KES1 provides an example where the dose on some objectives stays the same or even increases while the dose to other organs decreases. Overall, the KES model finds a weight vector that puts an emphasis only on a single objective - whichever was used for normalization - and the resulting treatment plans seem to vary significantly.

Table 2 shows the results from applying the IOPr model, LIOPr model, and SLP algorithm to three example patients. For the SLP algorithm, the solution time is the total running time through all iterations. As expected, the component-wise ratio obtained by the IOPr model is equivalent to $\epsilon^{*}$ unless the corresponding weight is zero in which case the ratio is at most $\epsilon^{*}$ (see Corollary 3.1). The practical interpretation of this result is that although the ratio $f_{k}(\mathbf{x}^{*})/f_{k}(\hat{\mathbf{x}})=\epsilon^{*}$ may be achievable, the objective value can be further decreased (i.e., improved) without any sacrifice. Note that preserving the trade-off using the IOPr model comes at a higher computational cost compared to the LIOPr model which approximately maintains the relative preservation. Also note that $\epsilon^{*}$ from LIOPr is less than $\epsilon^{*}$ from IOPr as expected, because LIOPr is an outer approximation to IOPr. Finally, the SLP algorithm strikes a middle ground between the LIOPr and IOPr models in terms of trade-off preservation and computational efficiency.

Figure 3 shows the performance of KES, IOPr, LIOPr and SLP across all 24 patients in terms of the solution time and trade-off preservation, the latter of which is quantified using the variance of the component-wise ratios $f_{k}(\mathbf{x}^{*})/f_{k}(\hat{\mathbf{x}})$ across $k$. Variances less than $2^{-14}(<0.0001)$ were considered to represent perfect trade-off preservation so their values were rounded to zero. For the KES results for each patient, we report the performance of the model with the normalization constraint $\alpha_{k^{\prime}}=1$, where $k^{\prime}$ is the structure with the highest inverse weight determined by the IOPr model. For example, in the case of patient #1, the IOPr model shows that the rectum weight is the highest, so we use the inverse weights derived from the KES model that employs the normalization constraint $\alpha_{2}=1$ (the rectum is the second objective). In Figure 3, we can see that the KES model can generally be solved quickly but exhibits high variance in the objective function ratios, while the IOPr model generally has the longest solution time but preserves the initial relative trade-off vector well. The three patients with variance in the range of $2^{-6}$ to $2^{-4}$ using the IOPr model had at least one zero-weighted objective. Finally, the LIOPr model and the SLP algorithm by extension strike a balance between the KES and IOPr models. Table 3 summarizes the numerical results across all patients, reinforcing the trade-offs between the models.

4.4 Impact of different inverse optimization models on clinical metrics

While Tables 1, 2, and 3 reinforce our theoretical results about the lack of trade-off preservation leading to inconsistent perturbations in objective space, we next comment on how these results translate into the clinical context for radiation therapy. In particular, we examine the impact of trade-off preservation on the resulting treatment plan’s ability to replicate the historical treatment plan with respect to key clinical metrics.

We first compared treatment plans generated by the KES weights to the clinical treatment plan for patient #1. We use a dose volume histogram (DVH), which shows for each plan the fractional volume of an organ that receives a certain dose or higher, to compare the KES and clinical plans. Figure 4 shows that while the solution using the KES1 weights (dashed lines) has significantly improved the bladder dose (dashed green), the rectum dose is generally worse, violating the clinical acceptability requirement of Princess Margaret Cancer Centre that less than 50% of the volume receive more than 50 Gy. On the other hand, Figure 5 shows that the treatment plan generated by the IOPr weights for patient #1 is much more similar to the clinical plan over all the four organs at risk. Given that the bladder, rectum, and the left/right femoral head objectives were to minimize radiation dose exceeding 50 Gy, 50 Gy, and 30 Gy, respectively, the DVHs of the IOPr plan show a similar dose reduction across the organs compared to the clinical DVHs above the dose thresholds. The results were similar for the remaining patients.

The ability to accurately replicate a clinical plan, including the relevant clinical trade-offs encoded in the plan, is an important feature of any knowledge-based, automated planning system. Such approaches are becoming increasingly popular as a means to efficiently create treatments for new patients. For example, given a new patient, the system can identify a “similar” historical patient (e.g., anatomically) and produce an appropriate weight vector using our trade-off preserving inverse optimization model. Alternatively, one can use a prediction model, trained using inversely optimized weight vectors and anatomy from historical patients, to generate an appropriate weight vector for a new patient (Wu et al. 2011, Moore et al. 2011, Chanyavanich et al. 2011, Lee et al. 2013, Boutilier et al. 2015). Such an approach is more likely to produce plans with acceptable clinical performance, assuming the historical plans were acceptable.

5 Conclusion

Trade-offs between objectives are critical in multi-objective optimization. In this paper, we developed a new approach to inverse convex multi-objective optimization that explicitly considers preserving trade-off encoded in a given initial solution. Our approach is general and maintains the complexity of the forward problem (i.e., is convex). The notion of trade-off preservation is generally applicable to multi-objective optimization problems where a decision maker’s intentions are encoded in the objective values of a given solution, and provides a connection with a relative or absolute optimality gap with respect to the given solution. We formulated a general inverse convex optimization framework that encompasses many of the inverse models in the literature and demonstrated how an inverse model can be adjusted to preserve the initial trade-off. We compared various inverse optimization models computationally using data from prostate cancer radiation therapy and demonstrated the ability of our inverse model to preserve the clinical quality of a historical treatment using imputed weights, which is important in the design of next generation automated treatment planning systems.

Acknowledgements

The authors would like to thank Dr. Michael B. Sharpe and Dr. Tim Craig at Princess Margaret Cancer Centre for providing the clinical data and expertise.