Optimization of Tree Ensembles

TL;DR

This paper develops a mixed-integer optimization framework for tuning controllable variables in tree ensemble models to maximize predictions, with scalable solution methods demonstrated on real-world case studies.

Contribution

It introduces a novel optimization formulation for tree ensembles, analyzes its strength, and proposes scalable algorithms for large-scale problems.

Findings

Efficiently solves large-scale tree ensemble optimization problems.

Outperforms heuristic methods in case studies.

Successfully applies to drug design and pricing scenarios.

Abstract

Tree ensemble models such as random forests and boosted trees are among the most widely used and practically successful predictive models in applied machine learning and business analytics. Although such models have been used to make predictions based on exogenous, uncontrollable independent variables, they are increasingly being used to make predictions where the independent variables are controllable and are also decision variables. In this paper, we study the problem of tree ensemble optimization: given a tree ensemble that predicts some dependent variable using controllable independent variables, how should we set these variables so as to maximize the predicted value? We formulate the problem as a mixed-integer optimization problem. We theoretically examine the strength of our formulation, provide a hierarchy of approximate formulations with bounds on approximation quality and exploit the structure of the problem to develop two large-scale solution methods, one based on Benders decomposition and one based on iteratively generating tree split constraints. We test our methodology on real data sets, including two case studies in drug design and customized pricing, and show that our methodology can efficiently solve large-scale instances to near or full optimality, and outperforms solutions obtained by heuristic approaches. In our drug design case, we show how our approach can identify compounds that efficiently trade-off predicted performance and novelty with respect to existing, known compounds. In our customized pricing case, we show how our approach can efficiently determine optimal store-level prices under a random forest model that delivers excellent predictive accuracy.

Tables

Table 1. Table 1: Summary of real data sets used in numerical experiments. Note: * = accessed via UCI Machine Learning Repository (Lichman 2013 ) ; ** = accessed via AppliedPredictiveModeling package in R (Kuhn and Johnson 2014 ) .

		Num.	Num.
Data set	Source	Vars.	Obs.	Description
winequalityred *	Cortez et al. (2009)	11	1599	Predict quality of (red) wine
concrete **	Yeh (1998)	8	1030	Predict strength of concrete
permeability **	Kansy et al. (1998)	1069	165	Predict permeability of compound
solubility **	Tetko et al. (2001),	228	951	Predict solubility of compound
	Huuskonen (2000)

Table 2. Table 2: Solution times for tree ensemble optimization experiment.

Data set	$T$	$N_{Levels}$	$N_{Leaves}$	${𝒯}_{MIO}$	${𝒯}_{StdLin,MIO}$	${𝒯}_{LS}$
solubility	10	1253	3157	0.1	215.2	0.2
	50	2844	15933	0.8	1800.3	1.8
	100	4129	31720	1.7	1801.8	8.8
	200	6016	63704	4.5	1877.8	33.7
	500	9646	159639	177.9	1800.3	147.0
permeability	10	2138	604	0.0	122.6	1.0
	50	2138	3056	0.1	1800.0	1.9
	100	2138	6108	0.2	1800.3	3.1
	200	2138	12214	0.5	1800.0	6.1
	500	2138	30443	2.7	1800.0	19.1
winequalityred	10	1370	3246	1.8	1800.1	0.0
	50	2490	16296	18.5	1800.1	0.6
	100	3000	32659	51.6	1800.1	2.5
	200	3495	65199	216.0	1800.2	11.4
	500	3981	162936	1159.7	1971.8	34.6
concrete	10	1924	2843	0.2	1800.8	0.1
	50	5614	14547	22.7	1800.1	1.3
	100	7851	29120	67.8	1800.1	4.3
	200	10459	58242	183.8	1800.2	20.2
	500	13988	145262	846.9	1809.4	81.6

Table 3. Table 3: Gaps for tree ensemble optimization experiment.

Data set	$T$	$G_{LO}$	$G_{StdLin,LO}$	$G_{StdLin,MIO}$	$G_{LS}$
solubility	10	0.0	485.5	0.0	18.6
	50	0.0	498.0	50.1	9.5
	100	0.0	481.2	70.5	0.3
	200	0.0	477.5	77.7	0.2
	500	0.0	501.3	103.2	0.2
permeability	10	0.0	589.5	0.0	6.1
	50	0.0	619.4	71.9	3.5
	100	0.0	614.1	75.0	1.8
	200	0.0	613.0	80.0	0.1
	500	0.0	610.4	85.9	0.0
winequalityred	10	1.5	11581.3	89.8	1.2
	50	3.4	11873.6	98.3	2.3
	100	4.3	12014.9	98.8	0.6
	200	4.3	12000.6	99.0	1.2
	500	4.5	12031.8	99.2	1.4
concrete	10	0.0	6210.6	72.5	0.0
	50	1.8	6657.1	95.0	0.0
	100	2.6	6706.6	98.3	0.0
	200	1.6	6622.2	98.5	0.0
	500	2.2	6652.6	98.8	0.0

Table 4. Table 4: Results of experiments comparing solution methods.

Data set	$T$	$N_{Levels}$	$N_{Leaves}$	${𝒯}_{Direct}$	${𝒯}_{Benders}$	${𝒯}_{SplitGen}$
solubility	100	4129	31720	1.7	1.1	0.8
	200	6016	63704	4.5	2.6	0.8
	500	9646	159639	177.9	28.6	13.7
permeability	100	2138	6108	0.2	0.0	0.0
	200	2138	12214	0.5	0.2	0.1
	500	2138	30443	2.7	0.4	0.7
winequalityred	100	3000	32659	51.6	41.1	56.5
	200	3495	65199	216.0	152.3	57.2
	500	3981	162936	1159.7	641.8	787.8
concrete	100	7851	29120	67.8	8.8	8.6
	200	10459	58242	183.8	12.5	15.4
	500	13988	145262	846.9	37.5	32.3

Table 5. Table 5: Summary of drug design data sets (see Ma et al. 2015 for further details).

Data set	Name	Num. Obs.	Num. Variables
1	3A4	37241	9491
2	CB1	8716	5877
3	DPP4	6148	5203
4	HIVINT	1815	4306
5	HIVPROT	3212	6274
6	LOGD	37388	8921
7	METAB	1569	4505
8	NK1	9965	5803
9	OX1	5351	4730
10	OX2	11151	5790
11	PGP	6399	5135
12	PPB	8651	5470
13	RAT_F	6105	5698
14	TDI	4165	5945
15	THROMBIN	5059	5552

Table 6. Table 6: Comparison of split generation and Benders decomposition for drug design data sets.

Data set	$N_{Levels}$	$N_{Leaves}$	${𝒯}_{Direct}$	${𝒯}_{SplitGen}$	${𝒯}_{Benders}$	${𝒯}_{LS}$	$G_{Direct}$	$G_{SplitGen}$	$G_{Benders}$	$G_{LS}$
1	27145	852533	151.2	97.3	7200.0	390.8	0.00	0.00	0.07	9.26
2	16480	289800	7201.6	6533.8	7200.0	132.6	0.13	0.00	1.30	8.64
3	13697	201265	7200.7	6252.2	7200.0	84.9	0.22	0.00	1.41	11.58
4	11790	59552	2.9	2.1	6.3	55.4	0.00	0.00	0.00	5.92
5	16426	109378	62.5	23.9	7200.0	108.9	0.00	0.00	0.07	6.99
6	26962	1307848	7203.5	7219.8	7200.1	409.3	2.40	2.15	12.82	23.71
7	12523	53934	16.5	12.3	1743.0	60.8	0.00	0.00	0.00	11.17
8	17319	328705	7202.5	7200.5	7200.1	146.9	0.37	0.27	2.67	5.22
9	12595	184841	370.2	55.0	7200.1	73.8	0.00	0.00	0.46	15.76
10	15780	379583	4101.1	1339.5	7200.2	124.7	0.00	0.00	4.71	12.03
11	15111	217395	281.1	81.0	7200.1	94.0	0.00	0.00	0.55	12.87
12	15737	291709	32.6	40.0	7200.0	94.4	0.00	0.00	0.02	12.17
13	17841	212926	7202.0	6731.5	7200.0	137.3	0.43	0.00	3.73	26.38
14	16272	145476	13.1	11.5	41.8	110.9	0.00	0.00	0.00	17.10
15	14863	169638	388.1	223.1	7200.0	111.7	0.00	0.00	0.89	13.10

Table 7. Table 7: Average and maximum proximity of split generation solutions for drug design data sets.

Data set	${\pi }_{\mathrm{avg}}$	${\pi }_{\max }$
1	0.00006	0.45
2	0.00033	0.37
3	0.00042	0.11
4	0.00157	0.71
5	0.00071	0.40
6	0.00005	0.15
7	0.00115	0.43
8	0.00029	0.42

Table 8. Table 8 : Comparison of out-of-sample profit prediction R 2 superscript 𝑅 2 R^{2} for the different models and brands, averaged over ten random splits of the data, for the first prediction task ( c i = 0 subscript 𝑐 𝑖 0 c_{i}=0 for all brands i 𝑖 i ). The best R 2 superscript 𝑅 2 R^{2} value for each brand is indicated in bold.

	Brand $R^2$
Model	1	2	3	4	5	6	7	8	9	10	11
HB-SemiLog	0.34	0.74	0.41	0.07	0.38	0.53	0.47	0.34	0.25	–0.43	0.55
HB-SemiLog-Own	0.48	0.80	0.52	0.35	0.43	0.62	0.53	0.42	0.28	0.42	0.65
HB-LogLog	0.39	0.75	0.47	0.08	0.38	0.55	0.61	0.38	0.45	–0.17	0.58
HB-LogLog-Own	0.50	0.80	0.59	0.36	0.45	0.63	0.66	0.45	0.49	0.50	0.66
RF, $T=20$	0.75	0.84	0.69	0.83	0.83	0.69	0.76	0.64	0.68	0.79	0.73
RF, $T=50$	0.77	0.85	0.70	0.84	0.84	0.70	0.77	0.66	0.69	0.79	0.74
RF, $T=100$	0.77	0.85	0.70	0.84	0.84	0.71	0.77	0.67	0.70	0.80	0.74
RF, $T=500$	0.78	0.85	0.70	0.84	0.84	0.71	0.77	0.67	0.70	0.80	0.74
RF-Own, $T=20$	0.73	0.84	0.68	0.81	0.79	0.70	0.75	0.64	0.68	0.74	0.72
RF-Own, $T=50$	0.74	0.84	0.69	0.82	0.80	0.70	0.77	0.65	0.68	0.76	0.73
RF-Own, $T=100$	0.74	0.84	0.69	0.82	0.80	0.71	0.77	0.66	0.68	0.76	0.73
RF-Own, $T=500$	0.75	0.84	0.70	0.82	0.80	0.71	0.77	0.66	0.69	0.77	0.74

Table 9. Table 9 : Comparison of out-of-sample profit prediction R 2 superscript 𝑅 2 R^{2} for the different models and brands, averaged over ten random splits of the data, for the second prediction task ( c i = 0.9 ​ p i , min subscript 𝑐 𝑖 0.9 subscript 𝑝 𝑖 c_{i}=0.9p_{i,\min} for all brands i 𝑖 i ). The best R 2 superscript 𝑅 2 R^{2} value for each brand is indicated in bold.

	Brand $R^2$
Model	1	2	3	4	5	6	7	8	9	10	11
HB-SemiLog	0.25	0.75	0.25	0.08	0.33	0.53	0.40	0.32	0.15	-0.62	0.55
HB-SemiLog-Own	0.39	0.81	0.39	0.31	0.38	0.63	0.53	0.44	0.30	0.32	0.66
HB-LogLog	0.32	0.76	0.30	0.13	0.36	0.55	0.49	0.34	0.31	-0.30	0.58
HB-LogLog-Own	0.44	0.81	0.43	0.32	0.42	0.64	0.59	0.44	0.43	0.42	0.67
RF, $T=20$	0.72	0.84	0.62	0.84	0.82	0.71	0.75	0.65	0.65	0.77	0.74
RF, $T=50$	0.74	0.84	0.63	0.84	0.83	0.72	0.76	0.66	0.67	0.79	0.75
RF, $T=100$	0.75	0.85	0.63	0.84	0.83	0.72	0.76	0.67	0.67	0.79	0.76
RF, $T=500$	0.75	0.85	0.63	0.84	0.83	0.72	0.76	0.67	0.67	0.79	0.76
RF-Own, $T=20$	0.69	0.83	0.61	0.81	0.76	0.71	0.75	0.64	0.64	0.73	0.73
RF-Own, $T=50$	0.70	0.84	0.62	0.81	0.77	0.72	0.76	0.66	0.66	0.74	0.74
RF-Own, $T=100$	0.71	0.84	0.63	0.82	0.78	0.72	0.76	0.67	0.66	0.75	0.75
RF-Own, $T=500$	0.72	0.84	0.63	0.82	0.78	0.72	0.77	0.67	0.66	0.75	0.75

Table 10. Table 10 : Gurobi parameters for RF profit optimization problem in Section  11.4 .

Parameter	Value
Heuristics	0
Cuts	0
VarBranch	1
InfUnbdInfo	1
PrePasses	1

Table 11. Table 11 : Average and maximum computation times (in seconds) of store-level optimal prices for the different models.

Price	RF		HB-SemiLog-Own		HB-LogLog-Own
Increment $\delta$	Avg.	Max.	Avg.	Max.	Avg.	Max.
0.05	11.4	26.3	1.9	2.5	2.0	2.6
0.10	7.3	16.5	1.2	1.6	1.2	4.9
0.20	1.8	2.3	0.7	0.9	0.7	1.0

Table 12. Table 12 : Multi-model comparison for optimal prices under the three different models.

Price		$\vartheta (m_1,m_2)$
Increment $\delta$		$m_2=$ HB-SemiLog-Own	$m_2=$ HB-LogLog-Own	$m_2=$ RF
0.05	$m_1=$ HB-SemiLog-Own	100.0	99.5	80.6
	$m_1=$ HB-LogLog-Own	99.5	100.0	83.8
	$m_1=$ RF	83.3	82.9	100.0
0.1	$m_1=$ HB-SemiLog-Own	100.0	99.5	83.1
	$m_1=$ HB-LogLog-Own	99.5	100.0	85.8
	$m_1=$ RF	84.0	83.5	100.0
0.2	$m_1=$ HB-SemiLog-Own	100.0	99.6	87.6
	$m_1=$ HB-LogLog-Own	99.7	100.0	89.5
	$m_1=$ RF	88.1	88.0	100.0

Equations

$$X_i\le a$$

$$X_i\in A$$

$$\sum _{t=1}^T{\lambda }_t{f}_t(\mathbf{X}).$$

$$\underset{\mathbf{X}\in \mathcal{X}}{\mathrm{maximize}}\sum _{t=1}^T{\lambda }_t{f}_t(\mathbf{X}).$$

$$x_{i,j}$$

$$\underset{\mathbf{x},\mathbf{y}}{\mathrm{maximize}}$$

$$\sum _{\ell \in \mathbf{leaves}(t)}{y}_{t,\ell }=1,\forall t\in \{1,\dots ,T\},$$

$$\sum _{\ell \in \mathbf{left}(s)}{y}_{t,\ell }\le \sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},$$

$$\forall t\in \{1,\dots ,T\},s\in \mathbf{splits}(t),$$

$$\sum _{\ell \in \mathbf{right}(s)}{y}_{t,\ell }\le 1-\sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},$$

$$\forall t\in \{1,\dots ,T\},s\in \mathbf{splits}(t),$$

$$\sum _{j=1}^{K_i}{x}_{i,j}=1,\forall i\in \mathcal{C},$$

$$x_{i,j}\le x_{i,j+1},\forall i\in \mathcal{N},j\in \{1,\dots ,K_i-1\},$$

$$x_{i,j}\in \{0,1\},\forall i\in \{1,\dots ,n\},j\in \{1,\dots ,K_i\}$$

$$y_{t,\ell }\ge 0,\forall t\in \{1,\dots ,T\},\ell \in \mathbf{leaves}(t).$$

$$x_{i,1}$$

$$x_{i,2}$$

$$x_{i,3}$$

$$⋮$$

$$x_{i,K_i}$$

$$x_{i,K_i+1}$$

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

$$\displaystyle\text{constraints~{}\eqref{prob:TEOMIO_categorical}-\eqref{prob:TEOMIO_xbinary}}.$$

$$\underset{\mathbf{x},\mathbf{y}}{\mathrm{maximize}}$$

$$y_{t,\ell }\le \sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},\forall t\in \{1,\dots ,T\},\ell \in \mathbf{leaves}(t),s\in \mathbf{LS}(\ell ),$$

$$y_{t,\ell }\le 1-\sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},\forall t\in \{1,\dots ,T\},\ell \in \mathbf{leaves}(t),s\in \mathbf{RS}(\ell ),$$

$$y_{t,\ell }\ge \sum _{s\in \mathbf{LS}(\ell )}\left(\sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j}\right)+\sum _{s\in \mathbf{RS}(\ell )}\left(1-\sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j}\right)-(|\mathbf{LS}(\ell )|+|\mathbf{RS}(\ell )|-1),$$

$$\forall t\in \{1,\dots ,T\},\ell \in \mathbf{leaves}(t),$$

$$\displaystyle\text{constraints~{}\eqref{prob:TEOMIO_categorical}-\eqref{prob:TEOMIO_ycontinuous}}.$$

$$\underset{\mathbf{x},\mathbf{y}}{\mathrm{maximize}}$$

$$\sum _{\ell \in \mathbf{left}(s)}{y}_{t,\ell }\le \sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},\forall (t,s)\in \bar{\Omega },$$

$$\sum _{\ell \in \mathbf{right}(s)}{y}_{t,\ell }\le 1-\sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},\forall (t,s)\in \bar{\Omega },$$

$$\displaystyle\text{constraints~{}\eqref{prob:TEOMIO_ysumtoone}, \eqref{prob:TEOMIO_categorical} - \eqref{prob:TEOMIO_ycontinuous}}.$$

$${\delta }_{t,s}=\max \left\{\underset{\ell \in \mathbf{left}(s)}{\max }{p}_{t,\ell }-\underset{\ell \in \mathbf{left}(s)}{\min }{p}_{t,\ell },\underset{\ell \in \mathbf{right}(s)}{\max }{p}_{t,\ell }-\underset{\ell \in \mathbf{right}(s)}{\min }{p}_{t,\ell }\right\}.$$

$${\Delta }_t=\underset{s\in \mathbf{splits}(t,d)}{\max }{\delta }_{t,s}.$$

$$Z_{MIO,d}^{\ast }-\sum _{t=1}^T{\lambda }_t{\Delta }_t\le Z_d\le Z_{MIO}^{\ast }\le Z_{MIO,d}^{\ast }.$$

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

$$\displaystyle\text{constraints~{}\eqref{prob:TEOMIO_categorical} - \eqref{prob:TEOMIO_xbinary}},$$

$$G_t(\mathbf{x})=$$

$$\sum _{\ell \in \mathbf{leaves}(t)}{y}_{t,\ell }=1,$$

$$\sum _{\ell \in \mathbf{left}(s)}{y}_{t,\ell }\le \sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},\forall s\in \mathbf{splits}(t),$$

$$\sum _{\ell \in \mathbf{right}(s)}{y}_{t,\ell }\le 1-\sum _{j\in \mathbf{C}(s)}{x}_{\mathbf{V}(s),j},\forall s\in \mathbf{splits}(t),$$

$$y_{t,\ell }\ge 0,\forall \ell \in \mathbf{leaves}(t).$$

$$y_{t,\ell}=\left\{\begin{array}[]{ll}1&\text{if}\ \ell=\ell^{*},\\ 0&\text{otherwise}.\end{array}\right.$$

Full text

\SingleSpacedXI\TheoremsNumberedThrough\ECRepeatTheorems\EquationsNumberedThrough

\RUNTITLE

Optimization of Tree Ensembles \RUNAUTHORMišić

\TITLE

Optimization of Tree Ensembles

\ARTICLEAUTHORS\AUTHOR

Velibor V. Mišić \AFFAnderson School of Management, University of California, Los Angeles, 110 Westwood Plaza, Los Angeles, CA, 90095, \EMAILvelibor.misic@anderson.ucla.edu

\ABSTRACT

Tree ensemble models such as random forests and boosted trees are among the most widely used and practically successful predictive models in applied machine learning and business analytics. Although such models have been used to make predictions based on exogenous, uncontrollable independent variables, they are increasingly being used to make predictions where the independent variables are controllable and are also decision variables. In this paper, we study the problem of tree ensemble optimization: given a tree ensemble that predicts some dependent variable using controllable independent variables, how should we set these variables so as to maximize the predicted value? We formulate the problem as a mixed-integer optimization problem. We theoretically examine the strength of our formulation, provide a hierarchy of approximate formulations with bounds on approximation quality and exploit the structure of the problem to develop two large-scale solution methods, one based on Benders decomposition and one based on iteratively generating tree split constraints. We test our methodology on real data sets, including two case studies in drug design and customized pricing, and show that our methodology can efficiently solve large-scale instances to near or full optimality, and outperforms solutions obtained by heuristic approaches.

\KEYWORDS

tree ensembles; random forests; mixed-integer optimization; drug design; customized pricing.

1 Introduction

A decision tree is a model used for predicting a dependent variable $Y$ using a collection of independent variables $\mathbf{X}=(X_{1},\dots,X_{n})$. To make a prediction, we start at the root of the tree, and check a query (e.g., “Is $X_{3}\leq 5.6$?”); we then proceed to the left child node if the query is true, and to the right child node if the query is false. We then check the new node’s query; the process continues until we reach a leaf node, where the tree outputs a prediction. A generalization of this type of model, called a tree ensemble model, involves making this type of prediction from each of a collection of trees and aggregating the individual predictions into a single prediction (for example, by taking a weighted average of the predictions for a regression setting, or by taking a majority vote of the predictions for a classification setting). An example of a decision tree and a prediction being made is given in Figure 1.

Many types of tree ensemble models have been proposed in the machine learning literature; the most notable examples are random forests and boosted trees. Tree ensemble models are extremely attractive due to their ability to model complex, nonlinear relationships between the independent variables $\mathbf{X}$ and the dependent variable $Y$. As a result, tree ensemble models have gained or are gaining widespread popularity in a number of application areas, including chemistry (Svetnik et al. 2003), genomics (Díaz-Uriarte and De Andres 2006), economics (Varian 2014, Bajari et al. 2015), marketing (Lemmens and Croux 2006) and operations management (Ferreira et al. 2015).

In many applications of tree ensemble models and predictive models in general, the independent variables that are used for prediction are exogenous and beyond our control as the decision maker. For example, one might build a random forest model to predict whether a patient is at risk of developing a disease based on the patient’s age, blood pressure, family history of the disease and so on. Clearly, features such as age and family history are not amenable to intervention. Such models are typically used for some form of post hoc action or prioritization. In the disease risk example, one can use the random forest model to rank patients by decreasing predicted risk of an imminent acute event, and this ranking can be used to determine which patients should receive some intervention, e.g., closer monitoring, treatment with a particular drug and so on.

In an increasing number of predictive modeling applications, however, some of the independent variables are controllable; that is to say, those independent variables are also decision variables. We provide a couple of examples:

Design of drug therapies. In a healthcare context, one might be interested in building a model to predict patient response given a particular drug therapy, and then using such a model to find the optimal drug therapy. The dependent variable is some metric of the patient’s health, and the independent variables may specify the drug therapy (which drugs and in what doses) and characteristics of the patient. A recent example of such an approach can be found in Bertsimas et al. (2016), which considers the design of clinical trials for combination drug chemotherapy for gastric cancer; the first step involves estimating a model (specifically, a ridge regression model) of patient survival and toxicity using information about the drug therapy, and the second step involves solving an optimization problem to find the drug therapy that maximizes the predicted survival of the given patient group subject to a constraint on the predicted toxicity. 2. 2.

Pricing/promotion planning. In marketing and operations management, a fundamental problem is that of deciding which products should be promoted when and at what price. In such a context, the data might consist of weekly sales of each product in a category, and the prices of the products in that category during that week and in previous weeks; one might use this data to build a predictive model of demand as a function of the prices, and then optimize such a model to decide on a promotion schedule (for a recent example see, e.g., Cohen et al. 2017).

In this paper, we seek to unlock the prescriptive potential of tree ensemble models by considering the problem of tree ensemble optimization. This problem is concerned with the following question: given a tree ensemble model that predicts some quantity $Y$ using a set of controllable independent variables $\mathbf{X}$, how should we set the independent variables $\mathbf{X}$ so as to maximize the predicted value of $Y$? This problem is of significant practical interest because it allows us to leverage the high accuracy afforded by tree ensemble models to obtain high quality decisions. At the same time, the problem is challenging, due to the highly nonlinear and large-scale nature of tree ensemble models.

We make the following contributions:

We propose the tree ensemble optimization problem and we show how to formulate this problem as a mixed-integer optimization (MIO) problem. The formulation can accommodate independent variables that are discrete, categorical variables as well as continuous, numeric variables. To the best of our knowledge, the problem of optimizing an objective function described as the prediction of a tree ensemble has not been previously proposed in either the machine learning or the operations research community. 2. 2.

From a theoretical standpoint, we develop a number of results that generally concern the tractability of the formulation. First, we prove that the tree ensemble optimization problem is in general NP-Hard. We then show that our proposed MIO formulation offers a tighter relaxation of the problem than an alternate MIO formulation, obtained by applying a standard linearization to a binary polynomial formulation of the problem. We develop a hierarchy of approximate formulations for the problem, obtained by truncating each tree in the ensemble to a depth $d$ from the root node. We prove that the objective value of such an approximate formulation is an upper bound that improves as $d$ increases, and show how to construct a complementary a priori lower bound that depends on the variability of each tree’s prediction below the truncation depth $d$. 3. 3.

From a solution methodology standpoint, we present two different strategies for tackling large-scale instances of our MIO formulation. The first is based on solving a Benders reformulation of the problem using constraint generation. Here, we analyze the structure of the Benders subproblem and show that it can be solved efficiently. The second is based on applying lazy constraint generation directly to our MIO formulation. For this approach, we propose an efficient procedure for identifying violated constraints, which involves simply traversing each tree in the ensemble. 4. 4.

We evaluate the effectiveness of our formulation and solution methods computationally using an assortment of real data sets. We show that the full MIO formulation can be solved to full optimality for small to medium sized instances within minutes, and that our formulation is significantly stronger in terms of relaxation bound and solution time than the aforementioned alternate formulation. We also show that our approach often significantly outperforms a simple local search heuristic that does not guarantee optimality. Lastly, we show that our customized solution methods can drastically reduce the solution time of our formulation. 5. 5.

We provide a deeper showcase of the utility of our approach in two applications. The first is a case study in drug design, using a publicly available data set from Merck Research Labs (Ma et al. 2015). Here, we show that our approach can optimize large-scale tree ensemble models with thousands of independent variables to full or near optimality within a two hour time limit, and can be used to construct a Pareto frontier of compounds that efficiently trade off predicted performance and similarity to existing, already-tested compounds. The second is a case study in customized pricing using a supermarket scanner data set (Montgomery 1997). Here, we show that a random forest model leads to considerable improvements in out-of-sample prediction accuracy over two state-of-the-art models based on hierarchical Bayesian regression, and that our optimization approach can find provably optimal prices at the individual store level within seconds.

The rest of the paper is organized as follows. In Section 2, we survey some of the related literature to this work. In Section 3, we present our formulation of the tree ensemble optimization problem as an MIO problem, and provide theoretical results on the structure of this problem. In Section 4, we present two solution approaches for large-scale instances of the tree ensemble optimization problem. In Section 5, we present the results of our computational experiments with real data sets. In Section 6 we present our case study in drug design and in Section 7, we summarize our case study in customized pricing (described in greater detail in Section 11 in the electronic companion). We give concluding remarks in Section 8.

2 Literature review

Decision tree models became popular in machine learning with the introduction of two algorithms, ID3 (iterative dichotomiser; see Quinlan 1986) and CART (classification and regression tree; see Breiman et al. 1984). Decision tree models gained popularity due to their interpretability, but were found to be generally less accurate than other models such as linear and logistic regression. A number of ideas were thus consequently proposed for improving the accuracy of tree models, which are all based on constructing an ensemble of tree models. Breiman (1996) proposed the idea of bootstrap aggregation, or bagging, where one builds a collection of predictive models, each trained with a bootstrapped sample of the original training set; the predictions of each model are then aggregated into a single prediction (by majority vote for classification and by averaging for regression). The motivation for bagging is that it reduces the prediction error for predictive models that are unstable/highly sensitive to the training data (such as CART); indeed, Breiman (1996) showed that bagged regression trees can be significantly better than ordinary regression trees in out-of-sample prediction error. Later, Breiman (2001) proposed the random forest model, where one builds a collection of bagged CART trees for which the subset of features selected for splitting at each node of each tree is randomly sampled from the set of all features (the so-called random subspace method; see Ho 1998). Concurrently, other research has considered the idea of boosting (Schapire and Freund 2012), wherein one iteratively builds a weighted collection of basic predictive models (such as CART trees), with the goal of reducing the prediction error with each iteration.

Tree ensembles occupy a central place in machine learning because they generally work very well in practice. In a systematic comparison of 179 different prediction methods on a broad set of benchmark data sets, Fernández-Delgado et al. (2014) found that random forests achieved best or near-best performance over all of these data sets. Boosted trees have been similarly successful: on the data science competition website Kaggle, one popular implementation of boosted trees, XGBoost, was used in more than half of the winning solutions in the year 2015 (Chen and Guestrin 2016). There exist robust and open source software implementations of many tree ensemble models. For boosted trees, the R package gbm (Ridgeway 2006) and XGBoost are widely used; for random forests, the R package randomForest (Liaw and Wiener 2002) is extremely popular.

There has also been a significant concurrent effort to develop a theoretical foundation for tree ensemble methods; we briefly survey some of this work for random forests. For random forests, the original paper of Breiman (2001) developed an upper bound on the generalization error of a random forest. Later research studied the consistency of simplified versions of the random forest model (e.g., Biau et al. 2008) as well as the original model (e.g., Scornet et al. 2015). For an overview of recent theoretical advances in random forests, the reader is referred to Biau and Scornet (2016).

At the same time, there is an increasing number of papers originating in operations research where a predictive model is used to represent the effect of the decision, and one solves an optimization problem to find the best decision with respect to this predictive model. Aside from the papers already mentioned in clinical trials and promotion planning, we mention two other examples. In pricing, data on historical prices and demand observed at those prices is often used to build a predictive model of demand as a function of price and to then determine the optimal price (e.g., Besbes et al. 2010). In assortment optimization, Mišić (2016) considers a two-step approach, where the first step involves estimating a ranking-based choice model from historical data, and the second involves solving an MIO model to optimize the revenue predicted under that model.

In the research literature where predictive models are used to understand and subsequently optimize decisions, the closest paper conceptually to this one is the paper of Ferreira et al. (2015). This paper considers the problem of pricing weekly sales for an online fashion retailer. To solve the problem, the paper builds a random forest model of the demand of each style included in a sale as a function of the style’s price and the average price of the other styles. The paper then formulates an MIO problem to determine the optimal prices for the styles to be included in the sale, where the revenue is based on the price and the demand (as predicted by the random forest) of each style. The MIO formulation does not explicitly model the random forest prediction mechanism – instead, one computes the random forest prediction for each style at each of its possible prices and at each possible average price of the other styles, and these predictions enter the MIO model as coefficients in the objective function. (The predictions could just as easily have come from a different form of predictive model, without changing the structure of the optimization problem.) In contrast, our MIO formulation explicitly represents the structure of each tree in the ensemble, allowing the prediction of each tree to be determined through the variables and constraints of the MIO. Although the modeling approach of Ferreira et al. (2015) is feasible for their pricing problem, it is difficult to extend this approach when there are many independent variables, as one would need to enumerate all possible combinations of values for them and compute the tree ensemble’s prediction for each combination of values. To the best of our knowledge, our paper is the first to conceptualize the problem of how to optimize an objective function that is given by a tree ensemble.

Methodologically, the present paper is most related to the paper of Bertsimas and Mišić (2019), which considers the problem of designing a product line or assortment to optimize revenue under a ranking-based choice model. The ranking-based model considered in Bertsimas and Mišić (2019) can be understood as a type of tree ensemble model; as such, the MIO formulation of Bertsimas and Mišić (2019) can be regarded as a special case of the more general formulation that we analyze here. Some of the theoretical results found in Bertsimas and Mišić (2019) – specifically, those results on the structure of the Benders cuts – are generalized in the present paper to tree ensemble models. Despite this similarity, the goals of the two papers are different. Bertsimas and Mišić (2019) consider an optimization approach specifically for product line decisions, whereas in the present paper, our methodology can be applied to any tree ensemble model, thus spanning a significantly broader range of application domains. Indeed, later in the paper we will present two different case studies – one on drug design (Section 6) and one on customized pricing (Section 11) – to illustrate the broad applicability of tree ensemble optimization.

Finally, we note that there is a growing literature on the use of mixed-integer optimization for the purpose of estimating decision tree models and other forms of statistical models. For example, Bertsimas and Dunn (2017) consider an exact MIO approach to constructing CART trees, while Bertsimas and King (2015) consider an MIO approach to model selection in linear regression. While the present paper is related to this previous work in that it also leverages the technology of MIO, the goal of the present paper is different. The above papers focus on the estimation of trees and other statistical models, whereas our paper is focused on optimization, namely, how to determine the optimal decision with respect to a given, fixed tree ensemble model.

3 Model

We begin by providing some background on tree ensemble models in Section 3.1 and defining the tree ensemble optimization problem. We then present our mixed-integer optimization model in Section 3.2. We provide results on the strength of our formulation in Section 3.3. Finally, in Section 3.4, we describe a hierarchy of approximate formulations based on depth truncation.

3.1 Background

In this section, we provide some background on tree ensemble models. We are given the task of predicting a dependent variable $Y$ using the independent variables $X_{1},\dots,X_{n}$; for convenience, we use $\mathbf{X}$ to denote the vector of independent variables. We let $\mathcal{X}_{i}$ denote the domain of independent variable $i$ and let $\mathcal{X}=\prod_{i=1}^{n}\mathcal{X}_{i}$ denote the domain of $\mathbf{X}$. An independent variable $i$ may be a numeric variable or a categorical variable.

A decision tree is a model for predicting the dependent variable $Y$ using the independent variable $\mathbf{X}$ by checking a collection of splits. A split is a condition or query on a single independent variable that is either true or false. More precisely, for a numeric variable $i$, a split is a query of the form

[TABLE]

for some $a\in\mathbb{R}$. For a categorical variable $i$, a split is a query of the form

[TABLE]

where $A\subseteq\mathcal{X}_{i}$ is a set of levels of the categorical variable. The splits are arranged in the form of a tree, with each split node having two child nodes. The left child corresponds to the split condition being true, while the right child corresponds to the condition being false. To make a prediction for an observation with the independent variable $\mathbf{X}$, we start at the root of the tree and check whether $\mathbf{X}$ satisfies the split condition; if it is true, we move to the left child, and if it is false, we move to the right child. At the new node, we check the split again, and move again to the corresponding node. This process continues until we reach a leaf of the tree. The prediction that we make is the value corresponding to the leaf we have reached.

In this paper, we will focus on predictive models that are ensembles or collections of decision trees. We assume that there are $T$ trees, where each tree is indexed from 1 to $T$. Each tree $t$ has a weight $\lambda_{t}$, and its prediction is denoted by the function $f_{t}$; for the independent variable $\mathbf{X}$, the prediction of tree $t$ is $f_{t}(\mathbf{X})$. For an observation with independent variable $\mathbf{X}$, the prediction of the ensemble of trees is given by

[TABLE]

The optimization problem that we would like to solve is to find the value of the independent variable $\mathbf{X}$ that maximizes the ensemble prediction:

[TABLE]

We shall make two key assumptions about the tree ensemble model $\sum_{t=1}^{T}\lambda_{t}f_{t}(\cdot)$ and our tree ensemble optimization problem (1):

First, we assume that we are only making a single, one-time decision and that the tree ensemble model is fixed. Extending our approach to the multistage setting is an interesting direction for future research. 2. 2.

Second, we assume that the tree ensemble model is an accurate representation of the outcome when we make the decision $\mathbf{X}$. In practice, some care must be taken here because depending on how the tree ensemble model is estimated and the nature of the data, the prediction $\sum_{t=1}^{T}\lambda_{t}f_{t}(\mathbf{X})$ may not necessarily be an accurate estimate of the causal effect of setting the independent variables to $\mathbf{X}$. This issue has been the focus of some recent work in prescriptive analytics (see Bertsimas and Kallus 2016, Kallus 2016). Our goal in this paper is to address only the question of optimization – how to efficiently and scalably optimize a tree ensemble function $\sum_{t=1}^{T}\lambda_{t}f_{t}(\cdot)$ – which is independent of such statistical questions. As such, we will assume that the tree ensemble model we are given at the outset is beyond suspicion.

Problem (1) is very general, and one question we may have is whether it is theoretically tractable or not. Our first theoretical result answers this question in the negative.

Proposition 3.1

The tree ensemble optimization problem (1) is NP-Hard.

The proof of Proposition 3.1, given in Section 9.2 of the e-companion, uses a reduction from the minimum vertex cover problem.

3.2 Optimization model

We now present an MIO formulation of (1). Before we present the model, we will require some additional notation. We let $\mathcal{N}$ denote the set of numeric variables and $\mathcal{C}$ denote the set of categorical variables; we have that $\mathcal{N}\cup\mathcal{C}=\{1,\dots,n\}$.

For each numeric variable $i\in\mathcal{N}$, let $\mathcal{A}_{i}$ denote the set of unique split points, that is, the set of values $a$ such that $X_{i}\leq a$ is a split condition in some tree in the ensemble $\{f_{t}\}_{t=1}^{T}$. Let $K_{i}=|\mathcal{A}_{i}|$ be the number of unique split points. Let $a_{i,j}$ denote the $j$th smallest split point of variable $i$, so that $a_{i,1}<a_{i,2}<\dots<a_{i,K_{i}}$.

For each categorical variable $i\in\mathcal{C}$, recall that $\mathcal{X}_{i}$ is the set of possible values of $i$. For convenience, we use $K_{i}$ in this case to denote the size of $\mathcal{X}_{i}$ (i.e., $K_{i}=|\mathcal{X}_{i}|$) and use the values $1,2,\dots,K_{i}$ to denote the possible levels of variable $i$.

Let $\mathbf{leaves}(t)$ be the set of leaves or terminal nodes of tree $t$. Let $\mathbf{splits}(t)$ denote the set of splits of tree $t$ (non-terminal nodes). Recall that the left branch of the split corresponds to “yes” or “true” to the split query, and the right branch corresponds to “no” or “false”. Therefore, for each split $s$ in $S_{t}$, we let $\mathbf{left}(s)$ be the set of leaves that are accessible from the left branch (all of the leaves for which the condition of split $s$ must be true), and $\mathbf{right}(s)$ be the set of leaves that are accessible from the right branch (all of the leaves for which the condition of split $s$ must be false). For each split $s$, we let $\mathbf{V}(s)\in\{1,\dots,n\}$ denote the variable that participates in split $s$, and let $\mathbf{C}(s)$ denote the set of values of variable $i$ that participate in the split query of $s$. Specifically, if $\mathbf{V}(s)$ is numeric, then $\mathbf{C}(s)=\{j\}$ for some $j\in\{1,\dots,K_{\mathbf{V}(s)}\}$, which corresponds to the split query $X_{i}\leq a_{i,j}$. If $\mathbf{V}(s)$ is categorical, then $\mathbf{C}(s)\subseteq\{1,\dots,K_{\mathbf{V}(s)}\}$, which corresponds to the query $X_{i}\in\mathbf{C}(s)$. For each leaf $\ell\in\mathbf{leaves}(t)$, we use $p_{t,\ell}$ to denote the prediction that tree $t$ makes when an observation reaches leaf $\ell$.

We now define the decision variables of the problem. There are two sets of decision variables. The first set is used to specify the independent variable value $\mathbf{X}$. For each categorical independent variable $i\in\mathcal{C}$ and each category/level $j\in\mathcal{X}_{i}$, we let $x_{i,j}$ be 1 if independent variable $i$ is set to level $j$, and 0 otherwise. For each numeric independent variable $i\in\mathcal{N}$ and each $j\in\{1,\dots,K_{i}\}$, we let $x_{i,j}$ be 1 if independent variable $i$ is set to a value less than or equal to the $j$th split point, and 0 otherwise. Mathematically,

[TABLE]

We use $\mathbf{x}$ to denote the vector of $x_{i,j}$ values.

The second set of decision variables is used to specify the prediction of each tree $t$. For each tree $t$ and each leaf $\ell\in\mathbf{leaves}(t)$, we let $y_{t,\ell}$ be a binary decision variable that is 1 if the observation encoded by $\mathbf{x}$ belongs to/falls into leaf $\ell$ of tree $t$, and 0 otherwise.

With these definitions, the MIO can be written as follows:

[TABLE]

The constraints have the following meaning. Constraint (2b) ensures that the observation falls in exactly one of the leaves of each tree $t$. Constraint (2c) ensures that, if $\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=0$, then $y_{t,\ell}$ is forced to zero for all $\ell\in\mathbf{left}(s)$; in words, the condition of the split is false, so the observation cannot fall into any leaf to the left of split $s$, as this would require the condition to be true. Similarly, constraint (2d) ensures that if the condition of split $s$ is satisfied, then $y_{t,\ell}$ is forced to zero for all $\ell\in\mathbf{right}(s)$; in words, the condition of the split is true, so the observation cannot fall into any leaf to the right of split $s$, as this would require the condition to be false. Constraint (2e) ensures that for each categorical variable $i\in\mathcal{C}$, exactly one of the $K_{i}$ levels is selected. Constraint (2f) requires that if the numeric independent variable $i$ is less than or equal to the $j$th lowest split point, then it must also be less than or equal to the $(j+1)$th lowest split point. Constraint (2g) defines each $x_{i,j}$ to be binary, while constraint (2h) defines each $y_{t,\ell}$ to be nonnegative. The objective represents the prediction of the ensemble of trees on the observation that is encoded by $\mathbf{x}$.

We now comment on several features of the model. The first is that the $y_{t,\ell}$ variables, despite having a binary meaning, are defined as continuous variables. The reason for this is that when $\mathbf{x}$ is binary, the constraints automatically force $\mathbf{y}$ to be binary. We will formally state this result later (Proposition 4.1 of Section 4.1). As a result, the only binary variables are those in $\mathbf{x}$, of which there are $\sum_{i=1}^{n}K_{i}$. Recall that for categorical independent variables, $K_{i}$ is the number of levels, whereas for numeric independent variables, $K_{i}$ is the number of unique split points found in the tree ensemble $\{f_{t}\}_{t=1}^{T}$.

The second is that our formulation does not model the exact value of each numeric independent variable $i\in\mathcal{N}$. In contrast, the formulation only models where the variable is in relation to the unique split points in $\mathcal{A}_{i}$ – for example, $x_{i,1}=1$ indicates that independent variable $i$ is set to be less than or equal to the first lowest split point. The reason for this is that each decision tree function $f_{t}(\cdot)$ is a piecewise constant function and therefore the tree ensemble function $\sum_{t=1}^{T}\lambda_{t}f_{t}(\mathbf{X})$ is also piecewise constant. Thus, for the purpose of optimizing the function $\sum_{t=1}^{T}\lambda_{t}f_{t}(\mathbf{X})$, it is not necessary to explicitly maintain the value $X_{i}$ of each numeric independent variable $i$.

The third is that numeric independent variables are modeled in terms of an inequality, that is, $x_{i,j}=\mathbb{I}\{X_{i}\leq a_{i,j}\}$. Alternatively, one could model numeric independent variables by using $x_{i,j}$ to represent whether $X_{i}$ is between two consecutive split points, e.g.,

[TABLE]

One would then re-define the set $\mathbf{C}(s)$ for each split involving the variable $i$ so as to include all of the relevant $j$ values under this new encoding. The advantage of our choice of encoding – using $x_{i,j}=\mathbb{I}\{X_{i}\leq a_{i,j}\}$ – is that the resulting formulation enhances the power of branching on fractional values of $x_{i,j}$ and leads to more balanced branch-and-bound trees (Vielma 2015). This type of encoding has been used successfully in scheduling and transportation applications (so-called “by” variables, representing an event happening by some period $t$; see for example Bertsimas et al. 2011); for further details, the reader is referred to Vielma (2015).

3.3 Theoretical properties

We now compare our formulation against an alternate MIO formulation of the tree ensemble optimization problem, which involves relating the $\mathbf{y}$ and $\mathbf{x}$ variables in a different way.

In particular, for any leaf $\ell$ of any tree $t$, let $\mathbf{LS}(\ell)$ be the set of splits for which leaf $\ell$ is on the left side (i.e., $s$ such that $\ell\in\mathbf{left}(s)$), and $\mathbf{RS}(\ell)$ be the set of splits for which leaf $\ell$ is on the right side (i.e., $s$ such that $\ell\in\mathbf{right}(s)$). The tree ensemble optimization problem can then be formulated as the following problem:

[TABLE]

The above problem is a binary polynomial problem. Note that the product term, $\prod_{s\in\mathbf{LS}(\ell)}\left(\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}\right)\cdot\prod_{s\in\mathbf{RS}(\ell)}\left(1-\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}\right)$, is exactly 1 if the observation is mapped to leaf $\ell$ of tree $t$, and 0 otherwise. The standard linearization of (3) (see Crama 1993) is the following MIO:

[TABLE]

Let $Z^{*}_{LO}$ be the optimal value of the linear optimization (LO) relaxation of (2) and let $Z^{*}_{LO,StdLin}$ be the optimal value of the LO relaxation of (4). The following result relates the two optimal values.

Proposition 3.2

$Z^{*}_{LO}\leq Z^{*}_{LO,StdLin}$.

The proof of Proposition 3.2 (see Section 9.3) consists of showing that an optimal solution of the relaxation of (2) is a feasible solution for the relaxation of (4) and achieves an objective value of exactly $Z^{*}_{LO}$. The significance of Proposition 3.2 is that it establishes that formulation (2) is a stronger formulation of the tree ensemble optimization problem than formulation (4). This is desirable from a practical perspective, as stronger MIO formulations are generally faster to solve than weaker MIO formulations. We shall see in Section 5.2 that the difference in relaxation bounds can be substantial and that formulation (4) is significantly less tractable than our formulation (2).

3.4 Depth $d$ approximation

In this section, we describe a hierarchy of relaxations of (2) that are based on approximating each tree in the ensemble up to a particular depth.

The motivation for this hierarchy of relaxations comes from the following observation regarding the size of (2). In particular, a key driver of the size of (2) is the number of left and right split constraints (2c) and (2d), respectively; these constraints are enforced for every single split in each tree in the ensemble. For a large number of trees that are deep (and thus have many splits), the resulting number of left and right split constraints will be large. At the same time, it may be reasonable to expect that if we do not represent each tree to its full depth, but instead only represent each tree up to some depth $d$ and only include splits that occur before (and including) depth $d$, then we might still obtain a reasonable solution to (2). In this section, we rigorously define this hierarchy of approximate formulations, and provide theoretical guarantees on how close such approximations are to the original formulation.

Let $\Omega=\{(t,s)\,|\,t\in\{1,\dots,T\},\ s\in\mathbf{splits}(t)\}$ be the set of tree-split pairs. Let $\bar{\Omega}\subseteq\Omega$ be a subset of all possible tree-split pairs. The $\bar{\Omega}$ tree ensemble problem is defined as problem (2) where constraints (2c) and (2d) are restricted to $\bar{\Omega}$:

[TABLE]

Problem (5) is obtained by removing constraints from (2), while retaining the same decision variables. Solving problem (5) for a fixed $\bar{\Omega}$ will result in a solution $(\mathbf{x},\mathbf{y})$ that is only guaranteed to satisfy constraints (2c) and (2d) for those $(t,s)$ pairs in $\bar{\Omega}$. As a result, for any $\bar{\Omega}\subseteq\Omega$, the objective value of (5) will be a valid upper bound on (2). We will now define a collection of subsets $\bar{\Omega}$, obtained by truncating each tree at a particular depth, for which we also have an easily-computable accompanying lower bound.

For any $d\in\mathbb{Z}_{+}$, let $\bar{\Omega}_{d}$ be the set of all tree-split pairs where the split is at a depth $d^{\prime}\leq d$ (a depth of 1 corresponds to the split at the root node). Let $Z^{*}_{MIO,d}$ denote the objective value of (5) with $\bar{\Omega}_{d}$, i.e., all splits up to and including depth $d$, and let $d_{\max}$ be the maximum depth of any split in any tree of the ensemble. Let $Z^{*}_{MIO}$ be the objective value of (2), where the split constraints are up to the full depth of each tree; note that $Z^{*}_{MIO}=Z^{*}_{MIO,d_{\max}}$. The following result establishes that $Z^{*}_{MIO,d}$ is an upper bound on $Z^{*}_{MIO}$ that becomes tighter as $d$ increases:

Proposition 3.3

$Z^{*}_{MIO,1}\geq Z^{*}_{MIO,2}\geq\dots\geq Z^{*}_{MIO,d_{\max}}=Z^{*}_{MIO}.$**

We now show how to construct a complementary lower bound.

Fix some depth $d$. Let $\mathbf{splits}(t,d)$ denote the set of splits at depth $d$; if the depth of the tree is strictly less than $d$, $\mathbf{splits}(t,d)$ is empty. Let us define the constant $\delta_{t,s}$ for each split $s\in\mathbf{splits}(t,d)$ of each tree $t$ as

[TABLE]

The constant $\delta_{t,s}$ is an upper bound on the maximum error possible (due to the depth $d$ truncation of the split constraints) in the prediction of tree $t$ for the observation encoded by $\mathbf{x}$, given that the observation reaches split $s$. Stated differently, if the observation $\mathbf{x}$ reaches split $s$, the error from omitting all splits below $s$ cannot be more than the difference between the highest and lowest leaf prediction on the left subtree or the right subtree at $s$.

We define $\Delta_{t}$ as the maximum of the $\delta_{t,s}$ values over all the depth $d$ splits of tree $t$:

[TABLE]

(In the case that $\mathbf{splits}(t,d)$ is empty, we set $\Delta_{t}=0$.)

Before stating our approximation guarantee, we note that given a solution $(\mathbf{x},\mathbf{y})$ that solves problem (5) with $\bar{\Omega}_{d}$, it is possible to find a solution $(\mathbf{x},\tilde{\mathbf{y}})$ that is a feasible solution for the full depth (2). Our approximation guarantee is given below.

Theorem 3.4

Suppose that $\lambda_{t}\geq 0$ for all $t\in\{1,\dots,T\}$ and $d\in\mathbb{Z}_{+}$. Let $(\mathbf{x},\mathbf{y})$ be an optimal solution of problem (5) with $\bar{\Omega}_{d}$. Let $Z_{d}$ be the true objective of $\mathbf{x}$ when embedded within the full-depth problem (5). We then have

[TABLE]

The above theorem, which we prove in Section 9.5, provides a guarantee on how suboptimal the $\bar{\mathbf{x}}$ solution, derived from the depth $d$ problem (5) with $\bar{\Omega}_{d}$, is for the true (full depth) problem (2). Note that the requirement of $\lambda_{t}$ being nonnegative is not particularly restrictive, as we can always make $\lambda_{t}$ of a given tree $t$ positive by negating the leaf predictions $p_{t,\ell}$ of that tree. This result is of practical relevance because it allows the decision maker to judiciously trade-off the complexity of the problem (represented by the depth $d$) against an a priori guarantee on the quality of the approximation. Moreover, the quantity $\sum_{t=1}^{T}\lambda_{t}\Delta_{t}$, which bounds the difference between $Z^{*}_{MIO,d}$ and $Z_{d}$, can be easily computed from each tree, allowing the bound to be readily implemented in practice. We shall see in Section 5.3 that although the lower bound can be rather conservative for small values of $d$, the true objective value of $\bar{\mathbf{x}}$ is often significantly better.

Theorem 3.4 suggests an interesting question regarding how the trees that are input to the tree ensemble optimization problem are estimated. In particular, for any given tree $t$, it may be possible to re-arrange the splits in the tree in a way that one obtains a new tree, $t^{\prime}$, such that $t^{\prime}$ gives identical predictions to $t$ for all $\mathbf{X}\in\mathcal{X}$. We may therefore ask: is there a way to re-arrange the splits of the trees so as to obtain a tighter upper bound $Z^{*}_{MIO,d}$, and a smaller approximation error bound $\sum_{t=1}^{T}\lambda_{t}\Delta_{t}$? This question constitutes an interesting direction for future work.

4 Solution methods

The optimization model in Section 3.2, although tractable for small to medium instances, can be difficult to solve directly for large instances. In this section, we present two solution approaches for tackling large-scale instances of (2). In Section 4.1, we present an approach based on Benders decomposition. In Section 4.2, we present an alternate approach based on iteratively generating the split constraints.

4.1 Benders decomposition

The first solution approach that we present is Benders decomposition. Recall that in (2), we have two sets of variables, $\mathbf{x}$ and $\mathbf{y}$; furthermore, $\mathbf{y}$ can be further partitioned as $\mathbf{y}=(\mathbf{y}_{1},\mathbf{y}_{2},\dots,\mathbf{y}_{T})$, where $\mathbf{y}_{t}$ is the collection of $y_{t,\ell}$ variables corresponding to tree $t$. For any two trees $t,t^{\prime}$ with $t\neq t^{\prime}$, notice that the variables $\mathbf{y}_{t}$ and $\mathbf{y}_{t^{\prime}}$ do not appear together in any constraints; they are only linked through the $\mathbf{x}$ variables.

The above observation suggests a Benders reformulation of (2). Re-write (2) as follows:

[TABLE]

where $G_{t}(\mathbf{x})$ is defined as the optimal value of the following subproblem:

[TABLE]

The first result we will prove is on the form of the optimal solution to (8). To do this, we first provide a procedure in Algorithm 1 for determining the leaf of the solution encoded by $\mathbf{x}$ for tree $t$. For ease of exposition, we will denote this procedure applied to a particular observation encoded by $\mathbf{x}$ and a given tree $t$ as $\textsc{GetLeaf}(\mathbf{x},t)$. We use $\mathbf{leftchild}(\nu)$ to denote the left child of a split node $\nu$, $\mathbf{rightchild}(\nu)$ to denote the right child of a split node $\nu$, and $\mathbf{root}(t)$ to denote the root split node of tree $t$. Note that if $d_{\max}$ is the maximum depth of the trees in the ensemble (i.e., the maximum depth of any split in any tree), then the time complexity of $\textsc{GetLeaf}(\mathbf{x},t)$ is $O(d_{\max})$, as we will traverse at most $d_{\max}$ split nodes to reach a leaf.

Having defined GetLeaf, we now present our first theoretical result (see Section 9.6 for the proof).

Proposition 4.1

Let $\mathbf{x}\in\{0,1\}^{\sum_{i=1}^{n}K_{i}}$ be a feasible solution of (7). Let $\ell^{*}=\textsc{GetLeaf}(\mathbf{x},t)$ be the leaf into which $\mathbf{x}$ falls, and let $\mathbf{y}_{t}$ be the solution to (8) defined as

[TABLE]

The solution $\mathbf{y}_{t}$ is the only feasible solution and therefore, the optimal solution of (8).

Since (8) is feasible and has a finite optimal value, then by LO strong duality the optimal value of (8) is equal to the optimal value of its dual. The dual of subproblem (8) is

[TABLE]

Letting $\mathcal{D}_{t}$ denote the set of dual feasible solutions $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ for subproblem $t$, we can re-write (7) as

[TABLE]

We can now solve (10) using constraint generation. In particular, we start with constraint (10b) enforced for a subset of dual solutions $\bar{\mathcal{D}}_{t}\subseteq\mathcal{D}_{t}$, and solve (10). This will yield a candidate integer solution $\mathbf{x}$. We then solve problem (9) for each tree $t$ to determine if there exists a solution $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})\in\mathcal{D}_{t}$ for which constraint (10b) is violated. If so, we add the constraint and solve the problem again. Otherwise, if no such $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ is found for any tree, then the current $\mathbf{x}$ solution is optimal.

In the above constraint generation scheme, the key step is to find a dual subproblem solution $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ whose constraint (10b) is violated. With this motivation, we now prove a result on the structure of an optimal solution to the dual problem (9) (see Section 9.7 for the proof).

Proposition 4.2

Let $\mathbf{x}\in\{0,1\}^{\sum_{i=1}^{n}K_{i}}$ be a feasible solution of problem (7). Let $\ell^{*}=\textsc{GetLeaf}(\mathbf{x},t)$. An optimal solution of dual subproblem (9) is then given as follows:

[TABLE]

The value of Proposition 4.2 is that we can check for violated constraints in problem (10) through a simple calculation, without invoking an LO solver.

We remark on a couple of important aspects of how Proposition 4.2 is used within our constraint generation approach. First, Proposition 4.2 requires that the candidate solution $\mathbf{x}$ to the master problem (10) be an integer solution. Thus, our constraint generation procedure can only be used to generate Benders cuts for candidate integer solutions of (10). It cannot be used to generate cuts for fractional solutions, such as those arising from the LO relaxation of (10).

Second, we solve (10) using branch-and-bound, with none of the Benders constraints (10b) enforced. At each integer solution $\mathbf{x}$ encountered in the branch-and-bound tree, we check whether it satisfies constraint (10b) for each tree $t$. This is accomplished by solving the subproblem (9) using Proposition 4.2 to produce a dual solution $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ and checking whether constraint (10b) is violated at that dual solution. This procedure is executed for each tree $t$, and any violated constraints are added to every node in the current branch-and-bound tree. (Note that each integer solution encountered in the branch-and-bound process may result in more than one constraint being added if constraint (10b) is violated for more than one tree in the ensemble.) In our numerical experiments, we implement this solution approach using lazy constraint generation, a computational paradigm implemented in modern solvers such as Gurobi that allows for constraints to be added as needed in a single branch-and-bound tree (Gurobi Optimization, Inc. 2015).

4.2 Split constraint generation

Recall from Section 3.4 that when there is a large number of trees and each tree is deep, the total number of splits will be large, and the number of left and right split constraints will be large. However, for a given encoding $\mathbf{x}$, observe that we do not need all of the left and right split constraints in order for $\mathbf{y}$ to be completely determined by $\mathbf{x}$. As an example, suppose for a tree $t$ that $s$ is the root split, and $\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=1$ (i.e., we take the left branch of the root split). In this case, the right split constraint (2d) will force all $y_{t,\ell}$ to zero for $\ell\in\mathbf{right}(s)$. It is clear that in this case, it is not necessary to include any left or right split constraint for any split node $s^{\prime}$ that is to the right of split $s$, because all of the $y_{t,\ell}$ values that could be affected by those constraints are already fixed to zero.

This suggests an alternate avenue to solving (2), based on iteratively generating the left and right split constraints. Rather than attempting to solve the full problem (2) with all of the left and right split constraints included in the model, start with a subset of left split constraints and a subset of right split constraints, and all of the decision variables ($x_{i,j}$ for all $i,j$, and $y_{t,\ell}$ for all $t,\ell$) of the full problem. This leads to a relaxation of (2), which we then solve. For the resulting solution $(\mathbf{x},\mathbf{y})$, determine whether there exist any tree-split pairs $(t,s)$ for which the left split constraint (2c) or right split constraint (2d) are violated. If a violated left or right split constraint is found, add the corresponding left or right constraint to the formulation, and solve it again. Repeat the procedure until no violated constraints are found, at which point we terminate with the current solution as the optimal solution.

The key implementation question in such a proposal is: how do we efficiently determine violated constraints? The answer to this question is given in the following proposition.

Proposition 4.3

Let $(\mathbf{x},\mathbf{y})\in\{0,1\}^{\sum_{i=1}^{n}K_{i}}\times\mathbb{R}^{\sum_{t=1}^{T}|\mathbf{leaves}(t)|}$ be a candidate solution to (2) that satisfies constraints (2b) and constraints (2e) to (2h). Let $t\in\{1,\dots,T\}$. The solution $(\mathbf{x},\mathbf{y})$ satisfies constraints (2c) and (2d) for all $s\in\mathbf{splits}(t)$ if and only if it satisfies constraint (2c) for $s\in\mathbf{RS}(\ell^{*})$ and constraint (2d) for $s\in\mathbf{LS}(\ell^{*})$, where $\ell^{*}=\textsc{GetLeaf}(\mathbf{x},t)$.

Proposition 4.3 (see Section 9.8 for the proof) states that, to check whether solution $(\mathbf{x},\mathbf{y})$ satisfies the split constraints for tree $t$, it is only necessary to check the split constraints for those splits that are traversed when the observation encoded by $\mathbf{x}$ is mapped to a leaf by the action of GetLeaf. This is a simple but extremely useful result, because it implies that we can check for violated constraints simply by traversing the tree, in the same way that we do to find the leaf of $\mathbf{x}$.

Algorithm 2 provides the pseudocode of this procedure. This algorithm involves taking the observation encoded by $\mathbf{x}$ and walking it down tree $t$, following the splits along the way. For each split we encounter, we determine whether we should proceed to the left child ($\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=1$) or to the right child ($\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=0$). If we are going to the left ($s\in\mathbf{LS}(\ell^{*})$ or equivalently, $\ell^{*}\in\mathbf{left}(s)$), then we check that $y_{t,\ell}$ is zero for all the leaves to the right of split $s$ (constraint (2d)). If we are going to the right ($s\in\mathbf{RS}(\ell^{*})$ or equivalently, $\ell^{*}\in\mathbf{right}(s)$), then we check that $y_{t,\ell}$ is zero for all the leaves to the left of split $s$ (constraint (2c)). In words, we are traversing the tree as we would to make a prediction, and we are simply checking that there is no positive $y_{t,\ell}$ that is on the “wrong” side of any left or right split that we take. If we reach a leaf node, we can conclude that the current solution $(\mathbf{x},\mathbf{y})$ does not violate any of the split constraints of tree $t$.

With regard to running time, Algorithm 2, like GetLeaf, has a complexity of $O(d_{\max})$: for a given tree, Algorithm 2 will check at most $d_{\max}$ left/right split constraints in the worst case. In contrast, the total number of split constraints for a given tree $t$ in (2) could be much larger than $O(d_{\max})$. For example, in a complete binary tree where every leaf is at depth $d_{\max}+1$, the total number of split constraints will be $O(2^{d_{\max}})$. Algorithm 2 thus provides a computationally tractable path to solving instances of (2) with deep trees and large numbers of split constraints.

As with our Benders approach, we comment on two important aspects of this solution approach. First, the separation procedure provided in Algorithm 2 can only be used for candidate integer solutions to (2); it cannot be used to separate fractional solutions, and thus cannot be used to solve the LO relaxation of (2). Extending Algorithm 2 and Proposition 4.3 to fractional solutions $(\mathbf{x},\mathbf{y})$ is an interesting direction for future research.

Second, we employ Algorithm 2 in a manner similar to our Benders approach in Section 4.1. In particular, we solve (2) using branch-and-bound, without any of the left/right split constraints enforced. Whenever a new candidate integer solution $(\mathbf{x},\mathbf{y})$ is found in the branch-and-bound process, we execute Algorithm 2 to determine if $(\mathbf{x},\mathbf{y})$ violates any left/right split constraints for any of the $T$ trees; if any violated constraints are found, they are added to every node in the branch-and bound tree. (This may result in more than one violated constraint being added, where each violated constraint will correspond to one tree.) This computational approach is also implemented using lazy constraint generation.

5 Computational experiments

In this section, we describe our first set of computational results. Section 5.1 provides the background of our experiments. Section 5.2 provides initial results on the full MIO formulation (2), while Section 5.3 provides results on the depth approximation scheme of Section 3.4. Finally, Section 5.4 compares the Benders and split generation solution methods against directly solving (2).

5.1 Background

We test our optimization formulation (2) and the associated solution methods from Section 4 using tree ensemble models estimated from real data sets, whose details are provided in Table 1. We wish to point out that in these data sets, the independent variables may in reality not be controllable. However, they are still useful in that they furnish us with real tree ensemble models for evaluating our optimization methodology.

We specifically focus on random forest models. Unless otherwise stated, all random forests are estimated in R using the randomForest package (Liaw and Wiener 2002), using the default parameters. All linear and mixed-integer optimization models are formulated in the Julia programming language (Bezanson et al. 2012), using the JuMP package (Julia for Mathematical Programming; see Lubin and Dunning 2015), and solved using Gurobi 6.5 (Gurobi Optimization, Inc. 2015). All experiments were executed on a late 2013 Apple Macbook Pro Retina laptop, with a quad-core 2.6GHz Intel i7 processor and 16GB of memory.

5.2 Full MIO formulation experiments

As part of our first experiment, we consider solving the unconstrained tree ensemble problem for each data set. For each data set, we consider optimizing the default random forest model estimated in R which uses 500 trees (the parameter ntree in randomForest is set to 500). For each data set, we also consider solving the tree ensemble problem using only the first $T$ trees of the complete forest, where $T$ ranges in $\{10,50,100,200\}$. For each data set and each value of $T$, we solve the MIO formulation (2), as well as its LO relaxation.

We compare our MIO formulation against two other approaches:

Local search: We solve the tree ensemble problem (1) using a local search heuristic. The details of this local search are provided in Section 10 of the e-companion; at a high level, it starts from a randomly chosen initial solution and iteratively improves the solution, one independent variable at a time, until a local optimum is reached. The heuristic is repeated from ten starting points, out of which we only retain the best (highest objective value) solution. We test such an approach to establish the value of our MIO-based approach, which obtains a globally optimal solution, as opposed to a locally optimal solution. 2. 2.

Standard linearization MIO: We solve the standard linearization MIO (4) and its relaxation, in order to obtain a relative sense of the strength of formulation (2). Due to this formulation being much harder to solve, we impose a 30 minute time limit on the solution time of the integer formulation.

We consider several metrics:

• •

$N_{Levels}$: the number of levels (i.e., dimension of $\mathbf{x}$), defined as $N_{Levels}=\sum_{i=1}^{n}K_{i}$.

• •

$N_{Leaves}$: the number of leaves (i.e., dimension of $\mathbf{y}$), defined as $N_{Leaves}=\sum_{t=1}^{T}|\mathbf{leaves}(t)|$.

• •

$\mathcal{T}_{MIO}$: the time (in seconds) to solve our MIO (2).

• •

$\mathcal{T}_{StdLin,MIO}$: the time (in seconds) to solve the standard linearization MIO (4).

• •

$\mathcal{T}_{LS}$: the time (in seconds) to run the local search procedure (value reported is the total for ten starting points).

• •

$G_{LS}$: the gap of the local search solution; if $Z_{LS}$ is the objective value of the local search solution and $Z^{*}$ is the optimal objective value of (2), then

[TABLE]

• •

$G_{LO}$: the gap of the LO relaxation of (2); if $Z_{LO}$ is the objective value of the LO relaxation and $Z^{*}$ is the optimal integer objective as before, then

[TABLE]

• •

$G_{StdLin,LO}$: the gap of the LO relaxation of the standard linearization MIO (4); if $Z_{StdLin,LO}$ is the optimal value of the relaxation, then

[TABLE]

• •

$G_{StdLin,MIO}$: the optimality gap of the standard linearization MIO (4); if $Z_{StdLin,UB}$ and $Z_{StdLin,LB}$ are the best upper and lower bounds, respectively, of problem (4) upon termination, then

[TABLE]

Note that all $\mathcal{T}$ metrics do not include the random forest estimation time in R.

Table 2 shows solution times and problem size metrics, while Table 3 shows the gap metrics. From these two tables, we can draw several conclusions. First, the time required to solve (2) is very reasonable; in the most extreme case (winequalityred, $T=500$), (2) can be solved to full optimality in about 20 minutes. (Note that no time limit was imposed on (2); all values of $\mathcal{T}_{MIO}$ correspond to the time required to solve (2) to full optimality.) In contrast, the standard linearization problem (4) was only solved to full optimality in two out of twenty cases within the 30 minute time limit. In addition, for those instances where the solver reached the time limit, the optimality gap of the final integer solution, $G_{StdLin,MIO}$, was quite poor, ranging from 50 to over 100%.

Second, the integrality gap $G_{LO}$ is quite small – on the order of a few percent in most cases. This suggests that the LO relaxation of (2) is quite tight. In contrast, the LO relaxation bound from (4) is weaker than that of (2), as predicted by Proposition 3.2, and strikingly so. The weakness of the relaxation explains why the corresponding integer problem cannot be solved to a low optimality gap within the 30 minute time limit. These results, together with the results above on the MIO solution times and the final optimality gaps of (4), show the advantages of our formulation (2) over the standard linearization formulation (4).

Third, although there are many cases where the local search solution performs quite well, there are many where it can be quite suboptimal, even when repeated with ten starting points. Moreover, while the local search time $T_{LS}$ is generally smaller than the MIO time $T_{MIO}$, in some cases it is not substantially lower (for example, solubility for $T=500$), and the additional time required by the MIO formulation (2) may therefore be justified for the guarantee of provable optimality.

5.3 Depth $d$ approximation experiments

In this section, we investigate the use of the depth $d$ tree problem (formulation (5) with $\bar{\Omega}_{d}$) for approximating the full depth problem (2). We focus on the same data sets as before with $T=100$. We solve (5) with $\bar{\Omega}_{d}$ and vary the depth $d$ of the approximation. We consider the upper bound $Z^{*}_{MIO,d}$ (denoted by “UB”), the actual value of the solution $Z_{d}$ (denoted by “Actual”) and the lower bound $Z^{*}_{MIO,d}-\sum_{t=1}^{T}\lambda_{t}\Delta_{t}$ (denoted by “LB”).

Figures 2 and 3 plot the above three metrics for the winequalityred and concrete data sets, respectively. From these plots, we can see that the upper bound is decreasing, while the lower bound and the actual objective are increasing. We can also see that the lower bound is quite loose, and the depth $d$ needs to increase significantly in order for the lower bound to be close to the upper bound. However, even when the depth $d$ is small and the lower bound is loose, the actual objective of the solution produced by the approximation is very good. In the case of winequalityred, the solution is essentially optimal ($Z*_{MIO,d}$ and $Z_{d}$ are close to or equal) after a depth of $d=15$ (compared to a maximum depth of 26); for concrete, this occurs for a depth of $d=9$ (compared to a maximum depth of 24).

To complement these results on the objective values, Figures 4 and 5 show the computation time of the depth approximation formulation as $d$ varies for the same two data sets. Here we can see that the solution time required to solve the depth approximation formulation initially increases in an exponential fashion as $d$ increases; this is to be expected, because with each additional layer of splits, the number of splits roughly doubles. Interestingly, though, the solution time seems to plateau after a certain depth, and no longer continues to increase. Together with Figures 2 and 3, these plots suggest the potential of the depth approximation approach to obtain near-optimal and optimal solutions with significantly reduced computation time relative to the full depth problem.

5.4 Solution method experiments

In this final set of experiments, we evaluate the effectiveness of the two solution methods from Section 4 – Benders decomposition and split constraint generation – on solving large instances of (2). We use the same data sets as before with $T=500$. For each instance, we consider $\mathcal{T}_{Benders}$, $\mathcal{T}_{SplitGen}$ and $\mathcal{T}_{Direct}$, which are the times to solve (2) to full optimality using the Benders, split constraint generation and direct solution approaches, respectively.

Table 4 shows the results from this comparison. Both approaches can lead to dramatic reductions in the solution time relative to solving (2) directly with all split constraints enforced at the start. In the most extreme case (concrete), we observe a reduction from about 800 seconds for the standard solution method to about 32 seconds for split constraint generation and about 37 seconds for the Benders approach – a reduction in solution time of over 95%. In some cases, Benders decomposition is slightly faster than split generation; for example, for winequalityred with $T=500$, the Benders approach requires just under 11 minutes whereas split generation requires just over 13 minutes. In other cases, split generation is faster (for example, solubility with $T=500$).

6 Case study 1: drug design

In this section, we describe our case study in drug design. Section 6.1 provides the background on the problem and the data. Section 6.2 shows results on the unconstrained optimization problem, while Section 6.3 shows results for when the similarity to existing compounds is constrained to be small.

6.1 Background

For this set of experiments, we use the data sets from Ma et al. (2015). These data sets were created for a competition sponsored by Merck and hosted by Kaggle. There are fifteen different data sets. In each data set, each observation corresponds to a compound/molecule. Each data set has a single dependent variable, which is different in each data set and represents some measure of “activity” (a property of the molecule or the performance of the molecule for some task). The independent variables in each data set are the so-called “atom pair” and “donor-acceptor pair” features, which describe the substructure of each molecule (see Ma et al. 2015 for further details). The goal of the competition was to develop a model to predict activity using the molecular substructure; such models are known as quantitative structure-activity relationship (QSAR) models.

The optimization problem that we will consider is to find the molecule that maximizes activity as predicted by a random forest model. Our interest in this problem is two-fold. First, this is a problem of significant practical interest, as new drugs are extremely costly to develop. Moreover, these costs are rising: the number of drugs approved by the FDA per billion dollars of pharmaceutical R&D spending has been decreasing by about 50% every 10 years since 1950 (a phenomenon known as “Eroom’s Law” – Moore’s Law backwards; see Scannell et al. 2012). As a result, there has been growing interest in using analytics to identify promising drug candidates in academia as well as industry (see for example Atomwise Inc. 2017). We note that random forests are widely used in this domain: the QSAR community was one of the first to adopt them (Svetnik et al. 2003) and they have been considered a “gold standard” in QSAR modeling (Ma et al. 2015).

Second, the problem is of a very large scale. The number of independent variables ranges from about 4000 to just under 10,000, while the number of observations ranges from about 1500 to just over 37,000; in terms of file size, the smallest data set is approximately 15MB, while the largest is just over 700MB. Table 5 summarizes the data sets. Estimating a random forest model using the conventional randomForest package in R on any one of these data sets is a daunting task; to give an example, a single tree on the largest data set requires more than 5 minutes to estimate, which extrapolates to a total computation time of over 8 hours for a realistic forest of 100 trees.

6.2 Unconstrained optimization results

In our first set of experiments, we proceed as follows. For each data set, we estimate a random forest model to predict the activity variable using all available independent variables. To reduce the computational burden posed by estimating random forest models from such large data sets, we deviate from our previous experiments by using the ranger package in R (Wright and Ziegler 2017), which is a faster implementation of the random forest algorithm suited for high dimensional data sets. In addition, we follow Ma et al. (2015) in restricting the number of trees to 100. For each such random forest model, we solve the (unconstrained) tree ensemble optimization problem (2) using the Benders approach of Section 4.1 and the split generation approach of Section 4.2, as well as directly using Gurobi. We impose a time limit of two hours. We also solve each tree ensemble optimization problem using local search with ten repetitions. We consider the following metrics:

• •

$G_{Direct},G_{SplitGen},G_{Benders}$: the optimality gap of the solution produced by solving the problem directly using Gurobi, solving it using the split generation method and solving it using the Benders methods, respectively. If $Z_{m,LB}$ is the lower bound of a method $m$ and $Z_{m,UB}$ is the upper bound, then $G_{m}$ are defined as

[TABLE]

• •

$G_{LS}$: the optimality gap of the local search solution, relative to the best split generation solution. If $Z_{LS}$ is the local search objective, it is defined as

[TABLE]

• •

$\mathcal{T}_{Direct}$ $\mathcal{T}_{SplitGen},\mathcal{T}_{Benders}$: the time (in seconds) to solve (2) directly using Gurobi, using the split generation method and using the Benders method, respectively. (A time that is below 7200 indicates that the problem was solved to full optimality.)

• •

$\mathcal{T}_{LS}$: the time (in seconds) to execute the local search procedure. (The time reported is the total of ten repetitions.)

• •

$N_{Levels}$ and $N_{Leaves}$: the number of levels and the number of leaves in the ensemble, respectively, defined as in Section 5.2.

Table 6 displays the results of this experiment. We first discuss the direct solution approach. For this approach, we can see that out of the fifteen data sets, nine were solved to full optimality within one hour, and one more data set was solved within the next hour. For the remaining five data sets (2, 3, 6, 8 and 13), the solver terminated after two hours with very low optimality gaps (four of the data sets having an optimality gap of below 0.5%, and one with an optimality gap of 2.4%). The high optimality gap for set 6 is to be expected, as this data set is among the two largest data sets in terms of the number of levels and the total number of leaves (which stems from the number of variables and the number of observations in that data set; see Table 5). With regard to the split generation approach, we can see that it improves on the direct solution method; with split generation, ten data sets are solved to full optimality within one hour, and another three are solved within the next hour. The greatest improvement is for data set 10, where the direct approach required over an hour, but split generation terminates in just over 20 minutes.

For the Benders approach, we can see that the performance is quite different. The optimality gap is substantially higher than that achieved by both the direct approach and split generation after two hours. The Benders approach is only able to solve two instances to full optimality within the two hour time limit and in both instances, the split generation approach is able to solve the same instance to full optimality more quickly.

The last important insight from Table 6 concerns the performance of the local search procedure. With regard to solution times, we can see that in some cases the total time required for the ten repetitions of the local search exceeds the time required for split generation (see data set 1 for example). In addition, and more importantly, the best solution obtained by local search in each data set is highly suboptimal, as evidenced by the high values of $G_{LS}$. In the best case (data set 8), $G_{LS}$ is about 5%, whereas in the worst case (data set 13), $G_{LS}$ is as high as 26%. These results suggest that local search is not adequate for this problem: our approaches, which are provably optimal and based on mixed-integer optimization, deliver significantly better solutions.

6.3 Controlling proximity

In the random forest literature, one concept that is useful for analyzing random forest models is that of proximity. The proximity of two observations $\mathbf{X},\mathbf{X}^{\prime}\in\mathcal{X}$ is defined as the proportion of trees for which $\mathbf{X}$ and $\mathbf{X}^{\prime}$ fall in the same leaf:

[TABLE]

where $\ell_{t}(\mathbf{X})$ is the leaf to which tree $t$ maps the observation $\mathbf{X}$. Proximity is valuable to consider because it allows one to use the random forest model to determine how similar two observations are, based on how each tree categorizes those two observations. In Svetnik et al. (2003), proximity was used as a similarity metric for clustering molecules, and was shown to be competitive with Tanimoto similarity (a different metric, commonly used in QSAR modeling; see Willett et al. 1998) in producing meaningful molecule clusters. In the remainder of this section, we first analyze the proximity of the solutions from Section 6.2, and verify that these solutions are distinct from the molecules in the data sets. Then, we consider solving (2) with constraints on the proximity. By imposing upper bounds on proximity, one can find molecules that maximize predicted activity while limiting how similar they are to existing molecules or equivalently, requiring them to be dissimilar to existing molecules; such a constrained optimization approach may potentially be useful in identifying novel molecules.

Table 7 displays, for each data set, the average proximity $\pi_{\mathrm{avg}}$ and maximum proximity $\pi_{\max}$ between the split generation solution and all of the molecules in the data set. We can see that the maximum proximity – the highest proximity between the solution and any point in the training data – is in general low. For example, for data set 3, $\pi_{\max}$ is 0.11: for any point in the data set that we consider, at most eleven trees out of the 100 trees in the forest will place both the training point and our solution in the same leaf. In addition, the average proximity is much lower than the maximum proximity, which suggests that for most training set points, the actual proximity to the solution is close to zero (i.e., there is very little similarity between the solution and the training point). Most importantly, all solutions have a $\pi_{\max}$ strictly lower than 1, indicating that all fifteen solutions are in fact different from all of the molecules of their respective training data sets.

We now turn our attention to solving (2) with an added constraint on the maximum proximity of the solution to the training points. Such a constraint is defined as follows: let $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(M)}$ be the set of observations in the training data, and define for each observation $m$ the vector $\mathbf{y}^{(m)}\in\mathbb{R}^{\sum_{t=1}^{T}|\mathbf{leaves}(t)|}$ as

[TABLE]

The proximity between the fixed observation $\mathbf{X}^{(m)}$ and the solution encoded by $(\mathbf{x},\mathbf{y})$ can then be written as an affine function of $\mathbf{y}$:

[TABLE]

(Note that $\mathbf{y}^{(1)},\dots\mathbf{y}^{(M)}$ are data, and not decision variables.) We can thus enforce a constraint on the proximity of the solution encoded by $\mathbf{x}$ to each observation $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(M)}$ to be at most $c\in[0,1]$ through the following family of linear constraints on $\mathbf{y}$:

[TABLE]

We solve (2) with constraint (11) and vary the parameter $c$ to generate a Pareto efficient frontier of solutions that optimally trade-off their maximum proximity to the training data and their predicted value under the random forest model. We impose a time limit of one hour on each solve of (2). We solve each constrained instance using the split generation approach.

To demonstrate, we focus on data sets 4 and 7. Figure 6 shows the proximity-objective value frontier (the points labeled “MIO”) for data sets 4 and 7. (Note that the objective value is expressed in terms of the maximum unconstrained objective value, i.e., the objective value attained when the proximity constraint is omitted.) We note that the right-most point in each frontier (maximum proximities of 0.71 for data set 4 and 0.43 for data set 7) corresponds to the original unconstrained solution. As a comparison, we also solve the (unconstrained) problem for each data set using local search with 100 repetitions, and plot the proximity and objective value of each solution from each repetition (the points labeled “LS”). We remark that both frontiers are approximate because the constrained optimization problem associated with each point is not necessarily solved to full optimality, due to the one hour time limit. However, upon closer examination, only six points out of a total of 112 across both plots did not solve to full optimality, with a maximum gap of only about 0.12%. Thus, although constraint (11) adds to the problem size, it does not drastically impact our ability to solve (2).

From both of these figures, we obtain several insights. First, in these two cases, we are able to push the proximity to its lowest; for both data sets, we can find solutions with maximum proximities of 0.01 (i.e., one tree out of the ensemble of 100). (A solution with maximum proximity of zero does not exist, because for any tree, the leaf $\ell_{t}(\mathbf{X})$ that the candidate solution $\mathbf{X}$ is mapped to will also be shared by at least one point in the training set.) Second, we can can see that the price of dissimilarity is low: as we decrease the maximum proximity, we can still obtain solutions with very good predicted performance. For example, for data set 4, we can see that if we lower the proximity to 0.01, the relative objective value decreases by only about 7%. Third, although the solutions obtained by local search have smaller maximum proximities than the unconstrained MIO solution, they are highly suboptimal with respect to objective value (best relative objectives for data sets 4 and 7 are roughly 94% and 90%, respectively) and are dominated in both maximum proximity and objective value by the MIO solutions. Overall, these results suggest that our MIO formulation can be used to systematically identify promising molecules that have good predicted performance and are sufficiently different from existing molecules in the training data.

7 Case study 2: customized pricing

In addition to our drug design case study, we also consider the problem of customized pricing. The results of this case study are fully reported in Section 11 of the e-companion. At a high level, we estimate random forest models of profits for eleven orange juice products as a function of their prices and store-level attributes, and then optimize these models to derive prices for each store. We also estimate log-linear and semi-log demand models in a hierarchical Bayesian (HB) framework, and heuristically optimize these models to obtain prices as well. We find that our random forest models can be optimized quickly (on the order of seconds). Moreover, the random forest-optimal prices perform well under the HB models, and are less extreme, in that fewer products are set to their highest or lowest allowable prices across the whole store chain. From a predictive standpoint, our random forest models also achieve significantly higher out-of-sample $R^{2}$ values than the HB models, which may be of independent interest.

8 Conclusion

In this paper, we developed a modern optimization approach to the problem of finding the decision that optimizes the prediction of a tree ensemble model. At the heart of our approach is a mixed-integer optimization formulation that models the action of each tree in the ensemble. We showed that this formulation is better than a general alternate formulation, that one can construct an hierarchy of approximations to the formulation with bounded approximation quality through depth-based truncation and that one can exploit the structure of the formulation to derive efficient solution methods, based on Benders decomposition and split constraint generation. We demonstrated the utility of our approach using real data sets, including two case studies in drug design and customized pricing. Given the prevalence of tree ensemble models, we believe that this methodology will become an important asset in the modern business analytics toolbox and is an exciting starting point for future research at the intersection of optimization and machine learning.

\ACKNOWLEDGMENT

The author sincerely thanks the area editor Marina Epelman, the associate editor and the two anonymous referees for their careful reading of the paper and thoughtful comments that have helped to greatly improve the quality of the paper. The author thanks Fernanda Bravo, Vishal Gupta and Auyon Siddiq for helpful conversations. The authors of Ma et al. (2015) and Merck Research Laboratories are gratefully acknowledged for the data set used in the drug design case study of Section 6. Dominick’s Finer Foods and Peter Rossi are gratefully acknowledged for the data set used in the customized pricing case study of Section 11.

9 Proofs

9.1 Auxiliary results

Before proceeding to the proofs, we first state two auxiliary results.

Lemma 9.1

Fix a tree $t$ and a leaf $\ell\in\mathbf{leaves}(t)$. Then

[TABLE]

Furthermore, define the collection $\mathcal{S}=\{\mathbf{right}(s)\,|\,s\in\mathbf{LS}(\ell)\}\cup\{\mathbf{left}(s)\,|\,s\in\mathbf{RS}(\ell)\}$. Then each pair of distinct sets $A,B\in\mathcal{S}$, $A\neq B$, is disjoint.

To gain some intuition for the right hand set in equation (12), recall that $\mathbf{LS}(\ell)$ is the set of splits for which we follow the left branch in order to reach $\ell$, and $\mathbf{RS}(\ell)$ is the set of splits for which we follow the right branch to reach $\ell$. For each $s\in\mathbf{LS}(\ell)$, $\mathbf{right}(s)$ is the set of leaves that is on the “wrong side” of split $s$ (we take the left branch to reach $\ell$, but each leaf in $\mathbf{right}(s)$ is only reachable by taking the right branch). Similarly, for each $s\in\mathbf{RS}(\ell)$, $\mathbf{left}(s)$ is the set of leaves that is on the wrong side of split $s$. The union of all leaves $\ell^{\prime}$ on the wrong side of each split $s\in\mathbf{RS}(\ell)\cup\mathbf{LS}(\ell)$ covers all leaves except $\ell$.

Proof 9.2

Proof of Lemma 9.1: We prove this by establishing two inclusions:

[TABLE]

For the first inclusion, let us fix $\ell^{\prime}\in\mathbf{leaves}(t)$ such that $\ell^{\prime}\neq\ell$. Then there must exist a split $\bar{s}\in\mathbf{splits}(t)$ such that either $\ell^{\prime}\in\mathbf{left}(\bar{s})$, $\ell\in\mathbf{right}(\bar{s})$, or $\ell^{\prime}\in\mathbf{right}(\bar{s})$, $\ell\in\mathbf{left}(\bar{s})$. In the former case, we have that $\bar{s}\in\mathbf{RS}(\ell)$, so that

[TABLE]

In the latter case, we have that $\bar{s}\in\mathbf{LS}(\ell)$, so that

[TABLE]

This proves the first inclusion.

For the second inclusion, we will argue the contrapositive. We have that

[TABLE]

where the first step follows by De Morgan’s law; the second follows by the fact that for any split $s\in\mathbf{splits}(t)$, $\mathbf{left}(s)$ and $\mathbf{right}(s)$ are disjoint; and the last by the definition of $\mathbf{LS}(\ell)$ and $\mathbf{RS}(\ell)$. This proves the second inclusion, and thus proves the equivalence.

Finally, to show that the sets in $\mathcal{S}$ are pairwise disjoint, order the splits in $\mathbf{LS}(\ell)\cup\mathbf{RS}(\ell)$ according to their depth:

[TABLE]

where $K$ is the total number of splits in $\mathbf{LS}(\ell)\cup\mathbf{RS}(\ell)$. Let us also define the sequence of sets $A_{1},A_{2},\dots,A_{K}$ as

[TABLE]

and the sequence of sets $B_{1},B_{2},\dots,B_{K}$ as

[TABLE]

We need to show that the collection $\{B_{1},\dots,B_{K}\}$ (this is the collection $\mathcal{S}$) is disjoint. Observe that by the definition of $\mathbf{LS}(\ell)$ and $\mathbf{RS}(\ell)$, the sets $A_{1},A_{2},\dots,A_{K}$ form a nested sequence, i.e.,

[TABLE]

Notice also that for each $i\in\{2,\dots,K\}$,

[TABLE]

and for each $i\in\{1,\dots,K\}$,

[TABLE]

It therefore follows that given $i,j\in\{1,\dots,K\}$ with $i<j$, that

[TABLE]

which establishes that $\{B_{1},\dots,B_{K}\}$ are pairwise disjoint. $\square$

In addition to Lemma 9.1, it will also be useful to state an analogous lemma for splits. With a slight abuse of notation, let us define $\mathbf{LS}(s)$ for a split $s\in\mathbf{splits}(t)$ as the sets of splits $s^{\prime}$ such that $s$ is on the left subtree of $s^{\prime}$; similarly, we define $\mathbf{RS}(s)$ for a split $s\in\mathbf{splits}(t)$ as the set of splits $s^{\prime}$ such that $s$ is on the right subtree of $s^{\prime}$. We then have the following lemma; the proof follows along similar lines to Lemma 9.1 and is omitted.

Lemma 9.3

For a given tree $t$, let $s\in\mathbf{splits}(t)$. We then have

[TABLE]

9.2 Proof of Proposition 3.1

To prove that problem (1) is NP-Hard, we will show that it can be used to solve the minimum vertex cover problem. An instance of the minimum vertex cover problem is defined by a graph $(V,E)$, where $V$ is a set of vertices and $E$ is the set of edges between these vertices. The minimum vertex cover problem is to find the smallest set of vertices $S$ from $V$ such that each edge in $E$ is incident to at least one vertex from $S$.

We now show how to cast this problem as a tree ensemble optimization problem. For convenience, let us index the vertices from $1$ to $|V|$ and the edges from $1$ to $|E|$. Also, let $e_{1}$ and $e_{2}$ be the nodes to which edge $e\in E$ is incident to. Suppose that our independent variable $\mathbf{X}$ is given by $\mathbf{X}=(X_{1},\dots,X_{|V|})$, where $X_{i}$ is a numeric variable that is 1 or 0. The tree ensemble we will consider will consist of the following two types of trees:

Type 1*: Trees $1$ to $|V|$, where for $i\in\{1,\dots,|V|\}$,*

[TABLE] 2. 2.

Type 2*: Trees $|V|+1$ to $|V|+|E|$, where for each edge $e\in E$,*

[TABLE]

These two types of trees are visualized in Figure 7.

We let the weight $\lambda_{t}$ of each tree $t$ be $-1$. The corresponding tree ensemble optimization problem is

[TABLE]

The above problem is identical to the minimum vertex cover problem. In particular, the independent variable $\mathbf{X}$ encodes the cover; $X_{i}=1$ indicates that vertex $i$ is part of the cover, and $X_{i}=0$ indicates that vertex $i$ is not in the cover. The type 1 trees count the size of the cover, while the type 2 trees penalize the solution if an edge is not covered by the set. More precisely, to understand the role of the type 2 trees, observe that:

• •

If the set of vertices encoded by $\mathbf{X}$ is a feasible cover, then $\sum_{t=|V|+1}^{|V|+|E|}f_{t}(\mathbf{X})=0$, and the objective only consists $-\sum_{t=1}^{|V|}f_{t}(\mathbf{X})$, which counts the number of vertices in the set encoded by $\mathbf{X}$.

• •

If the set of vertices encoded by $\mathbf{X}$ is not a feasible cover, then $f_{t}(\mathbf{X})=|V|+1$ for at least one $t\in\{|V|+1,\dots,|V|+|E|\}$, and therefore the objective satisfies the bound

[TABLE]

Observe that the bound in inequality (14) is strictly worse than selecting all of the vertices, that is, setting $X_{1}=X_{2}=\dots=X_{|V|}=1$; using all of the vertices corresponds to an objective value of $-|V|$. Therefore, at optimality, the set of vertices encoded by $\mathbf{X}$ must be a feasible cover. As stated above, the objective value of $\mathbf{X}$ when it corresponds to a feasible cover reduces to

[TABLE]

which is (the negative of) the size of the set of vertices encoded by $\mathbf{X}$. Maximizing this quantity is equivalent to minimizing its negative, which is the same as minimizing the size of the set of vertices that covers $E$. Therefore, solving (13) is equivalent to solving the minimum vertex cover problem for $(V,E)$.

*Since the minimum vertex cover problem is NP-Complete [Garey and Johnson, 1979], it follows that the tree ensemble optimization problem (1) is NP-Hard. *$\square$

9.3 Proof of Proposition 3.2

To prove the proposition, we will show that any optimal solution $(\mathbf{x},\mathbf{y})$ of the relaxation of (2) is a feasible solution of the relaxation of the standard linearization problem (4). Since the objective functions of problems (2) and (4) are the same, it will follow that the objective of $(\mathbf{x},\mathbf{y})$, which is $Z^{*}_{LO}$, is less than or equal to $Z^{*}_{LO,StdLin}$, which is the optimal objective value of the relaxation of (4).

Let $(\mathbf{x},\mathbf{y})$ be an optimal solution of the relaxation of problem (2). To show it is feasible for the relaxation of (4), we need to show that it satisfies the constraints of that formulation. We only need to show that constraints (4b) – (4d) are satisfied, since the other constraints of (4) are the same as in (2).

To verify constraint (4b), observe that $(\mathbf{x},\mathbf{y})$ satisfies constraint (2c). For any $t\in\{1,\dots,T\}$, $\ell\in\mathbf{leaves}(t)$ and $s\in\mathbf{LS}(\ell)$, we have

[TABLE]

where the first inequality is exactly constraint (2c), and the second inequality follows because $\ell\in\mathbf{left}(s)$ (this is because $s\in\mathbf{LS}(\ell)$) and all $y_{t,\ell^{\prime}}$’s are nonnegative (by constraint (2h)). Therefore, $(\mathbf{x},\mathbf{y})$ satisfies constraint (4b). Similar reasoning can be used to establish constraint (4c).

To verify constraint (4d), observe that $(\mathbf{x},\mathbf{y})$ satisfies

[TABLE]

for any tree $t$, by virtue of constraint (2b). Fix a tree $t$ and a leaf $\ell\in\mathbf{leaves}(t)$, and re-arrange the above to obtain

[TABLE]

We then have

[TABLE]

*where the first inequality is just inequality (15) from earlier, and the first equality follows from Lemma 9.1; the second inequality follows by constraints (2c) and (2d); and the last two equalities follow by simple algebra. This establishes that $(\mathbf{x},\mathbf{y})$ satisfies constraint (4d). This establishes that $(\mathbf{x},\mathbf{y})$ is feasible for (4), which concludes the proof. *$\square$

9.4 Proof of Proposition 3.3

Observe that the sets of tree-split pairs are nested in the following way:

[TABLE]

As a consequence, the feasible region of (5) at depth $d$ is a superset of the feasible region of (5) at depth $d+1$ and so we have

[TABLE]

*The equality $Z^{*}_{MIO,d_{\max}}=Z^{*}_{MIO}$ follows by the definition of $d_{\max}$ as the maximum depth of any tree in the ensemble. Combining this equality with the above sequence of inequalities establishes the result. *$\square$

9.5 Proof of Theorem 3.4

Let $(\mathbf{x},\mathbf{y})$ be an optimal solution of the depth $d$ problem (i.e., problem (5) with $\bar{\Omega}_{d}$). Let $Z^{*}_{MIO,d}$ be the objective value of $(\mathbf{x},\mathbf{y})$ within problem (5) with $\bar{\Omega}_{d}$.

Let $(\mathbf{x},\tilde{\mathbf{y}})$ be the solution for (2) (the full-depth problem), obtained by finding the unique value of $\tilde{\mathbf{y}}$ such that $(\mathbf{x},\tilde{\mathbf{y}})$ is feasible for (2). The existence and uniqueness of such a $\tilde{\mathbf{y}}$ is guaranteed by Proposition 4.1. Let $Z_{d}$ be the objective value of $(\bar{\mathbf{x}},\tilde{\mathbf{y}})$ within (2).

For a given tree $t\in\{1,\dots,T\}$, let us consider the difference of the prediction of tree $t$ for $(\mathbf{x},\mathbf{y})$ and $(\mathbf{x},\tilde{\mathbf{y}})$:

[TABLE]

In order to understand this quantity, we need to understand which part of the tree $\bar{\mathbf{x}}$ will get mapped to. We will do this through the following procedure:

Initialize $\nu$ to the root node of the tree. 2. 2.

If $\nu\in\mathbf{leaves}(t)$ or $\nu\in\mathbf{splits}(t,d)$, stop; otherwise:

• •

If $\sum_{j\in\mathbf{C}(\nu)}x_{\mathbf{V}(\nu),j}=1$, set $\nu$ to its left child node;

• •

If $\sum_{j\in\mathbf{C}(\nu)}x_{\mathbf{V}(\nu),j}=0$, set $\nu$ to its right child node;

and repeat step 2.

Upon termination, $\nu$ is some node in the tree – either a leaf or a depth $d$ split node. Regardless of the type of node, we know that for any $s\in\mathbf{LS}(\nu)$, the depth of $s$ is in $\{1,\dots,d\}$, and so $(\mathbf{x},\mathbf{y})$ satisfies the right split constraint (5c) of the depth $d$ problem for $s$. By the definition of the procedure, it also must be that $\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=1$, which implies that

[TABLE]

For $s\in\mathbf{RS}(\nu)$, similar reasoning using the left split constraint (5b) allows us to assert that for any $s\in\mathbf{RS}(\nu)$,

[TABLE]

We thus know the following about $\mathbf{y}$:

[TABLE]

We can also assert the same about $\tilde{\mathbf{y}}$, since $(\mathbf{x},\tilde{\mathbf{y}})$ satisfies the left and right split constraints (2c) and (2d) at all depths:

[TABLE]

We now consider three possible cases for the type of node $\nu$ is:

Case 1*: $\nu\in\mathbf{leaves}(t)$. In this case, by Lemma 9.1, we can assert that*

[TABLE]

[TABLE]

Since both $\mathbf{y}$ and $\tilde{\mathbf{y}}$ are nonnegative and sum to one, it follows that $y_{t,\nu}=\tilde{y}_{t,\nu}=1$. We therefore have that the prediction difference (16) is simply

[TABLE] 2. 2.

Case 2*: $\nu\in\mathbf{splits}(t,d)$ and $\sum_{j\in\mathbf{C}(\nu)}x_{\mathbf{V}(\nu),j}=1$. In this case, by Lemma 9.3, we have that*

[TABLE]

In addition, since $\sum_{j\in\mathbf{C}(\nu)}x_{\mathbf{V}(\nu),j}=1$, then by the right split constraints, we additionally have $y_{t,\ell}=0$ and $\tilde{y}_{t,\ell}=0$ for all $\ell\in\mathbf{right}(\nu)$, which implies that

[TABLE]

We can use the above properties of $\mathbf{y}_{t}$ and $\tilde{\mathbf{y}}_{t}$ to bound the prediction difference (16) as follows:

[TABLE] 3. 3.

Case 3*: $\nu\in\mathbf{splits}(t,d)$ and $\sum_{j\in\mathbf{C}(\nu)}x_{\mathbf{V}(\nu),j}=0$. Similar reasoning as in case 2 can be used to establish that*

[TABLE]

and to bound the prediction difference (16) as

[TABLE]

Given cases 2 and 3, observe that if we know that $\nu\in\mathbf{splits}(t,d)$, then a valid upper bound on the prediction difference is simply

[TABLE]

Now, if we do not know what type of node $\nu$ is – whether it is a leaf, or which split in $\mathbf{splits}(t,d)$ it is – then we can construct an upper bound, based on cases 1, 2 and 3 above, for the prediction difference as

[TABLE]

where the maximum is defined to be zero if $\mathbf{splits}(t,d)$ is empty. (Note that in the case $\nu$ is a leaf, the above bound is valid, since all $\delta_{t,\nu}$ values are nonnegative by definition.)

Let us now unfix the tree $t$. Applying the above bound to bound the prediction difference of all trees $t$, and using the fact that $\lambda_{t}\geq 0$ for all $t$, it follows that the difference between the objective $Z^{*}_{MIO,d}$ of $(\mathbf{x},\mathbf{y})$ and the objective $Z_{d}$ of $(\mathbf{x},\tilde{\mathbf{y}})$ can be written as

