A Stochastic Primal-Dual Method for Optimization with Conditional Value at Risk Constraints

TL;DR

This paper introduces a stochastic primal-dual method for optimizing problems with CVaR constraints, providing convergence guarantees and analyzing the impact of risk parameters on computational complexity.

Contribution

It proposes a simple modification to existing primal-dual algorithms that avoids complex bounds, enabling efficient risk-constrained optimization with theoretical guarantees.

Findings

Achieves $rac{ ext{constant}}{\sqrt{K}}$-approximate feasibility and optimality within $K$ iterations.

Provides bounds on step sizes and characterizes the computational cost related to risk parameters.

Draws parallels between sample complexity of CVaR optimization and chance-constrained programming.

Abstract

We study a first-order primal-dual subgradient method to optimize risk-constrained risk-penalized optimization problems, where risk is modeled via the popular conditional value at risk (CVaR) measure. The algorithm processes independent and identically distributed samples from the underlying uncertainty in an online fashion, and produces an $\eta/\sqrt{K}$-approximately feasible and $\eta/\sqrt{K}$-approximately optimal point within $K$ iterations with constant step-size, where $\eta$ increases with tunable risk-parameters of CVaR. We find optimized step sizes using our bounds and precisely characterize the computational cost of risk aversion as revealed by the growth in $\eta$. Our proposed algorithm makes a simple modification to a typical primal-dual stochastic subgradient algorithm. With this mild change, our analysis surprisingly obviates the need for a priori bounds or complex adaptive bounding schemes for dual variables assumed in many prior works. We also draw interesting parallels in sample complexity with that for chance-constrained programs derived in the literature with a very different solution architecture.

Equations

$$\begin{array}{rlrl}{\mathcal{P}}^{\mathrm{CVaR}}:& \underset{\mathbf{x}\in \mathbb{X}}{\mathrm{minimize}}& & F(\mathbf{x}):={\mathrm{CVaR}}_{\alpha }[f_{\omega }(\mathbf{x})],\\ & \text{subject to }& & G^i(\mathbf{x}):={\mathrm{CVaR}}_{{\beta }^i}[g_{\omega }^i(\mathbf{x})]\le 0,i=1,\dots ,m,\end{array}$$

$${\mathrm{CVaR}}_{\delta }[y_{\omega }]=\underset{u\in \mathbb{R}}{\min }\left\{u+\frac{1}{1-\delta }\text{mathds}E[y_{\omega }-u{]}^{+}\right\},$$

$${\mathrm{CVaR}}_{\delta }[0]=0,{\mathrm{CVaR}}_{\delta }[y_{\omega }^1+y_{\omega }^2]\le {\mathrm{CVaR}}_{\delta }[y_{\omega }^1]+{\mathrm{CVaR}}_{\delta }[y_{\omega }^2],$$

$${\mathrm{CVaR}}_{\delta }[t{y}_{\omega }]=t{\mathrm{CVaR}}_{\delta }[y_{\omega }],{\mathrm{CVaR}}_{\delta }[y_{\omega }+t^{\prime }]={\mathrm{CVaR}}_{\delta }[y_{\omega }]+t^{\prime }$$

$$\begin{array}{rlrl}{\mathcal{P}}^{\mathrm{E}}:& \underset{\mathbf{x}\in \mathbb{X}}{\mathrm{minimize}}& & F(\mathbf{x}):=\text{mathds}E[f_{\omega }(\mathbf{x})],\\ & \text{subject to }& & G^i(\mathbf{x}):=\text{mathds}E[g_{\omega }^i(\mathbf{x})]\le 0,i=1,\dots ,m.\end{array}$$

$${\mathbf{x}}_{k+1}$$

$${\mathbf{z}}_{k+1}$$

$$\mathcal{L}(\mathbf{x},\mathbf{z}):=F(\mathbf{x})+{\mathbf{z}}^{⊺}\mathbf{G}(\mathbf{x})=\text{mathds}E[{\mathcal{L}}_{\omega }(\mathbf{x},\mathbf{z})],$$

$$p_{\star }^{\mathrm{E}}:=\underset{\mathbf{x}\in \mathbb{X}}{\min }\underset{\mathbf{z}\in {\mathbb{R}}_{+}^m}{\max }\mathcal{L}(\mathbf{x},\mathbf{z}).$$

$$d_{\star }^{\mathrm{E}}:=\underset{\mathbf{z}\in {\mathbb{R}}_{+}^m}{\max }\underset{\mathbf{x}\in \mathbb{X}}{\min }\mathcal{L}(\mathbf{x},\mathbf{z}).$$

$$\mathcal{L}({\mathbf{x}}_{\star },\mathbf{z})\le \mathcal{L}({\mathbf{x}}_{\star },{\mathbf{z}}_{\star })\le \mathcal{L}(\mathbf{x},{\mathbf{z}}_{\star })$$

$$\begin{array}{c}0\in {\partial }_x\mathcal{L}({\mathbf{x}}_{\star },{\mathbf{z}}_{\star })+{\mathcal{N}}_{\mathbb{X}}({\mathbf{x}}_{\star }),G^i({\mathbf{x}}_{\star })\le 0,z_{\star }^i{G}^i({\mathbf{x}}_{\star })=0\end{array}$$

$$\nabla F({\mathbf{x}}_{\star })+\sum _{i=1}^m{z}_{\star }^i\nabla G^i({\mathbf{x}}_{\star })+\mathbf{n}=0.$$

$$\underset{\le F(\mathbf{x})-F({\mathbf{x}}_{\star })}{\underbrace{⟨\nabla F({\mathbf{x}}_{\star }),\mathbf{x}-{\mathbf{x}}_{\star }⟩}}+\sum _{i=1}^m\underset{\le z_{\star }^i[G^i(\mathbf{x})-G^i({\mathbf{x}}_{\star })]}{\underbrace{z_{\star }^i⟨\nabla G^i({\mathbf{x}}_{\star }),\mathbf{x}-{\mathbf{x}}_{\star }⟩}}+\underset{\le 0}{\underbrace{⟨\mathbf{n},\mathbf{x}-{\mathbf{x}}_{\star }⟩}}=0.$$

$$\mathcal{L}({\mathbf{x}}_{\star },{\mathbf{z}}_{\star })-\mathcal{L}({\mathbf{x}}_{\star },\mathbf{z})={\mathbf{z}}_{\star }^{⊺}\mathbf{G}({\mathbf{x}}_{\star })-{\mathbf{z}}^{⊺}\mathbf{G}({\mathbf{x}}_{\star })\ge 0,$$

$$\begin{array}{c}{P}_1:=2\Vert {\mathbf{x}}_1-{\mathbf{x}}_{\star }{\Vert }^2+4\Vert \text{mathds}1+{\mathbf{z}}_{\star }{\Vert }^2,\\ P_2:=8(4C_F^2+{\sigma }_F^2)+2\Vert {\mathbf{D}}_G{\Vert }^2,P_3:=8m(4\Vert {\mathbf{C}}_G{\Vert }^2+\Vert {\mathbf{\sigma }}_G{\Vert }^2).\end{array}$$

$$\text{mathds}E[F({\bar{\mathbf{x}}}_{K+1})]-p_{\star }^{\mathrm{E}}$$

$$\text{mathds}E[G^i({\bar{\mathbf{x}}}_{K+1})]$$

$$\text{mathds}E[F({\bar{\mathbf{x}}}_{K+1})]-p_{\star }^{\mathrm{E}}\le \frac{\eta }{\sqrt{K}},\text{mathds}E[G^i({\bar{\mathbf{x}}}_{K+1})]\le \frac{\eta }{\sqrt{K}}$$

$${\mathbf{x}}_{k+1}$$

$${\mathbf{z}}_{k+1}$$

