Risk Bounds and Calibration for a Smart Predict-then-Optimize Method
Year: 2021

Abstract:
The predict-then-optimize framework is fundamental in practical stochastic decision-making problems: first predict unknown parameters of an optimization model, then solve the problem using the predicted values. A natural loss function in this setting is defined by measuring the decision error induced by the predicted parameters, which was named the Smart Predict-then-Optimize (SPO) loss by Elmachtoub and Grigas [2021]. Since the SPO loss is typically nonconvex and possibly discontinuous, Elmachtoub and Grigas [2021] introduced a convex surrogate, called the SPO+ loss, that importantly accounts for the underlying structure of the optimization model. In this paper, we greatly expand upon the consistency results for the SPO+ loss provided by Elmachtoub and Grigas [2021]. We develop risk bounds and uniform calibration results for the SPO+ loss relative to the SPO loss, which provide a quantitative way to transfer the excess surrogate risk to excess true risk. By combining our risk bounds with generalization bounds, we show that the empirical minimizer of the SPO+ loss achieves low excess true risk with high probability. We first demonstrate these results in the case when the feasible region of the underlying optimization problem is a polyhedron, and then we show that the results can be strengthened substantially when the feasible region is a level set of a strongly convex function. We perform experiments to empirically demonstrate the strength of the SPO+ surrogate, as compared to standard `1 and squared `2 prediction error losses, on portfolio allocation and cost-sensitive multi-class classification problems.

text: The predict-then-optimize framework is fundamental in practical stochastic decision-making problems: first predict unknown parameters of an optimization model, then solve the problem using the predicted values. A natural loss function in this setting is defined by measuring the decision error induced by the predicted parameters, which was named the Smart Predict-then-Optimize (SPO) loss by Elmachtoub and Grigas [2021]. Since the SPO loss is typically nonconvex and possibly discontinuous, Elmachtoub and Grigas [2021] introduced a convex surrogate, called the SPO+ loss, that importantly accounts for the underlying structure of the optimization model. In this paper, we greatly expand upon the consistency results for the SPO+ loss provided by Elmachtoub and Grigas [2021]. We develop risk bounds and uniform calibration results for the SPO+ loss relative to the SPO loss, which provide a quantitative way to transfer the excess surrogate risk to excess true risk. By combining our risk bounds with generalization bounds, we show that the empirical minimizer of the SPO+ loss achieves low excess true risk with high probability. We first demonstrate these results in the case when the feasible region of the underlying optimization problem is a polyhedron, and then we show that the results can be strengthened substantially when the feasible region is a level set of a strongly convex function. We perform experiments to empirically demonstrate the strength of the SPO+ surrogate, as compared to standard `1 and squared `2 prediction error losses, on portfolio allocation and cost-sensitive multi-class classification problems.

