<figure><img src="image_1.png" /></figure>

HAL
open science

# Risk-averse optimization: Linear Programming implementation of CVaR

Pierre Haessig

► To cite this version:

Pierre Haessig. Risk-averse optimization: Linear Programming implementation of CVaR. Centrale-Supélec. 2025. hal-04966655v2

HAL Id: hal-04966655

https://hal.science/hal-04966655v2

Submitted on 28 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Risk-averse optimization: Linear Programming implementation of CVaR

Pierre Haessig, February 2025

**Objective:** This technical note describes the practical implementation of a risk-averse optimization based on CVaR. The CVaR (Conditional Value-at-Risk) of a random cost function can be minimized using a Linear Programming (LP) implementation. It is a simplification of the ideas introduced in the original CVaR optimization article by Rockafellar & Uryasev in 2000.

**License:** this document is freely available under the CC BY 4.0.

## References:

* R. T. Rockafellar and S. Uryasev, "Optimization of conditional value-at-risk," The Journal of Risk, vol. 2, no. 3, pp. 21-41, 2000, doi: 10.21314/jor.2000.038.

* Original article proposing the alternative "optimization-ready" definition of CVaR

* See also slides S. Uryasev, "Conditional Value-at-Risk (CVaR): Algorithms and Applications." Aug. 2001. http://www-iam.mathematik.hu-berlin.de/romisch/SP01/Uryasev.pdf

* S. Boyd and L. Vandenberghe, Convex Optimization. Cambridge University Press, 2004. Freely available at https://web.stanford.edu/~boyd/cvxbook/

## Notations and definitions

Let be $J(x, \theta)$ a random cost function:

* $x$: decision variable (Here and Now, i.e. deterministic)

* $\theta$: random Variable with known distribution

* represented in practice using $n$ samples $\theta_i$ with probability $p_i$ for $i = 1, n$

## Reminder of optimization under uncertainty

In order to be "optimization-ready", the random cost needs to be **mapped to a scalar** (mapping $\mathbb{R} \cdot V \rightarrow \mathbb{R}$). The simplest statistical mapping is the expectation:

$$ \mathbb{E}_{\theta}(J(x, \theta)) = \sum_{i} p_i \cdot J(x, \theta_i) $$

and minimizing the expectation is said *risk neutral* in the sense that values below and above are supposed to compensate in a symmetrical fashion by the law of large numbers.

Here we are instead interested in the *risk-averse* mapping CVaR (Conditional Value at Risk)

## "Base" definition of CVaR (statistical definition)

Definition: CVaR$_{\alpha}$ with risk level $\alpha \in [0, 1]$ (but likely close to one, e.g. 90%) is:

$$ \text{CVaR}_{\alpha}(J(x, \theta)) = \mathbb{E}_{\theta}[J(x, \theta)|J(x, \theta) \geq \text{VaR}_{\alpha}(J(x, \theta))] $$
that is, CVaR<sub>α</sub> is the conditional expectation of values above the VaR<sub>α</sub> with same risk level α, where VaR is the Value at Risk, which a synonymous to the α-quantile, i.e. the value such that $P_{\theta}[J(x, \theta) \le \text{VaR}_{\alpha}(J(x, \theta))] = \alpha$. Here is a graphical illustration using the probability density function of a random variable X (which plays the role of the random cost J here, sorry for the notational mismatch):

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Remark 1: See (Rockafellar & Uryasev, 2000) and Uryasev's slides for some subtleties on the VaR definition, when the cumulative distribution is *not continuous* (in the presence of discrete probability masses).

Remark 2: In the equations above and the subsequent ones, we always write explicitly the dependency in the decision variable $x$ to recall that the final goal is to optimization a scalar function of $x$ (Expectation: $x \mapsto E_{\theta}(J(x, \theta))$ or CVaR: $x \mapsto \text{CVaR}_{\alpha}(J(x, \theta))$). However, in all these changes the $x$-dependency is transparently propagated through the mappings, so that you are free to regard the random cost function $J(x, \theta)$ as a simpler random variable $J(\theta)$, computed for "some value" of $x$.

# Linear Programming implementation of CVaR

## Alternative CVaR definition (optimization-based definition)

The basic definition of CVaR may look easy to implement in an optimization setting for the expectation part, but is difficult because the condition depends on the VaR which is itself to be computed. However, this can be bypassed using an alternative definition of CVaR which, oddly enough is not explicit but instead the outcome of the following *partial* optimization problem:

$$ \text{CVaR}_{\alpha}(J(x, \theta)) = \min_{v} \Phi_{\alpha}^{J}(x, v) $$

with minimization variable $v \in \mathbb{R}$ and with the function $\Phi_{\alpha}^{J}$ to be minimized:

$$ \Phi_{\alpha}^{J}(x, v) = v + \frac{1}{1 - \alpha} E_{\theta}|J(x, \theta) - v|_{+} $$

where $|\cdot|_{+} = \max(\cdot, 0)$ is shorthand notation of the *positive part*. Using the $n$ uncertainty samples $\theta_i$, this function is implemented by:

$$ \Phi_{\alpha}^{J}(x, v) = v + \frac{1}{1 - \alpha} \sum_{i} p_{i} \cdot |J(x, \theta_{i}) - v|_{+} $$

Interpretation of the argmin $v^*$
Property: the argmin $v^*$ of the partial optimization of $v \mapsto \Phi_{\alpha}^{J}(x, v)$ is in fact $\text{VaR}_{\alpha}(J(x, \theta))$:

$$ \text{VaR}_{\alpha}(J(x, \theta)) = v^* $$

Remark: See again (Rockafellar & Uryasev, 2000) for a more rigorous analysis (when the argmin is not unique).

## Why it works

While this fact is not self-evident (at least for me), once known it is then quite easy to see that $\Phi_{\alpha}^{J}(x, v = \text{VaR}_{\alpha})$ falls back to the CVaR (hereafter we shorten the notation $\text{VaR}_{\alpha} = \text{VaR}_{\alpha}(J(x, \theta))$):

$$ \Phi_{\alpha}^{J}(x, \text{VaR}_{\alpha}) = \text{VaR}_{\alpha} + \frac{1}{1-\alpha} \mathbb{E}_{\theta} |J(x, \theta) - \text{VaR}_{\alpha}|_{+} $$

Indeed, for computing $\mathbb{E}$ in the second term, there are only two possibilities for the value of the positive part $|\cdot|_{+}$, depending on the sign of the difference $J(x, \theta) - \text{VaR}_{\alpha}$. By definition of $\text{VaR}_{\alpha}$:

* with probability $\alpha$, we have $J(x, \theta) < \text{VaR}_{\alpha}$ and thus the positive part is zero

* with probability $1 - \alpha$, we have the converse and thus the positive part is the value $J(x, \theta) - \text{VaR}_{\alpha}$ itself

Thus by the law of total expectation:

$$ \mathbb{E}_{\theta} [J(x, \theta) - \text{VaR}_{\alpha}]_{+} = \alpha \times \mathbb{E}_{\theta} [0 | J(x, \theta) < \text{VaR}_{\alpha}] + (1 - \alpha) \times \mathbb{E}_{\theta} [J(x, \theta) - \text{VaR}_{\alpha} | J(x, \theta) \geq \mathbb{V}] $$

and the quantity $\mathbb{E}_{\theta} [J(x, \theta) - \text{VaR}_{\alpha} | J(x, \theta) \geq \text{VaR}_{\alpha}]$ is the initial definition of the CVaR, but translated by $-\text{VaR}_{\alpha}$. Thus:

$$ \mathbb{E}_{\theta} [J(x, \theta) - \text{VaR}_{\alpha}]_{+} = (1 - \alpha)(\text{CVaR}_{\alpha} - \text{VaR}_{\alpha}) $$

In the end

$$ \Phi_{\alpha}^{J}(x, \text{VaR}_{\alpha}) = \text{VaR}_{\alpha} + \frac{1-\alpha}{1-\alpha}(\text{CVaR}_{\alpha} - \text{VaR}_{\alpha}) = \text{CVaR}_{\alpha} $$

For the complete proof that Var and CVaR are respectively the argmin and min of $v \mapsto \Phi_{\alpha}^{J}(x, v)$, see (Rockafellar & Uryasev, 2000) where they first prove the convexity of this function and then compute the derivative with respect to $v$ and solve for $(\partial/\partial v)\Phi_{\alpha}^{J}(x, \text{VaR}_{\alpha}) = 0$.

## CVaR minimization

We called the alternative CVaR definition a *partial optimization problem* in the sense that the result of minimizing $\Phi_{\alpha}^{J}(x, v)$ over $v$ is still a function of decision $x$. This *implicit function* $\text{CVaR}(x)$ is left to be further optimized. However, the good news is that the full optimization problem $\min_{x} \text{CVaR}_{\alpha}(J(x, \theta))$ can be cast as the natural *joint optimization* of $\Phi_{\alpha}^{J}$:

$$ \min_{x} \text{CVaR}_{\alpha}(J(x, \theta)) = \min_{x,v} \Phi_{\alpha}^{J}(x, v) $$

Proof: See Theorem 2 of (Rockafellar & Uryasev, 2000)

## Linear Programming implementation
The alternative definition of CVaR (as the minimum of $\Phi_{\alpha}^{J}(x, v)$) can be implemented as a Linear Program (LP) thanks to the "epigraph reformulation" which is a "classical" LP modeling trick that can transform any convex piecewise linear (PWL) expression into a set of inequality (see section 4.1.3 Equivalent problems in the book (Boyd & Vandenberghe, 2004)). Notice that it only works when PWL expression is meant to be minimized.

<figure><img src="image_5.png" /></figure>

Examples:

**Absolute value minization:**

The nonlinear (but convex piecewise linear) problem:

$$ \min_{x} |x| $$

is reformulated as:

$$ \min_{t,x} t $$

with the two inequality constraints $t \ge x$ and $t \ge -x$. Then, *at the optimum*, because $t$ is being minimized, one of the two inequalities will be tight, i.e. $t^* = +x^*$ or $t^* = -x^*$ depending on which one is larger, i.e. $t^* = |x^*|$

**Positive part minimization:**

The nonlinear (but convex piecewise linear) problem:

$$ \min_{x} |x|_{+} $$

is reformulated as:

$$ \min_{t,x} t $$

with the two inequality constraints $t \ge x$ and $t \ge 0$. Then, *at the optimum*, because $t$ is being minimized, one of the two inequalities will be tight, i.e. $t^* = x^*$ or $t = 0$ depending on which one is larger, i.e. $t^* = \max(x^*, 0) = |x^*|_+$

And this last template can be applied to reformulate $\Phi_{\alpha}^{J}(x, v)$ (which is to be minimized to get the CVaR):

$$ \Phi_{\alpha}^{J}(x, v) = v + \frac{1}{1-\alpha} \sum_{i} p_{i} \cdot |J(x, \theta_{i}) - v|_{+} $$

but since there are $i = 1...n$ instances of the positive part $|\cdot|_+$ = $\max(0, 0)$, we need to create $n$
extra variables $t_i$

$$ \text{CVaR}_{\alpha}(J(x, \theta)) = \min_{v} v + \frac{1}{1-\alpha} \sum_{i} p_i \cdot t_i $$

with constraints $t_i \ge 0$ and:

$$ t_i \ge J(x, \theta_i) - v \text{ for } i = 1 \dots n $$

so at the end of the road, it's only hidden in this last set of constraint that we encode which random quantity (here $J(x, \theta)$) we would like to compute the CVaR on...

## Notation mapping to (Rockafellar & Uryasev, 2000)

This document follows the reasoning of (Rockafellar & Uryasev, 2000), albeit with different notations to be consistent with the other lectures:

* decision variable(s): $x$ for both!

* random parameter(s): $\theta$ here, $y$ for R&U

* random cost function: $J(x, \theta)$ here, $f(x, y)$ for R&U (like manu optimization books, but I keep $f$ for the PDFs)

* risk level to define the VaR and CVaR: $\alpha$ here, $\beta$ for R&U

* VaR: $VaR_{\alpha}$ here, $\alpha_{\beta}$ for R&U (I find this notation very confusing...)

* and I use $v$ when it is a free variable which is VaR at the optimum

* CVaR: $CVaR_{\alpha}$ here, $\phi_{\beta}$ for R&U

* Function which minimum is the CVaR: $\Phi_{\alpha}^{J}(x, v)$ here, $F_{\beta}(x, \alpha)$ for R&U

## Numerical implementation in Matlab/YALMIP

In Matlab/Octave using YALMIP, but any Linear Programming modeling language could do it (JuMP in Julia, Pyomo in Python...):

### Simplest test case $J = \theta$

i.e. compute $CVaR_{\alpha}(\theta)$, with no decision variable $x$

In [1]: % Add YALMIP folder and subfolders to the MATLAB/Octave path
addpath(genpath('/home/pierre/Programmes/others/YALMIP-2021-03-31'))

Try to adjust graphics toolkit (see https://github.com/Calysto/octave_kernel and https://octave.discourse.group/t/error-graphics-toolkit-qt-toolkit-is-not-available/2146/6)

### Uncertainty samples

Uncertainty samples for Random Variable $\theta$: Uniform law for which we know

* VaR = 900, for risk level $\alpha$ = 0.9

* CVaR = 950

In [2]: n = 51 % number of samples
```matlab
% bounds of Uniform law:
theta_min = 0;
theta_max = 1000;

theta Samples = linspace theta_min, theta_max, n);
theta proba = ones(1,n)/n;
disp("Samples subset (start...end):")
disp([theta samples(1:3) theta samples(end-2:end)])
```
%stem theta samples, theta proba); % broken inline Octave plot

```
n = 51
Samples subset (start...end):
0      20     40     960    980   1000
```

## Auxiliary variables

auxiliary variable $v$ (to be the VaR at the optimum)

In [3]: $v = sdpvar(1)$

Linear scalar (real, 1 variable)
Coefficients range: 1 to 1

auxiliary variables $t_i$ for $i = 1...n$ to implement the epigraph form of the positive part:

In [4]: $tvec = sdpvar(1,n)$

Linear matrix variable 1x51 (full, real, 51 variables)
Coefficient range: 1 to 1

## Objective function & constraints

Desired CVaR probability $\alpha \in [0, 1]$:

In [5]: $alpha = 0.90$

$alpha = 0.9000$

Then the objective function is

$$ \Phi_{\alpha}^{J}(x, v) = v + \frac{1}{1 - \alpha} \sum_{i} p_{i} \cdot |J(x, \theta_{i}) - v|_{+} $$

which is relaxed as $v + \frac{1}{1-\alpha} \sum_{i} p_{i} \cdot t_{i}$

In [6]: $Phi\_alpha = v + 1/(1-\alpha) \cdot \sum(\theta\_proba \cdot tvec)$

Linear scalar (real, 52 variables)
Coefficients range: 0.19608 to 1

Constraints: It's only at the level of the $t \ge \cdot$ constraints that we formulate the quantity for which we want to calculate the CVaR.

As a simplistic case, we use $J(x, \theta) = \theta$ (i.e. no decision variable)

In [7]: $Jvec = theta\_samples$

In [8]: $Cons = [tvec>=0]$; $Cons = [Cons, tvec >= Jvec-v]$;
```matlab
% can also be written as
% for i = 1:n
%     Cons = [Cons, tvec(i) >= Jvec(i)-v];
% end
Cons
```

<table><thead><tr><th>ID</th><th>Constraint</th><th>Coefficient range</th></tr></thead><tbody><tr><td>#1</td><td>Element-wise inequality 51x1</td><td>1 to 1</td></tr><tr><td>#2</td><td>Element-wise inequality 51x1</td><td>1 to 1000</td></tr></tbody></table>

In [9]: optimize(Cons, Phi alpha)

ans =

scalar structure containing the fields:

    yalmipversion = 20210331
    matlabversion = 6.4.0
    yalmiptime = 1.3265
    solvetime = 4.7750e-03
    info = Successfully solved (GLPK-GLPKMEX)
    problem = 0

At optimum, v should be the VaRα (i.e. quantile)

In [10]: value(v) % 900 expected

ans = 900

while, $\Phi_{\alpha}^{J}(v^*)$ should be the CVaR$_\alpha$ (up to the discretisation of the uniform law)

In [12]: value(Phi alpha) % 950 if infinitely many samples

ans = 958.82

at the optimum, $t_i^*$ is indeed equal to $|\theta_i - \text{VaR}(\theta)|_+$:

In [13]: [theta_sample(\*end-7:\*end); value(tvec(\*end-7:\*end))]
ans =

ans =

860    880    900    920    940    960    980   1000
    0     0     0     20     40     60     80    100

## More interesting variant with decision variables

We consider two products we can buy in quantities $x_k \ge 0$ at a random unit prices $c_k$, for $k = 1, 2$.
Total cost is:

$$J = c_1 \cdot x_1 + c_2 \cdot x_2$$

