[Week 4] Robust Optimization [2]

ESC·2023년 12월 13일

2023-Spring

목록 보기

7/9

Robust SOCPs

이번에는 원래의 constraint 자체가 second-order cone 형태로 주어진 robust SOCP에 대해 알아보자. SOCP의 목적함수는 LP와 동일하며, constraint는 일반적으로 다음과 같이 주어진다. $A$ 의 각 행벡터를 $a_i\in \mathbb{R}^n$ 이라 하자.

$||Ax + b||_2 \leq c^Tx + d,$
$\text{where}\,\,A\in \mathbb{R}^{m \times n},\,b\in \mathbb{R}^{m},\,c\in \mathbb{R}^{n},\,d\in \mathbb{R}$

SOCPs with interval uncertainty

가장 간단하게 생각할 수 있는 uncertainty set은 행렬 $A$ 에 원소 크기가 $\delta$ 로 제한된 $\Delta$ 가 더해진 형태이다. 이를 수식으로 나타내면 아래와 같다. 이전에 LP에 conic uncertainty를 더하면 SOCP로 변환됨을 확인했기 때문에, 우변에 다음과 같이 자유롭게 더할 수 있다.

$||(A+\Delta)x + b||_2 \leq (c+u)^Tx + d$
$\text{for all}\,\,\Delta\,\text{s.t.}\,||\Delta||_{\infty} \leq \delta,$
$u \in \mathcal{U}=\left\{u:Fu +g \succeq_K 0\right\}$

이때 좌변과 우변을 $t \geq 0$ 에 대해 $||(A+\Delta)x + b||_2 \leq t$ 와 $t \leq (c+u)^Tx + d$ 로 분리해 생각하자. 그러면 우변은 이전과 동일하게 $t \leq c^Tx + d - \lambda^Tg,\,\,F^T\lambda + x = 0,\,\,\lambda \succeq_{K^*} 0$ 의 세 부등식으로 나타난다. 좌변은 계산을 통해 확인해 보자.
$\sup ||(A+\Delta)x + b||_2 = \sup\left(\sum_{i=1}^m [(a_i + \Delta_i)^T x + b_i]^2 \right)^{1/2} = \sup \left\{||z||_2\,|\,z_i = a_i^T x + \Delta_i^T x + b_i,\,||\Delta_i||_{\infty} \leq \delta \right\} = \inf \left\{||z||_2\,|\,z_i \geq |a_i^T x + b_i| + \delta||x||_1 \right\}$

따라서 처음의 무한 가지 경우의 수가 있는 부등식은, 아래의 $m+3$ 개의 부등식으로 축약된다.

$||z||_2 \leq c^Tx + d - \lambda^Tg,\,\,F^T\lambda + x = 0,\,\,\lambda \succeq_{K^*} 0$
$z_i \geq |a_i^T x + b_i| + \delta||x||_1\,\,\text{for}\,\,i=1,2,...,m$

SOCPs with ellipsoidal uncertainty

Interval uncertainty를 약간 변형한 형태로, 이전과 유사하게 constraint를 도출할 수 있다. Uncertainty를 반영한 inequality는 다음과 같다.

$\left(\sum_{i=1}^m [(a_i + P_iu)^T x + b_i]^2 \right)^{1/2} \leq t$
$\text{for all}\,\,||u||_2 \leq 1$

변수 $z$ 를 도입해 $z_i \geq \sup_u |a_i^Tx + b_i + u^TP_i^Tx|$ 으로 쓰면 이는 $z_i \geq |a_i^Tx + b_i| + ||P_i^Tx||_2$ 와 동치이고, 따라서 처음의 부등식은 아래와 같이 축약된다.

$||z||_2 \leq t,\,\,z_i \geq |a_i^Tx + b_i| + ||P_i^Tx||_2$
$\text{for}\,\,i=1,2,...,m$

SOCPs with matrix uncertainty

이번엔 $\Delta$ 의 l2-induced norm으로 주어진 형태의 uncertainty set에 대해 살펴보자. l2-induced norm은 $||\Delta|| = \sup_{||u||_2=||v||_2=1} u^T\Delta v$ 로 정의되며, $A$ 의 maximum singular value와 같다. 이런 경우에는 좀 더 복잡한 증명이 필요하다. $P \in \mathbb{R}^{m \times n}$ , $\Delta \in \mathbb{R}^{n \times n}$ 에 대해 matrix uncertainty를 반영한 inequality는 아래와 같다.

$||(A + P\Delta)x + b ||_2 \leq t,\,\,||\Delta|| \leq \delta$

이제 다음 두 사실을 이용해 이 부등식을 풀 수 있는 형태로 만들어 보자. 증명은 생략하겠다.

1. Schur complement representation
$||x||_2 \leq t$ 와 $\begin{bmatrix} t & x^T \\ x & tI_n \end{bmatrix} \succeq 0$ 은 동치이다.

2. The(homogeneous) S-lemma
$x^TAx \geq 0 \rightarrow x^TBx \geq 0\,\,\text{iff}\,\,\exists \lambda \geq 0\,\,\text{s.t.}\,\,B \succeq \lambda A$

Schur complement를 이용해 이 부등식을 $\begin{bmatrix} t & ((A + P\Delta)x + b)^T \\ (A + P\Delta)x + b & tI_m \end{bmatrix} \succeq 0,\,||\Delta|| \leq 1$
로 바꿀 수 있다. 이는 임의의 스칼라 $s\in \mathbb{R}$ 과 벡터 $v \in \mathbb{R}^m$ 에 대해
$ts^2 + 2s((A + P\Delta)x + b)^Tv + t||v||_2^2 \geq 0$ 임과 동치이다.
이 식에서 $||\Delta|| \leq 1$ 에 대해 infimum을 취하면 결국
$ts^2 + 2s(Ax + b)^Tv + t||v||_2^2 - 2||sx||_2||P^T v||_2 \geq 0$ 과 같다.
이때 $\sup \left\{u^Tx\,|\,||u||_2^2 \leq ||P^T v||_2^2 \right\} = ||P^T v||_2||x||_2$ 임을 이용해, 이 부등식을
$ts^2 + 2s(Ax + b)^Tv + t||v||_2^2 + 2su^Tx \geq 0$
$\text{for all}\,\,||u||_2^2 \leq ||P^T v||_2^2$
로 나타낸 뒤, 행렬 형태로 나타내면

$\begin{bmatrix} s \\ v \\ u \end{bmatrix}^T \begin{bmatrix} 0 & 0 & 0 \\ 0 & PP^T & 0 \\ 0 & 0 & -I_n \end{bmatrix}\begin{bmatrix} s \\ v \\ u \end{bmatrix} \succeq 0 \rightarrow \begin{bmatrix} s \\ v \\ u \end{bmatrix}^T \begin{bmatrix} t & (Ax+b)^T & x^T \\ Ax+b & t & 0 \\ x & 0 & 0 \end{bmatrix}\begin{bmatrix} s \\ v \\ u \end{bmatrix} \succeq 0$
그리고 S-lemma를 이용하면, 최종적으로 matrix uncertainty inequality는 다음을 만족하는 $x\in\mathbb{R}^n,\,\lambda$ 를 찾는 문제와 같다.

$\begin{bmatrix} t & (Ax+b)^T & x^T \\ Ax+b & t - \lambda PP^T & 0 \\ x & 0 & \lambda I_n \end{bmatrix} \succeq 0$

Example: robust regression

예시를 통해 위에서 살펴본 이론을 적용해 보자. Regression problem은 우리가 잘 알고 있듯, 주어진 데이터를 잘 설명하는 직선의 계수를 찾는 문제로 데이터 행렬 $A \in \mathbb{R}^{m \times n},\,b \in \mathbb{R}^{m}$ 가 주어져 있을 때 $||Ax - b||_2$ 를 최소화하는 $x \in \mathbb{R}^{n}$ 을 찾는 것이 목표이다. 이때 데이터 자체에 오염이 발생했다고 하고, 그 오염은 $\Delta \in \mathbb{R}^{m \times n}$ 이 $A$ 에 더해진 형태로 나타난다. 이때 $\Delta$ 의 범위를 위의 세 경우와 같이 분류해 보자. 즉 $|\Delta_{ij}|$ 또는 $||\Delta_i||_2$ 또는 $||\Delta||$ 가 높은 확률로 제한된 경우이다. $\Delta$ 를 어떤 확률 분포를 따르는 random noise라고 하면 이 가정은 현실적임을 알 수 있다. 가령 gaussian noise의 경우 noise의 범위가 평균(0)에서 그리 멀지 않은 부분에 집중될 것이기 때문이다.
다음 theorem을 이용해 $\Delta$ 에 관한 probability bound를 도출해 보자.

