[Week 4] Robust Optimization [3]

이번에는 이전에 살펴보았던 예시를 일반화해서, constraint 만족을 높은 확률로 보장하는 기회 제한 최적화(chance-constrained optimization)에 대해 알아보자. Robust optimization에서 solution이 너무 보수적(too conservative)인지 혹은 리스크가 너무 큰지, 문제 상황의 불확실성을 잘 반영하는지 여부는 문제 표현의 실용성을 평가하는 중요한 기준이 된다. 기회 제한 최적화에서 우리는 어떤 분포(알려져 있을 수도 있고, 아닐 수도 있다.)를 따르는 확률변수 $U$에 대해 아래의 문제를 풀려고 한다.

$\text{minimize}\,\,f_0(x)$
$\text{subject to}\,\,\mathbb{P}(f_i(x, U) > 0) \leq \epsilon,\,\,i=1,2,...,m$

여기에서 $\epsilon$은 constraint를 만족하지 못할 확률로, $0.05,\, 0.01$ 등의 작은 값을 보통 선택한다. 이런 형태의 문제는 일반적으로 non-convex이기 때문에, convex 문제로 safe approximation을 시도한다. 다시 말해 어떤 convex function $g_i(x) \leq 0$이면 $\mathbb{P}(f_i(x, U) > 0) \leq \epsilon$이 보장되도록 하는 $g_i$를 찾는 것이다. 가장 직접적인 해결책은 적절한 uncertainty set $\mathcal{U}$를 찾는 것이다. 구체적으로,
$\mathbb{P}(U \in \mathcal{U}) \geq 1-\epsilon$이면 $f_i(x, u) \leq 0\,\,\text{for}\,\,u \in \mathcal{U}$가 $\mathbb{P}(f_i(x, U) > 0) \leq \epsilon$을 보장한다. 이 때 확률 조건을 만족하는 $\mathcal{U}$를 최대한 작게 설정하는 것이 문제가 된다. 또는 uncertainty set을 직접 선택하는 대신 chance constraint의 확률을 approximation하는 식의 간접적인 접근도 가능하며, 지금부터 차근차근 살펴보자.

Value at Risk (VaR)

확률변수 $Z$에 대해 value at risk $\epsilon$은 다음과 같이 계산된다.

$\textbf{VaR}(Z;\epsilon) =\inf\left\{\gamma\,|\,\mathbb{P}(Z \leq \gamma) \geq 1-\epsilon\right\}=\inf\left\{\gamma\,|\,\mathbb{P}(Z > \gamma) \leq \epsilon\right\}.$

즉, $Z$가 $\gamma$보다 $1-\epsilon$의 확률로 작거나 같음이 보장되는 가장 작은 $\gamma$값이다. 금융 분야에서는 보통 $1-\epsilon =\eta$로 반대로 표기하는데, 이 때 $Z$를 자산의 수익률이라고 하면 $\textbf{VaR}(Z;\eta)$는 자산 수익률의 lower $\eta$ quantile의 lower bound, 즉 최대화 문제의 objective가 된다. 컨벡스 최적화 문제에서, $Z = f(x, U)$, 즉 높은 확률로 0보다 작거나 같길(만족되길) 원하는 constraint라고 하면
$\textbf{VaR}(Z;\epsilon) \leq 0\,\,\text{iff}\,\,\mathbb{P}(Z > 0) \leq \epsilon$
이다. 다만 일부 특수한 경우를 제외하고 value at risk 역시 다루기 어려운 non-convex constraint를 도출하기 때문에, 조금 다른 접근이 필요하다.

Example. normal distribution and value at risk

Value at risk가 convex constraint를 도출하는 특수한 예시이다. $U \in \mathbb{R}^n$이 평균이 $\mu$, covariance matrix가 $\sigma^2I$인 정규분포를 따르고 $f(x,U)$가 $x$와 $U$에 대한 쌍선형 형식(bilinear form)일 때, 즉 $f(x,u) = \sum_{i,k=1}^n a_{ik}x_i u_k$를 만족할 때 value at risk는 convex constraint를 도출한다. $Z = x^TU$라고 하자. 그러면
$\mathbb{P}(Z \leq \gamma) = \Phi\left(\frac{\gamma - \mu^Tx}{\sigma||x||_2}\right)$이므로
$\textbf{VaR}(U^Tx - \gamma;\epsilon) \leq 0$은 $\gamma \geq \mu^T x + \sigma \Phi^{-1}(1-\epsilon)||x||_2$의 second-order cone constraint와 동치이다. $(\epsilon \leq 1/2)$

Safe convex approximations for chance constraints

Chance constraint를 $1 -\epsilon$의 확률로 보장할 수 있는 convex function을 construct하는 방법으로 돌아가 보자. 확률변수 $Z$에 대해 $\phi: \mathbb{R} \rightarrow \mathbb{R}_+$가 non-negative, non-decreasing이면 $1(z \geq 0) \leq \phi(z)$임이 성립한다고 하자. 여기로부터
$\mathbb{P}(Z \geq 0) \leq E\phi(\alpha^{-1}Z)\,\,\text{for}\,\,\alpha >0$이 성립함을 확인할 수 있다. (Indicator function의 기댓값은 확률과 같고, $z \geq 0$일 때 $\phi(z) \geq 1$인데 $\alpha^{-1} z \geq 0$이기 때문이다.) 따라서 $E\phi(\alpha^{-1}Z) \leq \epsilon$임을 확인하면 $\mathbb{P}(Z \geq 0) \leq \epsilon$을 보장할 수 있다.
또한, $\alpha$를 적절히 조정하여 더욱 가까운 근사를 해낼 수도 있다. $\phi$가 nonnegative, nondecreasing convex function일 때, $\alpha E\phi\left(\frac{f(x,U)}{\alpha}\right) \leq \alpha \epsilon$에서 perspective function $(w, \alpha) \rightarrow \phi(\frac{w}{\alpha})$는 $w,\,\alpha$에 대해 jointly convex이므로, 우리는 $x$와 $\alpha$를 optimization variable로 놓고 풀 수 있다. 이 때 구해진 $\alpha$는
$\underset{\alpha \geq 0}{\inf}\left\{\alpha E\phi\left(\frac{f(x,U)}{\alpha}\right) - \alpha \epsilon \right\} \leq 0$, 즉 가장 가까운(tightest) convex approximation이다.

Tightest convex bounds and conditional value at risk

convex function $\phi(z) = [1+z]_+ = \max\left\{0, 1+z \right\}$로 설정하자. 그러면 $1(z \geq 0) \leq [1+z]_+$을 만족한다. 위에서 얻은 tightest convex approximation을 적용하면 $\underset{\alpha \geq 0}{\inf}\left\{\alpha E\left[\frac{f(x,U)}{\alpha} + 1\right]_+ - \alpha \epsilon \right\} = \underset{\alpha \geq 0}{\inf}\left\{E\left[f(x,U) + \alpha \right]_+ - \alpha \epsilon \right\} \leq 0$을 얻는다. 여기에서 $\alpha$대신 $-\alpha$를 대입하면 $E\left[f(x,U) - \alpha \right]_+ + \alpha \epsilon \leq 0$를 얻는다. 이것은 conditional value at risk의 표현식으로, 아래와 같다.

$\textbf{CVaR}(Z;\epsilon) = \underset{\alpha }{\inf}\left\{\frac{1}{\epsilon}E[Z-\alpha]_+ +\alpha \right\}$

Conditional value at risk는 $x$에 대한 convex function이므로 $\textbf{CVaR}(f(x,Z);\epsilon) \leq 0$은 convex constraint가 된다.

CVaR의 정의로부터 $\alpha$를 minimize out 해보자. 그러면 $\frac{\alpha}{\partial \alpha} \left\{\alpha + \frac{1}{\epsilon}E[Z-\alpha]_+ \right\} = 1 - \frac{1}{\epsilon}E1(Z\geq \alpha) = 1-\frac{1}{\epsilon}\mathbb{P}(Z\geq\alpha) = 0$이 되게 하기 위해, $\epsilon = \mathbb{P}(Z \geq\alpha^*)$가 되도록 $\alpha^*$를 정하면 아래의 관계식을 얻는다.

$\textbf{CVaR}(Z;\epsilon) = \frac{1}{\epsilon}E[Z-\alpha^*]_+ +\alpha^*$
$=\frac{1}{\epsilon}E[Z-\alpha^*] + \textbf{VaR}(Z;\epsilon)$

즉 conditional value at risk는 항상 value at risk보다 크거나 같고, 추가된 term은 $\alpha^*$로부터의 deviation에 대한 penelty term이다. 정확한 approximation은 필연적으로 높은 deviation을 동반할 수밖에 없기 때문이다. 따라서 CVaR을 minimize하는 것은 $Z$의 upper $\epsilon$ quantile과 deviation penelty를 동시에 minimize하는 것과 같고, VaR에 비해 보수적인 결과를 낸다.

또한 CVaR은 VaR 밖의 tail에 대한 expection으로 해석할 수 있다. 왜냐하면 $E[Z\,|\,Z \geq \alpha^*] = E[\alpha^* + (Z-\alpha^*)\,|\,Z \geq \alpha^*] = \alpha^* + \frac{E[Z-\alpha^*]_+}{\mathbb{P}(Z \geq \alpha^*)} = \alpha^* + \frac{1}{\epsilon}E[Z-\alpha^*]_+ = \textbf{CVaR}(Z;\epsilon)$이기 때문이다.

CVaR은 convex constraint를 도출하지만, $E[f(x,U)-\alpha^*]_+$를 analytic하게 계산할 수 없는 경우가 많아 보통 monte carlo simulation을 이용해 계산하게 된다. 따라서 계산 비용이 너무 많이 들 수 있어, 규모가 큰 문제에서는 다른 방법을 활용해야 한다.

Analytic approximation using moment generating functions

계산 비용을 줄이기 위한 수단으로, mgf를 이용해 constraint를 analytic하게 도출하는 방법에 대해 살펴보자. Safe convex approximation으로 다시 돌아가, $\phi(z) = e^z$로 설정하자. 이 경우에도 $1(z \geq 0) \leq e^z$가 만족된다. 그러면 Chernoff bound로부터 모든 $\lambda \geq 0$에 대해 $\mathbb{P}(Z \geq t) \leq E\exp(\lambda Z - \lambda t)$가 성립한다. (Chernoff bound는 https://en.wikipedia.org/wiki/Chernoff_bound)를 참고하자. Markov inequality에 $u(z) = \exp(\lambda z),\,c=M(\lambda)$를 대입하면 얻어낼 수 있다.) 따라서 $E\exp(\alpha^{-1}f(x,U)) \leq \epsilon$은 $\mathbb{P}(Z > 0) \leq \epsilon$에 대한 safe approximation임을 알 수 있다.
만약 random function $f(x, U)$가 모든 $U$에서 $x$에 대한 convex function이면, $\text{log-sum-exp function}$인 $\log E \exp(f(x,U))$ 역시 $x$에 대해 convex임이 성립한다. 그 증명은 아래와 같다.
$\lambda \in [0,1]$에 대해
$\log E \exp(f(\lambda x + (1-\lambda)y,U)) \leq \log E\exp(\lambda f(x, U) + (1-\lambda)f(y, U)) = \log E\exp(f(x, U)^\lambda + f(y, U)^{(1-\lambda)}) \leq \log\left[(E\exp(f(x, U)))^\lambda (E\exp(f(y, U)))^{(1-\lambda)} \right]= \lambda \log E e^{f(x, U)} + (1-\lambda)\log E e^{f(y, U)}$
3번째에서 4번째 줄은 Hölder's inequality로, $\frac{1}{p} + \frac{1}{q} = 1$일 때 $a^Tb \leq ||a||_p||b||_q$가 성립한다.

따라서 $\log E\exp(\alpha^{-1}f(x,U)) \leq \log \epsilon$은 safe convex approximation이다. 아까의 perspective function 개념을 적용해 $\alpha$에 대해 minimization하면 아래와 같은 conservative approximiation을 얻을 수 있다.

$\underset{\alpha \geq 0}{\inf}\left\{\alpha \log E\exp\left(\frac{f(x,U)}{\alpha}\right) + \alpha\log \frac{1}{\epsilon} \right\} \leq 0$

만약 moment generating function $E \exp (f(x,U))$를 계산하기 어렵다면 대체 수단을 찾아야 하는데, 좋은 방법은 적절한 mgf의 적절한 convex upper bound $E \exp (f(x,U)) \leq \psi(x)$를 찾아 대체하는 것이다. (less conservative하지만 여전히 좋은 approximation일 수 있다.)

예시로 random variable이 sub-Gaussian일 때, 그리고 $f(x,U) = U^Tx$로 linear 형태일 때를 살펴보자. Sub-gaussian 확률 분포는 정규분포에 비해 평균 근처에 집중되어 있고, variability가 제한되어 있어 극단적인 값이 나올 확률이 극히 적은 분포를 말한다. Sub-gaussian 확률 분포를 따르는 평균이 0인 independent random vector $U=(U_1,\,U_2,\,...,\,U_n)$에 대하여, $\lambda \in \mathbb{R}$과 $\sigma_i \geq 0$에 대해 $E \exp(\lambda U_i) \leq \exp\left(\frac{\lambda^2 \sigma_i^2}{2} \right)$이 성립한다. 또한 이 부등식은 gaussian일 경우$(U_i \sim N(0, \sigma_i^2))$, bounded random variable $(U_i \in [-\sigma_i, \sigma_i]\,\, a.s.)$인 경우에도 성립한다. (이 경우 이전의 Hoeffding's inequality와 동일하다.)
따라서 $f(x, t;U) = U^Tx -t$라 할 때 $\log E \exp(\frac{U^Tx}{\alpha} - \frac{t}{\alpha}) \leq \sum_{i=1}^n \frac{x_i^2\sigma_i^2}{2\alpha^2} - \frac{t}{\alpha}$이다. 따라서 위의 safe convex approximation은 아래의 조건으로 대체할 수 있다.

$\underset{\alpha \geq 0}{\inf}\left\{\frac{1}{2\alpha}\sum_{i=1}^n x_i^2 \sigma_i^2 -t + \alpha \log \frac{1}{\epsilon} \right\} = \sqrt{2\log \frac{1}{\epsilon}}||\text{diag}(\sigma)x||_2 -t \leq 0$

Example: normal random variable

위의 normal random variable의 예시를 다시 살펴보자. $U \in \mathbb{R}^n$이 평균이 $\mu$, covariance matrix가 $\sigma^2I$인 정규분포를 따르고 $f(x,U)$가 $x$와 $U$에 대한 쌍선형 형식(bilinear form)일 때,
$\mathbb{P}(U^Tx > 0) \leq \epsilon \iff \textbf{VaR}(U^Tx;\epsilon) \leq 0 \iff 0 \geq \mu^T x + \sigma \Phi^{-1}(1-\epsilon)||x||_2$ 라고 했다. 위와 달리 mgf approximation을 이용해 구한 constraint는
$\mu^Tx + \sqrt{2\log \frac{1}{\epsilon}}||\text{diag}(\sigma)x||_2 \leq 0$이다. 밑의 그림으로부터 mgf approximation이 VaR에 비해 더 보수적임을 알 수 있다.

Example: analytic versus monte carlo approximations

이전의 포트폴리오 최적화를 다시 살펴보자. 이전에는 자산의 평균 수익률을 $\mu_i = 1.05 + 0.3/i$로, 수익률을 $R_i \in [\mu_i - 0.6/i,\,\mu_i + 0.6/i]$의 bounded random variable로 설정한 뒤 robust optimization을 통해 결과를 확인해 보았다. 사실 두 번째 objective인 $\mu^Tx - \sqrt{2\log \frac{1}{\epsilon}}||\text{diag}(u)x||_2$는 lower bound 관점에서의 $\textbf{VaR}(R^Tx;\epsilon)$과 같다. 이번에는 convex constraint의 upper bound 관점에서 CVaR과 mgf approximation을 비교해 보자. 이 최적화 문제를 수식으로 나타내면

$\text{maximize}\,\,t$
$\text{subject to}\,\,\mathbf{1}^Tx = 1,\,x\succeq 0,$
$\underset{\alpha \geq 0}{\inf}\,E[t-R^Tx + \alpha]_+ - \alpha \epsilon \leq 0$
$\text{OR:}\,\,t-\mu^Tx + \sqrt{\frac{1}{2}\log \frac{1}{\epsilon}}||\text{diag}(u)x||_2 \leq 0$

이다. 다음 코드를 통해 결과를 확인할 수 있다.

import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp
%config InlineBackend.figure_format='retina'

np.random.seed(1)

n = 15  # Dimensions
m = 2 * 10**3  # Monte carlo simulations
mu = 1.05 + 0.3 / np.arange(1, n+1)  # Means for random variables
epsilon = 0.1  # Probability of failing to satisfy constraint
ranges = 0.6 / np.arange(1, n+1)

# Construct objective for analytic and MC problem
t_analytic = cp.Variable(1)
x_analytic = cp.Variable(n)
obj_analytic = cp.Maximize(t_analytic)
constraints_analytic = [
    t_analytic - cp.sum(mu @ x_analytic) +
    cp.sqrt(np.log(1 / epsilon) / 2) * cp.norm(2 * cp.multiply(ranges, x_analytic), 'fro') <= 0,
    x_analytic >= 0,
    cp.sum(x_analytic) == 1
]
prob_analytic = cp.Problem(obj_analytic, constraints_analytic)

# Random problem
t_mc = cp.Variable(1)
x_mc = cp.Variable(n)
alpha_mc = cp.Variable(1)
obj_mc = cp.Maximize(t_mc)
Uvals = mu[:, np.newaxis] + ranges[:, np.newaxis] * np.sign(1 - 2 * np.random.rand(n, m))
constraints_mc = [
    (1 / m) * cp.sum(cp.maximum(t_mc - Uvals.T @ x_mc + alpha_mc, 0)) - alpha_mc * epsilon <= 0,
    x_mc >= 0,
    cp.sum(x_mc) == 1
]
prob_mc = cp.Problem(obj_mc, constraints_mc)

prob_analytic.solve(solver=cp.ECOS, max_iters=1000)
prob_mc.solve(solver=cp.ECOS, max_iters=1000)

# Monte carlo simulation
np.random.seed(5)

m = 5 * 10**4
Rs = mu[:, np.newaxis] + ranges[:, np.newaxis] * np.sign(1 - 2 * np.random.rand(n, m))

ta = t_analytic.value
values_analytic = Rs.T @ x_analytic.value
num_analytic_violations = np.sum(ta > values_analytic)

tmc = t_mc.value
values_mc = Rs.T @ x_mc.value
num_mc_violations = np.sum(tmc > values_mc)

plt.figure()

plt.hist([values_analytic, values_mc], bins=15, color=([0, 0.5, 1], [1, 0.5, 0]), label=("analytic", "monte carlo"))

ytick_vals = plt.yticks()[0]
proportions = ytick_vals / m
plt.yticks(ytick_vals, proportions, fontsize=14)

miny, maxy = plt.axis()[2:4]
plt.plot([ta, ta], [miny, maxy], linewidth=2, color=[0, 0.25, 0.5])
plt.plot([tmc, tmc], [miny, maxy], linewidth=2, color=[0.5, 0.25, 0])
plt.legend()

print("Fraction of violations under analytic: ", num_analytic_violations / m)
print("Fraction of violations under Monte Carlo: ", num_mc_violations / m)

print("Mean return under analytic: ", np.sum(mu @ x_analytic.value))
print("Mean return under Monte Carlo: ", np.sum(mu @ x_mc.value))

수직선 왼쪽은 chance constraint를 만족하지 못한 경우. Analytic approximation이 CVaR보다 좀 더 보수적인 결과를 냄을 확인할 수 있다.