and are total amount of allocation is fixed: $x_1 + x_2 = 1$ so that the optimization problem $\min_J J$ is
to select the least expensive mix.

The prices random prices are fully correlated so that only one random parameter is needed.
* $\theta \sim U(0, 1000)$ like above

* $c_1 = \theta$, i.e 500 € in expectation, 1000 € in the worst case

* $c_2 = 400 + \theta * 2/5$, i.e 600 € in expectation, 800 € in the worst case

* VaR90 = 400+0.9*400 = 760 €

* CVaR90 = 760+20 = 780 €

Expected optimal allocation:

* if risk level $\alpha$ is close to 0, then CVaR falls back to the regular expectation and we should be a full allocation on product 1 (least expected price): $x_1 = 1, x_2 = 0$

* but if risk level $\alpha$ is close to 1, we should get a full allocation on product 2 (least worst price)
$x_1 = 0, x_2 = 1$

* TODO: compute analytically the tipping point $\alpha_{threshold}$ at we switch the allocation

In [16]:
c1 Samples = theta Samples;
c2 Samples = theta Samples * 2/5 + 400;
disp("Min/Mean/Max of c1 and c2:")
disp([min(c1 Samples), mean(c1 Samples), max(c1 Samples)
      min(c2 Samples), mean(c2 Samples), max(c2 Samples)])

Min/Mean/Max of c1 and c2:
   0    500   1000
  400    600    800

In [17]:
x = sdpvar(1,2)

Linear matrix variable 1x2 (full, real, 2 variables)
Coefficient range: 1 to 1

In [18]:
Jvec = c1 Samples * x(1) + c2 Samples * x(2)

Linear matrix variable 1x51 (full, real, 2 variables)
Coefficient range: 20 to 1000

auxiliary variables for CVaR

In [19]:
v = sdpvar(1)
tvec = sdpvar(1,n)

Linear scalar (real, 1 variable)
Coefficients range: 1 to 1
Linear matrix variable 1x51 (full, real, 51 variables)
Coefficient range: 1 to 1

Desired CVaR probability $\alpha \in [0, 1]$:

In [20]:
alpha = 0.90
alpha = 0.9000

Objective function

$$\Phi_{\alpha}^{J}(x, v) = v + \frac{1}{1 - \alpha} \sum_{i} p_{i} \cdot |J(x, \theta_{i}) - v|_{+}$$

which is relaxed as $v + \frac{1}{1-\alpha} \sum_{i} p_{i} \cdot t_{i}$

In [21]:
Phi alpha = v + 1/(1-alpha) * sum theta proba * tvec
Linear scalar (real, 52 variables)
Coefficients range: 0.19608 to 1

In [22]:
Cons = [x>=0];
Cons = [Cons, sum(x)==1];
Cons = [Cons, tvec>=0, tvec >= Jvec-v];
Cons

<table><thead><tr><th>ID</th><th>Constraint</th><th>Coefficient range</th></tr></thead><tbody><tr><td>#1</td><td>Element-wise inequality 2x1</td><td>1 to 1</td></tr><tr><td>#2</td><td>Equality constraint 1x1</td><td>1 to 1</td></tr><tr><td>#3</td><td>Element-wise inequality 51x1</td><td>1 to 1</td></tr><tr><td>#4</td><td>Element-wise inequality 51x1</td><td>1 to 1000</td></tr></tbody></table>

In [23]: optimize(Cons, Phi alpha)

ans =

scalar structure containing the fields:

yalmipversion = 20210331
matlabversion = 6.4.0
yalmiptime = 1.3426
solvertime = 8.5998e-04
info = Successfully solved (GLPK-GLPKMEX)
problem = 0

Optimal allocation:

In [24]: value(x)

ans =
0 1.0000

Optimal v (should be VaR)

In [25]: value(v)

ans = 760

Optimal CVaR (theoretically 780 for α = 0.9)

In [26]: value(Phi alpha)

ans = 783.53

at the optimum, $t_i^*$ is indeed equal to $|J_i - VaR(\theta)|_+$:

In [27]: [theta_sample (end-6:end); value(tvec(end-6:end))]/10

ans =

88.0000   90.0000   92.0000   94.0000   96.0000   98.0000   100.0000
    0       0       0.8000   1.6000   2.4000   3.2000   4.0000