Theorem.
$t \geq 0$ 에 대해
1. $\mathbb{P}(|\Delta_{ij}| \geq t) \leq 2\exp(-\frac{t^2}{2})$
2. $\mathbb{P}(||\Delta_{i}||_2 \geq \sqrt{n} + t) \leq \exp(-\frac{t^2}{2})$
3. $\mathbb{P}(||\Delta_{i}|| \geq \sqrt{m} + \sqrt{n} + t) \leq \exp(-\frac{t^2}{2})$

이를 이용해 각각 $t_{\infty}^2(\delta) = 2\log(2mn/\delta),\,t_{2}^2(\delta) = 2\log(m/\delta),\, \,t_{op}^2(\delta) = 2\log \frac{1}{\delta}$ 로 설정하면 위의 upper bound가 $\delta$ 가 되도록 조정할 수 있다. 따라서 uncertainty set의 bound를 각각 위의 세 값으로 설정하면 각각 다른 robust least squares solution을 얻을 수 있다.

mport cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(1)

# Generate data
m = 50
n = 20
A = 10 * np.random.randn(m, n)
u = np.random.randn(n)
u = u / np.linalg.norm(u)
b = A @ u + np.random.randn(m)

# Probabilities of exceeding the bounds on robustness
num_deltas = 23
deltas = np.logspace(-11, -1, num_deltas)

x_inf_sols = np.zeros((n, num_deltas))
x_2_sols = np.zeros((n, num_deltas))
x_op_sols = np.zeros((n, num_deltas))

nom_vals_inf = np.zeros(num_deltas)
nom_vals_2 = np.zeros(num_deltas)
nom_vals_op = np.zeros(num_deltas)

for ind, delta in enumerate(deltas):
    t_op = np.sqrt(n) + np.sqrt(m) + np.sqrt(2 * np.log(1 / delta))
    t_2 = np.sqrt(n) + np.sqrt(2 * np.log(m / delta))
    t_inf = np.sqrt(2 * np.log(2 * m * n / delta))
    
    print(f"Solving iteration {ind+1} of {num_deltas} with delta = {delta}")
    
    # First with infinity norm-bounded perturbations
    x_inf = cp.Variable(n)
    z_inf = cp.Variable(m)
    objective_inf = cp.Minimize(cp.norm(z_inf, 2))
    constraints_inf = [z_inf[i] >= cp.abs(A[i, :] @ x_inf - b[i]) + t_inf * cp.norm(x_inf, 1) for i in range(m)]
    problem_inf = cp.Problem(objective_inf, constraints_inf)
    problem_inf.solve(solver=cp.ECOS, max_iters=1000, verbose=False)
    x_inf_sols[:, ind] = x_inf.value
    
    # Now with ellipsoid (l2) perturbations
    x_2 = cp.Variable(n)
    z_2 = cp.Variable(m)
    objective_2 = cp.Minimize(cp.norm(z_2, 2))
    constraints_2 = [z_2[i] >= cp.abs(A[i, :] @ x_2 - b[i]) + t_2 * cp.norm(x_2, 2) for i in range(m)]
    problem_2 = cp.Problem(objective_2, constraints_2)
    problem_2.solve(solver=cp.ECOS, max_iters=1000, verbose=False)
    x_2_sols[:, ind] = x_2.value
    
    # Now with the operator norm perturbation
    x_op = cp.Variable(n)
    objective_op = cp.Minimize(cp.norm(A @ x_op - b, 2) + t_op * cp.norm(x_op, 2))
    problem_op = cp.Problem(objective_op)
    problem_op.solve(solver=cp.ECOS, max_iters=1000, verbose=False)
    x_op_sols[:, ind] = x_op.value

for ind, delta in enumerate(deltas):
    nom_vals_inf[ind] = np.linalg.norm(A @ x_inf_sols[:, ind] - b, 2)
    nom_vals_2[ind] = np.linalg.norm(A @ x_2_sols[:, ind] - b, 2)
    nom_vals_op[ind] = np.linalg.norm(A @ x_op_sols[:, ind] - b, 2)

x_nominal = np.linalg.lstsq(A, b, rcond=None)[0]
nominal_value = np.linalg.norm(A @ x_nominal - b)

# ... Additional code for Monte Carlo simulations and plotting ...

위의 코드를 통해 각각의 상황을 시뮬레이션할 수 있다. 결과는 아래와 같다.