Computing Monetary Risk Measures in Linear Time


Abstract

Monetary risk measures have gained popularity for expressing decision-makers’ risk aversion. Value-at-Risk (VaR) and Conditional-Value-at-Risk (CVaR), in particular, are used commonly for this purpose. This paper proposes new efficient algorithms to compute these risk measures for a discrete random variable in expected linear time with respect to the size of its domain. First, we propose a QuickVaR algorithm that computes the VaR of a discrete random variable. Then, we leverage QuickVaR to propose QuickDivergence, an algorithm for computing a class of \(\varphi\)-divergence risk measures, including the popular CVaR risk measure. The QuickVaR algorithm adapts the well-known Quickselect algorithm, while QuickDivergence builds on polymatroid optimization algorithms. Numerical results show that our new algorithms offer an order-of-magnitude speedup for large domains, and a library implementation of the algorithms is available at https://github.com/RiskAverseRL/RiskMeasures.jl.

1 Introduction↩︎

Monetary measures of risk generalize the expectation operator to be able to represent the risk aversion of decision-makers in domains that range from robotics [1][3] to finance [4] to infrastructure maintenance [5] to disaster response [6]. Value-at-Risk (VaR) and Conditional-Value-at-Risk (CVaR) are the most popular risk measures and have made inroads in reinforcement learning, machine learning, and decision making [7]. Machine learning and reinforcement learning algorithms typically must compute the risk measures quickly for large discrete random variables. For example, in risk-averse reinforcement learning, one must evaluate the risk measure in every iteration of the algorithm, which can represent a significant computational bottleneck [8][12].

In this paper, we propose linear-time algorithms for computing risk measures of discrete random variables. First, we develop QuickVaR, an algorithm that computes VaR in linear time and generalizes the well-known quickselect algorithm. Second, we develop QuickDivergence, an algorithm that computes, in linear time, \(\varphi\)-divergence risk measures [13] that satisfy a specific piecewise linearity property. This algorithm reduces \(\varphi\)-divergence risk measures to linear minimization problems over polymatroids [4]. We illustrate this algorithm for the well-known CVaR risk measure and a related TVaR risk measure. Both of these risk measures are in the \(\varphi\)-divergence family and satisfy the required linearity property.

Our algorithms’ time complexity improves on the state of the art by a logarithmic factor. Virtually all implementations of VaR and CVaR algorithms require sorting the input arrays. Unbounded sorting requires \(\Omega(n \log n)\) operations where \(n\) is the size of the probability space [4]. Our numerical results show that the logarithmic acceleration can lead to an order-of-magnitude speedup for realistic large problems without incurring any penalty for small problems. In addition, our linear-time algorithms have simple implementations that add minimal complexity over standard algorithms.

Linear-time algorithms similar to our QuickVaR and QuickDivergence have been studied in many domains, and we summarize some of the most relevant. Some examples of such structured optimization problems are the projection onto the simplex [14], [15], projection onto the \(k\)-capped simplex [16], [17], linear optimization on the intersection of \(L_p\) or \(\varphi\)-divergence balls and the probability simplex in robust optimization [10], optimization for \(\varphi\)-divergence risk measures [13], and on polymatroids [18]. While some of these papers discuss the possibility of designing linear-time algorithms, they usually do not study them, with a few exceptions [17], [19]. Finally, QuickVaR is closely related to the computation of the weighted quantile but returns the same value only when the quantile is unique [20].

The remainder of the paper is organized as follows. 2 describes the necessary background including risk measures, \(\varphi\)-divergence, and polymatroid optimization. Then, 3 proposes QuickVaR to compute VaR in linear time. 4 builds on QuickVaR to develop QuickDivergence that computes certain \(\varphi\)-divergence risk measures in linear time. We evaluate the algorithms numerically on synthetic and realistic problems in 5.

2 Preliminaries↩︎

Notation: We restrict our attention to random variables defined on finite probability spaces \((\Omega, 2^{\Omega}, \boldsymbol{p})\) where \(\Omega := \left\{ 1, \dots , n \right\}\), \(\boldsymbol{p} \in \Delta_n\), and \(\Delta_n\) is the \(n\)-dimensional probability simplex. We use a tilde to indicate a random variable, such as \(\tilde{x}\colon \Omega \to \Real\). Vectors are denoted with a lower-case bold font, such as \(\boldsymbol{x} \in \Real^n\) and indexed as \(\boldsymbol{x}_i \text{ for } i = 1, \dots , n\). We use the function notation \(\tilde{x}\) and vector notation \(\boldsymbol{x}\) interchangeably, such that \(x_\omega = x(\omega), \forall \, \omega \in \Omega\). Sets are denoted with calligraphic letters, \(\bar{\Real} := \Real \cup \left\{ -\infty, \infty \right\}\) denotes the extended real line, and \(\Real_{\ge 0}\) are non-negative reals. The symbol \(\mathbb{N}\) denotes all 0-based natural numbers and \(a{:}b := a, a+1, \dots , b\) for \(a \le b \in \mathbb{N}\) and the empty set otherwise. We use \(\boldsymbol{x}_{a{:}b}\) to denote a subset of elements of \(\boldsymbol{x}\) and comparisons such as \(\boldsymbol{x}_{a{:}b} < y\) are interpreted element-wise. Finally, the symbol \(\boldsymbol{1} \in \Real^n\) denotes a vector of all ones and \(\boldsymbol{1}_{\mathcal{O}} \in \Real^n\) for \(\mathcal{O} \subseteq \Omega\) is an indicator vector of the set \(\mathcal{O}\).

Monetary risk measures \(\rho\colon \Delta_n \times \Real^n \to \bar{\Real}\) are functions of a probability measure and a random variable that are monotone and translation invariant [4]. Risk measures generalize the expectation operator to account for uncertainty in a random variable’s distribution. We apply risk measures to random variables representing rewards and suppose that higher risk values are preferable. For the sake of streamlining the notation, we omit the probability measure when it is obvious from the context and only write \(\rho[\tilde{x}]\) to compute the risk value.

Value-at-Risk (VaR) is a popular risk measure and is defined for \(\alpha\in [0,1)\) and \(\tilde{x}\in \mathbb{R}^n\) as [4], [21] \[\label{eq:var-definition} \var{\alpha}{ \tilde{x} } \;:=\; \max\; \left\{\tau \in \bar{\Real} \mid \P{\tilde{x} < \tau} \le \alpha\right\} = \max\; \left\{\tau \in \bar{\Real} \mid \P{\tilde{x} \ge \tau} \ge 1 - \alpha\right\},\tag{1}\] where the maximum exists because \(\tau \mapsto \P{\tilde{x} < \tau}\) is lower semi-continuous. One can readily see that \(\varo_0[\tilde{x}] = \mathop{\mathrm{ess}}\inf [\tilde{x}]\), and typically the risk measure is extended as \(\varo_{1}[ \tilde{x} ] = \infty\).

Table 1: Parameters for relevant \(\varphi\)-divergence risk measures for \(\alpha\in (0,1)\).
Risk measure \(\varphi(y)\) \(\varphi\opt(x)\) \(\beta\)
\(\mathbb{E}\) \(\begin{cases} 0 &\text{if } y = 1 \\ \infty & \text{otherwise} \end{cases}\) \(x\) \(0\)
\(\cvaro_{\alpha}\) \(\begin{cases} 0 &\text{if } y \in[0, \frac{1}{\alpha}] \\ \infty & \text{otherwise} \end{cases}\) \(\frac{1}{\alpha} \cdot \max \left\{ 0, x \right\}\) \(0\)
\(\evaro_{\alpha}\) \(\begin{cases} y \log y &\text{if } y > 0 \\ 0 &\text{if } y = 0 \\ \infty &\text{if } y < 0 \end{cases}\) \(\exp(x - 1)\) \(\log \frac{1}{\alpha}\)
\(\tvaro_{\alpha}\) \(|y - 1|\) \(\begin{cases} x &\text{if } |x| \le 1\\ \infty &\text{otherwise} \end{cases}\) \(\min \left\{ \sqrt{ 2 \log \frac{1}{\alpha} }, 2\right\}\)

A \(\varphi\)-divergence risk measure \(\mathop{\mathrm{DR}}_{\varphi, \beta }\colon \mathbb{R}^n \to \bar{\Real}\) is defined for \(\varphi\colon \bar{\Real} \to \bar{\Real}\) and \(\beta\in \Real_{\ge 0}\) as [13] \[\begin{align} \mathop{\mathrm{DR}}_{\varphi, \beta}[ \tilde{x} ] \tag{2} &:= \inf \left\{\boldsymbol{q}\tr \boldsymbol{x} \mid \boldsymbol{q}\in \Delta_n, H_{\varphi}(\boldsymbol{q} \| \boldsymbol{p}) \le \beta, \boldsymbol{q} \ll \boldsymbol{p} \right\} \\ \tag{3} &= \sup_{z\in \Real} \left( z - \inf_{t>0} t\cdot \E { \varphi\opt \left( \frac{z - \tilde{x}}{t} + \beta \right)} \right), \end{align}\] where \(\varphi\opt(x)\) is the convex conjugate of \(\varphi(y)\), \(\varphi\) is a closed, convex function that satisfies \(\varphi(1) = 0\), and \(H_{\varphi}\) is a \(\varphi\)-divergence measure [22]: \[H_\varphi(\boldsymbol{q} \| \boldsymbol{p}) := \sum_{\omega \in \Omega} p_{\omega} \cdot \varphi\left(\frac{q_{\omega}}{p_{\omega}}\right),\] defined for \(\boldsymbol{q} \ll \boldsymbol{p}\). Here, \(\boldsymbol{q} \ll \boldsymbol{p}\) denotes that \(\boldsymbol{q}\) is absolutely continuous with respect to \(\boldsymbol{p}\), meaning that \(p_{\omega} = 0 \implies q_{\omega} = 0\) for each \(\omega\in 1{:}n\). We also adopt the common convention that \(1 / 0 = \infty\) and \(0 \cdot \infty = 0\) [23]. It is important to note that \(\varphi\)-divergence risk measures are coherent, which means they are concave and positively homogeneous. Any coherent risk measure \(\rho\colon \mathbb{R}^n\to \bar{\Real}\) satisfies that \(\rho[\tilde{x}] = \inf_{\boldsymbol{q}\in\mathcal{Q}} \boldsymbol{q}\tr\boldsymbol{x}\) for some closed convex non-empty set \(\mathcal{Q} \subseteq \Delta_n\) [4].

