Gradient boosting for convex cone predict and optimize problems

Many problems in engineering and statistics involve both predictive forecasting and decision-based optimization. Traditionally, predictive models are optimized independently from the final decision-based optimization problem. In contrast, a `smart, predict then optimize' (SPO) framework optimizes prediction models to explicitly minimize the final downstream decision loss. In this paper we present dboost, a gradient boosting algorithm for training prediction model ensembles to minimize decision regret. The dboost framework supports any convex optimization program that can be cast as convex quadratic cone program and gradient boosting is performed by implicit differentiation of a custom fixed-point mapping. To our knowledge, the dboost framework is the first general purpose implementation of gradient boosting to predict and optimize problems. Experimental results comparing with state-of-the-art SPO methods show that dboost can further reduce out-of-sample decision regret.


Introduction
Recently there has been a growing body of research on decision-aware predictive modelling (see for example [5,7,6,20,21,24,29,37]). A naive 'predict, then optimize' framework treats the prediction estimation problem independently from its use in the downstream decision optimization problem. As such, an 'objective mismatch' [28] can occur whereby improved prediction accuracy does not necessarily result in smaller decision error.
Conversely, the Smart 'Predict, then Optimize' (SPO) [20] framework optimizes prediction models in order to minimize the final downstream decision regret. To date, the SPO framework has been studied in a general setting for linear and decision-tree regression models [20,21]. In this paper we present dboost, a general purpose framework that combines the strength of gradient boosting ensembles with the SPO framework. Previous work [27] considers gradient boosting for integrated prediction and optimization problems but only considers a small subset of optimization problems with linear inequality constraints. In contrast, the dboost framework is capable of supporting any optimization problem that can be cast as a convex quadratic cone program (CQP); and thus supports linear, quadratic and second-order cone programming with general convex cone constraints. As a secondary technical contribution, we present a novel fixed-point implicit differentiation algorithm for efficiently computing the gradient of the SPO loss with respect to all of the cone program variables. Thirdly, the dboost framework is provided as an open-source R package, available here: https://github.com/adsb85/dboost The remainder of the paper is outlined as follows. We begin with a brief overview of related work on integrated prediction and optimization. In Section 2 we present convex quadratic cone programming in the context of the the SPO framework and provide the fixed-point implicit differentiation algorithm. We then present the dboost framework as a general extension of the gradient boosting algorithm proposed by Friedman [22]. Experimental results are provided in Section 3 and demonstrate that training prediction models with dboost can reduce decision regret by anywhere from 15%-90% in comparison to existing solutions. We conclude with a summary of results and a discussion of future work.

Related work
The dboost framework applies gradient boosting, as proposed by Friedman [22] to the SPO loss function, proposed by Elmachtoub and Grigas [20], and described in more detail in Section 2. In general, optimizing the SPO loss by gradient descent methods is challenging as it requires computing the gradient of the optimal solution with respect to predicted costs. To date, several local gradient based methods [29,33] and convex approximations [16,20] have been proposed and have proven to be effective in comparison to naive 'predict, then optimize'.
Related, the work of Agrawal et al. [2] and Amos and Kolter [3] provide general frameworks for learning convex optimization programs in an endto-end trainable neural network. Backpropagation is performed by implicit differentiation of specialized convex optimization layers. The development of efficient and modular convex optimization layers has been fundamental to the advancement of end-to-end integration of prediction and optimization [1,9]. Growing empirical evidence advocates for a fully integrated approach (see for example [4,14,16,19,29,33,36,37]).
Most relevant to our work is the work Elmachtoub et al. [21] who present SPO trees (SPOT) for optimizing classification and regression trees (CART) [11] to minimize downstream decision regret. The authors also consider a random forest [10] implementation, but do not consider a gradient boosting ensemble approach. Konishi and Fukunaga [27] consider gradient boosting problems under the SPO framework for a subset of optimization problems; namely with linear inequality constraints, but do not consider more general convex optimization problems. To our knowledge, dboost is the first 'smart' gradient boosting implementation that can support any convex optimization program that can be cast as convex quadratic cone program. A necessary requirement of our framework is the ability to compute gradients of optimal decisions with respect to predicted costs. In the following section we present the dboost framework and provide a novel fixed-point implicit differentiation algorithm for computing said gradients.

Methodology
We consider convex quadratic cone programs (QCP) with the following primal-dual form [32]: where z ∈ R dz , y ∈ R dy and s ∈ R dy denote the primal, dual and slack variables, respectively. The objective input variables are a symmetric positive semi-definite matrix P ∈ R dz×dz and cost c ∈ R dz . The feasible region is defined by A ∈ R dy×dz , b ∈ R dy and K, a nonempty convex cone with associated dual cone K * . We denote the optimal solution to Program (1) as ζ * = (z * , y * , s * ) ∈ Z ⊂ R dz ×K×K * with feasible region: In practice, the true cost c is not directly observable at decision time, and rather an associated feature vector x ∈ R dx is observed. Given a training data sets D = {(x (i) , c (i) )} m i=1 we seek to estimate a prediction model f : R dx × R d θ → R dz such that: with optimal prediction model parameter θ * given by: In this paper we consider 'additive' prediction models of the form: with parameter θ = {(β n , α n )} N n=1 and nonnegative weight β n ≥ 0∀n . In particular, we focus on the case where each of the functions h(x, α m ) is a regression tree and therefore the parameters α n encode the feature component and splitting information.
A naive 'predict, then optimize' approach would first estimate θ * by minimizing a prediction loss (such as least-squares) and then simply 'plug-in' the estimateĉ (i) to Program (1) in order to retrieve the 'optimal' decision. In contrast, a 'smart predict, then optimize' (SPO) approach estimates the prediction model parameters by minimizing the decision regret: where z * (ĉ) = argmin We conclude this section by noting that the SPO estimation process can be posed as a bi-linear optimization: Solving Program (8) to global optimality is challenging as the problem is generally non-convex. As such, we choose to solve Program (8) to local optimality by applying functional gradient descent [22]. Specifically, at each iteration of gradient descent we must first solve m QCPs to optimality. The gradient, ∇ĉ (i) QSPO , is then obtained by specialized argmin differentiation [1, 2, 3], described in more detail below.

Fixed-point argmin differentiation
Program (1) is solved by applying a Douglas-Rachford splitting to a homogeneous embedding of the QCP as described in O'Donoghue [32]. Specifically, we define the convex cone C = R dz × K * and denote: where I z denotes the identity matrix of same dimension as z. O'Donoghue [32] demonstrate that a direct application of operator splitting produces the following procedure; from any initial u 0 and v 0 then the following iterations converge to the optimal ζ * (if it exists): where Π C denotes the euclidean projection operator onto the set C.
We follow Butler and Kwon [13] and recast the iterative procedure (10) as a fixed-point iteration over variable w k =ũ k+1 − v k . We then apply the implicit function theorem [18] to the fixed-point iteration in order to derive an efficient scheme for computing the required gradients. We begin with the following proposition. Note that all proofs are available in the Appendix.
Then the iterations in Equation (10) can be cast as a fixed-point iteration of the form F (w, θ) = w given by: We follow Busseti et al. [12] and denote the derivative of the projection operator as DΠ C . Note for ease of notation we drop the index k. The Jacobian, ∇ w F , is therefore defined as: (12) Direct application of the implicit function theorem gives the desired Jacboian ∇ θ w(θ), with respect to the parameter θ at optimality w * : From the definition of w * we have that u * = Π C (w * ) and therefore ∇ w * u * = DΠ C (w * ). In practice, it is inefficient to from the Jacobian ∇ θ u * (θ) directly and instead we compute the left matrix-vector product of the Jacobian with respect to the relevant gradient, ∂ /∂ z * , as outlined below.
Proposition 2. Denote the matrix G as: Letd z andd y be defined as: Then the gradients of the loss function, , with respect to problem variables P, c, A and b are given by: Proposition 2 provides a framework for computing the gradient of any arbitrary loss function with respect to the cost variable c by implicit differentiation of the fixed-point mapping of the transformed operator splitting iterations. We now present the dboost algorithm which applies gradient boosting to generate a boosted prediction model ensemble optimized for decision regret.

Motivating example
We are motivated by the work of [20,21] who demonstrate that optimizing prediction model parameters to minimize decision regret produces prediction models with lower complexity and improved out-of-sample performance. Below, we provide a motivating example to help illustrate the behaviour of dboost in comparison to a naive gradient boosting model trained to minimize prediction meansquare error (MSE). We seek to maximize the return (minimize the cost) of a portfolio of two assets subject to the linear constraints: where z 1 and z 2 denote the proportion of capital allocated to asset 1 and 2, respectively. Here, the norm penalty z 2 2 enforces portfolio diversification when differences in returns are small (see for example [15,17]). We generate a dataset of 500 observations where the feature, x, is drawn from the standard uniform distribution and the return of each asset is given by r 1 = x + and r 2 = x +sin(3 x) + where ∼ N (0, 0.10).    (18)) at each iteration of gradient boosting. Recall, that the naive model is unaware as to how the predictions will ultimately be used in the context of the final downstream optimization model, and as such continues to add trees to the ensemble in an effort to minimize the prediction MSE. In contrast, the dboost model explicitly minimizes decision regret; the boosting algorithm terminates after 5 iterations with equal decision accuracy. Figure 2 plots the predictions model forecasts as a function of the feature variable x. In this example there are two decision boundaries: x ≈ 1 and x ≈ 2, where the optimal decision changes from favouring one asset over the other. The dboost prediction model correctly identifies these approximate decision boundaries with substantially fewer trees and therefore avoids having to overfit the training dataset.

Experiments
We present several experiments that compare the out-of-sample decision cost of dboost to 5 alternative methods: 1. CART: naive 'predict, then optimize' with prediction model given by the classification and regression tree as proposed by Breiman et al. [11]. Trees are optimized to minimize prediction MSE.
2. CART Forest: A random forest [10] implementation of CART.
4. SPOT Forest: random forest implementation of SPOT.
5. MSE Boosting: naive 'predict, then optimize' with prediction model given by gradient boosting as proposed by Friedman [22].
Experiments are conducted on an Apple Mac Pro computer (2.7 GHz 12-Core Intel Xeon E5,128 GB 1066 MHz DDR3 RAM) running macOS 'Catalina'. All computations are run on an unloaded, singlethreaded CPU. The software was written using the R programming language (version 4.0.0). All ensemble based approaches use a maximum of 100 trees. Random forest implementations perform bagging across feature variables and training data observations with a 50% sampling rate. All models consider maximum tree depths: {0, 1, 2}. All experiments are evaluated over 10 independent trials. Each trial consists of a randomly generated training dataset D = {(x (i) , c (i) )} m i=1 with m = 1000 observations and an out-of-sample training dataset with m = 1000. Performance is evaluated with respect to the total excess decision cost, given by: Synthetic data is generated as follows. Feature variables, x, are generated by randomly drawing from the Uniform distribution U (−1, 1). The cost values, c, are generated according to the polynomial model: with intercept H 0 and regression coefficients H j ∈ R dz×dx . The parameter p controls the polynomial degree. Regression coefficients are sparse with each element of H j having a 50% probability of being 0 and 50% probability of being nonzero U(−1, 1).
We let ∼ N (0, 1) and the scalar value τ controls the amount of noise in the data. All experiments consider three noise levels: τ ∈ {0.0, 0.5, 1.0}.We present three optimization problems: 1. Noisy network-flow. We consider a continuous network flow problem over a directed graph with 5 nodes. Network edges are randomly generated with the probability that node i flows to node j given by P r(i → j) = 0.75 |i−j−1| . The number of decision variables, d z , is determined by the number of edges. Cost values are generated according to (19) with H 0 ∼ N (− 1, 1), p = 3 and d x = 5. We augment the standard network-flow linear program with L 2 norm regularization to regulate flow across the network. The optimization problem is given by: 2. Noisy quadratic program. We consider an equality constrained quadratic program of the form: with number of decision variables d z = 25 and randomly generated constraint matrix A ∈ R 3×dz where P r(A jk = 0) = P r(A jk = 1) = 0.50. The vector b is chosen to guarantee that the problem is non-empty. Cost values are generated according to (19) with H 0 = 0, p = 3 and d x = 5. The positive definite matrix, P, is subject to estimation error: P = P +0.1 Ξ where Ξ = 1 n L T L, L ∈ R 10dz×dz and entries of L ∼ N (0, 1). Program (22) occurs in many applications to statistics [34,35], machinelearning [23,25], signal-processing [26] and finance [8,31]. 3. Noisy portfolio optimization. We consider the Markowitz [30] mean-variance optimization program: where z j denotes the proportion of capital allocated to asset j. The cost values, c, are the negative asset returns and are generated according to (19) with H 0 =∼ N (0, 1), p = 3 and d x = 5. The covariance matrix, V, is generated according to a linear factor model V = L T L +ε I z , with ε = 0.01 and factor matrix L ∈ R 4×dz with entries L ∼ U (−1, 1). The second-order cone constraint limits the maximum risk of the portfolio. In each trial we set Figure 3 reports the out-of-sample excess cost for the noisy network flow problem with noise level τ ∈ {0.0, 0.5, 1.0}. First, observe that for all noise levels, the SPOT models produce a smaller excess cost in comparison to the CART counterparts. In general, the random forest implementations provide a further marginal reduction in excess cost. Gradient boosting models produce excess costs that are on average 50% lower than the corresponding CART and SPOT costs. In all cases, however, the dboost model most effectively minimizes the decision cost, further reducing the excess cost of the MSE boosting model by 30%-60%, on average. We observe that the reduction in excess cost is greatest when the noise level is low and in general, deeper tree models tend to produce smaller excess costs. Figure 4 reports the out-of-sample excess cost for the noisy quadratic program problem with noise level τ ∈ {0.0, 0.5, 1.0}. Here we note that the CART models produce smaller excess cost in comparison to the SPOT counterparts. Among the CART and SPOT models, the CART Forest produces the smallest excess cost across all noise levels and tree depths. Again, we observe that both gradient boosted models produce the smallest outof-sample costs, and in most cases, the excess costs are 25%-75% smaller than the corresponding CART and SPOT costs. When τ = 0, the average excess cost of the MSE Boosting model is approximately 22.2% whereas the dboost model produces an average excess cost of approximately 5.8%; a further 75% reduction in excess cost over the 'naive' gradient boosting model. In the majority of cases, the dboost model produces the smallest excess cost, however, as the noise level increases the difference in excess costs in comparison to MSE Boosting is less substantial. In fact, when τ = 1 and tree depth is greater than zero, the MSE Boosting model pro-duces marginally lower excess costs. Finally, observe that when the noise level is high (τ = 1) both gradient boosted models with tree depth 0 produce the smallest excess cost; an argument in favour of low model complexity when the training dataset is noisy.   Lastly, Figure 5 reports the out-of-sample excess cost for the noisy portfolio optimization program problem with noise level τ ∈ {0.0, 0.5, 1.0}. As in the network flow problem, we note that the SPOT models produce smaller excess cost in comparison to the CART counterparts. Among the CART and SPOT models, the SPOT Forest pro-duces the smallest excess cost across all noise levels and tree depths. Once again, we observe that both gradient boosted models produce excess costs that, in the most cases, are anywhere from 15%-90% smaller than the corresponding CART and SPOT costs. When τ = 0, we observe that the MSE Boosting model produces almost 0% excess cost;  the smallest excess cost across all methods. In contrast to the prior experiments, in this experiment, dboost provides no marginal benefit over 'naive' gradient boosting and we observe that as the noise level increases, the MSE Boosting and dboost produce very similar excess costs. However, across all noise levels, the MSE Boosting with tree depth 0 produces the smallest excess cost; again an argument in favour of prediction models with low model complexity.

Conclusion
In this paper, we presented dboost, a general purpose framework for building gradient boosting prediction models within a 'smart predict, then optimize' (SPO) framework. In contrast to prior work, dboost is capable of supporting any optimization problem that can be cast as a convex quadratic cone program (CQP). Prediction models are estimated by functional gradient descent, which ultimately requires computing the action of the Jacobian of optimal CQP solutions with respect to the prediction model parameters. We present a novel fixed-point implicit differentiation algorithm for computing the gradient of the SPO loss with respect to all of the CQP program variables. Our experimental results should be interpreted as a proof-of-concept and we acknowledge that further testing on alternative datasets or with different problem variable assumptions is required in order to better determine the efficacy of a 'smart' gradi-ent boosting approach. In this paper we considered three experiments: a noisy network flow, a noisy quadratic program and noisy portfolio optimization problem, with a variety of problem specifications. We compared dboost with CART and SPOT tree models, as well as a naive gradient boosting model optimized for prediction MSE. In the majority of cases, both gradient boosting prediction models produce anywhere from a 15% − 90% reduction in the out-of-sample excess costs in comparison to CART and SPOT models. Moreover, for the network flow and quadratic programming experiment, dboost produces the smallest decision costs, with excess costs that are 30%-75% smaller than the those produced by the naive MSE boosting model. In the portfolio optimization experiment, dboost provides no marginal benefit over 'naive' gradient boosting and in fact across all noise levels, the MSE Boosting with tree depth 0 produces the small excess cost. Nonetheless, in this experiment dboost remains competitive with MSE Boosting and outperforms both the CART and SPOT models.
We note that training gradient boosting models to minimize downstream decision regret in general is computationally demanding. Indeed, at each iteration of gradient boosting, we must solve m convex quadratic cone programs to obtain optimal decisions: {z * (ĉ (i) )} m i=1 , and then compute the resulting gradient of the SPO loss at optimality, both of which can present computational bottlenecks. The gradient, ∇ĉ (i) QSPO , is com-puted by implicit argmin differentiation and therefore requires solving m linear systems of dimension d z + d y , which for large scale problems can be computational expensive. All experiments herein were performed with relatively small training datasets (m = 1000) and small decision optimization problems (d z ≤ 25). Improving the performance and scalability of dboost to support large scale problems is therefore an important area of future research. Further experimentation on alternative optimization problems and datasets to better understand the character of 'smart' gradient boosting is another area of active research.

A.1 Proof of Proposition 1
The iterations in Equation (10) can be cast as a fixed-point iteration as follows. Let w k =ũ k+1 − v k , then: and : Substituting Equations (23) and (24) into Equation (10a) gives the desired fixed-point iteration: A.2 Proof of Proposition 2 We begin with the following definitions.
Definition 1. Let F : R dw × R d θ → R dw be a continuously differentiable function with variable w and parameter θ. We define w * as a fixed-point of F at (w * , θ) if: Definition 2. The residual map, G : R dw × R d θ → R dw of a fixed point, (w * , θ), of F is given by: G(w * , θ) = F (w * , θ) − w * = 0.
The implicit function theorem, as defined by Dontchev and Rockafellar [18], then provides the conditions on G for which the Jacobian of the solution mapping with respect to θ is well defined. Theorem 1. Let G : R dw × R d θ → R dw be a continuously differentiable function in a neighborhood of (w * , θ) such that G(w * , θ) = 0. Denote the non-singular partial Jacobian of G with respect to w as ∇ w G(w * , θ). Then w(θ) is an implicit function of θ and is continuously differentiable in a neighborhood, Θ, of θ with Jacobian: Corollary 1. Let F : R dw × R d θ → R dw be a continuously differentiable function with fixed-point (w * , θ). Then w(θ) is an implicit function of θ and is continuously differentiable in a neighborhood, Θ, of θ with Jacobian: We define F : R dz+dy × R d θ → R dz+dy as: The Jacobian, ∇ w F , is therefore defined as: