A Discretization Framework for Robust Contextual Stochastic Optimization

Thank you for your comments. We have done our best to address them in the revised version of the paper. In terms of terminology, indeed our approach is not solving a form of the traditional robust optimization problem. However, the end result is still to make decisions which are robust to uncertainty.

We apologize for missing the proof of Lemma A.1 (now Lemma B.2). We have included it in the revised version of the paper. For your question, we assume the hypothesis stability holds true for each $\hat{p}_k(\cdot)$ separately. We have revised the proof of the theorem to be clearer and easier to follow as well. Please let us know if you have any questions.

We also agree we should compare against other methods of solving the minimum-violation objective we propose. We have added the approach of solving this using the method of Bertsimas and Kallus (2020) as you suggested.

For the ML methods such as Donti et al. (2017) and the ``policy optimizer method'' we noticed this approach unfortunately does not work well. We noticed the following pattern. Whenever the objective for a datapoint is below the threshold $\phi$, the gradient is zero and the weight update does not take these datapoints into account. Rather the update minimizes cost for any datapoints that previously resulted in cost above $\phi$. However, this can result in the datapoints which had cost below $\phi$ before the gradient update step, to now have a cost above $\phi$ after the update. Ultimately, this only creates a mask for which datapoints are updated at each step. In the experiments, we saw that $\phi$ had very little effect (to no effect) on changing the results at any of the quantiles or the mean. We are happy to add some discussion of this and experimental results in the appendix in the following days.

We still hope you find it valuable to compare against other objectives such as minimizing the mean and the robust optimization approach. From the comments of another reviewer, we also now included comparisons with CVaR. The goal is to validate the choice itself of using this minimum-violation objective.

We made this choice of objective because it very naturally aligns with the construction of the sets $H^\epsilon_k$. We have added a discussion on this in the revised version of the paper at the end of Section 3. It is perhaps less clear whether there is such a clean connection when considering an approach such as CVaR instead. Nevertheless, we also agree that this is an interesting question and a possible extension to the paper (thank you!).

Thank you very much for your comments and review. Taking your comments into account, we hope that the new section we added regarding the impact of $\epsilon$ upon each of the terms $c \cdot \mathcal{E}$ and $c \cdot OBJ({w})$ addresses your questions and concerns. Ultimately, both $c \cdot \mathcal{E}$ and $c \cdot OBJ({w})$ are always bounded by a constant for any $\epsilon$. The only term left is the generalization term $c \cdot \sqrt{\log(1 / \delta) / N}$. From here we give a procedure on how to choose $\epsilon$ based on the amount of data available. This gives us a guide to understand the magnitude needed for $\epsilon$, but ultimately something like cross-validation can be used to make a final choice for $\epsilon$ in practice. Please see the discussion following the theorem statement that contains more details on this.

We have also added experimental results for a broader range of quantiles, as well as comparing our approach to CVaR as you have suggested. We notice that the CVaR approach produces decisions that change discretely as the robustness parameter changes (the quantile being targeted), whereas the minimum-violation objective produces continuously changing decisions. This continuity might be explained, as you pointed out, by the additional probability term being multiplied. From the computational experiments we also see that our proposed approach outperforms CVaR (and all approaches we considered) by the following metric: for the same average cost, the cost at the $q^{th}$ quantile is generally lower for our approach. This holds in particular for higher quantiles. Please see the plots in Figure 2 that provides many comparisons.

This is an interesting question about considering the online learning setting. In our case, we decided to use regret because it essentially normalized the objective function --- so that when applying the optimal decision we obtain a regret of $0$. However, if one only uses the objective function (and not consider regret), the optimal decision can have a different cost for every realization of ${\nu}$. Of course, the method we propose applies just as easily when considering only the objective, and not regret.

Thank you for the thoughtful questions and suggestions. We agree especially that the learnability of the problem and choice of $\epsilon$ are very important. To address these issues, we have added a discussion on the effect of $\epsilon$ on each term of the bound. Indeed, you are correct, in the worst case one would need an exponential amount of data for small $\epsilon$. Nevertheless, instead, one can increase the value of $\epsilon$ so that the number of positive labels increases. Of course, in the worst case this may require the choice of $\epsilon$ to be a constant fraction of the diameter of the feasible region.

Primarily, the choice of $\epsilon$ will depend on the distribution of $w^*(\nu)$. There are two factors that can reduce the dimensionality of this. (1) The $w^*({\nu})$ often belong to a lower-dimensional subspace of $\mathbb{R}^d$ since they must satisfy constraints. For example, if $\mathcal{P}$ is the set of {$\{ w : Aw = b, w \geq 0 \}$} then the effective dimension of the problem is equal to the dimension of the null space of $A$. For instance, if $Aw = b$ defines a 1-dimensional line, we can view the problem as 1-dimensional regardless of the ambient dimension of the decisions $w \in \mathbb{R}^d$. (2) the distribution further depends on the data itself. The distribution of ${\nu}$, and as a result of $w^*({\nu})$, does not pathologically cover the entire space when presented with real-world data.

This dimensionality issue is not unique to our approach. The same problem is encountered for example, for the methods of Bertsimas and Kallus (2020) and more generally for methods such as KNN. That being said, we agree with you that it is important to address this issue in the paper.

Finally, that is an interesting extension to consider if the feasible region also depends on the contextual data and uncertainty. That is something we have been looking into and can add in the future.

We really appreciate your comments and review. We have done our best to incorporate your suggestions and updated the paper accordingly. We would also like to take the opportunity and answer some of your questions in this response as well.

We agree that it would be very helpful to include a more thorough discussion about how the model is trained, and how the parameters are chosen. In particular, please see the new discussion we added in the revised paper following the theorem about the role of $\epsilon$ on the bound and on learning the $\hat{p}^{k}({x})$ (the learned probabilities of $w^*({\nu_x}) \in H^\epsilon_k$).

In what follows we try to address some of the questions you raised.

Indeed, the union of $H_k^\epsilon$ need not cover the entire feasible region. This is part of the bound in Theorem 4 and captured in c which depends on $\mathbb{P}(R_{\nu_y}(w^*(\nu_x) \leq \epsilon)$. This is the probability that the objective of a solution for $\nu_x$ is within $\epsilon$ of the optimal solution with respect to some other ${\nu}_{{y}}$.

The values of $\hat{p}_k({x})$ learned need not sum to 1. This is because the sets $H^\epsilon_k$ are not disjoint and do not necessarily cover the entire feasible region. There is a lot of freedom in how the $\hat{p}_k({x})$ can be learned. In the first experiment on the inventory problem, each $\hat{p}_k({x})$ comes from a separate knn predictor. For the second experiment on the generator scheduling problem, $\hat{p}_k({x})$ are trained jointly and are the output of a neural network with $N$ outputs (one for each $k =1, \dots, N$).

As also pointed out by reviewer 9FbH, the method indeed has similarities to that of (Berstimas and Kallus (2020)). In particular, one can also view the objective function as some new $f_{{\nu}}({w}) = \max \{ R_{{\nu}}({w}) - \phi, 0 \}$ and directly apply their method on this objective.

To the best of our knowledge, we are the first to propose this new objective function. More importantly, the significant difference lies in how the weights $\hat{p}_{k}({x})$ are defined. In our case, these are explicitly informed by the optimization problem since $H_k^\epsilon$ are directly dependent on the data and objective function.

From a computational viewpoint, indeed the problem increases linearly with the amount of data $N$. However, there are three possible approaches to improve the runtime. One is to simply choose a random subset of the data to use to generate sets $H^\epsilon_k$. The second is, as we run the optimization problem for various out-of-sample points ${x}$, to keep track of the ${x}$ and the corresponding solutions $\hat w_{\epsilon, \gamma}(x)$ found. When given a new ${x}'$ we can for example find the nearest neighbor, say ${x}''$ and use $\hat{w}_{\epsilon, \gamma}({x}'')$ as a warm start when solving the optimization problem for ${x}'$. Finally, once calculating all values of $\hat{p}_k({x})$, many of these will have very low value close to zero. One approximation is to simply remove these from the optimization problem. Or alternatively, only keep the largest $k$ values of $\hat{p}_k({x})$ and only optimize with respect to those. We are happy to add more discussion on this in the appendix of the revised paper.