Sensitivity analysis of Wasserstein distributionally robust optimization problems

We start with the simple example of risk-neutral pricing of a call option written on an underlying asset $(S_{t})_{t\leq T}$. Here, $T,K>0$ are the maturity and the strike respectively, $f(x,a)=(S_{0}x-K)^{+}$ and $\mu$ is the distribution of $S_{T}/S_{0}$. We set interest rates and dividends to zero for simplicity. In the Black and Scholes, (1973) model, $\mu$ is a log-normal distribution, i.e., $\log(S_{T}/S_{0})\sim\mathcal{N}(-\sigma^{2}T/2,\sigma^{2}T)$ is Gaussian with mean $-\sigma^{2}T/2$ and variance $\sigma^{2}T$. In this case, $V(0)$ is given by the celebrated Black-Scholes formula. Note that this example is particularly simple since $f$ is independent of $a$. However, to ensure risk-neutral pricing, we have to impose a linear constraint on the measures in $B_{\delta}(\mu)$, giving

To compute its sensitivity we encode the constraint using a Lagrangian and apply Theorem 2, see (Bartl et al., 2021a, , Rem. 11, Thm. 15). For $p=2$, letting $k=K/S_{0}$ and $\mu_{k}=\mu([k,\infty))$, the resulting formula, see (Bartl et al., 2021a, , Example 18), is given by

Let us specialise to the log-normal distribution of the Black-Scholes model above and denote the quantity in (3) as $\mathcal{R}BS(\delta)$. It may be computed exactly using methods from Bartl et al., (2020) and Figure 1 compares the exact value and the first order approximation.

We have $\Upsilon=S_{0}\sqrt{\Phi(d_{-})(1-\Phi(d_{-}))}$, where $d_{-}=\frac{\log(S_{0}/K)-\sigma^{2}T/2}{\sigma\sqrt{T}}$ and $\Phi$ is the cdf of $\mathcal{N}(0,1)$ distribution. It is also insightful to compare $\Upsilon$ with a parametric sensitivity. If instead of Wasserstein balls, we consider $\{\mathcal{N}(-\tilde{\sigma}^{2}T/2,\tilde{\sigma}^{2}T):|\sigma-\tilde{\sigma}|\leq\delta\}$ the resulting sensitivity is known as the Black-Scholes Vega and given by $\mathcal{V}=S_{0}\Phi^{\prime}(d_{-}+\sigma\sqrt{T})$. We plot the two sensitivities in Figure 2. It is remarkable how, for the range of strikes of interest, the non-parametric model sensitivity $\Upsilon$ traces out the usual shape of $\mathcal{V}$ but shifted upwards to account for the idiosyncratic risk of departure from the log-normal family. More generally, given a book of options with payoff $f=f^{+}-f^{-}$ at time $T$, with $f^{+},f^{-}\geq 0$, we could say that the book is $\Upsilon$-neutral if the sensitivity $\Upsilon$ was the same for $f^{+}$ and for $f^{-}$. In analogy to Delta-Vega hedging standard, one could develop a non-parametric model-agnostic Delta-Upsilon hedging. We believe these ideas offer potential for exciting industrial applications and we leave them to further research.

We turn now to the classical notion of optimized Certainty Equivalent (OCE) of Ben Tal and Teboulle, (1986). It is a decision theoretic criterion designed to split a liability between today and tomorrow’s payments. It is also a convex risk measure in the sense of Artzner et al., (1999) and covers many of the popular risk measures such as expected shortfall or entropic risk, see Ben Tal and Teboulle, (2007). We fix a convex monotone function $l\colon\mathbb{R}\to\mathbb{R}$ which is bounded from below and $g\colon\mathbb{R}^{d}\to\mathbb{R}$. Here $g$ represent the payoff of a financial position and $l$ is the negative of a utility function, or a loss function. We take $\|\cdot\|=|\cdot|$ and refer to (Bartl et al., 2021a, , Lemma 22) for generic sufficient conditions for Assumptions 1 and 4 to hold in this setup. The OCE corresponds to $V$ in (1) for $f(x,a)=l(g(x)-a)+a$ and $\mathcal{A}=\mathbb{R}$, $\mathcal{S}=\mathbb{R}^{d}$. Theorems 2 and 5 yield the sensitivities

where, for simplicity, we took $p=q=2$ for the latter.
A related problem considers hedging strategies which minimise the expected loss of the hedged position, i.e., $f(x,a)=l\left(g(x)+\langle a,x-x_{0}\rangle\right)$, where $\mathcal{A}=\mathbb{R}^{k}$ and $(x_{0},x)$ represent today and tomorrow’s traded prices. We compute $\Upsilon$ as

Further, we can combine loss minimization with OCE and consider $a=(H,m)\in\mathbb{R}^{k}\times\mathbb{R}$, $f(x,(h,m))=l(g(x)+\langle H,x-x_{0}\rangle+m)-m$. This gives $V^{\prime}(0)$ as the infimum over $(H^{\star},m^{\star})\in\mathcal{A}^{\star}_{0}$ of

The above formulae capture non-parametric sensitivity to model uncertainty for examples of key risk measurements in financial economics. To the best of our knowledge this has not been achieved before.

Finally, we consider briefly the classical mean-variance optimization of Markowitz, (1952). Here $\mu$ represents the loss distribution across the assets and $a\in\mathbb{R}^{d}$, $\sum_{i=1}^{d}a_{i}=1$ are the relative investment weights. The original problem is to minimise the sum of the expectation and $\gamma$ standard deviations of returns $\langle a,X\rangle$, with $X\sim\mu$. Using the ideas in (Pflug et al.,, 2012, Example 2) and considering measures on $\mathbb{R}^{d}\times\mathbb{R}^{d}$, we can recast the problem as (1). Whilst Pflug et al., (2012) focused on the asymptotic regime $\delta\to\infty$, their non-asymptotic statements are related to our Theorem 2 and either result could be used here to obtain that $V(\delta)\approx V(0)+\sqrt{1-\gamma^{2}}\delta$.

We specialise now to quantifying robustness of neural networks (NN) to adversarial examples. This has been an important topic in machine learning since Szegedy et al., (2013) observed that NN consistently misclassify inputs formed by applying small worst-case perturbations to a data set. This produced a number of works offering either explanations for these effects or algorithms to create such adversarial examples, e.g. Goodfellow et al., (2014); Li et al., (2019); Carlini and Wagner, (2017); Wong and Kolter, (2017); Weng et al., (2018); Araujo et al., (2019); Mangal et al., (2019) to name just a few. The main focus of research works in this area, see Bastani et al., (2016), has been on faster algorithms for finding adversarial examples, typically leading to an overfit to these examples without any significant generalisation properties. The viewpoint has been mainly pointwise, e.g., Szegedy et al., (2013), with some generalisations to probabilistic robustness, e.g., Mangal et al., (2019).

In contrast, we propose a simple metric for measuring robustness of neural networks which is independent of the architecture employed and the algorithms for identifying adversarial examples. In fact, Theorem 2 offers a simple and intuitive way to formalise robustness of neural networks: for simplicity consider a $1$-layer neural network trained on a given distribution $\mu$ of pairs $(x,y)$, i.e. $(A^{\star}_{1},A_{2}^{\star},b_{1}^{\star},b_{2}^{\star})$ solve

where the $\inf$ is taken over $a=(A_{1},A_{2},b_{1},b_{2})\in\mathcal{A}=\mathbb{R}^{k\times d}\times\mathbb{R}^{d\times k}\times\mathbb{R}^{k}\times\mathbb{R}^{d}$, for a given activation function $\sigma:\mathbb{R}\to\mathbb{R}$, where the composition above is understood componentwise. Set $f(x,y;A,b):=|y-(A_{2}(\cdot)+b_{2})\circ\sigma\circ(A_{1}(\cdot)+b_{1})(x)|^{p}$. Data perturbations are captured by $\nu\in B_{\delta}^{p}(\mu)$ and (2) offers a robust training procedure. The first order quantification of the NN sensitivity to adversarial data is then given by

A similar viewpoint, capturing robustness to adversarial examples through optimal transport lens, has been recently adopted by other authors. The dual formulation of (2) was used by Shafieezadeh-Abadeh et al., (2019) to reduce the training of neural networks to tractable linear programs. Sinha et al., (2020) modified (2) to consider a penalised problem $\inf_{a\in\mathcal{A}}\sup_{\nu\in\mathcal{P}(\mathcal{S})}\int_{\mathcal{S}}f\left(x,a\right)\,\nu(dx)-\gamma W_{p}(\mu,\nu)$ to propose new stochastic gradient descent algorithms with inbuilt robustness to adversarial data.

In the context of UQ the measure $\mu$ represents input parameters of a (possibly complicated) operation $G$ in a physical, engineering or economic system. We consider the so-called reliability or certification problem: for a given set $E$ of undesirable outcomes, one wants to control $\sup_{\nu\in\mathcal{P}}\nu(G(x)\in E)$, for a set of probability measures $\mathcal{P}$. The distributionally robust adversarial classification problem considered recently by Ho-Nguyen and Wright, (2020) is also of this form, with Wasserstein balls $\mathcal{P}$ around an empirical measure of $N$ samples. Using the dual formulation of Blanchet and Murthy, (2019), they linked the problem to minimization of the conditional value-at-risk and proposed a reformulation, and numerical methods, in the case of linear classification. We propose instead a regularised version of the problem and look for

