Predictive stochastic programming

Abstract

Several emerging applications call for a fusion of statistical learning and stochastic programming (SP). We introduce a new class of models which we refer to as Predictive Stochastic Programming (PSP). Unlike ordinary SP, PSP models work with datasets which represent random covariates, often refered to as predictors (or features) and responses (or labels) in the machine learning literature. As a result, these PSP models call for methodologies which borrow relevant concepts from both learning and optimization. We refer to such a methodology as Learning Enabled Optimization (LEO). This paper sets forth the foundation for such a framework by introducing several novel concepts such as statistical optimality, hypothesis tests for model-fidelity, generalization error of PSP, and finally, a non-parametric methodology for model selection. These new concepts, which are collectively referred to as LEO, provide a formal framework for modeling, solving, validating, and reporting solutions for PSP models. We illustrate the LEO framework by applying it to a production-marketing coordination model based on combining a pedagogical production planning model with an advertising dataset intended for sales prediction.

Appendices

Appendix I: Details of LEO-Wyndor and computational results

We consider the following statistical models in the computational results. The errors used for training (T) as well as validation (V) will depend on the model q, which in this case $q \in \{1, \ldots ,4 \}$ as shown below.

• 1. (DF) For Deterministic Forecasts (DF) we simply use the sales given by ${\hat{m}}_1(z) = {\hat{\beta _0}} + {\hat{\beta }}_1 z_1 + {\hat{\beta }}_2 z_2$. Thus, for deterministic predictions, we do not account for errors in forecast, i.e., training errors for DF are assumed to be zero. The remaining cases below will use errors based on the data.

• 2. (NDU) One approximation of the sales forecast is where the correlation between the coefficients are ignored. The resulting model takes the form $m_2(z, \xi ) := ({\hat{\beta }}_0 +\xi _0) + ({\hat{\beta }}_1 + \xi _1)z_1 + ({\hat{\beta }}_2 + \xi _2) z_2$, $z=(z_1, z_2)$ represents a generic outcome vector, and $(\xi _0, \xi _1, \xi _2)$ are generic error outcomes associated with ${{\hbox { normally distributed uncorrelated (NDU)}}}$ random variables $({\tilde{\xi }}_0, {\tilde{\xi }}_1, {\tilde{\xi }}_2)$ with mean zero, and a diagonal variance matrix $\varSigma _\beta$ as computed by the linear regression. The outcomes of errors associated with specific outcomes $(Z_i, W_i)$ are denoted $\xi _i = (\xi _{i0}, \xi _{i1}, \xi _{i2})$ and can be obtained by solving the following problem.

  $$\begin{aligned}&\underset{\xi _i}{\min } \left( \xi _i\right) ^\top \varSigma _\beta ^{-1}(\xi _i) \end{aligned}$$

  (19a)
  $$\begin{aligned}&\left( {\hat{\beta _0}}+\xi _{i0}\right) +\left( {\hat{\beta _1}} +\xi _{i1}\right) Z_{i1} + \left( {\hat{\beta _2}} +\xi _{i2}\right) Z_{i2} = W_i \end{aligned}$$

  (19b)

• 3. (NDC) This model of sales forecast ${\underline{\hbox {allows correlations between coefficients}}}$ of $m_3(z,\xi ) := ({\hat{\beta }}_0 + \xi _0) + ({\hat{\beta }}_1 + \xi _1)z_1 + ({\hat{\beta }}_2 + \xi _2) z_2$. Again, $(\xi _0, \xi _1, \xi _2)$ represent generic errors associated with generic outcomes z, whereas, identifying specific errors associated with each specific outcome $(Z_i, W_i)$ requires determining $(\xi _{i0}, \xi _{i1}, \xi _{i2})$ by solving (19), although the variance-covariance matrix $\varSigma _{\beta }$ is not diagonal in the case of NDC.

• 4. (EAE) This is the empirical additive error model, where $m_4(z,\xi _0)= ({\hat{\beta _0}} + \xi _{0}) + {\hat{\beta }}_1 z_{1} + {\hat{\beta }}_2 z_{2}$, and $\xi _0$ denotes a generic outcome of the error random variable. The specific values of these outcomes are $\xi _{i0} = W_i - m_4(Z_i, 0)$. We refer to this model as empirical additive errors.

As in the set up for (5), the index q refers to alternative error models (DF, NDU, NDC and EAE). Since all decision variables x (in this example) belong to the same space as outcomes of random variable Z, we remind the reader to explicitly use the recommended decision $x_q$ in defining the collection of validation errors by substituting $Z=z=x_q$ for the model indexed by q.

${{\hbox {LEO-Wyndor Model:}}}$ Index Sets and Variables

$i \equiv$ index of product, $i\in \{A, B\}.$ $y_i \equiv$ number of batches of product i produced.

$$\begin{aligned} f(x)&=-0.1x_1-0.5x_2+{\mathbb {E}}\left[ H\left( x, ~ {\tilde{\xi }}_q ~|~Z=z=x\right) \right] \end{aligned}$$

(20a)

$$\begin{aligned}&\text {s.t.} \quad x_1 + x_2 \le 200 \end{aligned}$$

(20b)

$$\begin{aligned}&x_1 - 0.5 x_2 \ge 0 \end{aligned}$$

(20c)

$$\begin{aligned}&L_1 := 0.7 \le x_1 \le U_1:=200, L_2 := 0 \le x_2 \le U_2 := 50 \end{aligned}$$

(20d)

$$\begin{aligned} h({x}, \xi _q ~|~ Z=z=x) \quad&= \quad \text {Max} \quad 3 y_A + 5 y_B \end{aligned}$$

(21a)

$$\begin{aligned} \text {s.t.} \quad y_A \quad \qquad&\le 8 \end{aligned}$$

(21b)

$$\begin{aligned} \quad \quad 2 y_B&\le 24 \end{aligned}$$

(21c)

$$\begin{aligned} \quad 3y_A + 2y_B&\le 36 \end{aligned}$$

(21d)

$$\begin{aligned} \quad y_A + y_B&\le m_q(z,\xi _q) \end{aligned}$$

(21e)

$$\begin{aligned} y_A, y_B&\ge 0. \end{aligned}$$

(21f)

As discussed in Sect. 5.1, $[L_1, U_1]$ and $[L_2, U_2]$ are available from the lower and upper limits of the regression data on covariates, and are assumed to allow decisions which are consistent with assumption A2. These values are same as upper and lower bounds for model (3a). Moreover, this instance is stated as a “maximization” model, whereas, our previous discussions were set in the context of “minimization”. When interpreting the results, it helps to keep this distinction in mind. The PSP models presented above are relatively general, allowing very general regression models such as kernel-based methods, projection pursuit, and others. However, our current computational infrastructure is limited to stochastic linear programming (SLP) and as a result the regressions are restricted to MLR models which are allowed to have one or more coefficients that are random, and are possibly correlated.

1.1 Computational results

We now present computations for the LEO-Wyndor data under alternative models. DF/LP represents predictive stochastic programming using deterministic forecasts, in which we use the expected value of the linear regression as the demand model. Because such a setup assume all coefficients to be deterministic, the resulting model is a linear program. Of course, one can also study other models where linear regression suggests alternative parameters: (a) the additive scalar error model (${\tilde{\xi }}_0$) using the empirical additive errors (EAE) and deterministic model coefficients ${\hat{\beta _0}}, {\hat{\beta _1}}, {\hat{\beta _2}}$ where ${\hat{\beta _1}}$ is the coefficient for TV expenditures, and ${\hat{\beta _2}}$ is the coefficient for radio expenditures; (b) a linear regression whose coefficients are random variables $\{{\tilde{\beta }}_j\}$, which are normally distributed and uncorrelated (NDU); (c) a linear regression whose coefficients are multi-variate random variables $\{{\tilde{\beta }}_\tau \}$ which and are possibly correlated (NDC). We reiterate that all three models EAE, NDU, NDC correspond to specific types of errors (which are indexed by q in the presentation in Sect. 3). Note that for models NDU and NDC, we have continuous random variables, and as a result we adopted SD as the solution methodology because it manages discrete and continuous random variables with equal ease. We refer to these results by NDU/SD and NDC/SD. However, note that for the case of EAE, the dataset is finite and reasonably manageable. Hence we will use both SAA and SD for this model, and refer to them by EAE/SAA and EAE/SD.

1.1.1 Results for alternative models (Error Terms)

The calculations begin with the first test as the top diamond block in Fig. 1b. Table 4 shows p-values and test results of $\chi ^2$ test for NDU/SD, NDC/SD and EAE. From values reported in Table 4, the fit appears to improve when a few of the data points near the boundary are eliminated. Since we retain more than 95% of the available training data, we accept the revised dataset as being appropriate for training.

Table 4 LEO-Wyndor: comparison of Chi-square test

1.1.2 Results for decisions and optimal value estimates

The decisions and various metrics discussed earlier are shown in Table 5. Although the prediction interval is listed in a symmetric way, the actual data points are asymmetric with respect to the estimated mean. The last two rows report the probability $\gamma$ and the corresponding tolerance level $\delta$, which are provided by SD algorithm based on the theorems in Sect. 4. We choose 1% of the mean value of validated objective to be $\varDelta$ in Theorem 2. Once again, notice that for both DF/LP and EAE/SAA, we do not report any probability because we use a deterministic solver as in (5).

The hypothesis test results for the recourse function objectives (the lowest diamond in Fig. 1b) for each model are reported in Table 6. The T-test rejected the DF/LP model as being appropriate for the given data, while the others were not rejected. The next two rows give the test results of variance based on F-statistic, and we conclude that none of the models can be rejected. We also performed a $\chi ^2$ test for the recourse function objectives using the training and validation sets. Again, the DF/LP model was rejected where as the others were not.

Remark 3

The concept of cross-validation ($\kappa$-fold) is uncommon in the stochastic programming literature.With $\kappa >1$, this tool is a computational embodiment of (15), and provides a prediction of error (Hastie et al. 2011). Without such $\kappa$-fold cross-validation ($\kappa >1$), it is often likely that model assessment can go awry. Although we have not explicitly reported on the case of $\kappa =1$ we have observed that the EAE/SAA model can get rejected, even though using $\kappa =5$ the EAE/SAA model is no longer rejected. This can be attributed to the fact that variance reduction due 5-fold cross-validation produced reduced Type I error (when compared with $\kappa =1$).

Table 7 reports the estimated optimization error, as well as the generalization error for all models. DF/LP shows the largest estimated optimization error, which indicates that it is not an appropriate model to recommend for this application. On the other hand, NDU/SD and NDC/SD have comparable and relatively small generalization errors. However the estimated optimization errors appear to be significant, therefore NDU/SD and NDC/SD do not yield the most profitable decisions. Had the criterion been simply the generalization error (as in SL), we might have chosen non-profit-maximizing decisions.

In Table 8 we present the pairwise comparisons due to the Kruskal-Wallis test which is a non-parametric method for testing similarities based on comparing counts in common bins for any pair of distributions. For the tests of DF/LP with other methodologies, the p-values are all smaller than 0.01, which implies that there are significant differences between DF/LP and each of the other four approaches. The stepped curve in Fig. 2 illustrates the ordering discovered by the Kruskal-Wallis test. Moreover, the curves for NDU/SD and NDC/SD are relatively close, whereas EAE/SAA and EAE/SD are indistinguishable. These similarities were quantified in Table 8 by the fact that the p-values for these comparisons are greater than 0.01. Finally, EAE/SAA and EAE/SD give the largest objective value, which is also reported in Table 5. The Kruskal-Wallis test suggests that the difference of EAE/SAA and EAE/SD is not significant, therefore both EAE/SAA and EAE/SD provide the most profitable decision.

Table 5 LEO-Wyndor: comparison of solutions for alternative models

Table 6 LEO-Wyndor: hypothesis test results for alternative models

Fig. 2

LEO-Wyndor: stochastic dominance of validated objectives under alternative models

Table 7 LEO-Wyndor: errors for alternative models

Table 8 LEO-Wyndor: Kruskal-Wallis test results ($p>0.01$)

Appendix II: algorithmic stochastic programming background

This Appendix is included for the sake of completeness, especially for those readers who are specialists in SL, but not SP.

1.1 Sample Average Approximation(SAA)

Sample Average Approximation is a standard sampling-based SP methodology, which involves replacing the expectation in the objective function by a sample average function of a finite number of data points. Suppose we have sample size of K, an SAA example is as follows:

$$\begin{aligned} \min _{x\in {\mathbf {X}}} F_K(x)=c^\top x+\frac{1}{K}\sum _{j=1}^K h(x,\xi ^j). \end{aligned}$$

(22)

As an overview, the SAA process may be summarized as follows.

Assumption SAA-a. The expectation function f(x) remains finite and well defined for all $x\in {\mathbf {X}}$. For $\delta > 0$ we denote by

$$\begin{aligned} S({\delta }) := \left\{ x\in {\mathbf {X}} : f(x)\le f^*+\delta \right\} ~~~\text {and}~~~ {\hat{S}}_K({\delta }) := \left\{ x\in {\mathbf {X}} : {\hat{f}}_K(x)\le {\hat{f}}_K^*+\delta \right\} , \end{aligned}$$

where $f^*$ denotes the true optimal objective, and $f_K^*$ denotes the optimal objective to the SAA problem with sample size K.

Assumption SAA-b. There exists a function $\kappa : \varXi \rightarrow {\mathbb {R}}_+$ such that its moment-generating function $M_\kappa (t)$ is finite valued for all t in a neighborhood of zero and

$$| F(x',\xi ) - F(x,\xi )| \le \kappa (\xi ) ||x'-x||$$

for a.e. $\xi \in \varXi$ and all $x', x \in {\mathbf {X}}$.

Assumption SAA-c. There exists a constant $\lambda > 0$ such that for any $x',x \in {\mathbf {X}}$ the moment-generating function $M_{x',x}(t)$ of random variable $[F(x',\xi )-f(x')]-[F(x,\xi )-f(x)]$, satisfies

$$\begin{aligned} M_{x',x}(t)\le \exp (\lambda ^2 ||x'-x||^2 t^2 /2), \forall t\in {\mathbb {R}}. \end{aligned}$$

(23)

From assumption SAA-b, $\Big |[F(x',\xi )-f(x')]-[F(x,\xi )-f(x)]\Big |\le 2L||x'-x||$ w.p. 1, and $\lambda =2L$.

Proposition 2

Suppose that assumptions SAA(a-c) hold, the feasible set ${\mathbf {X}}$ has a finite diameter D, and let $\delta >0, \varDelta \in [0,\delta ), \epsilon \in (0,1)$, and $L={\mathbb {E}}[\kappa (\xi )]$. Then for the sample size K satisfying

$$\begin{aligned} K(\epsilon , \delta )\ge \frac{8\lambda ^2D^2}{(\delta - \varDelta )^2 }\left[ n \ln \left( \frac{8 LD}{\delta -\varDelta }\right) +\ln \left( \frac{1}{1-\epsilon }\right) \right] , \end{aligned}$$

(24)

we have

$$\begin{aligned} \text {Pr}\left( {\hat{S}}_K(\varDelta ) \subset S(\delta ) \right) \ge \epsilon . \end{aligned}$$

Proof

This is Corollary 5.19 of Shapiro et al. (2009) with the assumption that the sample size K is larger than that required by large deviations theory (see 5.122 of Shapiro et al. (2009)). $\square$

1.1.1 Stochastic Decomposition (SD)

Unlike SAA which separates sampling from optimization, SD is based on sequential sampling and the method discovers sample sizes “on-the-fly” (Higle and Sen 1991, 1994). Because of sampling, any stochastic algorithm must contend with both variance reduction in objective values as well as solutions. SD uses M independent replications of value functions denoted $f^\nu$, $\nu = 1,\ldots ,M$. Each of these functions is a max-function whose affine pieces represent some Sample Average Subgradient Approximations (SASA). Because SD uses proximal point iterations, Higle and Sen (1994) shows that the maximum number of affine pieces is $n_1 + 3$, where $n_1$ is the number of first stage decisions, and these pieces are automatically identified using positive Lagrange multiplier estimates during the iterations. In theory, one can control the number of affine pieces to be smaller, but that can also be chosen “on-the-fly” depending on the size of the first stage decision variables $n_1$. When the number of affine functions reduces to only 1, the SD method reduces to a proximal stochastic variance reduction method (prox-SVRG) (Xiao and Zhang 2014). Among other strengths such as parallelizability and variance reduction, one of the main usability issues which SD overcomes is that when a model has a large number of random variables (as in the case of multi-dimensional random variables) it does not require users to choose a sample size because the sequential process automatically discovers an appropriate sample. Together with the notion of Statistical Optimality as set forth in Sect. 4, SD provides a very convenient optimization tool for PSP models.

For SLP models, Sen and Liu (2016) have already presented significant computational evidence of the advantage of SD over plain SAA. The reduced computational effort also facilitates replications for variance reduction (VR). VR in SD is achieved by creating the so-called compromise decision, denoted ${\mathbf {x}}^c$, which minimizes a grand-mean approximation $\bar{F}_M (x) := \frac{1}{M}\sum _{\nu =1}^M f^\nu (x)$ , where $\{f^\nu \}_{\nu =1}^M$ denotes a terminal value function approximation for each replication $\nu$. Suppose that solutions $x^\nu (\varepsilon )\in (\varepsilon -\text {argmin}~\{f^\nu (x) ~|~ x \in {\mathbf {X}} \})$ and ${\mathbf {x}}^c(\delta )\in (\delta -\text {argmin} \{{\bar{F}}_M(x) ~|~ x \in {\mathbf {X}}\})$. Then, Sen and Liu (2016) has proved consistency in the sense that $\lim _{\delta \rightarrow 0}\text {Pr}({\bar{F}}_M({\mathbf {x}}^c(\delta ))-f^*)\rightarrow 0$. Here are the critical assumptions of SD (Higle and Sen 1996b).

Assumption SD-a. The objective functions in the first and second stage models are either linear or quadratic, and all constraints are linear. Moreover, the set of first stage solutions is compact.

Assumption SD-b. The second stage optimization problem is feasible, and possesses a finite optimal value for all  $x \in {\mathbf {X}}$, and outcomes $\xi$ (i.e., the relatively complete recourse assumption holds).

Assumption SD-c. The second stage constraint functions are deterministic (i.e., fixed recourse), although the right hand side can be governed by random variables. The set of outcomes of the random variables is compact.

Assumption SD-d. The recourse function h is non-negative. So long as a lower bound on the optimal value is known, we can relax this assumption. (Higle and Sen 1996b)

A high-level structure of SD algorithm can be summarized as follows. For greater specifics regarding two-stage SD, we refer to section 3.1 of Sen and Liu (2016), and for the multi-stage version we refer to Sen and Zhou (2014).

The value function approximation for replication $\nu$ is denoted $f^\nu$ and the terminal solution for that replication is ${x}^\nu$. Note that we generate SASA functions using $K_\nu (\varepsilon )$ observations. Since these observations are i.i.d, the in-sample stopping rule ensures an unbiased estimate of the second stage objective is used for the objective function estimate at ${x}^\nu$. Hence, the Central Limit Theorem (CLT) implies that $[K_\nu (\varepsilon )]^{\frac{1}{2}} [f({x}^\nu ) - f^\nu ({x}^\nu )]$ is asymptotically normal ${\mathbf {N}}(0, \sigma _{\nu }^2)$, where $\sigma _{\nu }^2 < \infty$ denotes the variance of $f^\nu ({x}^\nu )$. Since

$$\begin{aligned} N = \min _\nu K_\nu (\varepsilon ), \end{aligned}$$

(25)

it follows that the error $[f({x}^\nu ) - f^\nu ({x}^\nu )]$ is no greater than $O_{p} (N^{-\frac{1}{2}})$. The following result provides the basis for compromise solutions $\mathbf{x }^c$ as proved in Sen and Liu (2016).

Proposition 3

Suppose that assumptions SD(a-d) stated in the Appendix II hold. Suppose $\bar{{\mathbf {x}}}$ is defined as in Theorem 1, and ${\mathbf {x}}^c = \bar{{\mathbf {x}}}$. Then,

(a) ${\mathbf {x}}^c$ solves

$$\begin{aligned} \underset{x \in {\mathbf {X}}}{\text {Min}} \quad \bar{F}_M (x) := \frac{1}{M} \sum _{\nu =1}^M f^\nu (x), \end{aligned}$$

(26)

(b)

$$\begin{aligned} f\left( {\mathbf {x}}^c\right) \le \bar{F}_M\left( {\mathbf {x}}^c\right) + O_p\left( (NM)^{-\frac{1}{2}}\right) , \end{aligned}$$

(27)

Proof

See Sen and Liu (2016). $\square$