$${\phi }_{\delta }(\mathbf{x},z^i):=\{\begin{array}{ll}{z}^i{g}^i(\mathbf{x})+\frac{\delta }{2}[g^i(\mathbf{x}){]}^2,& \text{if }\delta g^i(\mathbf{x})+z^i\ge 0,\\ -\frac{(z^i{)}^2}{2\delta },& \text{otherwise}\end{array}$$

$${\mathbf{z}}_{k+1}$$

$$\begin{array}{rl}& {\gamma }_k\text{mathds}E[\mathcal{L}({\mathbf{x}}_{k+1},\mathbf{z})-\mathcal{L}(\mathbf{x},{\mathbf{z}}_k)]+\frac{1}{2}\text{mathds}E\Vert {\mathbf{x}}_{k+1}-\mathbf{x}{\Vert }^2+\frac{1}{2}\text{mathds}E\Vert {\mathbf{z}}_{k+1}-\mathbf{z}{\Vert }^2\\ & \le \frac{1}{2}\text{mathds}E\Vert {\mathbf{x}}_k-\mathbf{x}{\Vert }^2+\frac{1}{2}\text{mathds}E\Vert {\mathbf{z}}_k-\mathbf{z}{\Vert }^2+\frac{1}{4}{P}_2{\gamma }_k^2+\frac{1}{4}{P}_3{\gamma }_k^2\text{mathds}E\Vert {\mathbf{z}}_k{\Vert }^2\end{array}$$

$$\text{mathds}E\Vert {\mathbf{z}}_k{\Vert }^2\le \frac{P_1+P_2{A}_K}{1-P_3{A}_K}$$

$$\begin{array}{rl}& \text{mathds}E[F({\mathbf{x}}_{k+1})-F(\mathbf{x})+{\mathbf{z}}_k^{⊺}\mathbf{G}({\mathbf{x}}_{k+1})-{\mathbf{z}}_k^{⊺}\mathbf{G}(\mathbf{x})|\text{mathpzc}{W}_k]+\frac{1}{2{\gamma }_k}\text{mathds}E[\Vert {\mathbf{x}}_{k+1}-\mathbf{x}{\Vert }^2|\text{mathpzc}{W}_k]\\ & \le \frac{1}{2{\gamma }_k}\Vert {\mathbf{x}}_k-\mathbf{x}{\Vert }^2+2{\gamma }_k(4C_F^2+{\sigma }_F^2)+2{\gamma }_{k}m(4\Vert {\mathbf{C}}_G{\Vert }^2+\Vert {\mathbf{\sigma }}_G{\Vert }^2)\Vert {\mathbf{z}}_k{\Vert }^2\end{array}$$

$$\begin{array}{rl}& \text{mathds}E[{\mathbf{z}}^{⊺}\mathbf{G}({\mathbf{x}}_{k+1})-{\mathbf{z}}_k^{⊺}\mathbf{G}({\mathbf{x}}_{k+1})|\text{mathpzc}{W}_k]+\frac{1}{2{\gamma }_k}\text{mathds}E[\Vert {\mathbf{z}}_{k+1}-\mathbf{z}{\Vert }^2|\text{mathpzc}{W}_k]\\ & \le \frac{1}{2{\gamma }_k}\Vert {\mathbf{z}}_k-\mathbf{z}{\Vert }^2+\frac{{\gamma }_k}{2}\Vert {\mathbf{D}}_G{\Vert }^2\end{array}$$

$$⟨{\mathbf{x}}_{k+1}-\mathbf{x},\nabla f_{\omega }({\mathbf{x}}_k)+\sum _{i=1}^m{z}_k^i\nabla {\mathbf{g}}_{\omega }^i({\mathbf{x}}_k)+\frac{1}{{\gamma }_k}({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)⟩\le 0.$$

$$\begin{array}{rl}⟨{\mathbf{x}}_{k+1}-\mathbf{x},\nabla f_{\omega }({\mathbf{x}}_k)⟩=& \underset{\ge F({\mathbf{x}}_{k+1})-F({\mathbf{x}}_k)}{\underbrace{⟨{\mathbf{x}}_{k+1}-{\mathbf{x}}_k,\nabla F({\mathbf{x}}_{k+1})⟩}}+⟨{\mathbf{x}}_k-\mathbf{x},\nabla f_{\omega }({\mathbf{x}}_k)⟩\\ & -⟨{\mathbf{x}}_{k+1}-{\mathbf{x}}_k,\nabla F({\mathbf{x}}_{k+1})-\nabla f_{\omega }({\mathbf{x}}_k)⟩,\end{array}$$

$$\text{mathds}E[⟨{\mathbf{x}}_k-\mathbf{x},\nabla f_{\omega }({\mathbf{x}}_k)⟩|\text{mathpzc}{W}_k]=⟨{\mathbf{x}}_k-\mathbf{x},\nabla F({\mathbf{x}}_k)⟩\ge F({\mathbf{x}}_k)-F(\mathbf{x}).$$

$$\begin{array}{rl}& \text{mathds}E[⟨{\mathbf{x}}_{k+1}-\mathbf{x},\nabla f_{\omega }({\mathbf{x}}_k)⟩|\text{mathpzc}{W}_k]\\ & \ge \text{mathds}E[F({\mathbf{x}}_{k+1})-F(\mathbf{x})-⟨{\mathbf{x}}_{k+1}-{\mathbf{x}}_k,\nabla F({\mathbf{x}}_{k+1})-\nabla f_{\omega }({\mathbf{x}}_k)⟩|\text{mathpzc}{W}_k].\end{array}$$

$$\begin{array}{rl}& \sum _{i=1}^m\text{mathds}E[⟨{\mathbf{x}}_{k+1}-\mathbf{x},z_k^i\nabla g_{\omega }^i({\mathbf{x}}_k)⟩|\text{mathpzc}{W}_k]\\ & \ge \sum _{i=1}^m\text{mathds}E[z_k^i{G}^i({\mathbf{x}}_{k+1})-z_k^i{G}^i(\mathbf{x})|\text{mathpzc}{W}_k]\\ & -\text{mathds}E[⟨{\mathbf{x}}_{k+1}-{\mathbf{x}}_k,z_k^i\nabla G^i({\mathbf{x}}_{k+1})-z_k^i\nabla g_{\omega }^i({\mathbf{x}}_k)⟩|\text{mathpzc}{W}_k]\\ & =\text{mathds}E[{\mathbf{z}}_k^{⊺}\mathbf{G}({\mathbf{x}}_{k+1})-{\mathbf{z}}_k^{⊺}\mathbf{G}(\mathbf{x})|\text{mathpzc}{W}_k]\\ & -\sum _{i=1}^m\text{mathds}E[⟨{\mathbf{x}}_{k+1}-{\mathbf{x}}_k,z_k^i\nabla G^i({\mathbf{x}}_{k+1})-{\mathbf{z}}_k^i\nabla g_{\omega }^i({\mathbf{x}}_k)⟩|\text{mathpzc}{W}_k].\end{array}$$

Full text

1 Introduction

We study iterative primal-dual stochastic subgradient algorithms to solve risk-sensitive optimization problems of the form

[TABLE]

where $\omega\in\Omega$ is random and $\alpha,\boldsymbol{\beta}:=(\beta^{1},\ldots,\beta^{m})$ in $[0,1)$ define risk-aversion parameters. The collection of real-valued functions $f_{\omega}$, $g^{1}_{\omega},\ldots,g^{m}_{\omega}$ are assumed convex but not necessarily differentiable, over the closed convex set $\mathbb{X}\subseteq\mathbb{R}^{n}$, where $\mathbb{R}$ and $\mathbb{R}_{+}$ stand for the set of real and nonnegative numbers, respectively. Denote by ${\bm{G}}$ and ${\bm{g}}_{\omega}$, the collection of $G^{i}$’s and $g^{i}_{\omega}$’s respectively, for $i=1,\ldots,m$. ${\rm CVaR}$ stands for conditional value at risk. For any $\delta\in[0,1)$, ${\rm CVaR}_{\delta}[y_{\omega}]$ of a scalar random variable $y_{\omega}$ with continuous distribution equals its expectation computed over the $1-\delta$ tail of the distribution of $y_{\omega}$. For $y_{\omega}$ with general distributions, ${\rm CVaR}$ is defined via the following variational characterization

[TABLE]

following [1]. For each ${\bm{x}}\in\mathbb{X}$, assume that $\mathds{E}[|f_{\omega}({\bm{x}})|]$ and $\mathds{E}[|g^{i}_{\omega}({\bm{x}})|]$ are finite, implying that $F$ and ${\bm{G}}$ are well-defined everywhere in $\mathbb{X}$.

${{{\cal P}}^{\rm CVaR}}$ offers a modeler the flexibility to indicate her risk preference in $\alpha,\boldsymbol{\beta}$. With $\alpha$ close to zero, she indicates risk-neutrality towards the uncertain cost associated with the decision. With $\alpha$ closer to one, she expresses her risk aversion towards the same and seeks a decision that limits the possibility of large random costs associated with the decision. Similarly, $\beta$’s express the risk tolerance in constraint violation. Choosing $\beta$’s close to zero indicates that constraints should be satisfied on average over $\Omega$ rather than on each sample. Driving $\beta$’s to unity amounts to requiring the constraints to be met almost surely. Said succinctly, ${{{\cal P}}^{\rm CVaR}}$ permits the modeler to customize risk preference between the risk-neutral choice of expected evaluations of functions to the conservative choice of robust evaluations.

There is a growing interest in solving risk-sensitive optimization problems with data. See [2, 3] for recent examples that tackle problems with generalized mean semi-deviation risk that equals $\mathds{E}[{y_{\omega}}]+c\mathds{E}[|{y_{\omega}}-\mathds{E}[{y_{\omega}}]|^{p}]^{1/p}$ for $p>1$ for a random variable $y_{\omega}$. There is a long literature on risk measures, e.g., see [1, 4, 5, 6, 7, 8]. We choose ${\rm CVaR}$ for three particular reasons. First, it is a coherent risk measure, meaning that it is normalized, sub-additive, positively homogeneous and translation invariant, i.e.,

[TABLE]

for random variables $y_{\omega},y^{1}_{\omega},y^{2}_{\omega}$, $t>0$ and $t^{\prime}\in\mathbb{R}$. An important consequence of coherence is that $F$ and ${\bm{G}}$ in ${{{\cal P}}^{\rm CVaR}}$ inherit the convexity of $f_{\omega}$ and ${\bm{g}}_{\omega}$. Convexity together with the variational characterization in (2) allow us to design sampling based primal-dual methods for ${{{\cal P}}^{\rm CVaR}}$ for which we are able to provide finite sample analysis of approximate optimality and feasibility. The popularity of the ${\rm CVaR}$ measure is our second reason to study ${{{\cal P}}^{\rm CVaR}}$. Following Rockafellar and Uryasev’s seminal work in [1], ${\rm CVaR}$ has found applications in various engineering domains, e.g., see [9, 10], and therefore we anticipate wide applications of our result. Our third and final reason to study ${{{\cal P}}^{\rm CVaR}}$ is its close relation to other optimization paradigms in the literature as we describe next.

${{{\cal P}}^{\rm CVaR}}$ without constraints and $\alpha=0$ reduces to the minimization of $\mathds{E}[f_{\omega}({\bm{x}})]$, the canonical stochastic optimization problem. With $\alpha\uparrow 1$, the problem description of ${{{\cal P}}^{\rm CVaR}}$ approaches that of a robust optimization problem (see [11]) of the form $\min_{{\bm{x}}\in\mathbb{X}}\text{ess}\sup_{\omega\in\Omega}f_{\omega}({\bm{x}})$, where $\text{ess}\sup$ denotes the essential supremum. Driving $\beta$’s to unity, ${{{\cal P}}^{\rm CVaR}}$ demands the constraints to be enforced almost surely. Such robust constraint enforcement is common in multi-stage stochastic optimization problems with recourse and discrete-time optimal control problems, e.g., in [12, 13, 14]. ${\rm CVaR}$-based constraints are closely related to chance-constraints introduced by Charnes and Cooper in [6], that enforce ${\sf Pr}\{g_{\omega}({\bm{x}})\leq 0\}>1-\varepsilon$ where ${\sf Pr}$ refers to the probability measure on $\Omega$. Even if $g_{\omega}$ is convex, chance-constraints typically describe a nonconvex feasible set. It is well-known that ${\rm CVaR}$-based constraints provide a convex inner approximation of chance-constraints. Restricting the probability of constraint violation does not limit the extent of any possible violation, while ${\rm CVaR}$-based enforcement does so in expectation. ${\rm CVaR}$ is also intimately related to the buffered probability of exceedence (bPOE) introduced and studied more recently in [7, 15]. In fact, bPOE is the inverse function of ${\rm CVaR}$ and hence, problems with bPOE-constraints can often be reformulated as instances of ${{{\cal P}}^{\rm CVaR}}$.

It can be challenging to compute ${\rm CVaR}$ of $f_{\omega}({\bm{x}})$ or ${\bm{g}}_{\omega}({\bm{x}})$ for a given decision variable ${\bm{x}}$ with respect to a general distribution on $\Omega$ for two reasons. First, if samples from $\Omega$ are obtained from a simulation tool, an explicit representation of the probability distribution on $\Omega$ may not be available. Second, even if such a distribution is available, computation of ${\rm CVaR}$ (or even the expectation) can be difficult. For example, with $f_{\omega}$ as the positive part of an affine function and $\omega$ being uniformly distributed over a unit hypercube, computation of $\mathds{E}[f_{\omega}]$ via a multivariate integral is #P-hard according to [16, Corollary 1]. Therefore, we do not assume knowledge of $F$ and ${\bm{G}}$ but rather study a sampling-based algorithm to solve ${{{\cal P}}^{\rm CVaR}}$.

Solution architectures for ${{{\cal P}}^{\rm CVaR}}$ via sampling come in two flavors. The first approach is sample average approximation (SAA) that replaces the expectation in (2) by an empirical average over $N$ samples. One can then solve the sampled problem as a deterministic convex program.111For the unconstrained problem, variance-reduced stochastic gradient descent methods can efficiently minimize the resulting finite sum as in [17, 18]. We take the second and alternate approach of stochastic approximation and process independent and identically distributed (i.i.d.) samples from $\Omega$ in an online fashion. Iterative stochastic approximation algorithms for the unconstrained problem have been studied since the early works by Robbins and Monro in [19] and by Kiefer and Wolfowitz in [20]. See [21] for a more recent survey. Zinkevich in [22] proposed a projected stochastic subgradient method to tackle constraints in such problems. Without directly knowing ${\bm{G}}$, we cannot easily project the iterates on the feasible set $\{{\bm{x}}\in\mathbb{X}\ |\ {\bm{G}}({\bm{x}})\leq 0\}$. We circumvent the challenge by associating Lagrange multipliers ${\bm{z}}\in\mathbb{R}^{m}_{+}$ to the constraints and iteratively updating ${\bm{x}},{\bm{z}}$ by using $f_{\omega},{\bm{g}}_{\omega}$ and their subgradients via a first-order stochastic primal-dual algorithm for ${{{\cal P}}^{\rm CVaR}}$ along the lines of [23, 24, 25].

In Section 2, we first design and analyze Algorithm 1 for ${{{\cal P}}^{\rm CVaR}}$ with $\alpha=0,\boldsymbol{\beta}=0$, i.e., the optimization problem

[TABLE]

First-order stochastic primal-dual algorithms have a long history, dating back almost forty years, including that in [26, 23, 27, 21, 28]. The analyses of these algorithms often require a bound on the possible growth of the dual variables. Borkar and Meyn in [29] stress the importance of compactness assumptions in their analysis of stochastic approximation algorithms. A priori bounds used in [28] are difficult to know in practice and techniques for iterative construction of such bounds as in [23] require extra computational effort. A regularization term in the dual update has been proposed in [30, 31] to circumvent this limitation. Instead, we propose a different modification to the classical primal-dual stochastic subgradient algorithm. We were surprised to find that, with this simple modification, the analysis allows us to completely bypass the need for explicit or implicit bounds on the dual variable in our analysis. While the classical primal-dual approach samples once for a single update of the primal and the dual variables, we sample twice–once to update the primal variable and then again to update the dual variable with the most recent primal iterate–thus, adopting a Gauss-Seidel approach in place of a Jacobi framework. For Algorithm 1, we bound the expected optimality gap and constraint violations at a suitably weighted average of the iterates by $\eta/\sqrt{K}$ for a constant $\eta$ with a constant step-size algorithm. Using these bounds, we then carefully optimize the step-size that allows us to reach within a given threshold of suboptimality and constraint violation with the minimum number of iterations. The additional sample required in our update aids in the analysis; however, it comes at the price of making the sample complexity double of the iteration complexity. We also provide stability analysis of our algorithm with decaying step-sizes that takes advantage of a dissipation inequality that we derive for our Gauss-Seidel approach.

In Section 3, we solve ${{{\cal P}}^{\rm CVaR}}$ with general risk aversion parameters $\alpha,\boldsymbol{\beta}$ using Algorithm 1 on an instance of ${{{\cal P}}^{\rm E}}$ obtained through a standard reformulation via the variational formula for ${\rm CVaR}$ in (2) from [1]. We then bound the expected suboptimality and constraint violation at a weighted average of the iterates for ${{{\cal P}}^{\rm CVaR}}$ by $\eta(\alpha,{\bm{\beta}})/\sqrt{K}$. Upon utilizing the optimized step-sizes from the analysis of ${{{\cal P}}^{\rm E}}$, we are then able to study the precise growth in the required iteration (and sample) complexity of ${{{\cal P}}^{\rm CVaR}}$ as a function of $\alpha,{\bm{\beta}}$. Not surprisingly, the more risk-averse a problem one aims to solve, the greater this complexity increases. A modeler chooses risk aversion parameters primarily driven by attitudes towards risk in specific applications. Our precise characterization of the growth in sample complexity with risk aversion will permit the modeler to balance between desired risk levels and computational challenges in handling that risk. We remark that the algorithmic architecture for the risk neutral problem may not directly apply to the risk-sensitive variant for general risk measures. For example, the algorithm described in [2] for general mean-semideviation-type risk measures is considerably more complex than that required for the risk-neutral problem. We are able to extend our algorithm and its analysis for ${{{\cal P}}^{\rm E}}$ to ${{{\cal P}}^{\rm CVaR}}$, thanks to the variational form in (2) that ${\rm CVaR}$ admits. See the discussion after the proof of Theorem 3.1 for a precise list of properties a risk measure must exhibit for us to apply the same trick. Using concentration inequalities, we also report an interesting connection of our results to that in [32, 33] on scenario approximations to chance-constrained programs. The resemblance in sample complexity is surprising, given that the approach in [32, 33] solves a deterministic convex program with sampled constraints, while we process samples in an online fashion.

We illustrate properties of our algorithm through a stylized example. Our experiments reveal that the optimized iteration count (and sample complexity) for even a simple example is quite high. This limitation is unfortunately common for subgradient algorithms and likely cannot be overcome in optimizing general nonsmooth functions that we study. While the bounds are order-optimal, our numerical experiments reveal that a solution with desired risk tolerance can be found in less iterations than obtained from the upper bound. This is an artifact of optimizing step-sizes based on upper bounds on suboptimality and constraint violation. We end the paper in Section 4 with discussions on possible extensions of our analysis.

Very recently, it was brought to our attention that the work in [34] done concurrently presents a related approach to tackle optimization of composite nonconvex functions under related but different assumptions. In fact, their work appeared at the same time as our early version and claims a similar result that does not require bounds on the dual variables. Our analysis does not require or analyze the case with strongly convex functions within our setup and therefore Nesterov-style acceleration remains untenable. As a result, our algorithm is different. Our focus on ${\rm CVaR}$ permits us to further analyze the growth in optimized sample complexity with risk aversion and its connection to chance-constrained optimization, that is quite different.

2 Algorithm for ${{{\cal P}}^{\rm E}}$ and its analysis

We present the primal-dual stochastic subgradient method to solve ${{{\cal P}}^{\rm E}}$ in Algorithm 1.

The notation $\langle\cdot,\cdot\rangle$ stands for the usual inner product in Euclidean space and $\|\cdot\|$ denotes the induced $\ell_{2}$-norm. Here, $\nabla h({\bm{x}})$ stands for a subgradient of an arbitrary convex function $h$ at ${\bm{x}}$. For our analysis, the subgradient in Algorithm 1 for functions $f_{\omega}({\bm{x}})$ and $g_{\omega}({\bm{x}})$ can be arbitrary elements of the closed convex subdifferential sets $\partial f_{\omega}({\bm{x}})$ and $\partial g_{\omega}({\bm{x}})$, respectively. We assume that these subdifferential sets are nonempty everywhere in $\mathbb{X}$.

The primal-dual method in Algorithm 1 leverages Lagrangian duality theory. Specifically, define the Lagrangian function for ${{{\cal P}}^{\rm E}}$ as

[TABLE]

for ${\bm{x}}\in\mathbb{X}$, ${\bm{z}}\in\mathbb{R}^{m}_{+}$, where {\cal L}_{\omega}({\bm{x}},{\bm{z}})\coloneqq f_{\omega}({\bm{x}})+{\bm{z}}^{{\mathchoice{\raisebox{0.75346pt}{\displaystyle\intercal}}{\raisebox{0.75346pt}{\textstyle\intercal}}{\raisebox{0.75346pt}{\scriptstyle\intercal}}{\raisebox{0.75346pt}{\scriptscriptstyle\intercal}}}}{\bm{g}}_{\omega}({\bm{x}}). Then, ${{{\cal P}}^{\rm E}}$ admits the standard reformulation as a min-max problem of the form

[TABLE]

Denote its optimal set by $\mathbb{X}_{\star}\subseteq\mathbb{X}$. Define the dual problem of ${{{\cal P}}^{\rm E}}$ as

[TABLE]

Denote its optimal set by $\mathbb{Z}_{\star}\subseteq\mathbb{R}^{m}_{+}$. Weak duality then guarantees ${{p}^{\rm E}_{\star}}\geq{{d}^{\rm E}_{\star}}$. When the inequality is met with an equality, the problem is said to satisfy strong duality. A point $({\bm{x}}_{\star},{\bm{z}}_{\star})\in\mathbb{X}\times\mathbb{R}^{m}_{+}$ is a saddle point of ${\cal L}$ if

[TABLE]

for all $({\bm{x}},{\bm{z}})\in\mathbb{X}\times\mathbb{R}^{m}_{+}$. The following well-known saddle point theorem (see [35, Theorem 2.156]) relates saddle points with primal-dual optimal solutions.

Theorem** (Saddle point theorem).**

A saddle point of ${\cal L}$ exists if and only if ${{{\cal P}}^{\rm E}}$ satisfies strong duality, i.e., ${{p}^{\rm E}_{\star}}={{d}^{\rm E}_{\star}}$. Moreover, the set of saddle points of ${\cal L}$ is given by $\mathbb{X}_{\star}\times\mathbb{Z}_{\star}$.

Our convergence analysis of Algorithm 1 requires ${{{\cal P}}^{\rm E}}$ to satisfy the following properties.

Assumption 2.1**.**

1. (a)

*Subgradients of $F$ and ${\bm{G}}$ are bounded, i.e.,

$\|\nabla F({\bm{x}})\|\leq C_{F}$, $\|\nabla G^{i}({\bm{x}})\|\leq C_{G}^{i}$ for each $i=1,\dots,m$ and all ${\bm{x}}\in\mathbb{X}$.* 2. (b)

$\nabla f_{\omega}$* and $\nabla g_{\omega}^{i}$ for $i=1,\dots,m$ have bounded variance, i.e.,

$\mathds{E}\|\nabla f_{\omega}({\bm{x}})-\mathds{E}[\nabla f_{\omega}({\bm{x}})]\|^{2}\leq\sigma_{F}^{2}$ and $\mathds{E}\|\nabla g_{\omega}^{i}({\bm{x}})-\mathds{E}[\nabla g_{\omega}^{i}({\bm{x}})]\|^{2}\leq[\sigma_{G}^{i}]^{2}$ for all ${\bm{x}}\in\mathbb{X}$.* 3. (c)

${\bm{g}}_{\omega}({\bm{x}})$* has a bounded second moment, i.e., $\mathds{E}\|g_{\omega}^{i}({\bm{x}})\|^{2}\leq[D_{G}^{i}]^{2}$ for all ${\bm{x}}\in\mathbb{X}$.* 4. (d)

${\cal L}$* for ${{{\cal P}}^{\rm E}}$ admits a saddle point $({\bm{x}}_{\star},{\bm{z}}_{\star})\in\mathbb{X}\times\mathbb{R}^{m}_{+}$.*

The subgradient of $F$ and the variance of its noisy estimate are assumed bounded. Such an assumption is standard in the convergence analysis of unconstrained stochastic subgradient methods. The assumptions regarding ${\bm{G}}$ are similar, but we additionally require the second moment of the noisy estimate of ${\bm{G}}$ to be bounded over $\mathbb{X}$. Boundedness of ${\bm{G}}$ in primal-dual subgradient methods has appeared in prior literature, e.g., in [24, 25]. The second moment remains bounded if ${g}^{i}_{\omega}$ is uniformly bounded over $\mathbb{X}$ and $\Omega$ for each $i$. It is also satisfied if ${\bm{G}}$ remains bounded over $\mathbb{X}$ and its noisy estimate has a bounded variance. Convergence analysis of unconstrained optimization problems typically assumes the existence of a finite optimal solution. We extend that requirement to the existence of a saddle point in the primal-dual setting, which by the saddle point theorem is equivalent to the existence of finite primal and dual optimal solutions. A variety of conditions imply the existence of such a point; the next result delineates two such sufficient conditions in (a) and (b), where (a) implies (b).

Lemma 2.1** (Sufficient conditions for existence of a saddle point).**

${\cal L}$* for ${{{\cal P}}^{\rm E}}$ admits a saddle point, if either of the following conditions hold:*

1. (a)

$\mathbb{X}_{\star}$* is nonempty, ${{p}^{\rm E}_{\star}}$ is finite and Slater’s constraint qualification holds, i.e., there exists ${\bm{x}}$ in the relative interior of $\mathbb{X}$ for which ${\bm{G}}({\bm{x}})<0$.

2. (b)

${{{\cal P}}^{\rm E}}$* admits a finite (${\bm{x}}_{\star},{\bm{z}}_{\star})\in\mathbb{X}\times\mathbb{R}_{+}^{m}$ that satisfies the Karush-Kuhn-Tucker (KKT) conditions given by*

[TABLE]

for $i=1,\dots,m$, where ${\cal N}_{\mathbb{X}}({\bm{x}}_{\star})$ denotes the normal cone of $\mathbb{X}$ at ${\bm{x}}_{\star}$.

Proof.

Part (a) is a direct consequence of [35, Theorem 1.265]. To prove part (b), notice that (11) ensures the existence of subgradients $\nabla F({\bm{x}}_{\star})\in\partial F({\bm{x}}_{\star})$, $\nabla{G}^{i}({\bm{x}}_{\star})\in\partial{G}^{i}({\bm{x}}_{\star})$, $i=1,\ldots,m$ and ${\bm{n}}\in{\cal N}_{\mathbb{X}}({\bm{x}}_{\star})$ for which

[TABLE]

Then, for any ${\bm{x}}\in\mathbb{X}$, we have

[TABLE]

The inequalities in the above relation follow from the convexity of $F$ and ${G}^{i}$’s, nonnegativity of ${\bm{z}}_{\star}$, and the definition of the normal cone. From the above inequalities, we conclude ${\cal L}({\bm{x}},{\bm{z}}_{\star})\geq{\cal L}({\bm{x}}_{\star},{\bm{z}}_{\star})$ for all ${\bm{x}}\in\mathbb{X}$. Furthermore, for any ${\bm{z}}\geq 0$, we have

[TABLE]

where the last step follows from the nonnegativity of ${\bm{z}}$ and (11), completing the proof. ∎

We now present our first main result that provides a bound on the distance to optimality and constraint violation at a weighted average of the iterates generated by the algorithm on ${{{\cal P}}^{\rm E}}$ under Assumption 2.1. Denote by ${\bm{C}}_{G}$, ${\bm{D}}_{G}$, and ${\bm{\sigma}}_{G}$ the collections of $C_{G}^{i}$, $D_{G}^{i}$, and $\sigma_{G}^{i}$ respectively. We make use of the following notation.

[TABLE]

Theorem 2.2** (Convergence result for ${{{\cal P}}^{\rm E}}$).**

Suppose Assumption 2.1 holds. For a positive sequence $\{\gamma_{k}\}_{k=1}^{K}$, if $P_{3}\sum_{k=1}^{K}\gamma_{k}^{2}<1$, then the iterates generated by Algorithm 1 satisfy

[TABLE]

for each $i=1,\dots,m$, where $\bar{{\bm{x}}}_{K+1}:=\frac{\sum_{k=1}^{K}\gamma_{k}{\bm{x}}_{k+1}}{\sum_{k=1}^{K}\gamma_{k}}$. Moreover, if $\gamma_{k}=\gamma/\sqrt{K}$ for $k=1,\ldots,K$ with $0<\gamma<P_{3}^{-1/2}$, then

[TABLE]

for $i=1,\ldots,m$, where $\eta:=\frac{P_{1}+P_{2}\gamma^{2}}{4\gamma(1-P_{3}\gamma^{2})}$.

A constant step-size of $\eta/\sqrt{K}$ over a fixed number of $K$ iterations yields the ${\cal O}(1/\sqrt{K})$ decay rate in the expected distance to optimality and constraint violation of Algorithm 1. This is indeed order optimal, as implied by Nesterov’s celebrated result in [26, Theorem 3.2.1].

Remark 2.3**.**

While we present the proof for an i.i.d. sequence of samples, we believe that the result can be extended to the case where $\omega$’s follow a Markov chain with geometric mixing rate following the technique in [36]. For such settings, the expectations in the definition of ${F},{\bm{G}}$ should be computed with respect to the stationary distribution of the chain. The results will then possibly apply to Markov decision processes with applications in stochastic control.

Given that the literature on primal-dual subgradient methods is extensive, it is important for us to relate and distinguish Algorithm 1 and Theorem 2.2 with prior work. Using the Lagrangian in (6), Algorithm 1 can be written as

[TABLE]

where $\operatorname{proj}_{\mathbb{A}}$ projects its argument on set $\mathbb{A}$. The vectors $\nabla_{x}{\cal L}_{\omega}$ and $\nabla_{z}{\cal L}_{\omega}$ are stochastic subgradients of the Lagrangian function with respect to ${\bm{x}}$ and ${\bm{z}}$, respectively. Therefore, Algorithm 1 is a projected stochastic subgradient algorithm that seeks to solve the saddle-point reformulation of ${{{\cal P}}^{\rm E}}$ in (7). Implicit in our algorithm is the assumption that projection on $\mathbb{X}$ is computationally easy. Any functional constraints describing $\mathbb{X}$ that makes such projection challenging should be included in ${\bm{G}}$.

Closest in spirit to our work on ${{{\cal P}}^{\rm E}}$ are the papers by Baes et al. in [37], Yu et al. in [25], Xu in [24], and Nedic and Ozdaglar in [23]. Stochastic mirror-prox algorithm in [37] and projected subgradient method in [23] are similar in their updates to ours except in two ways. First, these algorithms in the context of ${{{\cal P}}^{\rm E}}$ update the dual variable ${\bm{z}}_{k}$ based on ${\bm{G}}$ or its noisy estimate evaluated at ${\bm{x}}_{k}$, while we update it based on the estimate at ${\bm{x}}_{k+1}$. Second, both project the dual variable on a compact subset of $\mathbb{R}^{m}_{+}$ that contains the optimal set of dual multipliers. While authors in [37] assume an a priori set to project on, authors in [23] compute such a set from a “Slater point” that satisfies ${\bm{G}}({\bm{x}})<0$. Specifically, Slater’s condition guarantees that the set of optimal dual solutions $\mathbb{Z}_{\star}$ is bounded (see [35, Theorem 1.265], [38]). Moreover, a Slater point can be used to construct a compact set that contains $\mathbb{Z}_{\star}$, e.g., using [23, Lemma 4.1]. While one can project dual variables on such a set in each iteration, execution of the algorithm then requires a priori knowledge of such a point. We do not assume knowledge of such a point (or any explicit bound on $\mathbb{Z}_{\star}$) to execute Algorithm 1. Rather, our proof provides an explicit bound on the growth of the dual variable sequence for Algorithm 1, much in line with Xu’s analysis in [24]. Much to our surprise, a minor modification of using a Gauss-Seidel style dual update as opposed to the popular Jacobi style dual update obviates the need for this crucial assumption in the literature for the proofs to work. Unfortunately, our Gauss-Seidel style dual update comes at an additional cost of an extra sample required per iteration of the primal-dual algorithm, making the sample complexity double of the iteration complexity. The constant factor of two, however, does not impact the order-wise complexity. We surmise that the additional sample and the Gauss-Seidel update of the dual variable helps to decouple the analysis of the primal and dual updates and points to a possible extension of our result to an asynchronous setting, often useful in engineering applications.

While sharing some parallels, our work has an important difference with that in [24]. Xu considers a collection of deterministic constraint functions, i.e., ${\bm{g}}^{\omega}$ is identical for all $\omega\in\Omega$, and considers a modified augmented Lagrangian function of the form $\tilde{{\cal L}}({\bm{x}},{\bm{z}}):=F({\bm{x}})+\frac{1}{m}\sum_{i=1}^{m}\varphi_{\delta}({\bm{x}},z^{i})$, where

[TABLE]

for $i=1,\ldots,m$ with a suitable time-varying sequence of $\delta$’s. His algorithm is similar to Algorithm 1 but performs a randomized coordinate update for the dual variable instead of (5). To the best of our knowledge, Xu’s analysis in [24] with such a Lagrangian function does not directly apply to our setting with stochastic constraints that is crucial for the subsequent analysis of the risk-sensitive problem ${{{\cal P}}^{\rm CVaR}}$.

Finally, Yu et al.’s work in [25] provides an analysis of the algorithm that updates its dual variables using

[TABLE]

where ${v}^{i}:={g}^{i}_{\omega_{k}}({\bm{x}}_{k})+\langle\nabla{g}^{i}_{\omega_{k}}({\bm{x}}_{k}),{\bm{x}}_{k+1}-{\bm{x}}_{k}\rangle$ for $i=1,\ldots,m$. In contrast, our ${\bm{z}}$-update in (5) samples $\omega_{k+1/2}$ and sets ${\bm{v}}_{k}:={\bm{g}}_{\omega_{k+1/2}}({\bm{x}}_{k+1})$ at the already computed point ${\bm{x}}_{k+1}$. We are able to recover the ${\cal O}(1/\sqrt{K})$ decay rate of suboptimality and constraint violation with a proof technique much closer to the classical analysis of subgradient methods in [39, 23]. Unlike [25], we provide a clean characterization of the constant $\eta$ in (20) that is crucial to study the growth in sample (and iteration) complexity of Algorithm 1 applied to a reformulation of ${{{\cal P}}^{\rm CVaR}}$.

2.1 Proof of Theorem 2.2

The proof proceeds in three steps.

1. (a)

We establish the following dissipation inequality that consecutive iterates of the algorithm satisfy.

[TABLE]

for any ${\bm{x}}\in\mathbb{X}$ and ${\bm{z}}\in\mathbb{R}_{+}^{m}$. 2. (b)

Next, we bound $\mathds{E}\|{\bm{z}}_{k}\|^{2}$ generated by our algorithm from above using step (a) as

[TABLE]

for $k=1,\ldots,K$, where $A_{K}:=\sum_{k=1}^{K}\gamma_{k}^{2}$. 3. (c)

We combine the results in steps (a) and (b) to complete the proof.

Define the filtration $\mathpzc{W}_{1}\subset\mathpzc{W}_{1+1/2}\subset\mathpzc{W}_{2}\subset\ldots$, where $\mathpzc{W}_{k}$ is the $\sigma$-algebra generated by the samples $\omega_{1},\ldots\omega_{k-1/2}$ for $k$ being multiples of $1/2$, starting from unity. Then, $\{{\bm{x}}_{1},{\bm{z}}_{1},\dots,{\bm{x}}_{k},{\bm{z}}_{k}\}$ becomes $\mathpzc{W}_{k}$-measurable, while $\{{\bm{x}}_{1},{\bm{z}}_{1},\dots,{\bm{x}}_{k},{\bm{z}}_{k},{\bm{x}}_{k+1}\}$ is $\mathpzc{W}_{k+1/2}$-measurable.

$\bullet$ Step (a) – Proof of (24): We first utilize the ${\bm{x}}$-update in (4) to prove

[TABLE]

for all ${\bm{x}}\in\mathbb{X}$. Then, we utilize the ${\bm{z}}$-update in (5) to prove

[TABLE]

for all ${\bm{z}}\in\mathbb{R}^{m}_{+}$. The law of total probability is then applied to the sum of (26) and (27) followed by a multiplication by $\gamma_{k}$ yielding the desired result in (24).

Proof of (26):

The ${\bm{x}}$-update in (4) yields

[TABLE]

We now simplify the inner product. The product with $\nabla f_{\omega}({\bm{x}}_{k})$ can be expressed as

[TABLE]

where $\nabla F({\bm{x}}_{k+1})$ denotes a subgradient of $F$ at ${\bm{x}}_{k+1}$. The inequality for the first term follows from the convexity of $F$. Since $\mathds{E}[\nabla f_{\omega}({\bm{x}}_{k})|\mathpzc{W}_{k}]\in\partial F({\bm{x}}_{k})$ from [40], the expectation of the second summand on the right hand side (RHS) of (29) satisfies

[TABLE]

Taking expectations in (29), the above relation implies

[TABLE]

Next, we bound the inner product with the second term on the RHS of (29). To that end, utilize the convexity of member functions in ${\bm{g}}_{\omega}$ and ${\bm{G}}$ along the above lines to infer

[TABLE]

To tackle the inner product with the third term in the RHS of (28), we use the identity

[TABLE]

The inequalities in (31), (32), and the equality in (33) together gives

[TABLE]

To simplify the above relation, apply Young’s inequality to obtain

[TABLE]

Recall that $\mathds{E}[\nabla f_{\omega}({\bm{x}}_{k})|\mathpzc{W}_{k}]\in\partial F({\bm{x}}_{k})$, subgradients of $F$ are bounded and $\nabla f_{\omega}$ has bounded variance. Therefore, we infer from the above inequality that

[TABLE]

Appealing to Young’s inequality $m$ times and a similar line of argument as above gives

[TABLE]

Leveraging the relations in (36) and (37) in (34), we get

[TABLE]

that upon simplification gives (26).

Proof of (27):

From the ${\bm{z}}$-update in (5), we obtain

[TABLE]

for all ${\bm{z}}\geq 0$. Again, we deal with the two summands in the second factor of the inner product of (39) separately. The expectation of the inner product with the first term yields

[TABLE]

In the above derivation, we have utilized Young’s inequality and the boundedness of the second moment of ${\bm{g}}_{\omega}$. Since $\mathpzc{W}_{k}\subset\mathpzc{W}_{k+{1}/{2}}$, the law of total probability can be used to condition (40) on $\mathpzc{W}_{k}$ rather than on $\mathpzc{W}_{k+{1}/{2}}$. To simplify the inner product with the second term in (39), we use the identity

[TABLE]

Utilizing (40) and (41) in (39) gives (27). Adding (26) and (27) followed by a multiplication by $\gamma_{k}$ yields

[TABLE]

Taking the expectation and applying the law of total probability completes the proof of (24).

$\bullet$ Step (b) – Proof of (25): Plugging $({\bm{x}},{\bm{z}})=({\bm{x}}_{\star},{\bm{z}}_{\star})$ in the inequality for the one-step update in (24) and summing it over $k=1,\ldots,\kappa$ for $\kappa\leq K$ gives

[TABLE]

for $\kappa=1,\ldots,K$. The above then yields

[TABLE]

Notice that $2\mathds{E}\|{\bm{z}}_{\kappa+1}-{\bm{z}}_{\star}\|^{2}+2\|{\bm{z}}_{\star}\|^{2}\geq\mathds{E}\|{\bm{z}}_{\kappa+1}\|^{2}$. This inequality and ${\bm{z}}_{1}=0$ in (44) gives

[TABLE]

We argue the bound on $\mathds{E}\|{\bm{z}}_{k}\|^{2}$ for $k=1,\ldots,K$ inductively. Since ${\bm{z}}_{1}=0$, the base case trivially holds. Assume that the bound holds for $k=1,\ldots,\kappa$ for $\kappa<K$. With the notation $A_{K}=\sum_{k=1}^{K}\gamma_{k}^{2}$, the relation in (45) implies

[TABLE]

completing the proof of step (b).

$\bullet$ Step (c) – Combining steps (a) and (b) to prove Theorem 2.2: For any ${\bm{z}}\geq 0$, the inequality in (24) with ${\bm{x}}={\bm{x}}_{\star}$ from step (a) summed over $k=1,\ldots,K$ gives

[TABLE]

Using ${\bm{z}}_{1}=0$ and an appeal to the saddle point property of $({\bm{x}}_{\star},{\bm{z}}_{\star})$ yields

[TABLE]

In deriving the above inequality, we have utilized the bound on $\mathds{E}\|{\bm{z}}_{k}\|^{2}$ from step (b) and the definition of $P_{1}$ and $A_{K}$. To further simplify the above inequality, notice that the saddle point property of $({\bm{x}}_{\star},{\bm{z}}_{\star})$ in (9) yields

[TABLE]

which implies {\bm{z}}_{\star}^{\mathchoice{\raisebox{0.75346pt}{\displaystyle\intercal}}{\raisebox{0.75346pt}{\textstyle\intercal}}{\raisebox{0.75346pt}{\scriptstyle\intercal}}{\raisebox{0.75346pt}{\scriptscriptstyle\intercal}}}{\bm{G}}({\bm{x}}_{\star})\geq 0. However, the saddle point theorem guarantees that ${\bm{x}}_{\star}$ is an optimizer of ${{{\cal P}}^{\rm E}}$, meaning that ${\bm{x}}_{\star}$ is feasible and ${\bm{G}}({\bm{x}}_{\star})\leq 0$, implying {\bm{z}}_{\star}^{\mathchoice{\raisebox{0.75346pt}{\displaystyle\intercal}}{\raisebox{0.75346pt}{\textstyle\intercal}}{\raisebox{0.75346pt}{\scriptstyle\intercal}}{\raisebox{0.75346pt}{\scriptscriptstyle\intercal}}}{\bm{G}}({\bm{x}}_{\star})\leq 0 as ${\bm{z}}_{\star}\in\mathbb{R}^{m}_{+}$. Taken together, we infer

[TABLE]

Since ${\cal L}({\bm{x}},{\bm{z}})$ is convex in ${\bm{x}}$, Jensen’s inequality and (50) implies

[TABLE]

where recall that $\bar{{\bm{x}}}_{K+1}$ is the $\gamma$-weighted average of the iterates. Utilizing (51) in (48), we get

[TABLE]

The above relation defines a bound on $\mathds{E}[{\cal L}(\bar{{\bm{x}}}_{K+1},{\bm{z}})]$ for every ${\bm{z}}\geq 0$. Choosing ${\bm{z}}=0$ and noting $\|\mathds{1}+{\bm{z}}_{\star}\|^{2}\geq 0$, we get the bound on expected suboptimality in (18). To derive the bound on expected constraint violation in (19), notice that the saddle point property in (9) and (50) imply

[TABLE]

where $\mathds{1}^{i}\in\mathbb{R}^{m}$ is a vector of all zeros except the $i$-th entry that is unity. Choosing ${\bm{z}}=\mathds{1}^{i}+{\bm{z}}_{\star}$ in (52) and the observation in (53) then gives

[TABLE]

for each $i=1,\ldots,m$. This completes the proof of (19). The bounds in (20) are immediate from that in (18)–(19). This completes the proof of Theorem 2.2.

Remark 2.4**.**

The bound in (19) can be sharpened to

[TABLE]

using ${\bm{z}}$ defined by ${\bm{z}}^{i}\coloneqq{\bm{z}}^{i}_{\star}+\operatorname{\mathbb{I}}_{\{{G}^{i}(\bar{{\bm{x}}}_{K+1})>0\}}$ for $i=1,\ldots,m$ in (52). Here, $\operatorname{\mathbb{I}}_{\{A\}}$ is the indicator function, evaluating to $1$ if $A$ holds and [math] otherwise. This improved bound was suggested to us by an anonymous reviewer. Notice that (55) suggests a much tighter bound on the expected constraint violation per constraint than (19) when $m$ is large.

In what follows, we offer insights into two specific aspects of our proof. First, we present our conjecture on where the Gauss-Seidel nature of our dual update obtained with an extra sample helps us circumvent the need for an a priori bound on the dual variable. Notice that our dual update allows us to derive the third line of (40) that ultimately yields the term -{\bm{z}}_{k}^{{\mathchoice{\raisebox{0.75346pt}{\displaystyle\intercal}}{\raisebox{0.75346pt}{\textstyle\intercal}}{\raisebox{0.75346pt}{\scriptstyle\intercal}}{\raisebox{0.75346pt}{\scriptscriptstyle\intercal}}}}{\bm{G}}({\bm{x}}_{k+1}) in (27). This term conveniently disappears when (27) is added to the inequality in (26) obtained from the primal update. We conjecture that this cancellation made possible by our dual update makes the theoretical analysis particularly easy. We anticipate that the classical Jacobi-style dual iteration derived with one sample shared within the primal and the dual steps will not lead to said cancellation and yield a term of the form {\bm{z}}_{k}^{{\mathchoice{\raisebox{0.75346pt}{\displaystyle\intercal}}{\raisebox{0.75346pt}{\textstyle\intercal}}{\raisebox{0.75346pt}{\scriptstyle\intercal}}{\raisebox{0.75346pt}{\scriptscriptstyle\intercal}}}}\left[{\bm{G}}({\bm{x}}_{k+1})-{\bm{G}}({\bm{x}}_{k})\right]. Bounding the growth of such a term might prove challenging without an available bound on $\|{\bm{z}}_{k}\|$ and will likely require a different argument. A detailed comparison between the proof techniques of the Jacobi and the Gauss-Seidel updates is left for future endeavors.

Second, we comment on the presence of a dimensionless constant $\mathds{1}$ in $P_{1}$ together with ${\bm{z}}_{\star}$. We use the inequality in (24) to establish (52) that is valid at all ${\bm{z}}\geq 0$. Inspired by arguments in [24], we then utilize (52) not only at the dual iterate ${\bm{z}}_{k}$–that is often the case with many prior analyses–but also at ${\bm{z}}=0$ and ${\bm{z}}=\mathds{1}^{i}+{\bm{z}}_{\star}$. Specifically, the nature of the Lagrangian function ${\cal L}({\bm{x}},{\bm{z}})$ in ${\bm{z}}$ permits us to relate these evaluations at ${\bm{z}}=0$ and ${\bm{z}}=\mathds{1}^{i}+{\bm{z}}_{\star}$ to the extents of suboptimality and constraint violation, respectively, using

[TABLE]

The deliberate inclusion of $\|{\mathds{1}+{\bm{z}}_{\star}}\|^{2}$ in the constant $P_{1}$ aids in drowning the effect of the term $\frac{1}{2}\|{{\bm{z}}}\|^{2}$ in (52) evaluated at ${\bm{z}}=\mathds{1}^{i}+{\bm{z}}_{\star}$ when deriving the bound on the extent of constraint violation, without impacting the same when evaluated at ${\bm{z}}=0$, used in deriving the bound on the extent of suboptimality.

2.2 Optimal step size selection

We exploit the bounds in Theorem 2.2 to select a step size that minimizes the iteration count to reach an $\varepsilon$-approximately feasible and optimal solution to ${{{\cal P}}^{\rm E}}$ and solve333The integrality of $K$ is ignored for notational convenience.

[TABLE]

The following characterization of optimal step sizes and the resulting iteration count from Proposition 2.5 will prove useful in studying the growth in iteration complexity in solving ${{{\cal P}}^{\rm CVaR}}$ with the risk-aversion parameters $\alpha,{\bm{\beta}}$ in the following section.

Proposition 2.5**.**

For any $\varepsilon>0$, the optimal solution of (57) satisfies

[TABLE]

where $y=1+\frac{P_{2}}{P_{1}P_{3}}$.

Proof.

It is evident from (58) that $\gamma_{\star}^{2}<P_{3}^{-1}$. Then, it suffices to show that $\gamma_{\star}$ from (58) minimizes

[TABLE]

over $\gamma>0$. To that end, notice that

[TABLE]

The above derivative is negative at $\gamma=0^{+}$ and vanishes only at $\gamma_{\star}$ over positive values of $\gamma$, certifying it as the global minimizer. ∎

Parameter $P_{1}$ is generally not known a priori. However, it is often possible to bound it from above. One can calculate $\gamma_{\star}$ and $K_{\star}$ using (58), replacing $P_{1}$ with its overestimate. Notice that

[TABLE]

It is straightforward to verify that $\frac{\partial K_{\star}}{\partial P_{1}}>0$, $\frac{dy}{dP_{1}}\leq 0$, and $\frac{\partial\gamma_{\star}}{\partial y}\leq 0$, and hence, overestimating $P_{1}$ results in a smaller $\gamma_{\star}$. Finally, $\frac{\partial K_{\star}}{\partial\gamma}>0$ for $\gamma>\gamma_{\star}$, implying that $K_{\star}$ calculated with an overestimate of $P_{1}$ is larger than the optimal iteration count–the computational burden we must bear for not knowing $P_{1}$. Our algorithm does require knowledge of $P_{3}$ to implement the algorithm, that in turn depend only on the nature of the functions defining the constraints and not a primal-dual optimizer.

2.3 Asymptotic almost sure convergence with decaying step-sizes

Subgradient methods are often studied with decaying non-summable square-summable step sizes, for which they converge to an optimizer in the unconstrained setting. The result holds even for distributed variants and for mirror descent methods (see [41]). Establishing convergence of Algorithm 1 to a primal-dual optimizer of ${{{\cal P}}^{\rm E}}$ is much more challenging without assumptions of strong convexity in the objective. With such step-sizes, we provide the following result to guarantee the stability of our algorithm, which is reminiscent of [42, Theorem 4].

Proposition 2.6**.**

Suppose Assumption 2.1 holds and $\{\gamma_{k}\}_{k=1}^{\infty}$ is a non-summable square-summable nonnegative sequence, i.e., $\sum_{k=1}^{\infty}\gamma_{k}=\infty,\sum_{k=1}^{\infty}\gamma_{k}^{2}<\infty$. Then, $({\bm{x}}_{k},{\bm{z}}_{k})$ generated by Algorithm 1 remains bounded and $\lim_{k\to\infty}{\cal L}({\bm{x}}_{k},{\bm{z}}_{\star})-{\cal L}({\bm{x}}_{\star},{\bm{z}}_{k})=0$ almost surely.

This ‘gap’ function ${\cal L}({\bm{x}},{\bm{z}}_{\star})-{\cal L}({\bm{x}}_{\star},{\bm{z}})$ looks notoriously similar to the duality gap at $({\bm{x}},{\bm{z}})$, but is not the same. We are unaware of any results on asymptotic almost sure convergence of primal-dual first-order algorithms to an optimizer for constrained convex programs with convex, but not necessarily strongly convex, objectives. A recent result in [43] establishes such a convergence in primal-dual dynamics in continuous time; our attempts at leveraging discretizations of the same have yet proven unsuccessful.

The proof of Proposition 2.6 takes advantage of the one-step update in (24), that makes it amenable to the well-studied almost supermartingale convergence result by Robbins and Siegmund in [19, Theorem 1].

Theorem** (Convergence of almost supermartingales).**

Let $m_{k},n_{k},r_{k},s_{k}$ be ${\cal F}_{k}$-measurable finite nonnegative random variables, where ${\cal F}_{1}\subseteq{\cal F}_{2}\subseteq\ldots$ describes a filtration. If $\sum_{k=1}^{\infty}s_{k}<\infty$, $\sum_{k=1}^{\infty}r_{k}<\infty$, and

[TABLE]

then $\lim_{k\to\infty}m_{k}$ exists and is finite and $\sum_{k=1}^{\infty}n_{k}<\infty$ almost surely.

Proof of Proposition 2.6.

Using notation from the proof of Theorem 2.2, (26) and (27) together yields

[TABLE]

We utilize the above to derive a similar inequality replacing $\mathds{E}[{\cal L}({\bm{x}}_{k+1},{\bm{z}}_{\star})|\mathpzc{W}_{k}]$ with ${\cal L}({\bm{x}}_{k},{\bm{z}}_{\star})$ by bounding the difference between them. Then, we apply the almost supermartingale convergence theorem to the result to conclude the proof. To bound said difference, the convexity of ${\cal L}$ in ${\bm{x}}$ and Young’s inequality together imply

[TABLE]

where ${\nabla}_{x}{\cal L}$ denotes a subgradient of ${\cal L}$ w.r.t. ${\bm{x}}$. To further bound the RHS of (64), Assumption 2.1 allows us to deduce

[TABLE]

for any ${\bm{x}}\in\mathbb{X}$. Furthermore, the ${\bm{x}}$-update in (21) and the non-expansive nature of the projection operator yield

[TABLE]

From Assumption 2.1, we get

[TABLE]

and along similar lines

[TABLE]

that together in (66) yield

[TABLE]

Combining the above with (65) in (64) gives

[TABLE]

Adding (70) to (63) and simplifying, we obtain

[TABLE]

The above inequality with

[TABLE]

becomes (62), where

[TABLE]

Each term is nonnegative, owing to (9), and ${\bm{\gamma}}$ defines a square summable sequence. Applying [19, Theorem 1], $m_{k}$ converges to a constant and $\sum_{k=1}^{\infty}n_{k}<\infty$. The latter combined with the non-summability of ${\bm{\gamma}}$ implies the result. ∎

3 Algorithm for ${{{\cal P}}^{\rm CVaR}}$ and its analysis

We now devote our attention to solving ${{{\cal P}}^{\rm CVaR}}$ via a primal-dual algorithm. To do so, we reformulate it as an instance of ${{{\cal P}}^{\rm E}}$ and utilize Algorithm 1 to solve that reformulation with constant step-sizes under a stronger set of assumptions given below.

Assumption 3.1**.**

1. (a)

*Subgradients of $F$ and ${\bm{G}}$ are bounded, i.e.,

$\|\nabla f_{\omega}({\bm{x}})\|\leq C_{F}$ and $\|\nabla g_{\omega}^{i}({\bm{x}})\|\leq C_{G}^{i}$ almost surely for all ${\bm{x}}\in\mathbb{X}$.* 2. (b)

${\bm{g}}_{\omega}(x)$* is bounded, i.e., $\|g^{i}_{\omega}({\bm{x}})\|\leq D_{G}^{i}$ for all ${\bm{x}}\in\mathbb{X}$, almost surely.* 3. (c)

${\cal L}$* for ${{{\cal P}}^{\rm CVaR}}$ admits a saddle point $({\bm{x}}_{\star},{\bm{z}}_{\star})\in\mathbb{X}\times\mathbb{R}^{m}_{+}$.444Lemma 2.1 provides sufficient conditions for the existence of such a saddle point.*

Using the variational characterization (2) of ${\rm CVaR}$, rewrite ${{{\cal P}}^{\rm CVaR}}$ as

[TABLE]

where ${\psi^{h}_{\omega}({\bm{x}},u;\delta)}\coloneqq{u+\frac{1}{1-\delta}[h_{\omega}({\bm{x}})-u]^{+}}$ for any collection of convex functions $h_{\omega}:\mathbb{R}^{n}\to\mathbb{R}$, $\omega\in\Omega$. Coupled with Assumption 3.1, we will show that we can bound $|u^{i}|\leq D_{G}^{i}$ for each $i=1,\dots,m$555${\rm CVaR}$ of any random variable can only vary between the mean and the maximum value that random variable can take., that allows us to rewrite ${{{\cal P}}^{\rm CVaR}}$ as

[TABLE]

where $|\cdot|$ denotes the element-wise absolute value. Call the optimal value of ${{{\cal P}}^{\rm CVaR}}$ as $p^{\rm CVaR}_{\star}$ in the sequel.

Theorem 3.1** (Convergence result for ${{{\cal P}}^{\rm CVaR}}$).**

Suppose Assumption 3.1 holds. The iterates generated by Algorithm 1 on ${{{\cal P}}^{\rm E}}^{\prime}$ for ${{{\cal P}}^{\rm CVaR}}$ with parameters $\alpha,{\bm{\beta}}$ satisfy

[TABLE]

*for $i=1,\ldots,m$ with step sizes $\gamma_{k}=\gamma/\sqrt{K}$ for $k=1,\dots,K$ with $0<\gamma<{P_{3}^{-1/2}(\alpha,{\bm{\beta}})}$, where $\eta(\alpha,{\bm{\beta}})\coloneqq\frac{P_{1}+\gamma^{2}P_{2}(\alpha,{\bm{\beta}})}{4\gamma(1-\gamma^{2}P_{3}(\alpha,{\bm{\beta}}))}$ and

[TABLE]

Proof.

We prove the result in the following steps.

1. (a)

Under Assumption 3.1, we revise $P_{2}$ and $P_{3}$ in Theorem 2.2 for ${{{\cal P}}^{\rm E}}$. 2. (b)

We show that if $f_{\omega},{\bm{g}}_{\omega}$ satisfy Assumption 3.1, then $\psi^{f}_{\omega}$ and $\psi^{g^{i}}_{\omega},i=1,\ldots,m$ satisfy Assumption 3.1, but with different bounds on the gradients and function values. Leveraging these bounds, we obtain $P_{2}(\alpha,{\bm{\beta}})$ and $P_{3}(\alpha,{\bm{\beta}})$ for ${{{\cal P}}^{\rm E}}^{\prime}$ using step (a). 3. (c)

Using Assumption 3.1, we prove that the Lagrangian function ${\cal L}^{\prime}:\mathbb{X}\times\mathbb{R}\times\mathbb{U}\times\mathbb{R}^{m}_{+}\to\mathbb{R}$ defined as

[TABLE]

admits a saddle point in $\mathbb{X}\times\mathbb{R}\times\mathbb{U}\times\mathbb{R}^{m}_{+}$, where $\mathbb{U}\coloneqq\left\{{\bm{u}}\in\mathbb{R}^{m}\mid|{\bm{u}}|\leq{\bm{D}}_{G}\right\}$. 4. (d)

We then apply Theorem 2.2 with $P_{2}(\alpha,{\bm{\beta}})$ and $P_{3}(\alpha,{\bm{\beta}})$ from step (b) on ${{{\cal P}}^{\rm E}}^{\prime}$ to derive the bounds in (78) and (79).

$\bullet$ Step (a) – Revising Theorem 2.2 with Assumption 3.1: Recall that in the derivation of (36) in the proof of Theorem 2.2, Assumption 2.1 yields

[TABLE]

Assumption 3.1 allows us to bound the same by $4C_{F}^{2}$, yielding $P_{2}=16C_{F}^{2}+2\|{\bm{D}}_{G}\|^{2}$. Along the same lines, we get $P_{3}=16m\|{\bm{C}}_{G}\|^{2}$.

$\bullet$ Step (b) – Deriving properties of $\psi_{\omega}$: Consider the stochastic subgradient of $\psi^{f}_{\omega}({\bm{x}},t;\alpha)$ given by

[TABLE]

where $\operatorname{\mathbb{I}}_{\{\cdot\}}$ is the indicator function. Recall that $\|\nabla f_{\omega}({\bm{x}})\|\leq C_{F}$ for all ${\bm{x}}\in\mathbb{X}$ almost surely. Therefore, we have

[TABLE]

Proceeding similarly, we obtain

[TABLE]

We also have

[TABLE]

Then, (80) follows from step (a) using (84), (85), and (86).

$\bullet$ Step (c) – Showing that ${\cal L}^{\prime}$ admits a saddle point: According to [1, Theorem 10], the minimizers of $\mathds{E}[\psi^{f}_{\omega}({\bm{x}},u^{0};\alpha)]$ over $u^{0}$ define a nonempty closed bounded interval (possibly a singleton). Thus, we have

[TABLE]

for some $u^{0}({\bm{x}})\in\mathbb{R}$ for each ${\bm{x}}\in\mathbb{X}$. Similarly, we infer

[TABLE]

for some $u^{i}({\bm{x}})\in\mathbb{R}$ for each ${\bm{x}}\in\mathbb{X}$. Moreover, for all $u^{i}>D_{G}^{i}$, we have

[TABLE]

and for $u^{i}<-D_{G}^{i}$, we have

[TABLE]

Thus $\mathds{E}[\psi^{g^{i}}_{\omega}({\bm{x}},u^{i};\beta^{i})]$ is non-increasing in $u^{i}$ below $-D_{G}^{i}$ and increasing in it beyond $D_{G}^{i}$. Hence, at least one among the minimizers of $\mathds{E}[\psi^{g^{i}}_{\omega}({\bm{x}},u^{i};\beta^{i})]$ must lie in $[-D_{G}^{i},D_{G}^{i}]$. In the sequel, let $u^{i}({\bm{x}})$ refer to such a minimizer.

Consider a saddle point $({\bm{x}}_{\star},{\bm{z}}_{\star})\in\mathbb{X}\times\mathbb{R}^{m}_{+}$ of ${{{\cal P}}^{\rm CVaR}}$. We argue that $({\bm{x}}_{\star},u^{0}({\bm{x}}_{\star}),{\bm{u}}({\bm{x}}_{\star}),{\bm{z}}_{\star})$ is a saddle point of ${\cal L}^{\prime}$. From the definitions of ${\cal L}$, ${\cal L}^{\prime}$, (87), (88), and the saddle point property of $({\bm{x}}_{\star},{\bm{z}}_{\star})$, we obtain

[TABLE]

for all $({\bm{x}},u^{0},{\bm{u}})\in\mathbb{X}\times\mathbb{R}\times\mathbb{U}$. Also, for all ${\bm{z}}\in\mathbb{R}^{m}_{+}$, we have

[TABLE]

$\bullet$ Step (d) – Proof of (78) and (79): By the saddle point theorem and (91), we have ${\cal L}({\bm{x}}_{\star},{\bm{z}}_{\star})=p_{\star}^{\rm CVaR}$, that also equals the optimal value of ${{{\cal P}}^{\rm E}}^{\prime}$. Applying Theorem 2.2 with revised $P_{2}$ and $P_{3}$ from step (b) to ${{{\cal P}}^{\rm E}}^{\prime}$ for which ${\bm{x}}_{0},\dots,{\bm{x}}_{K+1}$ and $u^{0}_{0},\dots,u^{0}_{K+1}$ are $\mathpzc{W}_{K+1/2}$-measurable, we obtain

[TABLE]

Following a similar argument for $i=1,\ldots,m$, we get

[TABLE]

completing the proof. ∎

Our proof architecture generalizes to problems with other risk measures as long as that measure preserves convexity of $f_{\omega},{\bm{g}}_{\omega}$, admits a variational characterization as in (2), and a subgradient for this modified objective can be easily computed and remains bounded over $\mathbb{X}$. We restrict our attention to ${\rm CVaR}$ to keep the exposition concrete.

Opposed to sample average approximation (SAA) algorithms, we neither compute nor estimate $F({\bm{x}})={\rm CVaR}[f_{\omega}({\bm{x}})]$, ${\bm{G}}({\bm{x}})={\rm CVaR}[{\bm{g}}_{\omega}({\bm{x}})]$ for any given decision ${\bm{x}}$ to run the algorithm. Yet, our analysis provides guarantees on the same at $\bar{{\bm{x}}}_{K+1}$ in expectation. If one needs to compute $F$ at any decision variable, e.g., at $\bar{{\bm{x}}}_{K+1}$, one can employ the variational characterization in (2). Such evaluation requires additional computational effort. Notice that Theorem 3.1 does not relate $F(\bar{{\bm{x}}}_{K+1})$ to $p_{\star}^{\rm CVaR}$ in an almost sure sense; it only relates the two in expectation according to (78), where the expectation is evaluated with respect to the stochastic sample path.

${\rm CVaR}$ of a random variable depends on the tail of its distribution. The higher the risk aversion, the further into the tail one needs to look, generally requiring more samples. Even if we do not explicitly compute the tail-dependent CVaR relevant to the objective or the constraints, it is natural to expect our sample complexity to grow with risk aversion, which the following result confirms.

Proposition 3.2**.**

Suppose Assumption 3.1 holds. For an $\varepsilon$-approximately feasible and optimal solution of ${{{\cal P}}^{\rm CVaR}}$ with risk aversion parameters $\alpha,{\bm{\beta}}$ using Algorithm 1 on ${{{\cal P}}^{\rm E}}^{\prime}$, then $\gamma_{\star}(\alpha,{\bm{\beta}})$ and $K_{\star}(\alpha,{\bm{\beta}})$ from Proposition 2.5, respectively decreases and increases with both $\alpha$ and ${\bm{\beta}}$.

Proof.

We borrow the notation from Proposition 2.5 and tackle the variation with $\alpha$ and ${\bm{\beta}}$ separately.

$\bullet$ Variation with $\alpha$: $P_{2}$ increases with $\alpha$, implying $\gamma_{\star}$ decreases with $\alpha$ because $\frac{d\gamma_{\star}^{2}}{dy}\leq 0$ and $\frac{dy}{dP_{2}}\geq 0$. Furthermore, using $\frac{\partial K_{\star}}{\partial\gamma_{\star}}<0$ for $\gamma<\gamma_{\star}$ and $\frac{\partial K_{\star}}{\partial P_{2}}\geq 0$ in

[TABLE]

we infer that $K_{\star}$ increases with $\alpha$.

$\bullet$ Variation with $\beta^{i}$: Both $P_{2}$ and $P_{3}$ increase with $\beta^{i}$ and

[TABLE]

Following an argument similar to that for the variation with $\alpha$, the first term on the RHS of the above equation can be shown to be nonpositive. Next, we show that the second term is nonpositive to conclude that $\gamma_{\star}$ decreases with $\beta^{i}$, where we use $\frac{dP_{3}}{d{\beta^{i}}}\geq 0$. Utilizing $\frac{P_{2}}{P_{1}P_{3}}=y-1$, we infer

[TABLE]

To characterize the variation of $K_{\star}$, notice that

[TABLE]

Again, the first term on the RHS of the above relation is nonnegative, owing to an argument similar to that used for the variation of $K_{\star}$ with $\alpha$. We show $\frac{\partial K_{\star}}{\partial P_{3}}\leq 0$ to conclude the proof. Treating $K_{\star}$ as a function of $P_{3}$ and $\gamma_{\star}$, we obtain

[TABLE]

It is straightforward to verify that the first summand is nonnegative. We have already argued that $\gamma_{\star}$ decreases with $P_{3}$, and $\frac{\partial K_{\star}}{\partial\gamma}<0$ for $\gamma<\gamma_{\star}$, implying that the second summand is nonnegative as well, completing the proof. $\square$ ∎

It is easy to compute the optimized iteration count $K_{\star}(\alpha,{\bm{\beta}})$ and the optimized constant step-size $\gamma_{\star}(\alpha,{\bm{\beta}})/\sqrt{K_{\star}(\alpha,{\bm{\beta}})}$ from Proposition 2.5. The formula is omitted for brevity. Instead, we derive additional insight by fixing ${\bm{\beta}}$ and driving $\alpha$ towards unity. For such an $\alpha,{\bm{\beta}}$, we have

[TABLE]

With $\alpha$ approaching unity, notice that ${{{\cal P}}^{\rm CVaR}}$ approaches a robust optimization problem. Thus, Algorithm 1 for ${{{\cal P}}^{\rm E}}^{\prime}$ is aiming to solve a robust optimization problem via sampling. Not surprisingly, the sample complexity exhibits unbounded growth with such robustness requirements, since we do not assume $\Omega$ to be finite. Also, this growth matches that of solving the SAA problem within $\varepsilon$-tolerance on the unconstrained problem to minimize $\widehat{F}({\bm{x}}):=\frac{1}{K}\sum_{j=1}^{K}\psi_{\omega^{j}}^{f}({\bm{x}},u;\alpha)$. To see this, apply Theorem 2.2 on $\widehat{F}({\bm{x}})$ with optimized step size from Proposition 2.5, where $P_{2}\sim\|\nabla\widehat{F}({\bm{x}})\|^{2}\sim(1-\alpha)^{-2}$ and $P_{3}=1$.

Parallelization can lead to stronger bounds. More precisely, run stochastic approximation in parallel on $N$ machines, each with $K$ samples and compute $\langle\bar{{\bm{x}}}\rangle_{K+1}:=\frac{1}{N}\sum_{j=1}^{N}\bar{{\bm{x}}}_{K+1}[j]$ using $\bar{{\bm{x}}}_{K+1}[1],\ldots,\bar{{\bm{x}}}_{K+1}[N]$ obtained from the $N$ separate runs. Then, we have

[TABLE]

for $i=1,\ldots,m$ and $\tau>0$. The steps combine coherence of ${\rm CVaR}$, convexity and uniform boundedness of $g^{i}_{\omega}$, Hoeffding’s inequality and Theorem 3.1. A similar bound can be derived for suboptimality. Thus, parallelized stochastic approximation produces a result whose ${\cal O}(1/\sqrt{K})$-violation occurs with a probability that decays exponentially with the degree of parallelization $N$.

The bound in (101) reveals an interesting connection with results for chance constrained programs. To describe the link, notice that ${\rm CVaR}_{\delta}[y_{\omega}]\leq 0$ implies ${\sf Pr}\{y_{\omega}\leq 0\}\geq 1-\delta$ for any random variable $y_{\omega}$ and $\delta\in[0,1)$. Therefore, (101) implies

[TABLE]

for constants $C,C^{\prime}$. Said differently, our stochastic approximation algorithm requires ${\cal O}(\log(1/{\nu}))$ samples to produce a solution that satisfies an ${\cal O}(1/\sqrt{\log(1/\sqrt{\nu})})$-approximate chance-constraint with a violation probability bounded by $\nu$. This result bears a striking similarity to that derived in [33], where the authors deterministically enforce ${\cal O}(\log(1/{\nu}))$ sampled constraints to produce a solution that satisfies the exact chance-constraint ${\sf Pr}\left\{g^{i}_{\omega}\left({\bm{x}}\right)\leq 0\right\}\geq 1-\beta^{i}$ with a violation probability bounded by $\nu$. This resemblance in order-wise sample complexity is intriguing, given the significant differences between the algorithms.

3.1 An illustrative example

We explore the use of our algorithm on the following example problem

[TABLE]

Let $\omega\sim\frac{1}{3}{\sf beta}(2,2)$ and consider the specific choice of risk parameters $\alpha=0.3,\beta=0.2$. To gain intuition into the optimal solution for this example, we numerically estimate $F(x)$ and $G^{1}(x)$ for each $x$ and plot them in Figure 1a. To that end, we first obtain a million samples of $\omega$. Then, for each value of the decision variable $x$, we sort the objective function value $f_{\omega}(x)$ and the constraint function value $g_{\omega}(x)$ with these samples. We then estimate $F$ and $G^{1}$ as the average of the highest $1-\alpha=70\%$ and $1-\beta=80\%$ among $f_{\omega}(x)$’s and $g_{\omega}(x)$’s, respectively, at each $x$ with those samples. The unique optimum for (103) is numerically evaluated as $x_{\star}\approx-0.1929$ for which $F(x_{\star})\approx 0.4042$ and $G^{1}(x_{\star})\approx 0$.

For this example, it is easy to show that $C_{F}=\frac{4}{3}$, $C_{G}=1$ and $D_{G}=\frac{5}{6}$ that yields $P_{2}(0.3,0.2)=\frac{8276}{93}$ and $P_{3}(0.3,0.2)=50$. To run Algorithm 1 on ${{{\cal P}}^{\rm E}}^{\prime}$ derived from (103), we can use constant step-size $\gamma_{k}=\gamma/\sqrt{K}$ with a pre-determined number of steps $K$ for any $0<\gamma<P_{3}^{-1/2}(0.3,0.2)=\frac{1}{5\sqrt{2}}$. With any given $K$, Theorem 3.1 guarantees that the expected distance to $F(x_{\star})$ and the expected constraint violation evaluated at $\bar{{x}}_{K+1}$ decays as $1/\sqrt{K}$. For a given $K$ and $\gamma<\frac{1}{5\sqrt{2}}$, calculating the precise bound $\eta(0.3,0.2)/\sqrt{K}$ requires the knowledge of $P_{1}$ or its overestimate. For this example, $|x_{\star}|\leq\frac{1}{2}$ and $|u^{1}_{\star}|\leq D_{G}=\frac{5}{6}$. Also, $|u^{0}_{\star}|$ is bounded above by the maximum value that $|f_{\omega}(x)|$ can take, that is given by $\frac{8}{9}$. Since we cannot determine $z_{\star}$ a priori, we assume $|z_{\star}|\leq 2$ (that will later be shown to be consistent with our result). Starting from $(x_{0},u^{0}_{0},u^{1}_{0},z_{0})=0$, we then obtain $P_{1}=\frac{3197}{81}$. To solve ${{{\cal P}}^{\rm CVaR}}$ (or equivalently ${{{\cal P}}^{\rm E}}^{\prime}$) with a tolerance of $\varepsilon=5\times 10^{-3}$, we require $\eta(0.3,0.2)/\sqrt{K}\leq 5\times 10^{-3}$. With this tolerance and the values of $P_{1},P_{2},P_{3}$, Proposition 2.5 yields an optimized $\gamma_{\star}=0.0808$ and $K_{\star}\approx 1.35\times 10^{9}$. We run Algorithm 1 on ${{{\cal P}}^{\rm E}}^{\prime}$ with constant step-size $\gamma_{\star}/\sqrt{K_{\star}}$ and plot $F$ and $G^{1}$ at the running ergodic mean of the iterates, i.e., at $\bar{x}_{k}:=\frac{1}{k}\sum_{j=1}^{k}x_{j}$ for each $k$. Again $F$ and $G^{1}$ are evaluated numerically using the ${\rm CVaR}$-estimation procedure we outlined above.

Notice that Theorem 3.1 only guarantees a bound on $F(\bar{x}_{K_{\star}+1})-F(\bar{x}_{\star})$ and $G^{1}(\bar{x}_{K_{\star}+1})$ in expectation. Thus, one would expect that only the average of the ${\rm CVaR}$ of $F$ and $G^{1}$ evaluated at $\bar{x}_{K_{\star}+1}$ over multiple sample paths to respect the $\varepsilon$-bound. However, our simulation yielded $\bar{x}_{K_{\star}+1}=-0.1926$ and $\bar{z}_{K_{\star}+1}=0.8976$, for which

[TABLE]

i.e., the ergodic mean after $K_{\star}$ iterations respects the $\varepsilon$-bound over the plotted sample path. The same behavior was observed over multiple sample paths. The ergodic mean of the dual iterate is indeed consistent with our assumption $|z_{\star}|\leq 2$ made in deriving $\eta(0.3,0.2)$. We point out that the ergodic mean in Figure 1b moves much more smoothly than our evaluation of $F$ and $G^{1}$ at those means, especially for large $k$. The noise in $F$ in $G^{1}$ emanate from the finitely many samples we use to evaluate $F$ and $G^{1}$. The errors appear much more pronounced at larger $k$, given the logarithmic scale of the plot.

The optimized iteration count $K_{\star}(\alpha,\beta)$ from Proposition 2.5 with a modest $\alpha=0.3,\beta=0.2$ is quite high even for this simple example. This iteration count only grows with increased risk aversion as Figure 2 reveals. Figure 1b suggests that the $\varepsilon=5\times 10^{-3}$ tolerance is met earlier than $K_{\star}$ iterations. This is the downside of optimizing upper bounds to decide step-sizes for subgradient methods. Carefully designed termination criteria may prove useful in practical implementations.

We end the numerical example with a remark about the comparison of Algorithm 1 that uses Gauss-Seidel-type dual update in (5) and another that uses the popular Jacobi-type dual update on ${{{\cal P}}^{\rm E}}^{\prime}$ for (103) with $\alpha=0.3,\beta=0.2$. This alternate dual update replaces ${\bm{g}}_{\omega_{k+1/2}}({\bm{x}}_{k+1})$ in (5) by ${\bm{g}}_{\omega_{k}}({\bm{x}}_{k})$. That is, the same sample $\omega_{k}$ is used for both the primal and the dual update. And, the primal iterate ${\bm{x}}_{k}$ is used instead of ${\bm{x}}_{k+1}$ to update the dual variable. We numerically compared this primal-dual algorithm with Algorithm 1 with various choices of step-sizes (consistent with the requirements of Theorem 3.1) and iteration count for our example and its variations. For each run, we found that the iterates from both these algorithms moved very similarly. The differences are too small to report. The Jacobi-type update requires half the number of samples compared to Algorithm 1. While the extra sample helps us in the theoretical analysis, our experience with this stylized example does not suggest any empirical advantage. A more thorough comparison between these algorithms, both theoretically and empirically, is left to future work.

4 Conclusions and future work

In this paper, we study a stochastic approximation algorithm for ${\rm CVaR}$-sensitive optimization problems. Such problems are remarkably rich in their modeling power and encompass a plethora of stochastic programming problems with broad applications. We study a primal-dual algorithm to solve that problem that processes samples in an online fashion, i.e., obtains samples and updates decision variables in each iteration. Such algorithms are useful when sampling is easy and intermediate approximate solutions, albeit inexact, are useful. The convergence analysis allows us to optimize the number of iterations required to reach a solution within a prescribed tolerance on expected suboptimality and constraint violation. The sample and iteration complexity predictably grows with risk-aversion. Our work affirms that a modeler must not only consider the attitude towards risk but also consider the computational burdens of risk in deciding the problem formulation.

Two possible extensions are of immediate interest. First, primal-dual algorithms find applications in multi-agent distributed optimization problems over a possibly time-varying communication network. We plan to extend our results to solve distributed risk-sensitive convex optimization problems over networks, borrowing techniques from [44, 45]. Second, the relationship to sample complexity for chance-constrained programs in [33] encourages us to pursue a possible exploration of stochastic approximation for such optimization problems.

5 Acknowledgements

We thank Eilyan Bitar, Rayadurgam Srikant, Tamer Başar and Stan Uryasev for helpful discussions. This work was partially supported by the International Institute of Carbon-Neutral Energy Research (I2CNER) and the Power System Engineering Research Center (PSERC).