分布的ロバスト最適化をロジスティク回帰に適用しcvxpyで実装

# はじめに
分布的ロバスト最適化をロジスティク回帰に適用した論文がNIPS2015に提案されました。
cvxpyで実装してみたので備忘録として残します。

分布的ロバスト最適化については、みるかさんの記事が参考になります。
最適輸送理論から導入して、分布的ロバスト最適を紹介されているのは、とても流れがわかりやすく素晴らしいなぁと思います。

https://mirucacule.hatenablog.com/entry/2022/12/24/000000

# ロジスティク回帰のモデル
標本の集合$\lbrace(x_1,y_1)\cdots(x_n,y_n)\rbrace$が与えられているとする。標本$i$の属性を$x\in\mathbb{R}^d$、標本$i$のラベルを$y_i\in\lbrace\pm1\rbrace$とする。このとき、$x$が与えられたときの$y$の確率は、

$$
\mathrm{Prob}(y|x)=[1+\exp(-y\langle\beta,x\rangle)]^{-1}.
$$

ここで、$\beta\in\mathbb{R}^d$は重みベクトルである。$\mathrm{Prob}(+1|x)+\mathrm{Prob}(-1|x)=1$となっているので確かに確率になっている。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/139678/912b7c3b-2f5f-c842-f864-8bbc822ae91d.png)

# 分布的ロバスト最適化による定式化
ロジスティク回帰の分布的ロバスト最適化による最適化問題は下式で表されれる。

```math
\begin{eqnarray}
\inf_{\lambda,s_i} &  &\lambda\epsilon+\dfrac{1}{N}\sum_{i=1}^N s_i\\
s.t. & &
l_{\beta}(\hat{x}_i,\hat{y}_i)\leq s_i \ \ (\forall i \leq N)\\
     & &
l_{\beta}(\hat{x}_i,-\hat{y}_i)-\lambda\kappa\leq s_i \ \ (\forall i \leq N)\\
     & &
\|\beta\|_{*}\leq\lambda     
\end{eqnarray}
```

ただし、$l_{\beta}(x,y)=\log(1+\exp(-y,\langle\beta ,x\rangle))$である。

# cvxpyによる実装

```python
import numpy as np
import cvxpy as cp
from sklearn.metrics import confusion_matrix

class dr_logistic_regression():
  def __init__(self,epsilon,kappa,pnorm):
    self.epsilon = epsilon
    self.kappa   = kappa
    self.pnorm   = pnorm

  def fit(self,x,y):
    self.N = x.shape[0]
    self.M = x.shape[1]
    self.x = x
    self.y = y
    beta_   = cp.Variable((self.M,1))
    lambda_ = cp.Variable((1,1))
    s       = cp.Variable((self.N,1))

    constraints = [
                 cp.logistic( cp.multiply(-self.y,self.x@beta_) ) <= s 
                ,cp.logistic( cp.multiply(+self.y,self.x@beta_) ) - lambda_*self.kappa <= s
                ,cp.norm(beta_, self.pnorm)<=lambda_
                  ]

    objective = lambda_*self.epsilon + 1/self.N*cp.sum(s)

    problem = cp.Problem(cp.Minimize(objective),constraints)
    problem.solve(verbose=False)
    self.beta_ = beta_.value
    return problem.value,beta_.value,lambda_.value,s.value
    
  def infer(self,x):
      y_est = []
      for i in 1/(1+np.exp(-x@self.beta_)):
        if i>=0.5:
          y_est.append(1)
        if i<0.5:
          y_est.append(-1)
      self.y_est = y_est
      return y_est
    
  @staticmethod
  def confusion_matrix(y,y_est):
    C = confusion_matrix(y, y_est)
    return C
```

コードは下記に置きました。よかったら使ってみてください。

https://github.com/Hideki105/robust-optimization


# 参照
- [Distributionally Robust Logistic Regression, Soroosh Shafieezadeh-Abadeh and Peyman Mohajerin Esfahani, NIPS2015](https://arxiv.org/pdf/1509.09259.pdf)
- [Daniel Kuhn: Data-driven and Distributionally Robust Optimization and Applications -- Part 1/2](https://www.youtube.com/watch?v=b4lJENGAeEA&list=PLKLR7D59yU0fgC_Ac6v0I0yGO8FcE98GH&index=7&t=550s)
- [Daniel Kuhn: Data-driven and Distributionally Robust Optimization and Applications -- Part 2/2](https://www.youtube.com/watch?v=M2NcdKwn824)
- [scikit-learn でクラス分類結果を評価する](https://pythondatascience.plavox.info/scikit-learn/%e5%88%86%e9%a1%9e%e7%b5%90%e6%9e%9c%e3%81%ae%e3%83%a2%e3%83%87%e3%83%ab%e8%a9%95%e4%be%a1)