[TABLE]

From here, it immediately follows that

[TABLE]

Combining this with the fact that $(\mathbf{x},\tilde{\mathbf{y}})$ is a feasible solution for (2) and Proposition 3.3 leads to the inequality,

[TABLE]

*as required. *$\square$

9.6 Proof of Proposition 4.1

Feasibility*. Let us first show that the proposed solution $\mathbf{y}_{t}$ is feasible. By construction, we have that $y_{t,\ell}\geq 0$ for all $\ell\in\mathbf{leaves}(t)$ and that $\sum_{\ell\in\mathbf{leaves}(t)}y_{t,\ell}=1$, so constraints (8b) and (8e) are satisfied. This leaves the left and right split constraints (8c) and (8d).*

For constraint (8c), let $s\in\mathbf{splits}(t)$. If $\ell^{*}\in\mathbf{left}(s)$, then it must be that $s\in\mathbf{LS}(\ell^{*})$ and so by the definition of GetLeaf, it must be that $\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=1$. Therefore, we have:

[TABLE]

Otherwise, if $\ell^{*}\notin\mathbf{left}(s)$, then

[TABLE]

which is automatically less than or equal to $\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}$ (the latter can only be 0 or 1).

For constraint (8d), let $s\in\mathbf{splits}(t)$. If $\ell^{*}\in\mathbf{right}(s)$, then $s\in\mathbf{RS}(\ell^{*})$ and by the definition of GetLeaf, it must be that $1-\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=1$. Therefore, applying similar reasoning as above, we get

[TABLE]

*Otherwise, if $\ell^{*}\notin\mathbf{right}(s)$, then again, we have $\sum_{\ell\in\mathbf{left}(s)}y_{t,\ell}=0$, which is automatically less than or equal to $1-\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}$ (again, it can only be 0 or 1). This establishes that $\mathbf{y}_{t}$ is a feasible solution to the subproblem (8). *

Unique feasible solution*. To establish that the proposed solution is the only feasible solution, we proceed as follows. We will show that if a solution $\mathbf{y}_{t}$ is a feasible solution of (8), then it must be equal to the solution of the statement of the proposition.*

Let $\ell\in\mathbf{leaves}(t)$ such that $\ell\neq\ell^{*}$. Then by Lemma 9.1, we have that

[TABLE]

Moreover, the collection of sets in the union above is disjoint. Therefore, either $\ell\in\mathbf{right}(\bar{s})$ for some $\bar{s}\in\mathbf{LS}(\ell^{*})$ or $\ell\in\mathbf{left}(\bar{s})$ for some $\bar{s}\in\mathbf{RS}(\ell^{*})$.

In the former case – that is, $\ell\in\mathbf{right}(\bar{s})$ for some $\bar{s}\in\mathbf{LS}(\ell^{*})$ – we have by constraint (8c) and constraint (8e) that

[TABLE]

Therefore, $y_{t,\ell}$ is upper bounded by $1-\sum_{j\in\mathbf{C}(\bar{s})}x_{\mathbf{V}(\bar{s}),j}$ and lower bounded by 0 (by constraint (8e)). Since $\bar{s}\in\mathbf{LS}(\ell^{*})$ and from the definition of GetLeaf, it must be that $\sum_{j\in\mathbf{C}(\bar{s})}x_{\mathbf{V}(\bar{s}),j}=1$, or equivalently, $1-\sum_{j\in\mathbf{C}(\bar{s})}x_{\mathbf{V}(\bar{s}),j}=0$. Therefore, $y_{t,\ell}$ must be equal to zero.

Similarly, if $\ell\in\mathbf{left}(\bar{s})$ for some $\bar{s}\in\mathbf{RS}(\ell^{*})$, then by constraints (8d) and (8e) we have that

[TABLE]

Therefore, $y_{t,\ell}$ is upper bounded by $\sum_{j\in\mathbf{C}(\bar{s})}x_{\mathbf{V}(\bar{s}),j}$ and lower bounded by 0 (by (8e)). Since $\bar{s}\in\mathbf{LS}(\ell^{*})$ and from the definition of GetLeaf, it must be that $\sum_{j\in\mathbf{C}(\bar{s})}x_{\mathbf{V}(\bar{s}),j}=0$. Therefore, $y_{t,\ell}$ must be equal to zero.

*From the above reasoning, we have shown that $y_{t,\ell}=0$ for every $\ell\neq\ell^{*}$. By constraint (8b), it must be that $y_{t,\ell^{*}}=1-\sum_{\ell\neq\ell^{*}}y_{t,\ell}=1-0=1$. The resulting solution is therefore exactly the same as the one proposed in the proposition; it follows that the proposed solution $\mathbf{y}_{t}$ is the only feasible solution to (8). *

Optimality*. Since $\mathbf{y}_{t}$ is the only feasible solution of (8), it must also be its optimal solution. This completes the proof. *$\square$

9.7 Proof of Proposition 4.2

First, let us check the dual objective of the proposed solution $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$. We have

[TABLE]

where the first step follows by the definition of $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$; the second step follows by the fact that $\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=0$ for $s\in\mathbf{RS}(\ell^{*})$ and $1-\sum_{j\in\mathbf{C}(s)}x_{\mathbf{V}(s),j}=0$ for $s\in\mathbf{LS}(\ell^{*})$; and the last two steps by algebra. The final value is exactly equal to the optimal primal objective value. If $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ is feasible for the dual problem, then the proposition will be proven.

*To verify feasibility, observe that by their definition, we have $\alpha_{t,s}\geq 0$ and $\boldsymbol{\beta}_{t,s}\geq 0$ for all $s\in\mathbf{splits}(t)$. Thus, we only need to check constraint (9b) for each $\ell\in\mathbf{leaves}(t)$. We consider two cases: *

Case 1*: $\ell=\ell^{*}$. In this case, proceeding from the left hand side of constraint (9b) for $\ell^{*}$, we have*

[TABLE]

*where the first equality follows because $\alpha_{t,s}=0$ for all $s\notin\mathbf{RS}(\ell^{*})$ and $\beta_{t,s}=0$ for all $s\notin\mathbf{LS}(\ell^{*})$ (note that $\mathbf{LS}(\ell^{*})\cap\mathbf{RS}(\ell^{*})=\emptyset$; a leaf cannot be both to the left of and to the right of the same split), and also because $\gamma_{t}=p_{t,\ell^{*}}$ by definition. *

Case 2*: $\ell\neq\ell^{*}$. In this case, by Lemma (9.1), we know that $\ell$ satisfies*

[TABLE]

Lemma 9.1 states that each set in the above union is disjoint. Therefore, we have that $\ell\in\mathbf{right}(\bar{s})$ for some $\bar{s}\in\mathbf{LS}(\ell^{*})$ or $\ell\in\mathbf{left}(\bar{s})$ for some $\bar{s}\in\mathbf{RS}(\ell^{*})$. We now show that the inequality holds in either of these two scenarios.

If $\ell\in\mathbf{right}(\bar{s})$ for some $\bar{s}\in\mathbf{LS}(\ell^{*})$, then we have

[TABLE]

where the first step follows because $\bar{s}$ must belong to $\mathbf{RS}(\ell)$ and because by definition, $\boldsymbol{\alpha}_{t}$ and $\boldsymbol{\beta}_{t}$ are nonnegative; the second step follows by the definition of $\beta_{t,s}$ for $s\in\mathbf{LS}(\ell^{*})$; the third step by the definition of the maximum; and the last step by algebra.

If $\ell\in\mathbf{left}(\bar{s})$ for some $\bar{s}\in\mathbf{RS}(\ell^{*})$, then we have

[TABLE]

which follows by logic similar to the first case ($\ell\in\mathbf{right}(\bar{s})$* for some $\bar{s}\in\mathbf{LS}(\ell^{*})$).*

*Thus, we have established that $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ is a feasible solution for the dual (9) and achieves the same objective as the primal optimal solution. By weak LO duality, the solution $(\boldsymbol{\alpha}_{t},\boldsymbol{\beta}_{t},\gamma_{t})$ must therefore be optimal. *$\square$

9.8 Proof of Proposition 4.3

For the $\Rightarrow$ direction, the implication is immediate, because if $(\mathbf{x},\mathbf{y})$ satisfies constraints (2c) and (2d) for all $s\in\mathbf{splits}(t)$, it will also satisfy them for arbitrary subsets of $\mathbf{splits}(t)$.

Thus, we only need to establish the $\Leftarrow$ direction. Let $s\in\mathbf{LS}(\ell^{*})$. For any tree $t$, any $s\in\mathbf{LS}(\ell^{*})$ and any $\ell\in\mathbf{right}(s)$,

[TABLE]

where the first step follows since $y_{t,\ell^{\prime}}$ is nonnegative for any $\ell^{\prime}$, the second step by the hypothesis of the implication, and the remaining two steps by algebra. Since $y_{t,\ell}$ is nonnegative, it must be that $y_{t,\ell}=0$. Similarly, for any $s\in\mathbf{RS}(\ell^{*})$, we can show that for each $\ell\in\mathbf{left}(s)$, $y_{t,\ell}=0$.

It therefore follows that for any $\ell\in\bigcup_{s\in\mathbf{RS}(\ell^{*})}\mathbf{left}(s)\cup\bigcup_{s\in\mathbf{LS}(\ell^{*})}\mathbf{right}(s)$, $y_{t,\ell}$. Invoking Lemma 9.1, we have that $y_{t,\ell}=0$ for all $\ell\in\mathbf{leaves}(t)\setminus\{\ell^{*}\}$. Since $\mathbf{y}$ is assumed to satisfy constraint (2b), it must be that $y_{t,\ell^{*}}=1$.

From here, we can see that $\mathbf{y}_{t}$, the collection of $y_{t,\ell}$ values corresponding to tree $t$, is defined exactly as in the statement of Proposition 4.1. Thus, invoking Proposition 4.1, we can assert that $\mathbf{y}$ satisfies

[TABLE]

*for all $s\in\mathbf{splits}(t)$. This concludes the proof. *$\square$

10 Local search procedure

We provide the pseudocode of our local search procedure for approximately solving (1) below as Algorithm 3. Before we define the procedure, we define the set $\bar{\mathcal{X}}_{i}$ as

[TABLE]

for categorical variables, and

[TABLE]

*for numeric variables, where $a_{i,j}$ is the *$j$th smallest split point of variable $i$ in the tree ensemble model. The $\bar{\mathcal{X}}_{i}$ are simply the domains of each independent variable defined in a way that will be helpful in defining our local search. For numeric variables, $\bar{\mathcal{X}}_{i}$ consists of the $K_{i}$ split points of variable $i$ and one extra point, $a_{i,K_{i}}+1$. The extra point $a_{i,K_{i}}+1$ is arbitrary. We can use any value here, as long as it is strictly larger than the largest split point of variable $i$, as this will allow us to choose to be on the right-hand side of a split with query $X_{i}\leq a_{i,K_{i}}$.

11 Case study 2: customized pricing

In this section, we apply our approach to customized pricing. Section 11.1 provides the background on the data, while Section 11.2 describes our random forest model as well as two alternative models based on hierarchical Bayesian regression. Section 11.3 compares the models in terms of out-of-sample predictions of profit. Finally, Section 11.4 formulates the profit optimization problem and compares the three models.

11.1 Background

We consider the data set from Montgomery [1997]**, which was accessed via the bayesm package in R [Rossi, 2012]. This data set contains price and sales data for eleven different refrigerated orange juice brands for the Dominick’s Finer Foods chain of grocery stores in the Chicago area.

In this data set, each observation corresponds to a given store in the chain at a given week. The data span 83 stores and a period of 121 weeks. Each observation consists of: the week $t$; the store $s$; the sales $q_{t,s,1},\dots,q_{t,s,11}$ of the eleven orange juice brands; the prices $p_{t,s,1},\dots,p_{t,s,11}$ of the eleven orange juice brands; dummy variables $d_{t,s,1},\dots,d_{t,s,11}$, where $d_{t,s,i}=1$ if orange juice brand $i$ had any in-store displays (such as in-store coupons) at store $s$ in week $t$; and dummy variables $f_{t,s,1},\dots,f_{t,s,11}$, where $f_{t,s,i}=1$ if brand $i$ was featured/advertised in store $s$ in week $t$. We use $\mathbf{p}$, $\mathbf{d}$ and $\mathbf{f}$ to denote vectors of prices, deal dummies and feature dummies, and subscripts to denote the observation; for example, $\mathbf{p}_{t,s}=(p_{t,s,1},\dots,p_{t,s,11})$ is the vector of brand prices at store $s$ in week $t$.

The data set also include 11 covariates for the 83 stores corresponding to demographic information of each store’s neighborhood, such as the percentage of the population with a college degree and the percentage of households with more than five members; we denote these covariates as $z_{s,1},\dots,z_{s,11}$. We denote the vector of covariates for store $s$ as $\mathbf{z}_{s}$; for notational convenience later, we assume that $z_{s,0}=1$ in $\mathbf{z}_{s}$. For more details, we refer the reader to Montgomery [1997]**.

11.2 Models

The standard approach for modeling this type of data in marketing is to posit a regression model within a hierarchical Bayes (HB) framework [Rossi et al., 2005]. We will consider two different HB specifications. In the first model, which we will refer to as HB-LogLog, we assume that the logarithm of sales of a given product is linear in the logarithm of prices of all eleven products, and linear in the deal and feature dummies. For a fixed focal brand $i$, the regression function is given below:

[TABLE]

where $\epsilon_{t,s,i}$ follows a univariate normal distribution with mean 0 and variance $\tau_{s}$:

[TABLE]

and the vector of regression coefficients $\boldsymbol{\beta}_{s,i}=(\beta_{s,i,0},\boldsymbol{\beta}_{s,i,\mathbf{p}},\boldsymbol{\beta}_{s,i,\mathbf{d}},\boldsymbol{\beta}_{s,i,\mathbf{f}})$ follows a multivariate normal distribution with mean $\Delta_{i}^{T}\mathbf{z}_{s}$ and covariance matrix $V_{\boldsymbol{\beta},i}$:

[TABLE]

where $\Delta_{i}$ is a 12-by-34 matrix. (Each row corresponds to one of the store level covariates, and each column corresponds to one of the $\beta$ coefficients.) This model assumes that the regression coefficients of each store $\boldsymbol{\beta}_{s,i}$ follow a normal distribution whose mean depends on the store-specific covariates $\mathbf{z}_{s}$. Specifically, the mean of $\boldsymbol{\beta}_{s,i}$ is a linear function of $\mathbf{z}_{s}$; recall that $z_{s,0}=1$, so that the first row of $\Delta$ specifies the intercept of the mean of $\boldsymbol{\beta}_{s,i}$. We note that log-log models like equation (18) are commonly used in demand modeling. The above distributions in equations (19) and (20) are sometimes referred to as first-stage priors.

For brand $i$’s sales, we model $\Delta_{i}$, $\tau_{1,i},\dots,\tau_{83,i}$ and $V_{\boldsymbol{\beta},i}$ as random variables with the following prior distributions:

[TABLE]

[TABLE]

[TABLE]

where $\mathrm{vec}(\bar{\Delta})$ is the elements of the matrix $\bar{\Delta}$ (also 12-by-34) stacked into a column vector; $\otimes$ denotes the Kronecker product; $IW(\nu,V)$ is the inverse Wishart distribution with degrees of freedom $\nu$ and scale matrix $V$; and $\chi^{2}_{\nu_{e}}$ is a chi-squared distributed random variable with $\nu_{e}$ degrees of freedom. The matrices $\bar{\Delta}$, $V$ and $A$ and the scalars $\nu$, $\nu_{\epsilon}$ and $C$ are the prior hyperparameters. The distributions (21), (22) and (23) are sometimes referred to as second-stage priors.

The second type of HB model that we will consider, which we will refer to as HB-SemiLog, we assume the same specification above except that instead of equation (18), where we assume the logarithm of sales is linear in the logarithm of prices, we assume it is linear in the actual prices themselves:

[TABLE]

This particular HB model is essentially the same model proposed in Montgomery [1997]**. One minor difference is that the above model allows for cross-brand promotion effects (e.g., brand 1 being on a deal or being featured can affect sales of brand 2). In our comparison of predictive accuracies in the next section, we will test both the above regression model, which allows for cross-brand promotion effects, and a simpler regression model with only own-brand promotion effects (i.e., equation (18) and (24) excludes $d_{t,s,i^{\prime}}$ and $f_{t,s,i^{\prime}}$ for brands $i^{\prime}$ different from $i$, as in Montgomery 1997).

In addition to these two HB models, we also consider a random forest model, which we denote by RF. The random forest model that we estimate will be different from the previous two HB models. Rather than predicting the sales $q_{t,s,i}$ or log sales $\log(q_{t,s,i})$, we will instead predict profit, that is,

[TABLE]

where $c_{i}$ is the unit cost of brand $i$. For each brand $i$, we estimate the random forest model using $p_{1},\dots,p_{11},d_{1},\dots,d_{11},f_{1},\dots,f_{11},z_{1},\dots z_{11}$ as the predictor variables, where $z_{1},\dots,z_{11}$ are the store-level covariates of the store at which the prediction will be made. Letting $F_{i}$ be the random forest prediction of profit for brand $i$, the overall profit that is predicted from the price and promotion decision $\mathbf{p},\mathbf{d},\mathbf{f}$ at a store with covariate vector $\mathbf{z}$ is simply the sum of these individual brand-level predictions:

[TABLE]

The choice to predict profit as opposed to sales or log of sales is motivated by tractability. By predicting the profit from each brand, the total profit function $F_{total}$ can be directly used within the tree ensemble formulation. In contrast, if $F_{i}(\cdot)$ was a prediction of sales, then our objective function would be

[TABLE]

which would require modifications to our MIO formulation – one would need to model the tree ensemble behavior of each $F_{i}$, and then model the product of the price $p_{i}$ and the predicted sales. Although our formulation can be suitably modified to do this, the resulting model is harder to solve than the basic problem (2).

11.3 Predictive accuracy results

Our first experiment does not consider optimization, but merely compares the out-of-sample predictive accuracy of the three models – HB-LogLog, HB-SemiLog and RF. To do so, we proceed as follows. We take the whole data set of store-week observations, and we split it randomly into a training set and a test set. We use the training set to estimate the HB-LogLog, HB-SemiLog and RF models for each brand. For each model, we predict the profit from brand $i$ with the store-week observations in the test set, and compute the $R^{2}$ of the test set profit predictions. We repeat this procedure with ten random splits of the observations into train and test sets. In each split, 80% of the observations are randomly chosen for the training set, and the remaining 20% are used for the test set.

For HB-LogLog, we estimate it using Markov chain Monte Carlo (MCMC) using the R package bayesm. We use default values for the prior hyperparameters provided by bayesm. We run MCMC for 2000 iterations to obtain samples from the posterior distribution of $\boldsymbol{\beta}_{s}$ for each store $s$, and samples from the posterior distribution of $\tau_{s}$ for each store $s$. Due to mild autocorrelation in the draws of $\boldsymbol{\beta}_{s}$ and $\tau_{s}$, we thin the draws by retaining every fifth draw. Of these thinned draws, we retain the last $J=100$ samples. We index the draws/posterior samples by $j=1,\dots,J$. For an arbitrary price and promotion decision $\mathbf{p},\mathbf{d},\mathbf{f}$, we compute the predicted sales of brand $i$ at store $s$ by computing an approximation to the posterior expectation of sales $q_{s,i}$ under HB-LogLog as

[TABLE]

*where $\log(\mathbf{p})$ is the component-wise logarithm of $\mathbf{p}$, and quantities with the superscript $(j)$ correspond to the *$j$th posterior sample from the appropriate posterior distribution (that of either $\boldsymbol{\beta}_{s,i}$ or $\tau_{s,i}$). With this prediction of sales, we predict the profit of brand $i$ as $(p_{i}-c_{i})\cdot\hat{q}_{s,i}$.

For HB-SemiLog, we proceed in the same way as for HB-LogLog, except that the predicted sales are computed as

[TABLE]

For RF, we use the R package ranger to estimate one model for each brand’s profit. We use default parameters, with the exception of the number of trees which we vary in $\{20,50,100,500\}$. For the RF model, we emphasize again that, unlike HB-LogLog and HB-SemiLog which first predict log sales and then translates this to revenue, the RF model directly predicts profits of each brand.

In addition to the above three models, we also consider slightly modified versions of the above three models where we do not allow for cross promotion effects (i.e., for each brand $i$’s sales or profit prediction, we leave out $d_{i^{\prime}}$ and $f_{i^{\prime}}$ for $i^{\prime}\neq i$ as independent variables). In our presentation below, we distinguish these models from the ones above by using the suffix “-Own” – thus, HB-LogLog-Own, HB-SemiLog-Own and RF-Own are the log-log, semi-log and random forest models with only own-brand promotion effects.

The prediction tasks we will consider involve predicting profit, which requires us to specify the unit cost $c_{i}$ of each brand. Since the data set does not include this information, we will test two different sets of values for $\mathbf{c}=(c_{1},\dots,c_{11})$. In the first set, we set $c_{i}=0$ for each brand $i$; we are thus effectively predicting the revenue from each brand. In the second set, we set $c_{i}=0.9\times p_{i,\min}$, where $p_{i,\min}$ is the lowest price at which brand $i$ was offered in the whole data set.

Table 8 displays the test set/out-of-sample $R^{2}$ values for the profit predictions of each brand in the first prediction task ($c_{i}=0$* for all brands $i$), averaged over the ten random splits of the data. From this table, we can see that both the HB-LogLog and HB-SemiLog models are very inaccurate for some brands, achieving $R^{2}$ values that are close to zero. In one case, namely brand 10, the out-of-sample $R^{2}$ is even negative, indicating that the model is worse than a naive model that just predicts the average training set profit. For the log-log model, when we remove cross-brand promotion effects and move from HB-LogLog to HB-LogLog-Own, the out-of-sample $R^{2}$ exhibits an absolute improvement ranging from 0.04 to 0.68, with an average over all brands of 0.15. For the semi-log model, when we remove cross-brand promotion effects and move from HB-SemiLog to HB-SemiLog-Own, the $R^{2}$ exhibits an absolute improvement ranging from 0.03 to 0.85, with an average of 0.17. With regard to the semi-log model, this finding is consistent with that of Montgomery [1997]**, where an own-brand promotion effect specification yielded a lower Schwartz information criterion value and higher out-of-sample accuracy than the cross-brand specification.*

Comparing the HB models to the RF model with 500 trees, we can see that RF provides a significant improvement in predictive accuracy. For example, for brand 8, the highest $R^{2}$ attained by HB-LogLog, HB-LogLog-Own, HB-SemiLog and HB-SemiLog-Own, is 0.45. In contrast, the $R^{2}$ attained by RF with 500 trees is 0.67, which is an absolute improvement of 0.25. Over all of the brands, the improvement of RF over the best HB model for each brand ranges from 0.05 (brand 2) to as much as 0.48 (brand 4).

Within the family of RF models, Table 8 gives us a sense of how the number of trees affects the predictive accuracy. In particular, while the out-of-sample accuracy decreases as the number of trees is decreased, we can see that the loss in accuracy is very modest. For example, for brand 1, the $R^{2}$ is 0.7778 with 500 trees, which is reduced to 0.7663 when we use 50 trees; note that this is still higher than any of the HB models. This suggests that RF can still achieve an improvement over the HB models even with a smaller number of trees.

We can also determine the impact of cross promotion effects within RF. Note that unlike HB-LogLog and HB-SemiLog, where the out-of-sample $R^{2}$ improves once cross-promotion effects ($d_{i^{\prime}}$* and $f_{i^{\prime}}$ for $i^{\prime}$ different to the focal brand $i$) are removed, the opposite happens with the random forest models: RF-Own has slightly lower $R^{2}$ values than RF.*

Table 9 presents the out-of-sample $R^{2}$ values for the profit predictions of each brand in the second prediction task ($c_{i}=0.9\times p_{i,\min}$* for all brands $i$), averaged over the ten random splits of the data. The same insights about the relative performance of the three different families of models derived from Table 8 hold for this case. Overall, these results provide evidence that a random forest model for profit predictions can outperform state-of-the-art models for this type of data. While a detailed comparison of random forests and hierarchical Bayesian models is beyond the scope of the present paper, we believe these result are encouraging and underscore the potential of tree ensemble models, such as random forests, to be used for profit/revenue prediction and for making marketing decisions.*

11.4 Optimization results

We now turn our attention to optimization. Using the complete data set, we estimate the HB-LogLog-Own model, HB-SemiLog-Own model and RF model with 50 trees per brand. We fix the price of each brand $c_{i}$ to $0.9\times p_{i,\min}$. We restrict the price vector $\mathbf{p}$ to a set $\mathcal{P}$. We will specifically consider the following choice of $\mathcal{P}$:

[TABLE]

where

[TABLE]

In the above definition of $\mathcal{P}_{i}$, $\delta$ is a discretization parameter – for example, $\delta=0.05$ indicates that prices go up in increments of $0.05 – and *$p_{i,q25}$* and *$p_{i,q75}$* are respectively the 25th percentile and 75th percentile of the prices observed for brand *$i$* in the whole data set. In words, the above expression restricts brand *$i$*’s price to go up in increments of *$\delta$*, starting at the smallest multiple of *$\delta$* above *$p_{i,q25}$* and ending at the largest multiple of *$\delta$* below *$p_{i,q75}$*; for example, if *$p_{i,q25}=1.44$* and *$p_{i,q75}=1.78$*, then brand *$i$*’s prices would be restricted to *${1.45,1.50,1.55,1.60,1.65,1.70,1.75}$*. We will consider values of the discretization parameter *$\delta\in{0.05,0.10,0.20}$*. The reason we restrict each product’s price to lie between the 25th and 75th percentile is to ensure the prices are sufficiently “inside” the range of previously offered prices, and to prevent them from taking extreme values. For the purpose of this experiment, we do not consider optimization of *$\mathbf{d}$* and *$\mathbf{f}$*, so we fix each *$d_{i}$* and *$f_{i}$ to zero.*

For HB-LogLog-Own, our profit optimization problem for store $s$ can thus be formulated as

[TABLE]

For HB-SemiLog-Own, our profit optimization problem for store $s$ is

[TABLE]

We solve both (28) and (29) using local search from ten different randomly chosen starting points. We remark here that local search is not guaranteed to obtain provably optimal solutions. However, for the coarsest price increment $\delta=0.20$, it is possible to solve both (28) and (29) by complete enumeration; in doing so, we find that the solutions returned by local search are provably optimal for both of these problems and for all 83 stores. Thus, it is reasonable to expect that at finer discretization levels ($\delta=0.05$* or 0.10), the solutions should remain close to optimal.*

For RF, our profit optimization problem for store $s$ is

[TABLE]

where $F_{i}$ is the random forest prediction function for the profit from brand $i$, $\mathbf{0}$ is a vector of zeros of the appropriate dimension (for the two arguments above, both are of length 11), and $\mathbf{z}_{s}$ is the vector of store-level covariates of store $s$. Regarding (30), we observe that:

The random forest defining each $F_{i}$ may contain not only splits on $\mathbf{p}$, but also splits on $\mathbf{d}$, $\mathbf{f}$ and $\mathbf{z}$. However, $\mathbf{d}$, $\mathbf{f}$ and $\mathbf{z}$ are fixed to $\mathbf{0}$, $\mathbf{0}$ and $\mathbf{z}_{s}$, and are not decision variables. 2. 2.

While the actual split points on $p_{1},\dots,p_{11}$ could take any value, each $p_{i}$ in the optimization problem (30) is restricted to values in $\mathcal{P}_{i}$.

These two observations are valuable because we can use them to simplify the tree model. In particular, the two observations can be used to identify splits/leaves that are unreachable in each tree, and to thus remove redundant splits. For example, if $z_{1}=10.5$ in $\mathbf{z}_{s}$ and a split in a tree has the query “Is $z_{1}\leq 10.1$?”, then all of the splits and leaves to the left of that split (the “yes” branch) can be removed because $z_{1}$ does not satisfy the query. As another example, suppose $\delta=0.05$ in the definition of $\mathcal{P}$, and we take the “yes” branch for a split with query “Is $p_{1}\leq 1.69$?”; if we then encounter the query “Is $p_{1}\leq 1.66$?”, we know that we cannot reach the node on the “no” branch because this would imply $1.66<p_{1}\leq 1.69$, but $p_{1}$ is restricted to multiples of $\delta=0.05$.

Once we have identified all such splits/leaves that cannot be reached, we can remove them and “collapse” the tree to obtain a much simpler, store-specific tree that is only in terms of $\mathbf{p}$. Using these collapsed trees, we formulate problem (30) using formulation (2) and solve it using the split generation approach. Due to the large number of MIO problems that need to be solved (83 stores by 3 discretization levels), we warm start each MIO by solving (30) using local search from ten random starting points, and supplying the best of those solutions as the initial solution to Gurobi. In addition, we also used slightly modified parameters for Gurobi, which we report in Table 10.

Table 11 displays the average and maximum time to solve problems (28), (29) and (30), for different values of $\delta$, where the average and maximum are taken over the 83 stores in the chain. (For RF, the time includes both the local search warm start time as well as the split generation MIO time.) From this table we can see that the two HB optimization problems are solved quite quickly, requiring no more than 5 seconds across all levels of discretization and all 83 stores. The RF problem (30) was solved in about 11.4 seconds on average per store for the finest discretization level $\delta=0.05$ and in all instances within 30 seconds. Although the RF problem is not solved quite as fast as the two HB-based problems, the increase is modest, and reasonable given the fact that the MIO approach provides a provably optimal solution, whereas the local search solution method needed for the two HB problems does not.

Aside from computation times, it is also interesting to compare the quality of the price prescriptions themselves. The challenge with comparing the prices produced by each model is that the ground truth model that maps prices to expected profit is unknown to us: we cannot simply plug the prices into a “true” model and see which method performs best. To deal with this, we take the following approach. Let $\phi(\mathbf{p};m,s)$ be the predicted profit of the price vector $\mathbf{p}$ at store $s$ according to model $m$, where $m\in\{\text{RF},\text{HB-LogLog-Own},\text{HB-SemiLog-Own}\}$. Let $\mathbf{p}_{m,s}$ be the optimal price vector if we assume model $m$ for store $s$. Lastly, let $\vartheta(m_{1},m_{2})$ denote the relative performance of the optimal price vector derived from model $m_{1}$ when evaluated according to $m_{2}$:

[TABLE]

The metric $\vartheta(m_{1},m_{2})$ measures how well our prices do if we assume reality behaves according to $m_{1}$, but reality actually behaves according to $m_{2}$; this performance is measured relative to the best possible profit under $m_{2}$ (obtained by using the price vector $\mathbf{p}_{m_{2},s}$), and is averaged over all of the stores. For example, $\vartheta(\text{RF},\text{HB-LogLog-Own})$ would mean we plug the random forest price vector into the HB-LogLog-Own model, and measure how much of the maximum profit under HB-LogLog-Own it recovers. Note that by definition, we always have $\vartheta(m,m)=100\%$. In essence, the idea of this approach is to be agnostic to the model, and to measure how much is lost if we are wrong about what the underlying model is. For example, suppose that there were only two models, $m_{1}$ and $m_{2}$. If it turned out that $\vartheta(m_{1},m_{2})=90\%$ but $\vartheta(m_{2},m_{1})=60\%$, then assuming $m_{1}$ is “safer” than assuming $m_{2}$: assuming $m_{1}$ and being wrong (i.e., reality follows $m_{2}$) means only missing out on 10% of the highest possible profit, whereas assuming $m_{2}$ and being wrong means missing out on 40% of the highest possible profit.

Table 12 shows the results of this multi-model comparison. From this table, we can see that the HB-SemiLog-Own price vectors perform very well under the HB-LogLog-Own model, and vice versa. Compared to RF, we can see that the RF price vectors are about as robust as the HB price vectors. For example, with $\delta=0.20$, assuming the ground truth is RF but reality behaves according to HB-SemiLog-Own means capturing only 88.1% of the best possible profit under HB-SemiLog-Own (a loss of 11.9%). Conversely, assuming the ground truth is HB-SemiLog-Own but reality behaving according to RF means capturing only 87.6% of the best possible profit under RF (a loss of 12.4%). Given the inherent uncertainty in how reality will behave, we can see that optimizing against the RF model is no more risky than optimizing under either of the HB models.

Finally, it is interesting to compare the prices themselves. Figure 8 shows the optimal prices under the three different models, over the 83 different stores, for each of the eleven brands; the solutions correspond to the increment $\delta=0.05$. From this plot, we can see that in the overwhelming majority of cases, the optimal price vector $\mathbf{p}$ under HB-LogLog-Own and HB-SemiLog-Own involves setting each brand’s price to the lowest or highest allowable price. In contrast, the price vectors obtained from RF are more interesting: for many stores, there are many brands whose prices are not set to the lowest or highest allowable price. Instead, the RF solutions cover a wider range of values within the price intervals of many brands, and are not as extreme as the prices found for HB-LogLog-Own and HB-SemiLog-Own. The extreme nature of prices under log-log and semi-log aggregate demand models has been highlighted previously in the marketing community [see, e.g., Reibstein and Gatignon, 1984, Anderson and Vilcassim, 2001].

Overall, from a prescriptive standpoint, we see that the RF profit optimization problem can be solved quickly to provable optimality (Table 11), leading to prices that are not extreme (Table 8) and still achieve high profits under HB models (Table 12). We believe that these results, combined with the strong out-of-sample predictive accuracy shown in Section 11.3, underscore the potential benefit of random forests and tree ensemble models in customized pricing, as well as the value of our optimization methodology in transforming such tree ensemble models into pricing decisions.