heading: 1 Introduction
text: The predict-then-optimize framework, where one predicts the unknown parameters of an optimization model and then plugs in the predictions before solving, is prevalent in applications of machine learning. Some typical examples include predicting future asset returns in portfolio allocation problems and predicting the travel time on each edge of a network in navigation problems. In most cases, there are many contextual features available, such as time of day, weather information, financial and business news headlines, and many others, that can be leveraged to predict the unknown parameters and reduce uncertainty in the decision making problem. Ultimately, the goal is to produce a high quality prediction model that leads to a good decisions when implemented, such as a position that leads to a large return or a route that induces a small realized travel time. There has been a fair amount of recent work examining this paradigm and other closely related problems in data-driven decision making, such as the works of Bertsimas and Kallus [2020], Donti et al. [2017], Elmachtoub and Grigas [2021], Kao et al. [2009], Estes and Richard [2019], Ho and Hanasusanto [2019], Notz and Pibernik [2019], Kotary et al. [2021], the references therein, and others.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In this work, we focus on the particular and important case where the optimization problem of interest has a linear objective with a known convex feasible region and where the contextual features are related to the coefficients of the linear objective. This case includes the aforementioned shortest path and portfolio allocation problems. In this context, Elmachtoub and Grigas [2021] developed the Smart Predict-then-Optimize (SPO) loss function, which directly measures the regret of a prediction against the best decision in hindsight (rather than just prediction error, such as squared error). After the introduction of the SPO loss, recent work has studied its statistical properties, including generalization bounds of the SPO loss function in El Balghiti et al. [2019] and generalization and regret convergence rates in Hu et al. [2020]. Moreover, since the SPO loss is not continuous nor convex in general [Elmachtoub and Grigas, 2021], which makes the training of a prediction model computationally intractable, Elmachtoub and Grigas [2021] introduced a novel convex surrogate loss, referred to as the SPO+ loss. Elmachtoub and Grigas [2021] highlight and prove several advantages of the SPO+ surrogate loss: (i) it still accounts for the downstream optimization problem when evaluating the quality of a prediction model (unlike prediction losses such as the squared `2 loss), (ii) it has a desirable Fisher consistency property with respect to the SPO loss under mild conditions, and (iii) it often performs better than commonly considered prediction losses in experimental results. Unfortunately, although a desirable property of any surrogate loss in this context, Fisher consistency is not directly applicable when one only has available a finite dataset, which is always the case in practice, because it relies on full knowledge of the underlying distribution. Motivated thusly, it is desirable to develop risk bounds that allow one to translate an approximate guarantee on the risk of a surrogate loss function to a corresponding guarantee on the SPO risk. That is, risk bounds (and the related notion of calibration functions) answer the question: to what tolerance δ should the surrogate excess risk be reduced to in order to ensure that the excess SPO risk is at most ? Note that, with enough data, it is possible in practice to ensure a (high probability) bound on the excess surrogate risk through generalization and optimization guarantees.
The main goal of this work is to provide risk bounds for the SPO+ surrogate loss function. Our results, to the best of our knowledge, are the first risk bounds of the SPO+ loss, besides the analysis of the 1-dimensional scenario in Ho-Nguyen and Kılınç-Karzan [2020]. Our results consider two cases for the structure of the feasible region of the optimization problem: (i) the case of a bounded polyhedron, and (ii) the case of a level set of a smooth and strongly convex function. In the polyhedral case, we prove that the risk bound of the SPO+ surrogate is O( 2), where is the desired accuracy for the excess SPO risk. Our results hold under mild distributional assumptions that extend those considered in Elmachtoub and Grigas [2021]. In the strongly convex level set case, we improve the risk bound of the SPO+ surrogate to O( ) by utilizing novel properties of such sets that we develop, namely stronger optimality guarantees and continuity properties. As a consequence of our analysis, we can leverage generalization guarantees for the SPO+ loss to obtain the first sample complexity bounds, with respect to the SPO risk, for the SPO+ surrogate under the two cases we consider. In Section 5, we present computational results that validate our theoretical findings. In particular we present results on entropy constrained portfolio allocation problems which, to the best of our knowledge, is the first computational study of predict-then-optimize problems for a strongly convex feasible region. Our results on portfolio allocation problems demonstrate the effectiveness of the SPO+ surrogate. We also present results for cost-sensitive multi-class classification that illustrate the benefits of faster convergence of the SPO risk in the case of strongly convex sets as compared to polyhedral ones.
Starting with binary classification, risk bounds and calibration have been previously studied in other machine learning settings. Pioneer works studying the properties of convex surrogate loss functions for the 0-1 loss include Zhang et al. [2004], Bartlett et al. [2006], Massart et al. [2006] and Steinwart [2007]. Works including Zhang [2004], Tewari and Bartlett [2007] and Osokin et al. [2017] have studied the consistency and calibration properties of multi-class classification problems, which can be considered as a special case of the predict-then-optimize framework [El Balghiti et al., 2019]. Most related to the results presented herein is the work of Ho-Nguyen and Kılınç-Karzan [2020], who study the uniform calibration properties of the squared `2 loss, and the related work of Hu et al. [2020], who also develop fast sample complexity results for the SPO loss when using a squared `2 surrogate.
Notation. Let represent element-wise multiplication between two vectors. For any positive integer m, let [m] denote the set {1, . . . ,m}. Let Ip denote the p by p identity matrix for any positive integer p. For c̄ ∈ Rd and a positive semi-definite matrix Σ ∈ Rd×d, let N (c̄,Σ) denote the normal distribution P(c) = e − 1 2
(c−c̄)TΣ−1(c−c̄)√ (2π)d det(Σ) . We will make use of a generic given norm ‖ · ‖ on w ∈ Rd,
as well as its dual norm ‖ · ‖∗ which is defined by ‖c‖∗ = maxw:‖w‖≤1 cTw. For a positive definite matrix A, we define the A-norm by ‖w‖A := √ wTAw. Also, we denote the diameter of the set S ⊆ Rd by DS := supw,w′∈S ‖w − w′‖2.

heading: 2 Predict-then-optimize framework and preliminaries
text: We now formally describe the predict-then-optimize framework, which is widely prevalent in stochastic decision making problems. We assume that the problem of interest has a linear objective, but the cost vector of the objective, c ∈ C ⊆ Rd, is not observed when the decision is made. Instead, we observe a feature vector x ∈ X ⊆ Rp, which provides contextual information associated with c. Let P denote the underlying joint distribution of the pair (x, c). Let w denote the decision variable and assume that we have full knowledge of the feasible region S ⊆ Rd, which is assumed to be non-empty, compact, and convex. When a feature vector x is provided, the goal of the decision maker is to solve the contextual stochastic optimization problem:
min w∈S Ec∼P(·|x)[cTw] = min w∈S Ec∼P(·|x)[c]Tw. (1)
As demonstrated by (1), for linear optimization problems the predict-then-optimize framework relies on a prediction of the conditional expectation of the cost vector, namely E[c|x] = Ec∼P(·|x)[c]. Let ĉ denote a prediction of the conditional expectation, then the next step in the predict-then-optimize setting is to solve the deterministic optimization problem with the cost vector ĉ, namely
P (ĉ) : min w∈S
ĉTw. (2)
Depending on the structure of the feasible region S, the optimization problem P (·) can represent linear programming, conic programming, and even (mixed) integer programming, for example. In any case, we assume that we can solve P (·) to any desired accuracy via either a closed-form solution or a solver. Let w∗(·) : Rd → S denote a particular optimization oracle for problem (2), whereby w∗(ĉ) is an optimal solution of P (ĉ). (We assume that the oracle is deterministic and ties are broken in an arbitrary pre-specified manner.)
In order to obtain a model for predicting cost vectors, namely a cost vector predictor function g : X → Rd, we may leverage machine learning methods to learn the underlying distribution P from observed data {(x1, c1), . . . , (xn, cn)}, which are assumed to be independent samples from P. Most importantly, following (1), we would like to learn the conditional expectation and thus g(x) can be thought of as an estimate of E[c|x]. We follow a standard recipe for learning a predictor function g where we specify a loss function to measure the quality of predictions relative to the observed realized cost vectors. In particular, for a loss function `, the value `(ĉ, c) represents the loss or error incurred when the cost vector prediction is ĉ (i.e., ĉ = g(x)) and the realized cost vector is c. Let R`(g;P) := E(x,c)∼P[`(g(x), c)] denote the the risk of given loss function ` and let R∗` (P) = infg′ R`(g′;P) denote the optimal `-risk over all measurable functions g′. Also, let R̂n` (g) := 1 n ∑n i=1 `(g(xi), ci) denote the empirical `-risk. Most commonly used loss functions are based on directly measuring the prediction error, including the (squared) `2 and the `1 losses. However, these losses do not take the downstream optimization task nor the structure of the feasible region S into consideration. Motivated thusly, one may consider a loss function that directly measures the decision error with respect to the optimization problem (2). Elmachtoub and Grigas [2021] formalize this notion in our context of linear optimization problems with their introduction of the SPO (Smart Predict-then-Optimize) loss function, which is defined by
`SPO(ĉ, c) := c Tw∗(ĉ)− cTw∗(c),
where ĉ ∈ Rd is the predicted cost vector and c ∈ C is the realized cost vector. Due to the possible non-convexity and possible discontinuities of the SPO loss, Elmachtoub and Grigas [2021] also propose a convex surrogate loss function, the SPO+ loss, which is defined as
`SPO+(ĉ, c) := max w∈S {(c− 2ĉ)Tw}+ 2ĉTw∗(c)− cTw∗(c).
Importantly, the SPO+ loss still accounts for the downstream optimization problem (2) and the structure of the feasible region S, in contrast to losses that focus only on prediction error. As discussed by Elmachtoub and Grigas [2021], the SPO+ loss can be efficiently optimized via linear/conic
optimization reformulations and with (stochastic) gradient methods for large datasets. Elmachtoub and Grigas [2021] provide theoretical and empirical justification for the use of the SPO+ loss function, including a derivation through duality theory, promising experimental results on shortest path and portfolio optimization instances, and the following theorem which provides the Fisher consistency of the SPO+ loss.
Theorem 2.1 (Elmachtoub and Grigas [2021], Theorem 1). Suppose that the feasible region S has a non-empty interior. For fixed x ∈ X , suppose that the conditional distribution P(·|x) is continuous on all of Rd, is centrally symmetric around its mean c̄ := Ec∼P(·|x)[c], and that there is a unique optimal solution of P (c̄). Then, for all ∆ ∈ Rp it holds that
Ec∼P(·|x) [`SPO+(c̄+ ∆, c)− `SPO+(c̄, c)] = Ec∼P(·|x) [ (c+ 2∆)T (w∗(c)− w∗(c+ 2∆)) ] ≥ 0.
Moreover, if ∆ 6= 0, then Ec∼P(·|x) [`SPO+(c̄+ ∆)− `SPO+(c̄)] > 0.
Notice that Theorem 2.1 holds for arbitrary x ∈ X , i.e., it employs a nonparametric analysis as is standard in consistency and calibration results, whereby there are no constraints on the predicted cost vector associated with x. Under the conditions of Theorem 2.1, given any x ∈ X , we know that the conditional expectation c̄ = Ec∼P(·|x)[c] is the unique minimizer of the SPO+ risk. Furthermore, since c̄ is also a minimizer of the SPO risk, it holds that the SPO+ loss function is Fisher consistent with respect to the SPO loss function, i.e., minimizing the SPO+ risk also minimizes the SPO risk. However, in practice, due to the fact that we have available only a finite dataset and not complete knowledge of the distribution P, we cannot directly minimize the true SPO+ risk. Instead, by employing the use of optimization and generalization guarantees, we are able to approximately minimize the SPO+ risk. A natural question is then: does a low excess SPO+ risk guarantee a low excess SPO risk? More formally, we are primarily interested in the following questions: (i) for any > 0, does there exist δ( ) > 0 such that RSPO+(g;P) − R∗SPO+(P) < δ( ) implies that RSPO(g;P)−R∗SPO(P) < ?, and (ii) what is the largest such value of δ( ) that guarantees the above?
Excess risk bounds via calibration. The notions of calibration and calibration functions provide a useful set of tools to answer the previous questions. We now review basic concepts concerning calibration when using a generic surrogate loss function `, although we are primarily interested in the aforementioned SPO+ surrogate. An excess risk bound allows one to transfer the conditional excess `-risk, E [`(ĉ, c)|x]− infc′ E [`(c′, c)|x], to the conditional excess `SPO-risk, E [`SPO(ĉ, c)|x]− infc′ E [`SPO(c′, c)|x]. Calibration, which we now briefly review, is a central tool in developing excess risk bounds. We adopt the definition of calibration presented by Steinwart [2007] and Ho-Nguyen and Kılınç-Karzan [2020], which is reviewed in Definition 2.1 below.
Definition 2.1. For a given surrogate loss function `, we say ` is `SPO-calibrated with respect to P if there exists a function δ`(·) : R+ → R+ such that for all x ∈ X , ĉ ∈ C, and > 0, it holds that
E [`(ĉ, c)|x]− inf c′ E [`(c′, c)|x] < δ`( )⇒ E [`SPO(ĉ, c)|x]− inf c′ E [`SPO(c′, c)|x] < . (3)
Additionally, if (3) holds for all P ∈ P , where P is a class of distributions on X × C, then we say that ` is uniformly calibrated with respect to the class of distributions P .
A direct approach to finding a feasible δ`(·) function and checking for uniform calibration is by computing the infimum of the excess surrogate loss subject to a constraint that the excess SPO loss is at least . This idea leads to the definition of the calibration function, which we review in Definition 2.2 below.
Definition 2.2. For a given surrogate loss function ` and true cost vector distribution Pc, the conditional calibration function δ̂`(·;Pc) is defined, for > 0, by
δ̂`( ;Pc) := inf ĉ∈Rd
{ E [`(ĉ, c)]− inf
c′ E [`(c′, c)] : E [`SPO(ĉ, c)]− inf c′ E [`SPO(c′, c)] ≥
} .
Moreover, given a class of joint distributions P , with a slight abuse of notation, the calibration function δ̂`(·;P) is defined, for > 0, by
δ̂`( ;P) := inf x∈X ,P∈P δ̂`( ;P(·|x)).
If the calibration function δ̂`(·;P) satisfies δ̂`( ;P) > 0 for all > 0, then the loss function ` is uniformly `SPO-calibrated with respect to the class of distributions P . To obtain an excess risk bound, we let δ∗∗` denote the biconjugate, the largest convex lower semi-continuous envelope, of δ`. Jensen’s inequality then readily yields δ∗∗` (RSPO(g,P)−R∗SPO(P)) ≤ R`(g,P)−R`(P), which implies that the excess surrogate risk R`(g,P)−R`(P) of a predictor g can be translated into an upper bound of the excess SPO risk RSPO(g,P)−R∗SPO(P). For example, the uniform calibration of the least squares (squared `2) loss, namely `LS(ĉ, c) = ‖ĉ − c‖22, was examined by Ho-Nguyen and Kılınç-Karzan [2020]. They proved that the calibration function is δ`LS( ) =
2/D2S , which implies an upper bound of the excess SPO risk by RSPO(g,P)−R∗SPO(P) ≤ DS(RLS(g,P)−R∗LS(P))1/2. In this paper, we derive the calibration function of the SPO+ loss and thus reveal the quantitative relationship between the excess SPO risk and the excess surrogate SPO+ risk under different circumstances.

heading: 3 Risk bounds and calibration for polyhedral sets
text: In this section, we consider the case when the feasible region S is a bounded polyhedron and derive the calibration function of the SPO+ loss function. As is shown in Theorem 2.1, the SPO+ loss is Fisher consistent when the conditional distribution P(·|x) is continuous on all of Rd and is centrally symmetric about its mean c̄. More formally, the joint distribution P lies in the distribution class Pcont, symm := {P : P(·|x) is continuous on all of Rd and is centrally symmetric about its mean, for all x ∈ X}. In Example 1, we later show that this distribution class is not restrictive enough to obtain a meaningful calibration function. Instead, we consider a more specific distribution class consisting of distributions whose density functions can be lower bounded by a normal distribution. More formally, for given parameters M ≥ 1 and α, β > 0, define PM,α,β := {P ∈ Pcont, symm : for all x ∈ X with c̄ = E[c|x], there exists σ ∈ [0,M ] satisfying ‖c̄‖2 ≤ βσ and P(c|x) ≥ α · N (c̄, σ2I) for all c ∈ Rd}. Intuitively, the assumptions on the distribution class PM,α,β ensure that we avoid a situation where the density of the cost vector concentrates around some “badly behaved points.” This intuition is further highlighted in Example 1. Theorem 3.1 is our main result in the polyhedral case and demonstrates that the previously defined distribution class is a sufficient class to obtain a positive calibration function. Recall that DS denotes the diameter of S and define a “width constant” associated with S by dS := minv∈Rd:‖v‖2=1 { maxw∈S v Tw −minw∈S vTw }
. Notice that dS > 0 whenever S has a non-empty interior. Theorem 3.1. Suppose that the feasible region S is a bounded polyhedron and define ΞS := (1 + 2
√ 3DS dS )1−d. Then the calibration function of the SPO+ loss satisfies
δ̂`SPO+( ;PM,α,β) ≥ αΞS
4 √ 2πe 3(1+β2) 2
·min { 2
DSM ,
} for all > 0. (4)
Theorem 3.1 yields an O( 2) uniform calibration result for the distribution class PM,α,β . The dependence on the constants is also natural as it matches the upper bound given by the cases with a `1-like unit ball feasible region S and standard multivariate normal distribution as the conditional probability P(·|x). Please refer to Example 2 in the Appendix for a detailed discussion. Let us now provide some more intuition on the parameters involved in the definition of the distribution class PM,α,β and their roles in Theorem 3.1. In the definition of PM,α,β , α is a lower bound on the ratio of the density of the distribution of the cost vector relative to a “reference” standard normal distribution. When α is larger, the distribution is behaved more like a normal distribution and it leads to a better lower bound on the calibration function (4) in Theorem 3.1. The parameter M is an upper bound on the standard deviation of the aforementioned reference normal distribution, and the lower bound (4) naturally becomes worse as M increases. The parameter β measures how the conditional mean deviates from zero relative to the standard deviation of the reference normal distribution. If this distance is larger then the predictions are larger on average and (4) becomes worse. The width constant dS measures the near-degeneracy of the polyhedron (dS = 0 is degenerate) and the bound becomes meaningless as dS → 0. When the feasible region S is near-degenerate, i.e., the ratio dSDS is close to zero, we tend to have a weaker lower bound on the calibration function, which is also natural.
We now state a remark concerning an extension of Theorem 3.1 and we describe the example that demonstrates that it is not sufficient to consider the more general distribution class Pcont, symm. Remark. In Theorem 3.1, we assume that the conditional distribution P(·|x) is lower bounded by a normal density on the entire space Rd. We can extend the result of Theorem 3.1 to the case when
P(·|x) is lower bounded by a normal density on a bounded set but the constant is more involved. For details, please refer to Theorem B.1 in the Appendix. Example 1. Let the feasible region be the `1 ball S = {w ∈ R2 : ‖w‖1 ≤ 1} and consider the distribution class Pcont, symm. For a fixed scalar > 0, let c1 = (9 , 0)T and c2 = (−7 , 0)T . Let the conditional distribution Pσ(c|x) be a mixture of Gaussians defined by Pσ(c|x) := 12 (N (c1, σ 2I) + N (c2, σ2I)) for any σ > 0, and we have Pσ(c|x) ∈ Pcont, symm. Let the predicted cost vector be ĉ = (0, )T , then the excess conditional SPO risk is E[`SPO(ĉ, c)− `SPO(c̄, c)|x] = . Then it holds that the excess conditional SPO+ risk E [`SPO+(ĉ, c)− `SPO+(c̄, c)|x] → 0 when σ → 0, and hence we have δ̂`( ;Pcont, symm) = 0.
The intuition of Example 1 is that the existence of a non-zero calibration function requires the conditional distribution of c given x to be “uniform” on the space Rd, but not concentrate near certain points. Example 1 highlights a situation that considers one such “badly behaved” case where a limiting distribution of a mixture of two Gaussians leads to a zero calibration function.
By combining Theorem 3.1 with a generalization bound for the SPO+ loss, we can develop a sample complexity bound with respect to the SPO loss. Corollary 3.1 below presents such a result for the SPO+ method with a polyhedral feasible region. The derivation of Corollary 3.1 relies on the notion of multivariate Rademacher complexity as well as the vector contraction inequality of Maurer [2016] in the `2-norm. In particular, for a hypothesis class H of cost vector predictor functions (functions from X to Rd), the multivariate Rademacher complexity is defined as Rn(H) = Eσ,x [ supg∈H 1 n ∑n i=1 σ T i g(xi) ] , where σi ∈ {−1,+1}d are Rademacher random vectors for i = 1, . . . , n. Please refer to Appendix A for a detailed discussion of multivariate Rademacher complexity and the derivation of Corollary 3.1. Corollary 3.1. Suppose that the feasible region S is a bounded polyhedron, the optimal predictor g∗(x) = E[c|x] is in the hypothesis classH, and there exists a constant C ′ such that Rn(H) ≤ C
′ √ n .
Let ĝnSPO+ denote the predictor which minimizes the empirical SPO+ risk R̂ n SPO+(·) over H. Then there exists a constant C such that for any P ∈ PM,α,β and δ ∈ (0, 12 ), with probability at least 1− δ, it holds that
RSPO(ĝ n SPO+;P)−R∗SPO(P) ≤
C √ log(1/δ)
n1/4 .
Notice that the rate in Corollary 3.1 is O(1/n1/4). However, the bound is with respect to the SPO loss which is generally non-convex and is the first such bound for the SPO+ surrogate. Hu et al. [2020] present a similar result for the squared `2 surrogate with a rate of O(1/ √ n), and an interesting open question concerns whether the rate can also be improved for the SPO+ surrogate.

heading: 4 Risk bounds and calibration for strongly convex level sets
text: In this section, we develop improved risk bounds for the SPO+ loss function under the assumption that the feasible region is the level set of a strongly convex and smooth function, formalized in Assumption 4.1 below. Assumption 4.1. For a given norm ‖ · ‖, let f : Rd → R be a µ-strongly convex and L-smooth function for some L ≥ µ > 0. Assume that the feasible region S is defined by S = {w ∈ Rd : f(w) ≤ r} for some constant r > fmin := minw f(w).
The results in this section actually hold in more general situations than Assumption 4.1. In Appendix C, we extend the results of this section to allow the domain of the strongly convex and smooth function in Assumption 4.1 to be any set defined by linear equalities and convex inequalities (Assumption C.1). Herein, we consider the simplified case where the domain is Rd for ease of exposition and conciseness. The results of this section and the extension developed in Appendix C allow for a broad choice of feasible regions, for instance, any bounded `q ball for any q ∈ (1, 2] and the probability simplex with entropy constraint. The latter example, which can also be thought of as portfolio allocation with an entropy constraint, is considered in the experiments in Section 5.
As in the polyhedral case, the distribution class Pcont, symm is not restrictive enough to derive a meaningful lower bound on the calibration function of the SPO+ loss. We instead consider two related classes of rotationally symmetric distributions with bounded conditional coefficient of variation.
These distribution classes are formally defined in Definition 4.1 below, and include the multi-variate Gaussian, Laplace, and Cauchy distributions as special cases. Definition 4.1. Let A be a given positive definite matrix. We define Prot symm,A as the class of distributions with conditional rotational symmetry in the norm ‖ · ‖A−1 , namely
Prot symm,A := {P : ∀x ∈ X ,∃q(·) : [0,∞]→ [0,∞] such that P(c|x) = q(‖c− c̄‖A−1)}.
Let c̄ denote the conditional expectation c̄ = E[c|x]. For constants α ∈ (0, 1] and β > 0, define Pβ,A := { P ∈ Prot symm,A : Ec|x[‖c− c̄‖2A−1 ] ≤ β 2 · ‖c̄‖2A−1 , ∀x ∈ X } ,
and Pα,β,A := { P ∈ Prot symm,A : Pc|x(‖c− c̄‖A−1 ≤ β · ‖c̄‖A−1) ≥ α, ∀x ∈ X } .
Under the above assumptions, Theorem 4.1 demonstrates that the calibration function of the SPO+ loss is O( ), significantly strengthening our result in the polyhedral case. Theorem C.1 in Appendix C extends the result of Theorem 4.1 to the aforementioned case where the domain of f(·) may be a subset of Rd, which includes the entropy constrained portfolio allocation problem for example. Theorem 4.1. Suppose that Assumption 4.1 holds with respect to the norm ‖ · ‖A for some positive definite matrix A. Then, for any > 0, it holds that δ̂`SPO+( ;Pβ,A) ≥ µ9/2 4(1+β2)L9/2 · and δ̂`SPO+( ;Pα,β,A) ≥ αµ9/2 4(1+β2)L9/2 · .
Let us now provide some more intuition on the parameters involved in the definitions of the distribution classes Pβ,A and Pα,β,A and their roles in Theorem 4.1. In both cases, β controls the concentration of the distribution of cost vector around the conditional mean. The more concentrated the distribution is, the better the bounds in Theorem 4.1 are. In the case of Pα,β,A, α relates to the probability that the cost vector is “relatively close” to the conditional mean. When α is larger, the cost vector is more likely to be close to the conditional mean and the bound will be better.
Our analysis for the calibration function (the proof of Theorem 4.1) relies on the following two lemmas, which utilize the special structure of strongly convex level sets to strengthen the first-order optimality guarantees and derive a “Lipschitz-like” continuity property of the optimization oracle w∗(·). The first such lemma strengthens the optimality guarantees of (2) and provides both upper and lower bounds of the SPO loss. Lemma 4.1. Suppose that Assumption 4.1 holds with respect to a generic norm ‖ · ‖. Then, for any c1, c2 ∈ Rd, it holds that
cT1 (w − w∗(c1)) ≥ µ 2 √ 2L(r − fmin) ‖c1‖∗‖w − w∗(c1)‖2, ∀w ∈ S,
and cT1 (w ∗(c2)− w∗(c1)) ≤ L 2 √ 2µ(r − fmin) ‖c1‖∗‖w∗(c1)− w∗(c2)‖2.
The following lemma builds on Lemma 4.1 to develop upper and lower bounds on the difference between two optimal decisions based on the difference between the two normalized cost vectors. Lemma 4.2. Suppose that Assumption 4.1 holds with respect to a generic norm ‖ · ‖. Let c1, c2 ∈ Rd be such that c1, c2 6= 0, then it holds that
‖w∗(c1)− w∗(c2)‖ ≥ √
2µ(r − fmin) L · ∥∥∥∥ c1‖c1‖∗ − c2‖c2‖∗ ∥∥∥∥ ∗ ,
and
‖w∗(c1)− w∗(c2)‖ ≤ √
2L(r − fmin) µ · ∥∥∥∥ c1‖c1‖∗ − c2‖c2‖∗ ∥∥∥∥ ∗ .
Note that the lower bound of cT1 (w−w∗(c1)) in Lemma 4.1 and the upper bound of ‖w∗(c1)−w∗(c2)‖ in Lemma 4.2 match bounds developed by El Balghiti et al. [2019]. Indeed, although El Balghiti et al. [2019] study the more general case of strongly convex sets, the constants are the same since Theorem 12 of Journée et al. [2010] demonstrates that our set S is a µ√
2L(r−fmin) -strongly convex
set. However, the upper bounds in Lemmas 4.1 and 4.2 appear to be novel and rely on the special properties of strongly convex level sets. It is important to emphasize that we generally do not expect all of the bounds in Lemmas 4.1 and 4.2 to holds for polyhedral sets. Indeed, for a polyhedron the optimization oracle w∗(·) is generally discontinuous at cost vectors that have multiple optimal solutions. The properties in Lemmas 4.1 and 4.2 drive the proof of Theorem 4.1 and hence lead to the improvement from O( 2) in the polyhedral case to O( ) in the strongly convex level set case.
By following similar arguments as in the derivation of Corlloary 3.1, Corollary 4.1 presents the sample complexity of the SPO+ method when the feasible region is a strongly convex level set. Corollary 4.1. Suppose that Assumption 4.1 holds with respect to the norm ‖ · ‖A for some positive definite matrix A. Suppose further that the optimal predictor g∗(x) = E[c|x] is in the hypothesis classH, and there exists a constant C ′ such that Rn(H) ≤ C
′ √ n . Let ĝnSPO+ denote the predictor which
minimizes the empirical SPO+ risk R̂nSPO+(·) over H. Then there exists a constant C such that for any P ∈ Pα,β ∪ Pβ and δ ∈ (0, 12 ), with probability at least 1− δ, it holds that
RSPO(ĝ n SPO+;P)−R∗SPO(P) ≤
C √ log(1/δ)
n1/2 .
Notice that Corollary 4.1 improves the rate of convergence toO(1/ √ n) as compared to theO(1/n1/4) rate of Corollary 3.1. This matches the rate for the squared `2 surrogate developed by Hu et al. [2020] (though their result is in the polyhedral case).

heading: 5 Computational experiments
text: In this section, we present computational results of synthetic dataset experiments wherein we empirically examine the performance of the SPO+ loss function for training prediction models, using portfolio allocation and cost-sensitive multi-class classification problems as our problem classes. We focus on two classes of prediction models: (i) linear models, and (ii) two-layer neural networks with 256 neurons in the hidden layer. We compare the performance of the empirical minimizer of the following four different loss function: (i) the previously described SPO loss function (when applicable), (ii) the previously described SPO+ loss function, (iii) the least squares (squared `2) loss function `(ĉ, c) = ‖ĉ− c‖22, and (iv) the absolute (`1) loss function `(ĉ, c) = ‖ĉ− c‖1. For all loss functions, we use the Adam method of Kingma and Ba [2015] to train the parameters of the prediction models. Note that the loss functions (iii) and (iv) do not utilize the structure of the feasible region S and can be viewed as purely learning the relationship between cost and feature vectors.
Entropy constrained portfolio allocation. First, we consider the portfolio allocation [Markowitz, 1952] problem with entropy constraint, where the goal is to pick an allocation of assets in order to maximize the expected return while enforcing a certain level of diversity through the use of an entropy constraint (see, e.g., Bera and Park [2008]). This application is an instance of our more general theory for strongly convex level sets on constrained domains, developed in Appendix C. Alternative formulations of portfolio allocation, including when S is a polyhedron or a polyhedron intersected with an ellipsoid, have been empirically studied in previous works (see, for example Elmachtoub and Grigas [2021], Hu et al. [2020]). The objective is to minimize cTw where c is the negative of the expected returns of d different assets, and the feasible region is S = {w ∈ Rd : w ≥ 0, ∑d i=1 wi =
1, ∑d i=1 wi logwi ≤ r} where r is a user-specified threshold of the entropy of portfolio w. Note that, due to the differentiability properties of of the optimization oracle w∗(·) in this case (see Lemma D.1 in the Appendix), it is possible to (at least locally) optimize the SPO loss using a gradient method even though SPO loss is not convex.
In our simulations, the relationship between the true cost vector c and its auxiliary feature vector x is given by c = φdeg(Bx) , where φdeg is a polynomial kernel mapping of degree deg, B is a fixed weight matrix, and is a multiplicative noise term. The features are generated from a standard multivariate normal distribution, we consider d = 50 assets, and further details of the synthetic data generation process are provided in Appendix D. To account for the differing distributions of the magnitude of the cost vectors, in order to evaluate the performance of each method we compute a “normalized” SPO loss on the test set. Specifically, let ĝ denote a trained prediction model and let {x̃i, c̃i}mi=1 denote the test set, then the normalized SPO loss is defined as ∑m i=1 `SPO(ĝ(x̃i),c̃i)∑m
i=1 z ∗(c̃i)
, where
z∗(c̃) := minw∈S c̃ Tw is the optimal cost in hindsight. We set the size of the test set to 10000. In all of our experiments, we run 50 independent trials for each setting of parameters.
Figure 1 displays the empirical performance of each method. We observe that with a linear hypothesis class, for smaller values of the degree parameters, i.e., deg ∈ {1, 2}, all four methods perform comparably, while the SPO and SPO+ methods dominate in cases with larger values of the degree parameters. With a neural net hypothesis class, we observe a similar pattern but, due to the added degree of flexibility, the SPO method dominates the cases with larger values of degree and SPO+ method is the best among all surrogate loss functions. The better results of the `1 loss as compared to the squared `2 loss might be explained by robustness against outliers. Appendix D also contains results showing the observed convergence of the excess SPO risk, in the case of polynomial degree one, for both this experiment and the cost-sensitive multi-class classification case.
Cost-sensitive multi-class classification. Here we consider the cost-sensitive multi-class classification problem. Since this is a multi-class classification problem, the feasible region is simply the unit simplex S = {w ∈ Rd : w ≥ 0, ∑d i=1 wi = 1} We consider an alternative model for generating the data, whereby the relationship between the true cost vector c and its auxiliary feature vector x is as follows: first we generate a score s = σ(φdeg(bTx) ), where φdeg is a degree-deg polynomial kernel, b is a fixed weight vector, is a multiplicative noise term, and σ(·) is the logistic function. Then the true label is given by lab = d10se ∈ {1, . . . , 10} and the cost vector c is given by ci = |i− lab| for i = 1, . . . , 10. The features are generated from a standard multivariate normal distribution, and further details of the synthetic data generation process are provided in Appendix D. Since the scale of the cost vectors do not change as we change different parameters, we simply compare the test set SPO loss for each method. We still set the size of the test set to 10000 and we run 50 independent trials for each setting of parameters. In addition to the regular SPO+, least squares, and absolute losses, we consider an alternative surrogate loss constructed by considering the SPO+ loss using a log barrier (strongly convex) approximation to the unit simplex. That is, we consider the SPO+ surrogate that arises from the set S̃ := {w ∈ Rd : w ≥ 0, ∑d i=1 wi = 1,− ∑d i=1 logwi ≤ r} for some r > 0. Details about how we chose the value of r are provided in Appendix D.
Herein we focus on the comparison between the standard SPO+ loss and the “SPO+ w/ Barrier” surrogate loss. (We include a more complete comparison of all the method akin to Figure 1 in Appendix D.) Figure 2 shows a detailed comparison between these alternative SPO+ surrogates as we vary the training set size. Note that the SPO loss is always measured with respect to the standard unit simplex and not the log barrier approximation. Interestingly, we observe that the “SPO+ w/ Barrier” surrogate tends to perform better than the regular SPO+ surrogate when the training set size is small, whereas the regular SPO+ surrogate gradually performs better as the training set size increases. These
results suggest that adding a barrier constraint to the feasible region has a type of regularization effect, which may also be explained by our theoretical results. Indeed, adding the barrier constraint makes the feasible region strongly convex, which improves the rate of convergence of the SPO risk. On the other hand, this results in an approximation to the actual feasible region of interest and, eventually for large enough training set sizes, the regularizing benefit of the barrier constraint is outweighed by the cost of this approximation.

heading: 6 Conclusions and future directions
text: Our work develops risk bounds and uniform calibration results for the SPO+ loss relative to the SPO loss, and our results provide a quantitative way to transfer the excess surrogate risk to excess true risk. We analyze the case with a polyhedral feasible region of the underlying optimization problem, and we strengthen the results when the feasible region is a level set of a strongly convex function. There are several intriguing future directions. In this work, we mainly focus on the worst case risk bounds of the SPO+ method under different types of feasible regions, and we consider the minimal possible conditions to guarantee the risk bounds. In many practical problems, it is reasonable to assume the true joint distribution of (x, c) satisfies certain “low-noise” conditions (see, for example, Bartlett et al. [2006], Massart et al. [2006], Hu et al. [2020]). In such conditions, one might be able to obtain improved risk bounds and sample complexities. Also, developing tractable surrogates and a corresponding calibration theory for non-linear objectives is very worthwhile.

heading: Acknowledgments
text: The authors are grateful to Othman El Balghiti, Adam N. Elmachtoub, and Ambuj Tewari for early discussions related to this work. PG acknowledges the support of NSF Awards CCF-1755705 and CMMI-1762744.

Review summaries:
The paper studies the excess SPO risk of the empirical SPO+ risk minimizers. When the feasible region of the downstream planning task is a bounded polyhedron, the excess risk is on the order of O(n^{-1/4}). When the feasible region is the level set of a strongly convex and smooth function, the excess risk is on the order of O(n^{-1/2}). Empirically, the authors show that the empirical minimizers of SPO+ perform better than the ones for \ell_1, \ell_2 for more complex hypothesis class. 

This paper derives uniform calibration bounds for the SPO+ loss relative to the SPO loss, when the feasible region is polyhedral or a level set of strongly convex functions. This enables us to translate excessive risk bound in terms of SPO+ loss, ie, the convex surrogate loss,  into bounds for the SPO loss. Notably, the authors show that the class of symmetric distributions considered in Elmachtoub and Grigas is not restrictive enough to derive calibration bounds, and instead stronger restrictions must be incorporated.  

The paper develops the risk bounds and uniform calibration results for the surrogate SPO+ loss relative to the SPO loss in the predict-then-optimize framework under the assumption of a linear loss function. The results are first obtained for a polyhedron feasible region of the optimization problem and further improved when the feasible region is a level set of a strongly convex function. Then, they are used to establish a sample complexity (with respect to the SPO loss) of the algorithm that minimizes the empirical SPO+ risk. The empirical performance of the SPO+ surrogate is demonstrated in portfolio optimization and cost-sensitive multi-class classification problems.