Mean robust optimization

Robust optimization is a tractable and expressive technique for decision-making under uncertainty, but it can lead to overly conservative decisions when pessimistic assumptions are made on the uncertain parameters. Wasserstein distributionally robust optimization can reduce conservatism by being data-driven, but it often leads to very large problems with prohibitive solution times. We introduce mean robust optimization, a general framework that combines the best of both worlds by providing a trade-off between computational effort and conservatism. We propose uncertainty sets constructed based on clustered data rather than on observed data points directly thereby significantly reducing problem size. By varying the number of clusters, our method bridges between robust and Wasserstein distributionally robust optimization. We show finite-sample performance guarantees and explicitly control the potential additional pessimism introduced by any clustering procedure. In addition, we prove conditions for which, when the uncertainty enters linearly in the constraints, clustering does not affect the optimal solution. We illustrate the efficiency and performance preservation of our method on several numerical examples, obtaining multiple orders of magnitude speedups in solution time with little-to-no effect on the solution quality.

The simulations presented in this article were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University. We would like to thank Daniel Kuhn for the useful feedback and for pointing us to the literature on scenario reduction techniques.

Authors and Affiliations

1. Department of Operations Research and Financial Engineering, Princeton University, Princeton, USA

   Irina Wang & Bartolomeo Stellato

2. Operations Research Center, Massachusetts Institute of Technology, Cambridge, USA

   Cole Becker

3. Stochastics Group, Centrum Wiskunde en Informatica, Amsterdam, The Netherlands

   Bart Van Parys

Corresponding author

Correspondence to Bartolomeo Stellato.

Conflict of interest

The authors have no relevant financial or non-financial interests to disclose.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Irina Wang and Bartolomeo Stellato are supported by the NSF CAREER Award ECCS 2239771.

Proof of the constraint reformulation in (15)

We give a prove for the general case of maximum-of-concave functions. When \(J=1\), we take \(\alpha _{1k} = w_k\), for all k. For a simpler proof for the case of single-concave functions, refer to [47, Appendix A].

To simplify notation, we define \(c_k(v_{jk}) =\Vert v_{jk} - \bar{d}_k \Vert ^p - \epsilon ^p\). Then, starting from the inner optimization problem of (14):

We applied the Lagrangian in the first equality. Then, as the summation is over upper-semicontinuous functions \(g_j(v_{jk},x)\) concave in \(v_{jk}\), we applied the Von Neumann-Fan minimax theorem [38] to interchange the inf and the sup. Next, we rewrite the formulation using an epigraph trick, and make a change of variables.

In the last step, we rewrote \(v_{jk} = (\alpha _{jk}v_{jk})/\alpha _{jk}\), and maximized over \({\alpha _{jk}v_{jk} \in \alpha _{jk}S}\). In the case \(\alpha _{ij} > 0\), the terms in the summation are unchanged, and maximizing over \(v_{jk}\) is equivalent to maximizing over \(\alpha _{jk}v_{jk}\). In the case \(\alpha _{jk}= 0\), we have in the transformed formulation \((\alpha _{jk}v_{jk})/\alpha _{jk} = 0/0 = 0\), and \(\alpha _{jk}(g_j(0, x) - \lambda c_k(0)) = 0(g_j(0, x) - \lambda c_k(0)) = 0\). In the original formulation, the term \(\alpha _{jk}(g_j(v_{jk}, x) - \lambda c_k(v_{jk}))\) is also equivalent to 0. As the terms in the summation are 0 regardless of the value of \(v_{jk}\), maximizing over \(v_{jk}\) is equivalent to maximizing over 0. Therefore, we note that the optimal value remains unchanged. Next, we make substitutions \({h_{jk}} = \alpha _{jk}v_{jk}\), and define functions \(g_j'(h_{jk},x) = \alpha _{jk}g_j(h_{jk}/\alpha _{jk},x)\), \(c_k'(h_{jk}) = \alpha _{jk}{c_k}(h_{jk}/\alpha _{jk})\).

For the new functions defined, we applied the definition of conjugate functions. We also define the characteristic function \(\chi _{S}(v)\) with \(\chi _{S}(v) = 0\) if \(v \in S\); \(=\infty \) otherwise. Now, using the conjugate form \(f^*(y) = \alpha g^*(y)\) of a right-scalar-multiplied function \(f(x) = \alpha g(x/\alpha )\), and noting that \(\chi _{\alpha _{jk}S}(h)\) takes the same value as \(\alpha _{jk} \chi _{S}(h/{\alpha _{jk}}) \), we can rewrite the constraint as

Next, borrowing results from Esfahani et al. [24, Theorem 4.2], Rockafellar and Wets [41, Theorem 11.23(a), p. 493], and Zhen et al. [51, Lemma B.8], with regards to the conjugate functions of infimal convolutions and p-norm balls, we note that:

Substituting this in, we have

Taking the supremum over \(\alpha \), noting that \(\sum _{j=1}^J \alpha _{jk} = w_k\) for all k, we arrive at

which is equivalent to (15).

Reformulation of the maximum-of-concave case for \({p = \infty }\)

We again give the general proof for \(J \ge 1\). For the case \(J=1\), refer to the simpler proof [47, Appendix B]. When \({p = \infty }\), we have

(17)

which has a reformulation where the constraint above is dualized,

with new variables \(s_k \in {\textbf{R}}\), \(z_{jk} \in {\textbf{R}}^m\), and \(y_{jk} \in {\textbf{R}}^m\). We prove this by starting from the inner optimization problem of (17):

We have again formulated the Lagrangian and applied the minmax theorem to the sum of concave functions in \(v_{jk}\). Now, we can rewrite this in epigraph form,

then use the definition of the dual norm to rewrite the constraints,

In the last step, we again rewrote \(v_{jk} = (\alpha _{jk}v_{jk})/\alpha _{jk}\) inside functions \(g_j\). Using similar logic as in Appendix A, we note that this does not change the optimal value for both \(\alpha _{jk} > 0\) and \(\alpha _{jk} = 0\), as taking the supremum over the transformed variables is equivalent to taking the supremum over the original variables. For conciseness, we omit the steps of creating the right-scalar-multiplied function and transformations. For details, please refer to Appendix A. We apply the definition of conjugate functions and arrive at the optimization problem

Now, substituting \(\lambda _k = \lambda _k w_k\), \(z_{jk} = -z_{jk}\), and substituting in the conjugate functions derived in Appendix A, we have

Note that rescaling \(\lambda _k\) did not affect value of the problem, as minimizing \(\lambda _k\) is equivalent to minimizing \(\lambda _k w_k\), as \(w_k >0.\) Lastly, taking the supremum over \(\alpha \), we arrive at

Proof of the dual problem reformulation as \(\mathbf {p \rightarrow \infty }\)

We prove the dual equivalence. For a proof of primal equivalence, see the appendix of [47].

Theorem 6

Let S be a bounded set. Define here

Then, \( \lim _{p\rightarrow \infty } \bar{g}^K(x;p)= \bar{g}^K(x;\infty ) \) for any \(x\in X\).

Proof

First, from Equation (5) we have for any \(p> 1\) that

where the first equality is established in Appendix A and the second equality follows from Lemma 1. Remark that the inequality in the second step simply follows as we introduce \(\lambda _k\) and do not impose that \(\lambda _k=\lambda \) for all \(k=1,\dots , K\). Hence, considering the limit for p tending to infinity gives us now \( {\liminf _{p\rightarrow \infty }} \bar{g}^K(x;p)\ge \bar{g}^K(x;\infty ). \) It remains to prove the reverse \( {\limsup _{p\rightarrow \infty }} \bar{g}^K(x;p) \le \bar{g}^K(x;\infty ). \)

Second, we have for any \(p>1\) with \(1/p+1/q=1\) that

which follows from the choice \(\lambda _k = \left[ \frac{q-1}{q} \epsilon ^{\frac{1}{1-q}} \max _{k'=1}^K\Vert z_{k'} \Vert _*\right] \) and by imposing the restrictions \((q-1)^{1/4} \le \Vert z_k \Vert _* \le (q-1)^{-1/4}\) for all \(k =1,\dots , K\). Next, using the identities \(\phi (q)\left[ \frac{q-1}{q}\right] ^{1-q} = (q-1)^{(q-1)}/q^q \left[ \frac{q-1}{q}\right] ^{1-q} = 1/q\) and \(p+\frac{1}{1-q} = \frac{1}{\frac{1}{p}} + \frac{1}{1-q} = \frac{1}{1-\frac{1}{q}}+\frac{1}{1-q} = \frac{-q}{-q+1}+\frac{1}{1-q}=\frac{1-q}{1-q}=1\), as well as pulling \(\epsilon \left\| z_{k} \right\| _*>0\) out of the last two terms, we can rewrite the first constraints as

Then, we note that \(\max _{k'=1}^K \Vert z_{k'} \Vert _*/\left\| z_{k} \right\| _* \ge \Vert z_{k} \Vert _*/\left\| z_{k} \right\| _* = 1\) and \(\max _{k'=1}^K \Vert z_{k'} \Vert _*/{\left\| z_{k} \right\| _*} \le (q-1)^{-1/2}\), and hence we can apply Lemma 2 to obtain

Let

Hence, as \(q=p/(p-1)\) and \(q-1=1/(p-1)\) we have \( \bar{g}^K(x;p) \le \bar{g}_u^K(x;p) \) for all \(p>1\). Hence, taking the limit \(p\rightarrow \infty \) we have \(\bar{g}^K(x;\infty )\le \limsup _{p\rightarrow \infty } \bar{g}^K(x;p)\). In fact, as the function \(D\left( \frac{p}{p-1}\right) \) defined in Lemma 2 is nonincreasing for all p sufficiently large this implies that \(\bar{g}^K_u(x;p)\) is nonincreasing for p sufficiently large and hence we have \({\bar{g}^K(x;\infty )\le \limsup _{p\rightarrow \infty } \bar{g}^K(x;p) }\le \liminf _{p\rightarrow \infty } \bar{g}_u^K(x;p) = \lim _{p\rightarrow \infty } \bar{g}_u^K(x;p)\). We now prove here that \(\lim _{p\rightarrow \infty } \bar{g}_u^K(x;p) = \liminf _{p\rightarrow \infty } \bar{g}_u^K(x;p)\le \bar{g}^K(x;\infty )\). Consider any feasible sequence \(\{(z^t_{k}, \,y^t_{k}, \,s^t_k=[-g]^*(z^t_{k} - y^t_{k}) + \sigma _{S}(y^t_{k}) - (z^t_{k})^T{d}_k +\epsilon \left\| z^t_{k} \right\| _* )\}_{t\ge 1}\) in the optimization problem characterizing \(\bar{g}^K(x;\infty )\) in Equation (19) so that \(\lim _{t\rightarrow \infty }\sum _{k=1}^K w_k s^t_k=\bar{g}^K(x;\infty )\). Let \(\tilde{z}^t_{k}\in \arg \max \{ \Vert z \Vert _* \mid z\in {\textbf{R}}^m, \, \Vert z-z_k^t \Vert _*\le 1/t \}\) for all \(t\ge 1\) and \(k=1,\dots , K\) and observe that \(\Vert \tilde{z}^t_k \Vert _*= 1/t+\Vert z_k^t \Vert _*\ge 1/t\). Consider now an increasing sequence \(\{p_t\}_{t\ge 1}\) so that \((p_t-1)^{1/4}\ge \max ^K_{k=1} \Vert \tilde{z}^t_k \Vert _*\) and \((p_t-1)^{-1/4}\le 1/t\). Finally observe that the auxiliary sequence \(\{(\tilde{z}^t_{k}, \,\tilde{y}^t_{k}=y^t_k + (\tilde{z}^t_{k}-z^t_{k}), \,\tilde{s}^t_k= [-g]^*(\tilde{z}^t_{k} - \tilde{y}^t_{k}) + \sigma _{S}(\tilde{y}^t_{k}) - (\tilde{z}^t_{k})^T{d}_k +\epsilon \left\| \tilde{z}^t_{k} \right\| _* D\left( p_t/(p_t-1)\right) )\}_{t\ge 1}\) is by construction feasible in the minimization problem characterizing the function \(\bar{g}_u^K(x;p_t)\) in Equation (20). Hence, finally, we have

To establish the third inequality observe first that \(-(\tilde{z}^t_{k})^T{d}_k = -(z^t_{k})^T{d}_k - (\tilde{z}^t_k-z^t_k)^Td_k\le -(z^t_{k})^T{d}_k + \Vert \tilde{z}^t_k-z^t_k \Vert _* \Vert d_k \Vert \le -(z^t_{k})^T{d}_k + \Vert d_k \Vert /t\). Second, remark that we have

as \(\Vert \tilde{z}_t-z_t \Vert \le 1/t\). Lemma 2 guarantees that \(\lim _{t\rightarrow \infty } D\left( p_t/(p_t-1)\right) =1\). Finally, \(\left\| z^t_{k} \right\| _*\le \left\| \tilde{z}^t_{k} \right\| _*\le (p_t-1)^{1/4}\) and

with \(1/p+1/q=1\) using again Lemma 2. \(\square \)

Lemma 1

We have

for any \(p>1\) and \(q>1\) for which \(1/p+1/q=1\), \(\phi (q)=(q-1)^{q-1}/q^q\) and \(\epsilon >0\).

Proof

Remark that as the objective function \(\lambda \mapsto \phi (q)\lambda \left\| z/\lambda \right\| ^q_* + \lambda \epsilon ^p\) is continuous and we have \(\lim _{\lambda \rightarrow 0} \phi (q)\lambda \left\| z/\lambda \right\| ^q_* + \lambda \epsilon ^p = \lim _{\lambda \rightarrow \infty } \phi (q)\lambda \left\| z/\lambda \right\| ^q_* + \lambda \epsilon ^p=\infty \) as \(\epsilon >0\) there must exist a minimizer \(\lambda ^\star \in \min _{\lambda \ge 0}\, \phi (q)\lambda \left\| z/\lambda \right\| ^q_* + \lambda \epsilon ^p\) with \(\lambda _\star >0\). The necessary and sufficient first-order convex optimality conditions of the minimization problem guarantee

where we exploit that \(1/p+1/q=1\) and \(\phi (q)=(q-1)^{q-1}/q^q\). Indeed, we have

Hence, we have

where we exploit that \(1/p+1/q=1\) and \(\phi (q)=(q-1)^{q-1}/q^q\). Indeed, we have

establishing the claim. \(\square \)

Lemma 2

Let \(q>1\) then

with \(\lim _{q\rightarrow 1} D(q)=1\) and \(\lim _{q\rightarrow 1} (q-1)^{1/4} \left( D(q)-1\right) =0\).

Proof

The objective function is convex in t. Convex functions attain their maximum on the extreme points of their domain. The limits can be verified using standard manipulations. \(\square \)

Proof of Theorem 4

We prove (i) \(\bar{g}^N(x) \le \bar{g}^K(x)\), (ii) \(\bar{g}^K(x) \le \bar{g}^{N*}(x) + (L/2) D(K)\), and (iii) when the support constraint does not affect the worst-case value, \({\bar{g}^K(x) \le \bar{g}^{N}(x) + (L/2) D(K)}\).

Proof of (i). We begin with a feasible solution \(v_1, \dots , v_{N}\) of (MRO-K) with \(K=N\). Then for \(K < N\), we set \(u_k = \sum _{i \in C_k}{v_i}/|C_k|\) for each of the K clusters. We see \(u_k\) with \(k = 1,\dots ,K\) satisfies the constraints of (MRO-K) for \(K < N\), as

where we have applied triangle inequality, Jensen’s inequality for the convex function \(f(x) = \Vert x\Vert ^p\), and the constraint of (MRO-K) with \(K=N\). In addition, since the support \(S\) is convex, for every k our constructed \(u_k\), as the average of select points \(v_i\) \(\in S\), must also be within \(S\). The same applies with respect to the domain of g.

Since we have shown that the \(u_k\)’s satisfies the constraints for (MRO-K) with \(K<N\), it is a feasible solution. We now show that for this pair of feasible solutions, in terms of the objective value, \(\bar{g}^K(x) \ge \bar{g}^N(x)\). By assumption, g is concave in the uncertain parameter, so by Jensen’s inequality,

Since this holds true for \(u_k\)’s constructed from any feasible solution \(v_i, \dots , v_N\), we must have \(\bar{g}^K(x) \ge \bar{g}^N(x)\).

Proof of (ii). Next, we prove \(\bar{g}^K(x) \le \bar{g}^{N*}(x) + (L/2) D(K)\) by making use of the L-smooth condition on \(-g\). We first solve (MRO-K) with \(K < N\) to obtain a feasible solution \(u_1, \dots , u_k\). We then set \(\Delta _k =u_k - \bar{d}_k\) for each \(k \le K\), and set \(v_i = d_i + \Delta _k \quad \forall i \in C_k, k = 1,\dots ,K\). These satisfy the constraint of (MRO-K) with \(K=N\) and no support constraints (which we will refer to as \(K=N^*\) from here on), as

Since the constraints are satisfied, the constructed \(v_i \dots v_N\) are a valid solution for (MRO-K) with \(K=N^*\). We note that these \(v_i\)’s are also in the domain of g, given that the uncertain data \(\mathcal {D}_N\) is in the domain of g. For monotonically increasing functions g, (e.g., \(\log (u)\), \(1/(1+u)\)), we must have \(\Delta _k = u_k - \bar{d_k} \ge 0\) in the solution of (MRO-K), as the maximization of g over \(u_k\) will lead to \(u_k \ge \bar{d_k}\). Therefore, \(v_i = d_i + \Delta _k\) is also in the domain, as the L-smooth and concave function g with only a potential lower bound will not have holes in its domain above the lower bound. For monotonically decreasing functions g, the same logic applies with a nonpositive \(\Delta _k\). We now make use of the convex and L-smooth conditions [2, Theorem 5.8] on \(-g: \forall v_1, v_2 \in S, \lambda \in [0,1],\)