1 summarizes several common \(\varphi\)-divergence risk measures; see 4.2 for the relevant definitions for \(\alpha\in [0,1]\). We focus on the Conditional-Value-at-Risk (CVaR), which is a \(\varphi\)-divergence risk measure and defined equivalently for \(\tilde{x} \in \mathbb{R}^n\) and \(\alpha \in (0,1]\) as [4] \[\label{eq:primal95cvar} \cvaro_\alpha[ \tilde{x} ] \; :=\; \max_{z \in \Real} \left(z - \frac{1}{\alpha} \cdot \E{z -\tilde{x}}_+\right) = \min \left\{ \boldsymbol{x}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_{n}, \, \boldsymbol{q} \le \alpha^{-1} \boldsymbol{p} \right\},\tag{4}\] and extended continuously to \(\cvaro_0[\tilde{x}] = \mathop{\mathrm{ess}}\inf[\tilde{x}]\).

We also study another \(\varphi\)-divergence risk measure, which we call Total-Value-at-Risk (TVaR) and which is defined for \(\alpha\in (0,1]\) as \[\label{eq:tvar-dual} \tvaro_{\alpha}[\tilde{x}] \;:=\; \min \left\{ \boldsymbol{x}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_{n}, \, \|\boldsymbol{p} - \boldsymbol{q}\|_1 \le \min \left\{ \sqrt{2 \log \frac{1}{\alpha }}, 2\right\}, \boldsymbol{q} \ll \boldsymbol{p} \right\},\tag{5}\] and extended continuously to \(\tvaro_{0}[\tilde{x}] = \mathop{\mathrm{ess}}\inf [\tilde{x}]\).

We also refer to the Entropic-Value-at-Risk (EVaR), which is the \(\varphi\)-divergence risk measure obtained from \(\varphi(y) := y \log y\) for \(y > 0\) and \(\beta = \log \frac{1}{\alpha }\) for \(\alpha > 0\) in 1. Because the corresponding \(\varphi\)-divergence is the Kullback–Leibler divergence \(H_{\varphi}(\boldsymbol{q} \| \boldsymbol{p}) = \mathop{\mathrm{KL}}(\boldsymbol{q} \| \boldsymbol{p}) = \sum_{\omega \in \Omega} q_{\omega} \log \frac{q_{\omega}}{p_{\omega}}\), the EVaR is equivalently defined for \(\alpha\in (0,1]\) and \(\tilde{x} \in \mathbb{R}^n\) as [13] \[\label{eq:evar-dual} \evaro_{\alpha}[\tilde{x}] \;:=\; \min \left\{ \boldsymbol{x}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_{n}, \, \mathop{\mathrm{KL}}(\boldsymbol{q} \| \boldsymbol{p}) \le \log \frac{1}{\alpha }, \boldsymbol{q} \ll \boldsymbol{p} \right\},\tag{6}\] and extended continuously to \(\evaro_{0}[\tilde{x}] = \mathop{\mathrm{ess}}\inf [\tilde{x}]\).

Figure 1: Comparison of risk measure values for a normal distribution with mean 0 and standard deviation 1 as a function of the risk level \alpha.

An important property of TVaR is that it bounds EVaR [11], [13], [24] from below, as we show next and illustrate in 1.

For each \(\tilde{x} \in \mathbb{R}^n\), the risk measures above extended to \(\alpha\in [0,1]\) satisfy that \[\var{\alpha}{\tilde{x}} \; \ge\; \cvar{\alpha}{\tilde{x}} \; \ge\; \evar{\alpha}{\tilde{x}} \; \ge\; \tvar{\alpha}{\tilde{x}}.\]

Proof. The inequality for VaR and CVaR is well known, and the inequality for CVaR and EVaR follows from [13]. Recall that Pinsker’s inequality [25] shows for \(\boldsymbol{p}, \boldsymbol{q}\in \Delta_n\) that \(\| \boldsymbol{p} - \boldsymbol{q}\|_1 \le \sqrt{2 \mathop{\mathrm{KL}}(\boldsymbol{q} \| \boldsymbol{p})}\). Then we prove that TVaR lower-bounds EVaR as \[\begin{align} \evaro_{\alpha}[\tilde{x}] &= \inf \left\{ \boldsymbol{q}\tr \boldsymbol{x} \mid \boldsymbol{q}\in \Delta_n, \mathop{\mathrm{KL}}(\boldsymbol{q} \| \boldsymbol{p}) \le \log \frac{1}{\alpha }, \boldsymbol{q} \ll \boldsymbol{p} \right\} \\ &\ge \inf \left\{ \boldsymbol{q}\tr \boldsymbol{x} \mid \boldsymbol{q}\in \Delta_n, \| \boldsymbol{p} - \boldsymbol{q} \|_1 \le \sqrt{ 2 \log \frac{1}{\alpha }} , \boldsymbol{q} \ll \boldsymbol{p} \right\} \\ &=\inf \left\{ \boldsymbol{q}\tr \boldsymbol{x} \mid \boldsymbol{q}\in \Delta_n, \| \boldsymbol{p} - \boldsymbol{q} \|_1 \le \min \left\{ \sqrt{ 2 \log \frac{1}{\alpha }}, 2\right\}, \boldsymbol{q} \ll \boldsymbol{p}\right\} = \tvaro_{\alpha }[\tilde{x}]. \end{align}\] We use the fact that \(\| \boldsymbol{p} - \boldsymbol{q}\|_1 \le 2\). ◻

3 Computing VaR in Linear Time↩︎

We now describe the QuickVaR algorithm, show that it computes VaR, and show that it runs in linear expected time.

The standard VaR algorithm for a discrete random variable \(\boldsymbol{x}\in \Real^n\) with distribution \(\boldsymbol{p}\in \Delta_n\) and risk level \(\alpha\) proceeds as follows (see listing 6.25 [26]). It sorts the entries of \(\boldsymbol{x}\) to get a permutation \(\sigma \colon 1{:}n \to 1{:}n\) such that \(i \mapsto x_{\sigma(i)}\) is monotone. Then, it finds an index \(i \in 1{:}n\) such that \[\sum_{j = 1}^{i-1} p_{\sigma(j)} \; <\; \alpha \; \le\; \sum_{j = 1}^i p_{\sigma(j)}.\] The value \(x_{\sigma(i)}\) is the optimal solution in 1 by 1 below. The computational complexity of this algorithm is \(\Omega(n \log n)\) because its runtime is dominated by the need to sort the entries of \(\boldsymbol{x}\). For the sake of completeness, [alg:slow-quantile] in the appendix summarizes the standard greedy approach to computing VaR.

QuickVaR, described in [alg:quick-quantile], computes VaR in \(O(n)\) time by avoiding the need to sort the elements of the input random variable \(\boldsymbol{x}\). In particular, QuickVaR partially sorts the indices of \(\boldsymbol{x}\), analogous to the quickselect algorithm [27].

[alg:quick-quantile] is recursive and applies to a random variable \(\boldsymbol{x}\) with a distribution \(\boldsymbol{p}\in \Delta_n\) as follows. In addition to the random variable, the input is \(f\) (front) and \(b\) (back), which delimit the indices of \(\boldsymbol{x}\) it manipulates. The algorithm first picks a pivot element \(v\) uniformly at random from \(f{:}b\). The random pivot selection is essential to achieve an average linear time complexity. Given a pivot, we use \(\operatorname{partition}\) to partially sort the entries of \(\boldsymbol{x}\) into three groups: 1. less than \(v\): \(\boldsymbol{x}_{f{:}(l-1)} < v\), 2. equal to \(v\): \(\boldsymbol{x}_{l{:}g} = v\), and 3. greater than \(v\): \(\boldsymbol{x}_{(g+1){:}b} > v\). Then, the algorithm checks for the conditions of 1 to determine which of the three groups contains the VaR value. Finally, if the VaR value is in group 1 or 3, the algorithm calls itself recursively with adjusted \(f\), \(b\), and \(\alpha\).

The \(\operatorname{partition}\) procedure in [alg:quick-quantile] uses the Dutch national flag algorithm to partition the input array into the three groups. The algorithm runs in linear time, and the properties of its output are given in 4 below. For the sake of completeness, we provide the reader with the algorithm in 5 in 7.

We present [alg:quick-quantile] recursively to simplify its analysis. However, a loop-based, non-recursive implementation, such as [alg:quickvar-nonrecursive] in 7, typically performs better in imperative languages.

\(v \gets x_{\operatorname{rand}(f{:}b)}\) \(l, g, \boldsymbol{x}', \boldsymbol{p}' \gets \operatorname{partition}(\boldsymbol{x}, \boldsymbol{p}, v, f, b)\) \(t \gets \sum_{j=f}^{l-1} p'_j\) \(e \gets \sum_{j=l}^{g} p'_j\) ()\(\alpha < t\) ()\(t + e \le \alpha\)

The following theorem states the correctness of QuickVaR.

Theorem 1. [alg:quick-quantile] executed with \(f = 1\) and \(b = n\) and \(\alpha\in [0,1)\) returns \(\varo_{\alpha}[ \tilde{x} ]\).

Before proving 1, we need to show several important lemmas. The following lemma, which holds only for discrete random variables, shows that the VaR value is the probability atom that satisfies specific probability properties.

Lemma 1. Assuming \(\alpha\in [0,1)\), there exists \(\omega\in\Omega\) with \(v := \var{\alpha }{\tilde{x}} = x_{\omega}\) and \[\label{eq:var-property-prob} \P{\tilde{x} < v} \le \alpha < \P{\tilde{x} \le v}.\qquad{(1)}\]

Proof. First, we show that any \(v\) that satisfies ?? solves 1 optimally. The implication \(\P{\tilde{x} < v} \le \alpha \implies v \in \left\{ t \in \Real \mid \P{ \tilde{x} < t } \le \alpha \right\}\) follows immediately. Second, we prove \(\alpha < \P{ \tilde{x} \le v } \implies \forall t \in \left\{ t \in \Real \mid \P{ \tilde{x} < t } \le \alpha \right\}, v \ge t\) by contradiction. Suppose that \(\exists t \in \left\{ t \in \Real \mid \P{ \tilde{x} < t } \le \alpha \right\}, v < t\). Then, by the monotonicity of the probability operator, we get a contradiction as \[\alpha < \P { \tilde{x} \le v } \le \P { \tilde{x} < t } \le \alpha.\] Then, the existence of \(\omega\) follows from 3 because ?? implies that \[\P{\tilde{x} = \var{\alpha }{\tilde{x}}} = \P{\tilde{x} \le \var{\alpha }{\tilde{x}}} - \P{\tilde{x} < \var{\alpha }{\tilde{x}}} > \alpha - \alpha = 0.\] ◻

The following lemma shows the recursive property of QuickVaR and is the main building block of 1.

Lemma 2. Suppose that the input conditions of [alg:quick-quantile] hold and \(\boldsymbol{x} \in \Real^n\) satisfies that \[\label{eq:qvar-precondition} \mathbf{x}_{1{:}(f-1)} < x_f, \qquad x_{f-1} < \boldsymbol{x}_{f{:}b} < x_{b+1}, \qquad x_b < \boldsymbol{x}_{(b+1){:}n},\qquad{(2)}\] where \(x_0 = -\infty\) and \(x_{n+1} = +\infty\). Then [alg:quick-quantile] returns \(x_{\omega}, \omega\in f{:}b\) such that \[\label{eq:qvar-postcondition} \P{\tilde{x} < x_{\omega}} \quad \le \quad \alpha + \sum_{i=1}^{f-1} p_i \quad < \quad \P{\tilde{x} \le x_{\omega}}.\qquad{(3)}\]

Proof. From ?? and 4, \(t\) and \(e\) on lines [ln:t] and [ln:e] satisfy that \[\label{eq:te-equivalence} t = \P{\tilde{x} < v} - \sum_{i=1}^{f-1} p_i , \qquad e = \P{\tilde{x} = v}.\tag{7}\] We proceed by induction on \(k = b - f\). Base case: \(k = 0\). Since \(b = f\), we have that \(v = x_b\) and \(l = g = b\), \(t = 0\), \(e = p_b\), and \(\alpha \ge t\). The input condition implies that \(t + e = p_b > \alpha\). Because \(t \le \alpha < p_b\), the algorithm returns \(v\) on line [ln:thirdif] that satisfies \[\P{\tilde{x} < v} \le \alpha + \sum_{i=1}^{f-1} p_i < \P{\tilde{x} \le v}.\] Above, the first inequality follows from \(v = x_f\) and \(0 \le \alpha\). The second inequality follows from the algorithm’s input condition.

Inductive step. Suppose that \(k>0\) and that the result holds for each \(k' < k\). Then:

(i) If \(\alpha < t\), the algorithm executes the first if statement on line [ln:firstif]. Because \(0 \le \alpha < t\), the recursive call satisfies the required input conditions: \(f \le l-1\) and \(\sum_{j=f}^{l-1} p_j = t > \alpha\). The recursive call also satisfies ?? by 4, and we can conclude that ?? follows from the inductive hypothesis because \(l-1 - f \le b-1 -f \le k-1 < k\).

(ii) If \(t + e \le \alpha\), the algorithm executes the second if statement on line [ln:secondif]. We show that the input conditions for the recursive call are met. Because \(t + e = \sum_{j = f}^g p_j \le \alpha < \sum_{j=f}^b p_j\) (the second inequality is the input condition), it must be that \(g < b\) and therefore \(g+1 \le b\). Also, by the algorithm’s input condition and algebraic manipulation, we get that \[\sum_{j=g+1}^b p_j = \sum_{j=f}^b p_j - t - e > \alpha - t - e,\] and the input conditions hold. Because ?? follows from 4 and \(b - g - 1 < b-f \le k\), the inductive hypothesis implies that there exists \(x_{\omega}\) such that \[\P{\tilde{x} < x_{\omega}} \quad \le \quad (\alpha - t - e) + \P{\tilde{x} \le v} \quad < \quad \P{\tilde{x} \le x_{\omega}}.\] Then, ?? follows by algebraic manipulation after substituting \(t,e\) from 7 .

(iii) If \(t \le \alpha < t + e\), the algorithm executes the else statement on line [ln:thirdif], and the inequality ?? holds by substituting \(t,e\) from 7 . ◻

We are now ready to prove the correctness of QuickVaR.

Proof of 1. With \(f = 1\) and \(b = n\), the condition in ?? holds vacuously, and the input condition \(\sum_{j=1}^n p_j = 1 > \alpha\) holds since \(\alpha\in [0,1)\). Therefore, 2 shows that QuickVaR returns \(x_{\omega}\) that satisfies \[\P{\tilde{x} < x_{\omega}} \quad \le \quad \alpha \quad < \quad \P{\tilde{x} \le x_{\omega}},\] which satisfies ?? and, therefore, it equals the VaR value by 1. ◻

With the correctness of [alg:quick-quantile] established in 2, we now show its computational complexity.

Theorem 2. [alg:quick-quantile] runs in expected \(O(n)\) time.

Proof. The proof is analogous to the standard quickselect algorithm, as in [27], except that we need to keep track of the risk level \(\alpha\). Fix some \(n \in \mathbb{N}\), \(\boldsymbol{x}\), and \(\boldsymbol{p}\), and we study the algorithm’s properties. Let \(T\colon 1{:}n \to [0,\infty]\) be the expected runtime of [alg:quick-quantile] such that \(T(k)\) is the expected runtime when executed with \(k := b - f + 1\). To simplify the notation, we assume without loss of generality in the remainder of the proof that \(f = 1\) and \(b = k\). Then, let \[l := \left| \left\{ i \in 1{:}k \mid x_i < \varo_{\alpha}[\tilde{x}]\right\} \right|.\] Note that by the assumption \(\sum_{j = f}^b p_j > \alpha\), we have that \(k \ge 1\) and \(\varo_{\alpha}[\tilde{x}] \le \max_{i\in f{:}b} x_i\). Therefore, \(l \in 0{:}(k-1)\). The algorithm terminates in constant time when \(k = 1\), and we upper-bound the runtime for the worst-case choice of \(l\) for \(k > 1\) as \[\begin{align} \nonumber T(k) &\le d\cdot k + \max_{l\in 0{:}(k-1)} \frac{1}{k} \sum_{i=0}^{k-1} \left( T(i) \cdot 1_{l < i} + T(k - i - 1) \cdot 1_{l > i} \right) \\ \nonumber &= d\cdot k + \max_{l\in 0{:}(k-1)} \frac{1}{k} \left( \sum_{i=l+1}^{k-1} T(i) + \sum_{i=0}^{l-1} T(k - i - 1) \right) \\ \label{eq:quickvar-complexity} &= d\cdot k + \max_{l\in 0{:}(k-1)} \frac{1}{k} \left( \sum_{i=l+1}^{k-1} T(i) + \sum_{i=k-l}^{k-1} T(i) \right), \end{align}\tag{8}\] where \(d\cdot k\) is the local runtime of the algorithm for some \(d > 0\) and \(T(i)\) and \(T(k - i - 1)\) are the average runtimes of the recursive calls. The average over \(i\) corresponds to the uniformly random selection of the pivot with \(i = \left| \left\{ j \in 1{:}k \mid x_j < v \right\} \right|\) in the first line of the inequality. The last two equalities follow by algebraic manipulation. Then, an argument identical to theorem 3.3 [27], since 8 is the same as equation (3.8) in [27], shows that \(T(k)\le 4\cdot d\cdot k\) for all \(k\ge 1\). ◻

An attentive reader may notice that QuickVaR resembles the algorithm for computing the weighted median [18]. However, the weighted median algorithm does not, in general, compute VaR. Define the quantile \({\mathfrak{q}}_{\alpha},\alpha \in (0,1)\) for \(\tilde{x}\) as the closed set \[{\mathfrak{q}}_{\alpha}[\tilde{x}] := \left\{ z\in \Real \mid \P{ \tilde{x} \le z } \ge \alpha, \P{ \tilde{x} \ge z } \ge 1-\alpha \right\},\] with bounds \({\mathfrak{q}}^+_{\alpha}[\tilde{x}] := \max {\mathfrak{q}}_{\alpha}[\tilde{x}]\) and \({\mathfrak{q}}^-_{\alpha}[\tilde{x}] := \min {\mathfrak{q}}_{\alpha}[\tilde{x}]\). One can readily show that \(\varo_{\alpha }[\tilde{x}] = {\mathfrak{q}}^+_{\alpha}[\tilde{x}]\). Although the precise definition of a weighted median varies, it is typically a value in \({\mathfrak{q}}_{\frac{1}{2}}[\tilde{x}]\). Because \({\mathfrak{q}}^+_{\alpha}[\tilde{x}] \neq {\mathfrak{q}}^{-}_{\alpha}[\tilde{x}]\), in general, VaR also does not equal the median. A careful analysis shows that the weighted median algorithm in [18] computes \({\mathfrak{q}}^-_{\frac{1}{2}}[\tilde{x}]\), which does not equal \(\varo_{\frac{1}{2}}[\tilde{x}]\).

4 Computing \(\varphi\)-divergence Risk Measures in Linear Time↩︎

In this section, we describe and analyze QuickDivergence, a new linear-time framework for computing a class of \(\varphi\)-divergence risk measures, including CVaR and TVaR. QuickDivergence treats risk evaluation as a linear optimization over a polymatroid. Standard polymatroid optimization algorithms must sort their input, which requires \(\Omega(n \log n)\) runtime. We identify a structural condition of polymatroid risk measures that allows us to avoid sorting and instead use the linear-time QuickVaR algorithm described in 3.

The remainder of the section is organized as follows. 4.1 describes a general linear-time algorithm suitable for \(\varphi\)-divergence risk measures that satisfy an EWS property. Then 4.2 shows that both CVaR and TVaR satisfy the EWS property and can be computed in linear time. Note that we use a fixed distribution \(\boldsymbol{p} \in \Delta_n\) and omit it in the notation when it is obvious from the context.

4.1 QuickDivergence↩︎

We first recall the definition of a submodular function [18] and a polymatroid [18].

Definition 1. A function \(f \colon 2^{1{:}n} \to \Real_{\ge 0}\) is submodular* if \[f(\mathcal{A}) + f(\mathcal{B}) \; \geq\; f(\mathcal{A} \cup \mathcal{B}) + f(\mathcal{A} \cap \mathcal{B}), \quad \forall \mathcal{A}, \mathcal{B} \subseteq 1{:}n.\] For a submodular \(f \colon 2^{1{:}n} \to \Real_{\ge 0}\), a polymatroid \(\mathcal{Q}_f\) is a polytope defined as \[\mathcal{Q}_f \;:=\; \left\{ \boldsymbol{q} \in \Real_{\ge 0}^n \mid \boldsymbol{1}\tr \boldsymbol{q}_{\mathcal{U}} \leq f(\mathcal{U}), \, \forall \, \mathcal{U} \subseteq 1{:}n \right\},\] and a normalized polymatroid is \(\widehat{\mathcal{Q}}_f := \mathcal{Q}_f \cap \Delta_n\).*

Polymatroids are useful in optimization because any coherent, comonotonic monetary risk measure can be formulated as linear optimization over a polymatroid [4]. A risk measure \(\rho\) is coherent if it is concave and positively homogeneous. A risk measure \(\rho\) is comonotonic if for any comonotone \(\tilde{x}, \tilde{y}\): \(\rho[ \tilde{x} + \tilde{y} ] = \rho[ \tilde{x} ] + \rho[ \tilde{y} ]\), where the random variables \(\tilde{x}, \tilde{y} \in \Real^n\) are comonotone if \(\tilde{x}(i) \ge \tilde{x}(j) \Leftrightarrow \tilde{y}(i) \ge \tilde{y}(j)\) for each \(i,j\in 1{:}n\).

A monetary risk measure \(\rho\colon \Delta_n \times \Real^n \to \Real\) is a coherent comonotonic risk measure if and only if \[\label{eq:comonotonic-optimization} \rho[\tilde{x}] = \min_{\boldsymbol{q}\in \widehat{\mathcal{Q}}_c} \boldsymbol{q}\tr \boldsymbol{x}, \qquad \text{where} \qquad c(\mathcal{A}) := -\rho[ -\tilde{1}_{\mathcal{A}} ], \quad \forall \mathcal{A} \subseteq 1{:}n,\tag{9}\] where \(\tilde{1}\) is an indicator random variable. Moreover, if the risk measure is law-invariant, then there exists \(g\colon \Real \to \Real\) such that \(c(\mathcal{A}) = g(\boldsymbol{1}\tr \boldsymbol{p}_{\mathcal{A}})\) for each \(\mathcal{A} \subseteq 1{:}n\).

[prop:coherent-polymatroid] follows from corollary 4.88 and theorem 4.94 in [4].

The standard solution to the polymatroid optimization problem in 9 is to use a greedy algorithm, such as 2. This algorithm adapts the standard polymatroid optimization, such as [18], to the normalized polymatroid \(\widehat{\mathcal{Q}}_f\). See [prop:standard-polymatroid-works] for the proof of correctness. Note that the greedy approach in 2 requires sorting and, therefore, runs in \(\Omega(n \log n)\) time.

Figure 2: Linear minimization over polymatroids.

We now define a special class of submodular functions, called EWS, that makes sorting unnecessary when optimizing a linear function over a polymatroid.

Definition 2. A function \(f\colon 2^{1{:}n} \to \Real_{\ge 0}\) is \(\mathop{\mathrm{EWS}}(\boldsymbol{\mu} ,m, c)\), or Element-Wise-Separable, for \(c \in \Real_{\ge 0}\), \(m\in \Real_{>0}\), such that \(c \le 1\), \(c + m \ge 1\), and \(\boldsymbol{\mu}\in \Delta_n\) if \(f(\mathcal{A}) := g(\boldsymbol{1}\tr \boldsymbol{\mu}_{\mathcal{A}})\) for all \(\mathcal{A}\subseteq 1{:}n\) and \[\label{eq:def-g} g(x) := \begin{cases} \min \{c + m \cdot x, 1\} &\text{if } x > 0, \\ 0 &\text{otherwise}. \end{cases}\qquad{(4)}\]

EWS functions satisfy the properties necessary to define a polymatroid as we show next.

An EWS function \(f\colon 2^{\Omega} \to \Real_{\ge 0}\) is submodular, monotone, and satisfies that \(f(\emptyset) = 0\), and \(f(\Omega) = 1\).

Proof. To show that \(f\) is submodular, it is sufficient, from 5, to show that \(g\) is concave on \([0,1]\). The function \(g\) is concave on \((0, 1]\) because it is the minimum of two linear functions [28]. The concavity of \(g\) on \([0, 1]\) follows from the definition. Since \(g(0) = 0\), for each \(\lambda \in (0,1]\) and \(b \in ( 0, 1]\) we get that \[\begin{align} (1 - \lambda)g(0) + \lambda \cdot g(b) = \lambda \cdot g(b) = \min\left\{\lambda \cdot c + m\cdot \lambda\cdot b, \lambda\right\} \leq \min \left\{ c + m \cdot \lambda \cdot b, 1 \right\} = g\left( (1 - \lambda)\cdot 0 + \lambda\cdot b\right). \end{align}\] Because the inequality also holds for \(\lambda = 0\), \(g\) is concave on \([0,1]\) by the definition. The monotonicity, non-negativity and the boundary value conditions for \(f\) follow by algebraic manipulation from the assumptions on \(c\) and \(m\). ◻

Figure 3: QuickDivergence: Linear minimization over EWS polymatroids.

3 describes QuickDivergence, an algorithm for polymatroid optimization. QuickDivergence differs from the standard linear polymatroid algorithm in 2 in that it does not require sorting the input vector. Instead, the algorithm uses QuickVaR to partially sort the input vector.

The following theorem states the correctness and runtime of 3.

Theorem 3. 3 solves 9 in expected \(O(n)\) runtime.

Proof. To prove the correctness of 3, we show that it returns \(\boldsymbol{q}\opt\) that equals \(\hat{\boldsymbol{q}}\) returned by the standard 2. Suppose that 2 chooses a permutation \(\sigma\colon 1{:}n \to 1{:}n\) such that \(i \le j \Rightarrow x_{\sigma(i)} \le x_{\sigma(j)}\) and \(x_{\sigma(i)} = x_{\sigma(j)} \wedge i < j \Rightarrow \sigma(i) < \sigma(j)\) for each \(i,j\in 1{:}n\). In other words, the permutation \(\sigma\) needs to be generated by a stable sort algorithm. We emphasize that the stability assumption is only needed for our proof technique; the actual implementation does not require it.

We now prove that \(q\opt_i = \hat{q}_i\) for each \(i\in 1{:}n\) by separately analyzing the following four cases.

(i) Suppose that \(i = k_{\min}\). Because \(p_i > 0\) from the construction of \(k_{\min}\) then \[q\opt_i = \min \left\{ m\cdot p_i + c, 1 \right\} = g(p_i) = g(p_i) - g(0) = f(\{ i \} \cup \mathcal{J}) - f(\mathcal{J}) = \hat{q}_i,\] where \(\mathcal{J} := \{ j \in 1{:}n \mid \sigma^{-1}(j) < \sigma^{-1}(i) \}\) and \(\boldsymbol{1}_{\mathcal{J}}\tr \boldsymbol{p} = 0\) by the construction of \(k_{\min}\).

(ii) Suppose that \(i \in \mathcal{K}_{<v} \wedge i \neq k_{\min}\). From 1, we get for each \(\mathcal{O} \subseteq \mathcal{K}_{<v}\) that \[\label{eq:lt-v-prob} \boldsymbol{1}\tr \boldsymbol{p}_{\mathcal{O}} \le \P{\tilde{x} < v} \le \frac{1-c}{m},\tag{10}\] and, hence, from ?? : \(f(\mathcal{O}) = g\left(\boldsymbol{1}\tr \boldsymbol{p}_{\mathcal{O}}\right) = c + m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\mathcal{O}}, \text{ when } \boldsymbol{1}\tr \boldsymbol{p}_{\mathcal{O}} > 0.\) Let \(k := \sigma^{-1}(i)\) and note that \(\sigma(1{:}k) \subseteq \mathcal{K}_{< v}\) because if \(l \in 1{:}k\), then \(x_{\sigma(l)} \le x_{\sigma(k)} = x_i < v\) and \(\sigma(l) \in \mathcal{K}_{< v}\). Then, since \(i \neq k_{\min}\): \[\hat{q}_i = \hat{q}_{\sigma(k)} = f(\sigma(1{:}k)) - f(\sigma(1{:}(k-1))) = m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}k) \setminus \sigma(1{:}(k-1)) } = m\cdot p_{\sigma(k)} = m\cdot p_i = q\opt_i.\] The third equality follows from \(i \neq k_{\min}\) and \(p_{k_{\min}} \neq 0\).

(iii) Suppose that \(i \in \mathcal{K}_{=v} \wedge i \neq k_{\min}\). Here, we use the assumption that \(\sigma\) is a stable sort and prove the claim by induction on \(l \in |\mathcal{K}_{<v} \cup \{k_{\min}\}|, \dots, |\mathcal{K}_{<v}\cup \{k_{\min}\}| + |\mathcal{K}_{=v} \setminus \{ k_{\min}\}|\) showing that \[\begin{align} \tag{11} \theta_l &= 1 - f(\sigma(1{:}l)), \\ \tag{12} \hat{q}_{\sigma(l+1)} &= q\opt_{\sigma(l+1)}. \end{align}\] To prove the base case, suppose that \(l = |\mathcal{K}_{<v} \cup \{k_{\min}\}|\). The equality in 11 holds by algebraic manipulation from ?? since \(\mathcal{K}_{<v} \cup \{ k_{\min} \} = \sigma(1{:}l)\). The equality in 12 follows from ?? and 10 . First, assume that \(\theta_l > 0\); then \[\label{eq:q-q-eq} \begin{align} \hat{q}_{i} = \hat{q}_{\sigma(l+1)} &= f ( \sigma(1{:}(l+1) )) - f ( \sigma(1{:}l) ) \\ &= \min \left\{ 1,\; c + m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}(l+1))} \right\} - c - m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}l)} \\ &= \min \left\{1 - c - m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}l)}, \; m \cdot p_{\sigma(l+1)}\right\} \\ &= \min \{ \theta_l, m \cdot p_i \} = q_i\opt. \end{align}\tag{13}\] When \(\theta_l = 0\), we get by algebraic manipulation that \(\hat{q}_i = 0 = q\opt_i\).

To prove the inductive case, suppose that 11 and 12 hold for \(l\), and we prove it for \(l+1\). The equality in 11 can be derived from 12 and the construction in 2 as \[\begin{align} \theta_{l+1} = \theta_l - q\opt_{\sigma(l+1)} = \theta_l - \hat{q}_{\sigma(l+1)} = 1 - f(\sigma(1{:}l)) - \left( f(\sigma(1{:}(l+1))) - f(\sigma(1{:}l)) \right) = 1 - f(\sigma(1{:}(l+1))). \end{align}\] The equality in 12 follows analogously to 13 if \(\theta_{l+1} > 0\): \[\label{eq:q-q-eq-2} \begin{align} \hat{q}_{i} = \hat{q}_{\sigma(l+2)} &= f ( \sigma(1{:}(l+2) )) - f(\sigma(1{:}(l+1))) \\ &= \min \left\{ 1,\; c + m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}(l+2))} \right\} - (1-\theta_{l+1}) \\ &= \min \left\{\theta_{l+1}, \;c + m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}(l+2))} - f(\sigma(1{:}(l+1))) \right\} \\ &= \min \left\{\theta_{l+1}, \; m \cdot p_{\sigma(l+2)} \right\} = q_i\opt. \end{align}\tag{14}\] If \(\theta_{l+1} = 0\), we get by algebraic manipulation that \(\hat{q}_i = 0 = q\opt_i\).

(iv) Suppose that \(i \in \mathcal{K}_{>v} := \{ j \in 1{:}n \mid x_j > v \}\) and let \(k = \sigma^{-1}(i)\). Then, using ?? we have that: \[\begin{align} 1 \ge f(\sigma(1{:}k)) \ge f(\sigma(1{:}(k-1))) \ge \min \{c + m \cdot \boldsymbol{1}\tr \boldsymbol{p}_{\sigma(1{:}(k-1))}, 1\} \ge 1, \end{align}\] since \(\boldsymbol{1}\tr\boldsymbol{p}_{\sigma(1{:}(k-1))} > \frac{1-c}{m}\) from 1 and the construction of \(v\) in 3. Then \[\hat{q}_i = \hat{q}_{\sigma(k)} = f(\sigma(1{:}k)) - f(\sigma(1{:}(k-1))) = 0 = q\opt_i.\] ◻

4.2 Solving CVaR and TVaR in Linear Time↩︎

We show that CVaR and TVaR are two examples of \(\varphi\)-divergence risk measures that can be formulated as linear optimization over a polymatroid with an EWS capacity function and, therefore, can be solved in linear time. To formulate the CVaR polymatroid optimization, we define the capacity function \(f_{\alpha}^{\cvaro}\colon 2^{\Omega } \to \Real\) for each \(\mathcal{O} \subseteq \Omega\) and \(\alpha \in [0,1]\) as \[\label{eq:f-cvar} f_{\alpha}^{\cvaro}(\mathcal{O}) \;:=\; \begin{cases} \min\left\{\frac{1}{\alpha} \cdot \P{\mathcal{O}}, 1\right\} & \text{if } \alpha \in (0,1], \\ 0 & \text{if } \alpha = 0 \wedge \P{\mathcal{O}} = 0, \\ 1 & \text{if } \alpha = 0 \wedge \P{\mathcal{O}} > 0. \end{cases}\tag{15}\] The Choquet capacity function \(f_{\alpha }^{\cvaro}\) for CVaR defined in 15 is well known, but is usually restricted to \(\alpha\in (0,1)\) [29]. The following proposition states its correctness.

For each \(\alpha\in [0,1]\), \(\tilde{x} \colon \Omega \to \Real\), and \(\boldsymbol{p} \in \Delta_n\): \[\label{eq:cvar-correct} \cvar{\alpha }{\tilde{x} } \;=\; \min_{\boldsymbol{q} \in \widehat{\mathcal{Q}}_{f_{\alpha }^{\cvaro}}} \boldsymbol{x}\tr \boldsymbol{q}.\tag{16}\] Moreover, \(f^{\cvaro}_{\alpha }\) is \(\mathop{\mathrm{EWS}}(\boldsymbol{p}, 1/\alpha, 0)\) when \(\alpha > 0\), and \(\mathop{\mathrm{EWS}}(\boldsymbol{p}, 1, 1)\) when \(\alpha = 0\).

To compute TVaR, we define the capacity function \(f_{\alpha}^{\tvaro}\colon 2^{\Omega } \to \Real\) for each \(\mathcal{O} \subseteq \Omega\) and \(\alpha \in [0,1]\) as \[\label{eq:f-tvar} f^{\tvaro}_{\alpha}(\mathcal{O}) \; :=\; \begin{cases} \min\left\{\P{\mathcal{O}} + \min \left\{ \sqrt{\frac{1}{2} \log \frac{1}{\alpha}}, 1\right\}, 1\right\} &\text{if } \alpha \in (0,1] \wedge \P{\mathcal{O}} > 0, \\ 0 & \text{if } \P{\mathcal{O}} = 0, \\ 1 & \text{if } \alpha = 0 \wedge \P{\mathcal{O}} > 0. \end{cases}\tag{17}\] TVaR was first proposed in the context of \(\varphi\)-divergence risk measures [30] as well as in the context of robust stochastic programming [31]. TVaR is also very closely related to the \(L_1\) ball used in robust Markov decision processes [8], [10]. However, we are not aware of prior work that defines or uses the capacity function in 17 . The following proposition shows that this function can be used to correctly formulate the TVaR evaluation using an EWS function.

For each \(\alpha\in [0,1]\), \(\tilde{x} \colon \Omega \to \Real\), and \(\boldsymbol{p} \in \Delta_n\): \[\label{eq:tvar-correct} \tvar{\alpha}{\tilde{x} } = \min_{\boldsymbol{q}\in \widehat{\mathcal{Q}}_{f_{\alpha}^{\tvaro}}} \boldsymbol{x}\tr \boldsymbol{q}.\tag{18}\] Moreover, \(f^{\tvaro}_{\alpha }\) is \(\mathop{\mathrm{EWS}}(\boldsymbol{p}, 1, \min \left\{ \sqrt{\frac{1}{2} \log \frac{1}{\alpha}}, 1 \right\})\) when \(\alpha > 0\), and \(\mathop{\mathrm{EWS}}(\boldsymbol{p}, 1, 1)\) when \(\alpha = 0\).

5 Experiments↩︎

a

b

Figure 4: Comparison of runtime of the algorithms for samples from a uniform random variable (left) and a sparse random variable (right). Shaded regions denote \(95\%\) confidence bands over the \(10\) runs..

We evaluate the proposed algorithms empirically to confirm that their theoretical performance guarantees translate into faster practical runtimes. We compare each linear-time algorithm against its standard sorting-based counterpart. Throughout this section, we prefix the name of our linear-time implementation with a “Q” to distinguish it from the sorting-based baseline, giving QVaR, QCVaR, and QTVaR. We also report the runtime of computing the standard expectation \(\mathbb{E}\) as a baseline.

Our implementation, together with the code and data needed to reproduce all results, is available at https://github.com/marekpetrik/paper_lineartime_risk, and the actual implementation of the risk measures can be found at https://github.com/RiskAverseRL/RiskMeasures.jl. All experiments were run on an AMD Ryzen Threadripper 3970X with 32 cores and 64 threads at \(4.55\) GHz with \(256\) GB of RAM using Julia 1.12.5. To reduce noise in the timing measurements, we disable the Julia garbage collector during each run.

We begin with synthetic data drawn from a uniform distribution with \(n\) samples, varying \(n\) from \(10^6\) to \(10^7\). For each value of \(n\) we run every algorithm \(10\) times on independently generated samples and report the mean runtime. 4 shows the results on a log-log scale for random variables with a uniform (left) and sparsely-supported (right) distributions. We fix the risk level to \(\alpha = 0.95\), though the results are similar for other levels except the degenerate cases \(\alpha = 0\) and \(\alpha = 1\). The shaded regions denote \(95\%\) confidence bands over the \(10\) runs; the bands are tight because the measured runtimes vary little between runs.

The results confirm that the quick algorithms are substantially faster than their sorting-based counterparts, with the gap widening as \(n\) grows. This behavior matches the theoretical analysis of [sec:var-linear] [sec:phi-divergence]: the sorting-based methods require \(O(n \log n)\) operations, whereas our algorithms run in \(O(n)\) expected time, so the logarithmic factor separating them grows with the size of the probability space. The acceleration is most consequential for VaR and CVaR, which are the most widely used risk measures in the applications discussed in 1 and are frequently evaluated on large random variables.

Table 2: Mean runtime of algorithms in the stock market scenario.
Method Mean (ms) St. dev. (ms)
VaR 0.37 0.00
QVaR 0.13 0.02
CVaR 0.39 0.00
QCVaR 0.22 0.03
TVaR 0.40 0.01
QTVaR 0.25 0.05
\(\mathbb{E}\) 0.01 0.00
Table 3: Mean runtime (ms) of algorithms on smaller sparse random domains.
n \(\mathbb{E}\) VaR QVaR CVaR QCVaR TVaR QTVaR
1000 0.00 0.03 0.01 0.03 0.02 0.03 0.02
2000 0.00 0.06 0.02 0.06 0.03 0.06 0.04
3000 0.01 0.09 0.04 0.10 0.05 0.10 0.07
4000 0.00 0.13 0.04 0.14 0.06 0.13 0.08
5000 0.01 0.16 0.05 0.17 0.10 0.18 0.10
6000 0.00 0.20 0.08 0.21 0.10 0.21 0.11
7000 0.00 0.23 0.08 0.25 0.12 0.24 0.13
8000 0.00 0.27 0.09 0.28 0.13 0.28 0.15
9000 0.01 0.30 0.09 0.32 0.12 0.31 0.17
10000 0.01 0.33 0.11 0.35 0.16 0.35 0.19

To confirm the speedups persist beyond large synthetic inputs, we repeat the comparison on two more realistic settings. The first one is a Stock Market scenario, in which the random variable models the returns of a financial portfolio. 2 reports the mean runtime and standard deviation of each method. The second one is small-scale problems, where the asymptotic advantage is least pronounced but can be practically important when the risk measures must be computed repeatedly. 3 reports the mean runtime on sparse random variables over probability spaces of size \(n\) ranging from \(1000\) to \(10\,000\), using the same protocol as above but with \(50\) trials per configuration in place of \(10\). Even at these sizes, the quick algorithms outperform their sorting-based counterparts for every method and every value of \(n\), confirming that the linear-time algorithms deliver consistent improvements across the full range of input sizes without incurring any penalty on small inputs.

6 Conclusion↩︎

We propose and analyze new algorithms for computing several common risk measures, including VaR, CVaR, and TVaR, in \(O(n)\) expected time, where \(n\) is the size of the probability space. This compares favorably with the \(O(n \log n)\) time required by sorting-based methods. This speedup is significant because risk measures are frequently evaluated on large random variables in applications such as robotics, finance, and infrastructure maintenance. An important direction for future work is to develop fast algorithms for other risk measures, such as EVaR and Expectile.

7 Supplemental Results for Section 3↩︎

In this section, we summarize the standard and supporting algorithms for solving VaR. We also include the selected technical proofs for the correctness of the proposed algorithm.

Lemma 3. For each \(t\in \Real\): \[\P{\tilde{x} = t} > 0 \implies \exists \omega \in \Omega, x(\omega) = t.\]

Proof. The result follows by induction on the size \(|\Omega|\), dropping the constraint \(\boldsymbol{1}\tr \boldsymbol{p} = 1\). ◻

Sort \(\boldsymbol{x}\) indices to get \(\sigma \colon 1{:}n \to 1{:}n\) such that \(i > j \implies x_{\sigma(i)} \ge x_{\sigma(j)}, \; \forall i, j \in 1{:}n\) \(q_0 \gets 0\) \(q_i \gets q_{i-1} + p_{\sigma(i)}\)

For each \(\boldsymbol{x} \in \Real^n\), \(\boldsymbol{p}\in \Delta_n\), and \(\alpha\in [0,1)\), [alg:slow-quantile] returns \(\var{\alpha}{\tilde{x}}\).

Proof. First, one can readily show by induction that for each \(i\in 1{:}n\) (note that values of \(\boldsymbol{x}\) need not be unique): \[\label{eq:proof-q-j} \P{\tilde{x} < x_{\sigma(i)}} \le q_{i-1} \le q_i \le \P{\tilde{x} \le x_{\sigma(i)}}.\tag{19}\] Since \(\alpha < 1 = q_n\), the algorithm is guaranteed to terminate with some \(i\opt \in 1{:}n\) and return \(x_{\sigma(i\opt)}\). Therefore, \(q_{i\opt-1} \le \alpha < q_{i\opt}\) and applying 19 gives us that \[\P{\tilde{x} < x_{\sigma(i\opt)}} \le q_{i\opt-1} \le \alpha < q_{i\opt} \le \P{\tilde{x} \le x_{\sigma(i\opt)}},\] and the result holds by 1. ◻

Figure 5: partition: Dutch national flag algorithm [32]

\(f_0 \gets 1\);    \(b_0 \gets n\);    \(k\gets 1\) \(v_k \gets x_{\operatorname{rand}(f_{k-1}{:}b_{k-1})}\) \(l_k, g_k, \boldsymbol{x}^k, \boldsymbol{p}^k \gets \operatorname{partition}(\boldsymbol{x}^{k-1}, \boldsymbol{p}^{k-1}, v_k, f_{k-1}, b_{k-1})\) \(t_k \gets \sum_{j=f_{k-1}}^{l_k-1} p^k_j\) \(e_k \gets \sum_{j = l_k}^{g_k} p^{k}_j\) \(f_k \gets f_{k-1}\)\(b_k \gets l_k-1\)\(\alpha_k \gets \alpha_{k-1}\) \(f_k \gets g_k+1\)\(b_k \gets b_{k-1}\)\(\alpha_k \gets \alpha_{k-1} - t_k - e_k\) \(k \gets k + 1\)

The following lemma proves the properties of the \(\operatorname{partition}\) operator defined in 5.

Lemma 4. If \((l, g, \boldsymbol{x}', \boldsymbol{p}') = \operatorname{partition}(\boldsymbol{x}, \boldsymbol{p}, v, f, b)\) for \(f \le b\), then: \[\begin{align} \label{eq:hoare-correspond} \boldsymbol{x}_{f{:}(l-1)}' < v, \qquad \boldsymbol{x}_{l{:}g}' = v, \qquad \boldsymbol{x}_{(g+1){:}b}' > v. \end{align}\qquad{(5)}\] Moreover, \(\boldsymbol{x}_{f{:}b} = \boldsymbol{x}'_{\sigma(f{:}b)}, \boldsymbol{p}_{f{:}b} = \boldsymbol{p}'_{\sigma(f{:}b)}\) for some permutation \(\sigma\colon f{:}b \to f{:}b\), and \(\boldsymbol{x}_{1{:}(f-1)} = \boldsymbol{x}'_{1{:}(f-1)}, \boldsymbol{p}_{1{:}(f-1)} = \boldsymbol{p}'_{1{:}(f-1)}\), and \(\boldsymbol{x}_{(b+1):n} = \boldsymbol{x}'_{(b+1):n}, \boldsymbol{p}_{(b+1):n} = \boldsymbol{p}'_{(b+1):n}\). The elements outside of the range \(f{:}b\) are not modified. The algorithm runs in \(O(b - f + 1)\) time.

Proof of 4. Let \(v\) be an input to 5 and consider the following invariants. \[\begin{align} \tag{20} x_j &< v, & j\in f_k{:}(l_k-1), \\ \tag{21} x_j &= v, & j \in l_k{:}(e_k-1), \\ \tag{22} x_j &> v, & j \in (g_k+1){:}b_k. \end{align}\]

Initialization. For \(k=0\), conditions 202122 are true because the pertinent intervals are empty.

Maintenance. At each iteration of the loop that starts on line [line:partition-while] of 5, there are 3 possible branches determined by the conditional logic: \(x_{e_k} < v\), \(x_{e_k} > v\), and \(x_{e_k} = v\). Since \(v, x_{e_k} \in \Real\) one of these branches must be executed for each iteration \(k\).

Each time the \(x_{e_k} < v\) or \(x_{e_k} > v\) branch is executed, conditions 20 and 22 , respectively, remain true since in each case we swap the necessary elements to preserve each condition. Notice also that \(e_k\) is incremented alongside \(l_k\) to maintain \(e_k - l_k = e_{k+1} - l_{k+1}\) as is necessary for condition 21 .

When \(x_{e_k} = v\), \(e_k\) is incremented; since \(e_0 = l_0\), this creates a difference between \(e_{k+1}\) and \(l_{k+1}\). The algorithm terminates when \(e_t = g_t + 1\), so \(g_t - l_t = e_t - l_t - 1 \ge 0\) and 21 remains true.

Termination. In the loop on line [line:partition-while] there are 3 possible paths of execution determined by the conditional logic \(x_{e_k} < v\), \(x_{e_k} > v\), and \(x_{e_k} = v\). Each of these paths of execution decreases \(g_k - e_k\) by 1. This occurs on lines [line:partition-e-1-1] and [line:partition-e-1-2] by increasing \(e_k\) by 1 and keeping \(g_k\) the same, or on line [line:partition-g-1] by decreasing \(g_k\) by 1 and keeping \(e_k\) the same.

Observe that \(g_0 - e_0 = b - f\) is non-negative since \(f\le b\). The difference \(g_k-e_k\) decreases by 1 at each iteration, and each iteration requires constant time to complete. Because the algorithm returns when \(g_t - e_t = -1\), it finishes in \(O(b-f)\) time.

On the last iteration \(g_t = e_t\), therefore by the same reasoning as above \(g_{t+1} - e_{t+1} = -1\) and condition 21 holds for \(g_{t+1} = e_{t+1} - 1\). ◻

8 Supplemental Results for Section 4↩︎

The following proposition states the correctness of the standard polymatroid optimization algorithm when applied to a normalized polymatroid.

Suppose that \(f\colon 2^{1{:}n} \to \Real\) is submodular, monotone, \(f(\emptyset) = 0\) and \(f(1{:}n) = 1\), and let \(\boldsymbol{x}\in \Real^n_{\ge 0}\). Then, 2 computes \(\hat{\boldsymbol{q}}\) such that \[\label{eq:normalized-poly-opt} \hat{\boldsymbol{q}} \in \mathop{\mathrm{arg\,min}}_{\boldsymbol{q}\in \widehat{\mathcal{Q}}_f} \boldsymbol{x}\tr \boldsymbol{q}.\tag{23}\]

Proof. First, we show that \(\hat{\boldsymbol{q}} \in \widehat{\mathcal{Q}}_f\). Because the construction of \(\hat{\boldsymbol{q}}\) in 2 matches the construction in proposition 14.10 in [18], since \(f(\emptyset) = 0\) and \(\hat{q}_{\sigma(1)} = f(\left\{ \sigma(1) \right\}) - f(\emptyset) = f(\left\{ \sigma(1) \right\})\), we have that \(\hat{\boldsymbol{q}} \in \mathcal{Q}_f\). In addition, \[\sum_{i=1}^n \hat{q}_i = \sum_{i=1}^n \hat{q}_{\sigma (i)} = f(1{:}n) - f(\emptyset) = 1,\] which implies that \(\hat{\boldsymbol{q}} \in \Delta_n \cap \mathcal{Q}_f = \widehat{\mathcal{Q}}_f\), \(\widehat{\mathcal{Q}}_f \neq \emptyset\), and the minimum 23 exists.

To show that \(\hat{\boldsymbol{q}}\) is optimal, we reduce the optimization in 23 to the standard linear polymatroid maximization problem with non-negative weights [18] as follows: \[\begin{align} \min_{\boldsymbol{q}\in \widehat{\mathcal{Q}}_f} \boldsymbol{x}\tr \boldsymbol{q} &= \min_{\boldsymbol{q}\in \widehat{\mathcal{Q}}_f} (\boldsymbol{x} - x_{\max}\cdot \boldsymbol{1}) \tr \boldsymbol{q} + x_{\max}\\ &= -\max_{\boldsymbol{q}\in \widehat{\mathcal{Q}}_f} (x_{\max}\cdot \boldsymbol{1} - \boldsymbol{x}) \tr \boldsymbol{q} + x_{\max}\\ &\ge -\max_{\boldsymbol{q}\in \mathcal{Q}_f} (x_{\max}\cdot \boldsymbol{1} - \boldsymbol{x}) \tr \boldsymbol{q} + x_{\max}\\ &= - (x_{\max}\cdot \boldsymbol{1} - \boldsymbol{x}) \tr \hat{\boldsymbol{q}} + x_{\max} \\ &= \boldsymbol{x}\tr \hat{\boldsymbol{q}}. \end{align}\] Here, \(x_{\max} := \max_{i\in 1{:}n} x_i\) and the equalities follow because \(\widehat{\mathcal{Q}}_f\) is a normalized polymatroid and \(\boldsymbol{1}\tr \hat{\boldsymbol{q}} = 1\). The second-to-last equality follows from theorem 14.11 in [18] and the construction in the greedy polymatroid algorithm. ◻

The following lemma is folklore in submodular optimization [33]. We state and prove it for completeness.

Lemma 5. Suppose that \(\boldsymbol{\mu}\in \Real^n_{\ge 0}\) and \(g \colon \Real_{\ge 0} \to \Real\) is concave and monotone. Then \(f(\mathcal{A}) := g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}})\) is a submodular monotone function.

Proof. The fact that \(f\) is monotone follows from the composition of two monotone functions. To prove that \(f\) is submodular, we use the following properties that hold for each \(\mathcal{A}, \mathcal{B} \subseteq 1{:}n\): \[\begin{align} \tag{24} \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}} &\le \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}} \le \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}, & \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}} &\le \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}} \le \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}, \\ \tag{25} \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}} + \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}} &= \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} + \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}. \end{align}\] The inequalities in 24 hold because \(\boldsymbol{\mu} \ge \boldsymbol{0}\) and the equality in 25 holds because \(\mathcal{A} \mapsto \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}}\) is modular. Suppose that \(\boldsymbol{1}\tr \boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}} < \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}\); otherwise the submodularity condition in 28 below holds with equality. Then, writing \(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}}\) and \(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}}\) as a convex combination of \(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}\cap \mathcal{B}}\) and \(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}\) using 24 and the concavity of \(g\) we get that \[\begin{align} \tag{26} g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}}) &\ge \frac{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}}}{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}} g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}) + \frac{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}}{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}} g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}), \\ \tag{27} g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}}) &\ge \frac{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}}}{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}} g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}) + \frac{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}}{\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}} - \boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}} g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}). \end{align}\] Then summing 26 and 27 and algebraic manipulation yields that \[\label{eq:submod-proof-cond} f(\mathcal{A}) + f(\mathcal{B}) = g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A}}) + g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{B}}) \ge g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cap \mathcal{B}}) + g(\boldsymbol{1}\tr\boldsymbol{\mu}_{\mathcal{A} \cup \mathcal{B}}) = f(\mathcal{A} \cap \mathcal{B}) + f(\mathcal{A} \cup \mathcal{B}),\tag{28}\] which proves that \(f\) is submodular. ◻

Proof of [prop:cvar-ews-correct]. The function \(f_{\alpha }^{\cvaro}\) is EWS by algebraic manipulation. The equality 16 for \(\alpha = 0\) follows by algebraic manipulation. To prove 16 for \(\alpha > 0\), we use [prop:coherent-polymatroid] and show for each \(\mathcal{O} \subseteq \Omega\) that \[\label{eq:cvar-eq-submodular} f_{\alpha}^{\cvaro}(\mathcal{O}) = \min\left\{ \alpha^{-1} \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p},\; 1 \right\} = \max \left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_{n}, \, \boldsymbol{q} \le \alpha^{-1} \boldsymbol{p} \right\} = -\cvar{\alpha }{- \tilde{1}_{\mathcal{O}}},\tag{29}\] using the dual CVaR formulation in 4 . We prove 29 as two inequalities.

(i) The inequality \[\min\left\{ \alpha^{-1} \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p},\; 1 \right\} \ge \max \left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_{n}, \, \boldsymbol{q} \le \alpha^{-1} \boldsymbol{p} \right\},\] holds because \(\boldsymbol{q}\) feasible in the right-hand side satisfies that \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \le \alpha^{-1} \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p}, \text{ and } \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \le \boldsymbol{1}\tr \boldsymbol{q} = 1.\)

(ii) To prove that \[\label{eq:cvar-proof-le} \min\left\{ \alpha^{-1} \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p},\; 1 \right\} \le \max \left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_{n}, \, \boldsymbol{q} \le \alpha^{-1} \boldsymbol{p} \right\},\tag{30}\] we analyze the following two cases.

  1. When \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} \ge \alpha\), let \(\boldsymbol{q} := (\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p})^{-1} \cdot \boldsymbol{p} \odot \boldsymbol{1}_{\mathcal{O}},\) which is feasible in the right-hand side of 30 and \(\odot\) is the Hadamard vector product.

  2. When \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} < \alpha\), let \(\boldsymbol{q} := \alpha^{-1} \cdot \boldsymbol{p} \odot \boldsymbol{1}_{\mathcal{O}} + d \cdot \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}),\) which is feasible in the right-hand side of 30 and \(d\in \Real\) solves that \(d \cdot \boldsymbol{1}\tr \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}) = 1 - \alpha^{-1} \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p}.\) Such \(d\) always exists because \(\boldsymbol{1}\tr \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}) = 1 - \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} > 1 - \alpha \ge 0.\) It is easy to verify that \(\boldsymbol{q} \ge \boldsymbol{0}\) and \(\boldsymbol{1}\tr \boldsymbol{q} = 1\). The condition \(\boldsymbol{q} \le \alpha^{-1} \boldsymbol{p}\) holds because \(d \le \alpha^{-1}\) is satisfied by multiplying both sides of the inequality by \(\boldsymbol{1}\tr \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}})\): \[d \cdot \boldsymbol{1}\tr \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}) \; \le\; 1 - \frac{1}{\alpha } \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} \; \le\; \frac{1}{\alpha } \cdot \boldsymbol{1}\tr \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}).\]

 ◻

Proof of [prop:tvar-ews-correct]. The function \(f_{\alpha }^{\tvaro}\) is EWS by algebraic manipulation. The equality 18 for \(\alpha\) such that \(\sqrt{\frac{1}{2} \log \frac{1}{\alpha} } \ge 1\), including \(\alpha = 0\), follows by algebraic manipulation.

It remains to prove 18 when \(\sqrt{\frac{1}{2} \log \frac{1}{\alpha} } < 1\). We use [prop:coherent-polymatroid] and show for each \(\mathcal{O} \subseteq \Omega\) that \[\label{eq:tvar-eq-submodular} \begin{align} f_{\alpha}^{\tvaro}(\mathcal{O}) &= \min\left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} + \sqrt{\frac{1}{2} \log \frac{1}{\alpha}},\; 1 \right\} \\ &= \max \left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_n, \, \|\boldsymbol{p} -\boldsymbol{q}\|_1 \le \sqrt{2 \log \frac{1}{\alpha}} \right\} = -\tvar{\alpha }{- \tilde{1}_{\mathcal{O}}}, \end{align}\tag{31}\] using the TVaR definition in 5 . We need to handle two special cases first. When \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} = 0\), then we have that \(0 = f^{\tvaro}_{\alpha}(\mathcal{O}) = g(0) = -\tvaro_{\alpha}[-\tilde{1}_{\mathcal{O}}]\) since \(\boldsymbol{q} \ll \boldsymbol{p}\) in the TVaR definition in 5 . When \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} = 1\), then we have that \(1 = f^{\tvaro}_{\alpha}(\mathcal{O}) = g(1) = -\tvaro_{\alpha}[-\tilde{1}_{\mathcal{O}}]\) since \(\boldsymbol{q} \ll \boldsymbol{p}\) and \(\boldsymbol{1}\tr \boldsymbol{q} = 1\) in the TVaR definition in 5 .

Finally, we prove 31 as two inequalities when \(0 < \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} < 1\). (i) To prove that \[\min\left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} + \sqrt{\frac{1}{2} \log \frac{1}{\alpha}},\; 1 \right\} \ge \max \left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_n, \, \|\boldsymbol{p} -\boldsymbol{q}\|_1 \le \sqrt{2 \log \frac{1}{\alpha}} \right\},\] observe that \(\boldsymbol{q}\) feasible in the right-hand side implies that \(\|\boldsymbol{p} - \boldsymbol{q}\|_1 = \boldsymbol{1}_{\mathcal{O}}\tr |\boldsymbol{p} - \boldsymbol{q}| + (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}})\tr |\boldsymbol{p} - \boldsymbol{q}| \leq \sqrt{2 \log \frac{1}{\alpha }}\) and \((\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}})\tr (\boldsymbol{p} - \boldsymbol{q}) = \boldsymbol{1}_{\mathcal{O}}\tr (\boldsymbol{q} - \boldsymbol{p})\) since \(\boldsymbol{p},\boldsymbol{q}\in \Delta_n\). Then \[\begin{align} \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} &= \boldsymbol{1}_{\mathcal{O}}\tr (\boldsymbol{q} - \boldsymbol{p}) + \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} = \frac{1}{2}\left( \boldsymbol{1}_{\mathcal{O}}\tr (\boldsymbol{q} - \boldsymbol{p}) + (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}})\tr (\boldsymbol{p} - \boldsymbol{q}) \right)+ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} \\ &\le \frac{1}{2}\left( \boldsymbol{1}_{\mathcal{O}}\tr |\boldsymbol{q} - \boldsymbol{p}| + (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}})\tr |\boldsymbol{p} - \boldsymbol{q}| \right)+ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} = \sqrt{\frac{1}{2} \log \frac{1}{\alpha }} + \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p}, \end{align}\] and \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \leq \boldsymbol{1}\tr \boldsymbol{q} = 1.\)

(ii) To prove that \[\min\left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} + \sqrt{\frac{1}{2} \log \frac{1}{\alpha }},\; 1 \right\} \le \max \left\{ \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{q} \mid \boldsymbol{q} \in \Delta_n, \, \|\boldsymbol{p} -\boldsymbol{q}\|_{1} \le \sqrt{2 \log \frac{1}{\alpha }} \right\},\] we analyze the following two cases:

  1. When \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} + \sqrt{\frac{1}{2} \log \frac{1}{\alpha }} \ge 1\), let \(\boldsymbol{q} \;:=\; (\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p})^{-1} \cdot \boldsymbol{p} \odot \boldsymbol{1}_{\mathcal{O}},\) where \(\odot\) is the Hadamard vector product. This \(\boldsymbol{q}\) is feasible in the right-hand side and satisfies the inequality.

  2. When \(\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} + \sqrt{\frac{1}{2} \log \frac{1}{\alpha }}< 1\), let \[\boldsymbol{q} \;:=\; \frac{\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} + \sqrt{\frac{1}{2} \log \frac{1}{\alpha }}}{\boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p}} \cdot \boldsymbol{p} \odot \boldsymbol{1}_{\mathcal{O}} + \frac{(\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}) \tr \boldsymbol{p} - \sqrt{\frac{1}{2} \log \frac{1}{\alpha }} }{(\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}) \tr \boldsymbol{p}}\cdot \boldsymbol{p} \odot (\boldsymbol{1} - \boldsymbol{1}_{\mathcal{O}}),\] which is feasible in the right-hand side and satisfies the inequality. Neither one of the divisions is by \(0\) since \(0 < \boldsymbol{1}_{\mathcal{O}}\tr \boldsymbol{p} < 1\).

 ◻

References↩︎

[1]
P. Akella et al., “Risk-aware robotics: Tail risk measures in planning, control, and verification,” arXiv preprint arXiv:2403.18972, 2024.
[2]
A. Majumdar and M. Pavone, “How should a robot assess risk? Towards an axiomatic theory of risk in robotics,” in Robotics research: The 18th international symposium ISRR, 2020, pp. 75–84.
[3]
M. Benrabah, C. Orou Mousse, E. Randriamiarintsoa, R. Chapuis, and R. Aufrère, “A review on traversability risk assessments for autonomous ground vehicles: Methods and metrics,” Sensors, vol. 24, no. 6, p. 1909, 2024.
[4]
H. Follmer and A. Schied, Stochastic finance: Introduction in discrete time, Fourth. De Gruyter Graduate, 2016.
[5]
A. Inzunza, R. Moreno, A. Bernales, and H. Rudnick, CVaR constrained planning of renewable generation with consideration of system inertial response, reserve services and demand participation,” Energy Economics, vol. 59, pp. 104–117, 2016.
[6]
F. Barahona, M. Ettl, M. Petrik, and P. M. Rimshnick, “Agile logistics simulation and optimization for managing disaster responses.” in Winter Simulation Conference, 2013, pp. 3340–3351.
[7]
M. Petrik, G. Tennenholtz, and M. Ghavamzadeh, “Bayesian regret minimization in offline bandits,” in Proceedings of the 41st International Conference on Machine Learning, 2024, pp. 40502–40522, Accessed: Apr. 28, 2025. [Online].
[8]
C. P. Ho, M. Petrik, and W. Wiesemann, “Partial policy iteration for L1-robust Markov decision processes,” Journal of Machine Learning Research, vol. 22, pp. 1–46, 2021.
[9]
C. P. Ho, M. Petrik, and W. Wiesemann, “Fast Bellman updates for robust MDPs,” in International Conference on Machine Learning (ICML), 2018, pp. 1979–1988.
[10]
C. P. Ho, M. Petrik, and W. Wiesemann, “Robust phi-divergence MDPs,” in Neural Information Processing Systems (NeurIPS), 2022, [Online]. Available: https://arxiv.org/abs/2205.14202.
[11]
J. L. Hau, E. Delage, M. Ghavamzadeh, and M. Petrik, “On dynamic programming decompositions of static risk measures in Markov decision processes,” in Neural Information Processing Systems (NeurIPS), 2023, Accessed: May 04, 2023. [Online]. Available: https://arxiv.org/abs/2304.12477.
[12]
J. L. Hau, E. Delage, E. Derman, M. Ghavamzadeh, and M. Petrik, “Q-learning for quantile MDPs: A decomposition, performance, and convergence analysis,” in Artificial Intelligence and Statistics (AISTATS), 2025, Accessed: Jul. 22, 2025. [Online]. Available: https://arxiv.org/abs/2410.24128.
[13]
A. Ahmadi-Javid, “Entropic Value-at-Risk: A new coherent risk measure,” Journal of Optimization Theory and Applications, vol. 155, no. 3, pp. 1105–1123, 2012, Accessed: Jan. 03, 2022. [Online].
[14]
L. Adam and V. Mácha, “Projections onto the canonical simplex with additional linear inequalities.” 2019, [Online]. Available: https://arxiv.org/abs/1905.03488.
[15]
W. Wang and M. Á. Carreira-Perpiñán, “Projection onto the probability simplex: An efficient algorithm with a simple proof, and an application.” 2013, [Online]. Available: https://arxiv.org/abs/1309.1541.
[16]
A. Ang, J. Ma, N. Liu, and K. Huang, “Fast projection onto the capped simplex with applications to sparse regression in bioinformatics,” in Conference on neural information processing systems, 2021, vol. 35.
[17]
C. H. Lim and S. J. Wright, “Efficient Bregman projections onto the permutahedron and related polytopes,” in Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, 2016, pp. 1205–1213, Accessed: Sep. 11, 2024. [Online].
[18]
B. Korte and J. Vygen, Combinatorial optimization: Theory and algorithms, Fifth. Springer, 2012.
[19]
L. Condat, “Fast projection onto the simplex and the L1 ball,” Mathematical Programming, vol. 158, no. 1–2, pp. 575–585, 2016.
[20]
T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction to algorithms, Third. Cambridge, MA: MIT Press, 2009.
[21]
A. Shapiro, D. Dentcheva, and A. Ruszczynski, Lectures on stochastic programming: Modeling and theory. SIAM, 2014.
[22]
G. Bayraksan and D. K. Love, “Data-driven stochastic programming using phi-divergences,” in The Operations Research Revolution, INFORMS, 2015, pp. 1–19.
[23]
R. T. Rockafellar and R. J.-B. Wets, Variational analysis. Springer, 1996.
[24]
X. Su, J. Grand-Clement, and M. Petrik, “Risk-averse total-reward MDPs with ERM and EVaR,” in AAAI Conference on Artificial Intelligence, 2025.
[25]
S. Boucheron, G. Lugosi, and P. Massart, Concentration inequalities: A nonasymptotic theory of independence. Oxford University Press, 2013.
[26]
N. Privault, “Notes on financial risk and analytics with R and Python,” Nanyang Technological University, 2026.
[27]
E. Horowitz, S. Sahni, and S. Rajasekaran, Computer algorithms. New York: Computer Science Press, 1998.
[28]
S. P. Boyd and L. Vandenberghe, Convex optimization. Cambridge University Press, 2004.
[29]
C. Acerbi, “Spectral measures of risk: A coherent representation of subjective risk aversion,” Journal of Banking & Finance, vol. 26, no. 7, pp. 1505–1518, 2002, Accessed: Apr. 29, 2026. [Online].
[30]
T. Kruse, J. C. Schneider, and N. Schweizer, “A toolkit for robust risk assessment using f-divergences,” Management Science, vol. 67, no. 10, pp. 6529–6552, 2021, Accessed: Dec. 03, 2024. [Online].
[31]
A. Shapiro, “Distributionally robust stochastic programming,” SIAM Journal on Optimization, vol. 27, no. 4, pp. 2258–2275, 2017.
[32]
E. W. Dijkstra, A discipline of programming. Englewood Cliffs, N.J: Prentice-Hall, 1976.
[33]
J. Bilmes, “Submodular functions, optimization, and applications to machine learning — spring quarter, lecture 4.” University of Washington, 2018.

  1. Equal contribution.↩︎