for a given safety level $\alpha$. We thus consider the average distance to the undesirable set, $d(G(x),E):=\inf_{e\in E}|G(x)-e|$, and not just its probability. The quantity $\delta(\alpha)$ could then be used to quantify the implicit uncertainty of the certification problem, where higher $\delta$ corresponds to less uncertainty. Taking statistical confidence bounds of the empirical measure in Wasserstein distance into account, see Fournier and Guillin, (2014), $\delta$ would then determine the minimum number of samples needed to estimate the empirical measure.

Assume that $E$ is convex. Then $x\mapsto d(x,E)$ differentiable everywhere except at the boundary of $E$ with $\nabla_{x}d(x,E)=0$ for $x\in{E}^{o}$ and $|\nabla_{x}d(x,E)|=1$ for all $x\in\bar{E}^{c}$. Further, assume $\mu$ is absolutely continuous w.r.t. Lebesgue measure on $\mathcal{S}$. Theorem 2, using (Bartl et al., 2021a, , Remark 11), gives a first-order expansion for the above problem

In the special case $\nabla_{x}G(x)=cI$ this simplifies to

and the minimal measure $\nu$ pushes every point $G(x)$ not contained in $E$ in the direction of the orthogonal projection. This recovers the intuition of (Chen et al.,, 2018, Theorem 1), which in turn rely on (Gao and Kleywegt,, 2016, Corollary 2, Example 7). Note however that our result holds for general measures $\mu$. We also note that such an approximation could provide an ansatz for dimension reduction, by identifying the dimensions for which the partial derivatives are negligible and then projecting $G$ on to the corresponding lower-dimensional subspace (thus providing a simpler surrogate for $G$). This would be an alternative to a basis expansion (e.g., orthogonal polynomials) used in UQ and would exploit the interplay of properties of $G$ and $\mu$ simultaneously.

We discuss two applications of our results in the realm of statistics. We start by highlighting the link between our results and the so-called influence curves (IC) in robust statistics. For a functional $\mu\mapsto T(\mu)$ its IC is defined as

Computing the $\mathrm{IC}$, if it exists, is in general hard and closed form solutions may be unachievable. However, for the so-called M-estimators, defined as optimizers for $V(0)$,

for some $f$ (e.g., $f(x,a)=|x-a|$ for the median), we have

under suitable assumptions on $f$, see (Huber and Ronchetti,, 1981, section 3.2.1). In comparison, writing $T^{\delta}$ for the optimizer for $V(\delta)$, Theorem 5 yields

under Assumption 4 and normalisation $\|\nabla_{x}f(x,T(\mu))\|_{L^{p}(\mu)}=1$. To investigate the connection let us Taylor-expand $\mathrm{IC}(y)$ around $x$ to obtain

Choosing $y=x+\delta\nabla f_{x}(x,T(\mu))$ and integrating both sides over $\mu$ and dividing by $\delta$, we obtain the asymptotic equality

for $\delta\to 0$ by (4). We conclude that considering the average directional derivative of IC in the direction of $\nabla f_{x}(x,T(\mu))$ gives our first-order sensitivity. For an interesting conjecture regarding the comparison of influence functions and sensitivities in KL-divergence we refer to (Lam,, 2018, Section 7.3) and (Lam,, 2016, Section 3.4.2).

Our second application in statistics exploits the representation of the LASSO/Ridge regressions as robust versions the standard linear regression. We consider $\mathcal{A}=\mathbb{R}^{k}$ and $\mathcal{S}=\mathbb{R}^{k+1}$. If instead of the Euclidean metric we take $\|(x,y)\|_{\ast}=|x|_{r}\mathbf{1}_{\{y=0\}}+\infty\mathbf{1}_{\{y\neq 0\}}$, for some $r>1$ and $(x,y)\in\mathbb{R}^{k}\times\mathbb{R}$, in the definition of the Wasserstein distance, then Blanchet et al., 2019a showed that

holds, where $1/r+1/s=1$. The $\delta=0$ case is the ordinary least squares regression. For $\delta>0$, the RHS for $s=2$ is directly related to the Ridge regression, while the limiting case $s=1$ is called the square-root LASSO regression, a regularised variant of linear regression well known for its good empirical performance. Closed-form solutions to (5) do not exist in general and it is a common practice to use numerical routines to solve it approximately. Theorem 5 offers instead an explicit first-order approximation of $a^{\star}_{\delta}$ for small $\delta$. We denote by $a^{\star}$ the ordinary least squares estimator and by $I$ the $k\times k$ identity matrix. Note that the first order condition on $a^{\star}$ implies that $\int(y-\langle a^{\star},x\rangle)x_{i}\mu(dx,dy)=0$ for all $1\leq i\leq k$. In particular, $V(0)=\int(y^{2}-\langle a^{\star},x\rangle y)\mu(dx,dy)$ and $a^{\star}=D^{-1}\int yx\mu(dx,dy)$, where we assume the system is overdetermined so that $D=\int xx^{T}\,\mu(dx,dy)$ is invertible. Letting $J=a^{\star}x^{T}+(Ia^{\star})(Ix)$ a direct computation, see Bartl et al., 2021a , yields

For $s=2$, $h(a^{\star})=a^{\star}/|a^{\star}|_{2}$ and for $s=1$, $h(a^{\star})=\text{sign}(a^{\star})$ and hence¹¹1The case $s=1$, inspecting the proof, we see that Theorem 5 still holds since $a^{\star}$ does not have zero components $\mu$-a.s., which are the only points of discontinuity of $h$. $a^{\star}_{\delta}$ is approximately

respectively. This corresponds to parameter shrinkage: proportional for square-root Ridge and a shift towards zero for square-root LASSO. To the best of our knowledge these are first such results and we stress that our formulae are valid in a general context and, in particular, parameter shrinkage depends on the direction through the $D^{-1}$ factor. Figure 3 compares the first order approximation with the actual results and shows a remarkable fit.

Furthermore, our results agree with what is known in the canonical test case for the (standard) Ridge and LASSO, see Tibshirani, (1996), when $\mu=\mu_{N}$ is the empirical measure of $N$ i.i.d. observations, the data is centred and the covariates are orthogonal, i.e., $D=\frac{1}{N}I$. In that case, (7) simplifies to

where $R^{2}$ is the usual coefficient of determination.

The case of $\mu_{N}$ is naturally of particular importance in statistics and data science and we continue to consider it in the next subsection. In particular, we characterise the asymptotic distribution of $\sqrt{N}(a^{\star}_{1/\sqrt{N}}-a^{\star})$, where $a^{\star}_{\delta}\in\mathcal{A}^{\star}_{\delta}(\mu_{N})$ and $a^{\star}\in\mathcal{A}^{\star}_{0}(\mu_{\infty})$ is the optimizer of the non-robust problem for the data-generating measure. This recovers the central limit theorem of Blanchet et al., 2019b , a link we explain further in section 44.2.

A benchmark of paramount importance in optimization is the so-called out-of-sample error, also known as the prediction error in statistical learning. Consider the setup above when $\mu_{N}$ is the empirical measure of $N$ i.i.d. observations sampled from the “true” distribution $\mu=\mu_{\infty}$ and take, for simplicity, $\|\cdot\|=|\cdot|_{s}$, with $s>1$. Our aim is to compute the optimal $a^{\star}$ which solves the original problem (1). However, we only have access to the training set, encoded via $\mu_{N}$. Suppose we solve the distributionally robust optimization problem (2) for $\mu_{N}$ and denote the robust optimizer $a^{{\star},N}_{\delta}$. Then the out-of-sample error

quantifies the error from using $a^{{\star},N}_{\delta}$ as opposed to the true optimizer $a^{\star}$.

While this expression seems to be hard to compute explicitly for finite samples, Theorem 5 offers a way to find the asymptotic distribution of a (suitably rescaled version of) the out-of-sample error. We suppose the assumptions in Theorem 5 are satisfied and note that the first order condition for $a^{\star}$ gives $\nabla_{a}V(0,a^{\star})=0$. Then, a second-order Taylor expansion gives

for some $\tilde{a}$, (coordinate-wise) between $a^{\star}$ and $a^{{\star},N}_{\delta}$. Now we write

where we define $a^{{\star},N}$ as the optimizer of the non-robust problem (1) with $\mu$ replaced by $\mu_{N}$. In particular the $\delta$-method for M-estimators implies that

where $H\sim\mathcal{N}\left(0,\int(\nabla_{a}f(x,a^{{\star}}))^{T}\nabla_{a}f(x,a^{{\star}})\,\mu(dx)\right)$ and $\Rightarrow$ denotes the convergence in distribution. On the other hand, for a fixed $N\in\mathbb{N}$, Theorem 5 applied to $\mu_{N}$ yields

where

Almost surely (w.r.t. sampling of $\mu_{N}$), we know that $\mu_{N}\to\mu$ in $W_{p}$ as $N\to\infty$, see Fournier and Guillin, (2014), and under the regularity and growth assumptions on $f$ in (Bartl et al., 2021a, , eq. (27)) we check that $\Delta_{N}\to 0$ a.s., see (Bartl et al., 2021a, , Example 28) for details. In particular, taking $\delta=1/\sqrt{N}$ and combining the above with (9) we obtain