which, we can apply iteratively, with the first iteration being

where \(\bar{v}_2 = \frac{1}{|C_k|-1}\sum _{i \in C_k, i\ne 1} v_i\). Note that \(v_1 - \bar{v}_2 = d_1 - \frac{1}{|C_k|-1}\sum _{i \in C_k, i\ne 1} d_i\), as they share the same \({\Delta _k}\). The next iteration will be applied to \(g(\bar{v}_2)\), and so on. For each cluster k, this results in:

where we used the equivalence \( \sum _{i=2} ^{|C_k|}((i-1)/{i}) \left\| d_i - {\sum _{j=1}^{i-1}d_j}/({i-1}) \right\| ^2_2 = \sum _{i \in C_k} \Vert d_i - \bar{d}_k\Vert ^2_2. \) Now, summing over all clusters, we have

Since this holds for any feasible solution of (MRO-K) with \(K<N\), we must have \( \bar{g}^K(x) \le \bar{g}^{N*}(x) + (L/2) D(K)\).

Proof of Theorem 5

Proof of the lower bound. We use the dual formulations of the MRO constraints. For the case \(p \ge 1\), we first solve (MRO-K-Dual) with \(K<N\) to obtain dual variables \(z_{jk}\), \(\gamma _{jk}\). For each data label i in cluster \(C_k\), for all clusters \(k = 1,\dots ,K\), and for all pieces \(j = 1,\dots , J\), if we set

we have obtained a valid solution for (MRO-K-Dual) with \(K=N\). The increase in the objective value from \(\bar{g}^K(x)\) to that of \(\bar{g}^N(x)\), \({i.e.}~\bar{g}^N(x)- \bar{g}^K(x)\), is

We note that \(\delta (K,z,\gamma ) \ge (|C_k|/N)\sum _{k = 1}^K (-z_{1k} - H^T\gamma _{1k})^T (1/|C_k|)\sum _{i\in C_k} (d_i - \bar{d}_k) = (|C_k|/N) \sum _{k = 1}^K (-z_{1k} - H^T\gamma _{1k})^T0 = 0\). The constructed feasible solution for (MRO-K-Dual) with \(K=N\) is an upper bound for its optimal solution, since it is a minimization problem. We then have \(\bar{g}^N(x) - \bar{g}^K(x) \le \delta (K,z,\gamma )\), which translates to \(\bar{g}^N(x) - \delta (K,z,\gamma ) \le \bar{g}^K(x)\).

Now, for \(p=\infty \), the same procedure can be applied; we obtain a solution to MRO-K-Dual-\(\infty \) with \(K < N\), and construct a feasible solution for MRO-K-Dual-\(\infty \) with \(K=N\). We modify the variables in the same manner, with the exception of \(\lambda \), which is now set by \(\lambda _i = \lambda _k\). The definition of \(\delta (K,z,\gamma )\) remains unchanged, so the same result follows.

Proof of the upper bound. We use the primal formulations of the MRO constraints. Similar to the single-concave proof, for \(p \ge 1\), we first solve the MRO problem with K clusters to obtain a feasible solution \(u_{11}, \dots , u_{JK}\), \(\alpha _{11}, \dots , \alpha _{JK}\). Note that the existence of \(\alpha \) removes the need to analyze which function \(g_j\) attains the maximum for each cluster. We then set \(\Delta _{jk} =u_{jk} - \bar{d}_k\) for each \(k \le K\), and set \(v_{ji} = d_i + \Delta _{jk}\), \(\alpha _{ji} = \alpha _{jk}/|C_k| \quad \forall i \in C_k, k = 1,\dots ,K\). These satisfy the constraint of the problem with N clusters and without support constraints, as

For \(p = \infty \), we repeat the process of obtaining and modifying a feasible solution, and observe, for \(i=1,\dots ,N\),

Therefore, we also have constraint satisfaction for \(p = \infty \). The constructed solutions for both cases remain in \(\mathop \textbf{dom}_u{g}\), following the arguments in the proof of (ii) in Appendix D.

Next, the objective functions for \(p\ge 1\) and \(p =\infty \) are identical, so the same analysis applies. For each cluster k and function \(g_j\), using the L-smooth condition on \(-g_j\), we observe

Then, summing over all the clusters and functions, we have

Since this holds for all feasible solutions of the problem with K clusters, we conclude that

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions