最適輸送距離に基づく分布的ロバスト最適化とその周辺

このエントリは「数理最適化 Advent Calendar 2022」の 24 日目の記事です．わたしの前後の記事は：

23 日目は @YamagenSakam さんによる『IIRフィルタの設計問題を焼きなまし法で解いてみる』
25 日目は @snowberryfield さんによる『整数計画問題のメタヒューリスティクス向け前処理手法について書きます』

です．Advent Calendar は 2020 年の「数理最適化 Advent Calendar 2020」以来 2 年ぶり 2 回目です．

前回の Advent Calendar では共役勾配法と呼ばれる最適化手法について扱いました．今回は，近年盛んに研究されている分布的ロバスト最適化というモデリング手法を紹介したいと思います．本エントリは，わたしが今年読んだ論文の中で特に印象的だった [Shafieezadeh-Abadeh, S., Kuhn, D., & Esfahani, P. M., 2019] が基になっています．本エントリを読んで興味を持った方は是非ご参照ください．それでは，やっていきましょう．

はじめに
機械学習における不確実性
不確実性を考慮した最適化手法
準備
最適輸送距離に基づく分布的ロバスト最適化
min-max 問題の凸計画問題への帰着
- 等価性の証明の概略
他の最適化手法との関係性
- ロバスト最適化問題との関係性
- 正則化手法との関係性
数値実験
おわりに
参考文献
補遺

はじめに

様々な現実のタスクに対して機械学習／深層学習モデルが人間に匹敵，あるいは，凌駕するような性能を発揮することがごく当たり前の世の中になっています．本エントリを執筆時点 (2022 年 12 月) では，OpenAI の発表した対話型 AI である ChatGPT が非常に大きな話題を集めています．試しに，本エントリの内容に即した質問として

What are some of the important points to keep in mind when using machine learning models and deep learning models for various real-world tasks?

を入力したときの ChatGPT の回答を補遺に掲載しておきます．

機械学習／深層学習の性能向上の一方で，機械学習／深層学習モデルは敵対的攻撃に対して脆弱であることが度々報告されており，性能の良いモデルであっても現実のタスクにそのまま適用する上で注意を要する場合があります．本エントリでは，数理最適化の立場から不確実性に対処するモデリング手法，特に近年盛んに研究されている分布的ロバスト最適化法について扱います．

記法．以下では，bold 体の文字 (e.g., $\boldsymbol{x}$) を適当な次元のベクトルを表します．また，$\|\cdot\|$ はユークリッド空間に定義される任意のノルム (e.g., $\ell_2$ ノルム，無限大ノルム) を表し，$\|\cdot\|_\ast$ でその双対ノルムを表します．双対ノルムは，ノルム $\|\boldsymbol{x}\|$ に対して $\|\boldsymbol{x}\|_{\ast}:= \sup \{ \left< \boldsymbol{x}, \boldsymbol{x}' \right> \mid \|\boldsymbol{x}'\| \leqslant 1 \}$ で定義されます．本エントリを読む上では，無限大ノルムに対する双対ノルムが $\ell_1$ ノルムである (i.e., $\|\cdot\|_{\ast} = \|\cdot\|_{1}$) こと，$\ell_2$ ノルムに対する双対ノルムが $\ell_2$ ノルムであることを抑えておけば十分です．その他の記号は以下の表を参照してください．

記号	意味	補足
$\mathcal{X}$	入力空間	$\mathcal{X}\subseteq\mathbb{R}^d$
$\mathcal{Y}$	出力空間	$\mathcal{Y}\subseteq\mathbb{R}$ or $\mathcal{Y}=\{-1,+1\}$
$\ell$	誤差関数	$\ell: \mathbb{R}^d \times \mathcal{Y} \to \mathbb{R}_{\geqslant 0}$
$L$	一変数誤差関数	$L: \mathbb{R} \to \mathbb{R}_{\geqslant 0}$
$\mathrm{lip}(f)$	$f$ のリプシッツ定数	$\mathrm{lip}(f) := \sup_{\boldsymbol{x},\, \boldsymbol{x}'} \{ \|f(\boldsymbol{x}) - f(\boldsymbol{x}')\| / \\| \boldsymbol{x} - \boldsymbol{x}\ \\| \mid \boldsymbol{x} \neq \boldsymbol{x}' \}$
$\mathbb{H}$	仮説集合	$\mathbb{H} \subseteq \mathbb{R}^{\mathcal{X}}$
$h$	仮説	$h \in \mathbb{H}$
$\boldsymbol{1}_{\mathcal{X}}$	定義関数	$\boldsymbol{1}_{\mathcal{X}} (\boldsymbol{x}) = 1 $ if $\boldsymbol{x} \in \mathcal{X}$; $\boldsymbol{1}_{\mathcal{X}} (\boldsymbol{x})=0$ otherwise.
$\delta_{\mathcal{X}}$	標示関数	$\delta_{\mathcal{X}} (\boldsymbol{x}) = 0 $ if $\boldsymbol{x} \in \mathcal{X}$; $\delta_{\mathcal{X}} (\boldsymbol{x})=+\infty$ otherwise.
$\hat{\mathbb{P}}_N$	経験分布	$\hat{\mathbb{P}}_N = \frac{1}{N}\sum_{n=1}^{N} \delta_{\boldsymbol{\xi}_n}$

本エントリでは Python 言語による実装を提供しており，以下のライブラリを読み込んでいることを前提とします：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# https://scipy.org/
from scipy.stats import wasserstein_distance
from scipy.spatial.distance import cdist
# https://scikit-learn.org/
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_breast_cancer
from sklearn.metrics import accuracy_score
# https://www.cvxpy.org/
import cvxpy as cp
# https://pythonot.github.io/
import ot
import ot.plot

また，本エントリで扱った数値計算については，こちらの Google Colaboratory に集約してありますので，手元で確認されたい方や本エントリで扱っていないパラメータでの実験を確認したい方はご参照ください．

機械学習における不確実性

機械学習モデル・深層学習モデルは様々な分野・多数の諸問題に対して高い性能を達成している一方，現実のタスクに対して直接適用する上で問題となる場合があります．これは，(ニューラルネットワークを含む) 機械学習モデルが敵対的な摂動に対して感度が高く脆弱性が高いことが要因として挙げられます．

例えば，[1] では医療分野における機械学習システムの脆弱性について取り上げています．具体的には，腫瘍の画像を入力としてその腫瘍が悪性か良性かを判定するモデルを取り上げ，高い確信度で良性と判定する画像に敵対的ノイズを加えることによって逆に高い確信度で悪性と判定してしまう例が紹介されています．また，敵対的攻撃に如何に対処すべきかについて論じています．他にも，

物体検出モデル (e.g., YOLO) に対して，人物にパッチを適用させることによって人物の検出を妨げる事例 [2]
道路標識の画像識別モデルに対して，画像に手を加えることによって標識の識別結果を変える事例 (e.g., 法定速度を 35 mph → 85 mph に誤認識させる) [3]
Speech-to-Text モデルで動作する音声認識 AI に対して微小なノイズを加えることで悪意のある命令 (e.g., 鍵の解錠) をさせる事例 [4]

など，多くの現実のタスクにおいて機械学習モデルが敵対的攻撃に脆弱であることを示す事例が報告されています．

機械学習モデルを通じて人々が意思決定を行う上でこのような不確実性を正しく扱うことは必要不可欠です．この不確実性を数理最適化の立場から扱うことを考えたいというのが本エントリのモチベーションの 1 つです．

不確実性を考慮した最適化手法

不確実性を考慮した最適化手法として，確率計画法，ロバスト最適化法，分布的ロバスト最適化法などがあります．それぞれの典型的な定式化は以下の表のようにまとめられます．本エントリではそれぞれの手法について詳ししく扱うことはしませんが，簡単にだけ触れておきます．

用語	定式化
確率計画法	$\displaystyle\min_{\boldsymbol{x}} \mathbb{E}_{\boldsymbol{\xi}\sim\mathbb{P}}[f(\boldsymbol{x}, \boldsymbol{\xi})]$
ロバスト最適化法	$\displaystyle\min_{\boldsymbol{x}} \max_{\boldsymbol{\xi}\in\Delta} f(\boldsymbol{x}, \boldsymbol{\xi})$
分布的ロバスト最適化法	$\displaystyle\min_{\boldsymbol{x}} \max_{\mathbb{Q}\in\mathcal{U}} \mathbb{E}_{\boldsymbol{\xi}\sim\mathbb{Q}}[f(\boldsymbol{x}, \boldsymbol{\xi})]$

確率計画法 (Stochastic Programming)

確率計画法では最適化すべき目的関数や制約条件を記述する関数に確率変数が含まれ，その確率変数が従う確率分布 (上記表における $\mathbb{P}$) が既知である (あるいは仮定した) 問題を扱います．

ロバスト最適化法 (Robust Optimization)

ロバスト最適化法では不確実性を有するデータの取り得る値の範囲 (上記表における $\Delta$) をあらかじめ設定し，その中の最悪ケースでのモデル化を行います．確率計画法と異なり，データに関する確率分布を考慮しない点が特徴です．

分布的ロバスト最適化法 (Distributionally Robust Optimization)

ロバスト最適化法では，不確実性を有するデータの取り得る最悪ケースでの最適化を行うため，過度に保守的な解が得られる可能性があります．確率計画法では，ロバスト最適化法では考慮しない不確実性を有するデータの確率分布を考慮することができる一方，サンプル数が十分に大きい場合でないと標本平均近似の精度が高くないため，本来解くべき最適化問題との乖離が大きくなる可能性があります．

この両者の良いところを取り入れたのが分布的ロバスト最適化法の立ち位置です．分布的ロバスト最適化法ではデータの従う分布に対する不確実性を考慮したモデリング手法で，最悪ケースの分布の下でのモデル化を行います．以下で詳しく見ていきましょう．

準備

分布的ロバスト最適化法について述べる準備をします．本エントリでは教師あり学習を例に説明するので問題の設定を明確にします．また，データの従う分布に対する不確実性を扱うため，2 つの分布が与えられたときにそれらの分布がどれくらい近いのかを測るための指標を導入します．

問題の設定

一般的な教師あり学習を考えます．すなわち，訓練データ集合 $\{(\hat{\boldsymbol{x}}_n, \hat{y}_n)\}_{n=1}^{N} \subset \mathcal{X} \times \mathcal{Y}$ が観測されたときに入出力関係を仮説 $h\in\mathbb{H}$ で説明することを考えます．ここで，$\mathbb{H}$ は仮説集合を表し，本エントリでは仮説集合 $\mathbb{H}$ を線形関数全体の空間，すなわち，

$$ \mathbb{H} = \left\{ \boldsymbol{x} \overset{h}{\mapsto} \left<\boldsymbol{w}, \boldsymbol{x} \right>\left| \, \boldsymbol{w} \in \mathbb{R}^d\right\}\right. $$

とします．

目標は，手元にある訓練データ集合 $\{(\hat{\boldsymbol{x}}_n, \hat{y}_n)\}_{n=1}^{N} \subset \mathcal{X} \times \mathcal{Y}$ を用いて記述される最適化問題：

$$ \min_{h\in\mathbb{H}} \mathbb{E}_{(\boldsymbol{x}, y)\sim\hat{\mathbb{P}}_N} [\ell (h(\boldsymbol{x}_n), y_n)] = \min_{\boldsymbol{w}} \frac{1}{N} \sum_{n=1}^{N} \ell (\left< \boldsymbol{w}, \hat{\boldsymbol{x}} \right>, \hat{y}_n) $$

を解くことによって未知のデータを予測する分類器 $h^\star \in \mathbb{H}$ を得ることです．上記の最適化問題を直接解くと訓練データに過適合してしまう問題があるため，目的関数に正則化項を加えた以下の問題を解くことがよく行われます：

$$ \min_{\boldsymbol{w}} \mathbb{E}_{(\boldsymbol{x}, y)\sim\hat{\mathbb{P}}_N} [\ell (\left<\boldsymbol{w}, \boldsymbol{x}\right>, y)] + c \,\Omega(\boldsymbol{w}). $$

ここで，$\Omega(\boldsymbol{w})$ は正則化項であり，$c > 0$ は正則化項の大きさを調整する正のパラメータです．なお，本エントリでは一変数の誤差関数 $L$ を導入します．例えば，回帰の問題であれば適当な凸関数 $L:\mathbb{R}\to\mathbb{R}_{\geqslant 0}$ に対して $\ell(h(\boldsymbol{x}), y)) = L(h(\boldsymbol{x}) - y)$ とします．回帰の問題であれば，同様に適当な凸関数 $L:\mathbb{R}\to\mathbb{R}_{\geqslant 0}$ に対して $\ell(h(\boldsymbol{x}), y)) = L(y h(\boldsymbol{x}))$ とします．凸関数 $L$ としては，区分線形な凸関数 $L(z)=\max_{j\in\{1,2,\dots,J\}} a_j z + b_j$ やヒンジロス関数 $L(z)=\max\{0, 1-z\}$ など多様な選択肢があります．

分布間の距離

次に，分布間の距離を測る Wasserstein 距離を導入します．分布間の距離の測り方には Kullback-Leibler divergence や f-divergence，Maximum Mean Discrepancy など様々な方法がありますが，本エントリでは最適輸送に基づく距離を考えます．最適輸送に関する日本語のドキュメントは非常に充実しており，ここで詳しく説明する必要はないと思われるため最小限の記述に留めます．詳しくは，[5, 6, 7] やそれらの中で参照されている文献 (例えば，[Peyré & Cuturi, 2019]) を参照してください．

本エントリでは，確率分布のサポートが有限である場合の最適輸送距が線形計画問題の最適値と一致することだけ確認します．分布 $\mathbb{Q}=(\mathbb{Q}_1,\dots,\mathbb{Q}_n)^\top$ と分布 $\mathbb{P} = (\mathbb{P}_1,\dots,\mathbb{P}_m)^\top$ の最適輸送距離は分布 $\mathbb{Q}$ を分布 $\mathbb{P}$ に輸送するための最小コストとして定義されます．すなわち，$\mathbb{Q}$ と $\mathbb{P}$ の最適輸送距離は各 $\mathbb{Q}_i$ を $\mathbb{P}_j$ にどれだけ輸送するのかを表す $\boldsymbol{\pi} = \{\pi_{ij}\}$ を決定変数とする以下の線形計画問題の最適値と一致します：

\begin{align} \begin{split} \mathcal{W}_d(\mathbb{Q},\mathbb{P}) = \left\{ \begin{array}{rl} \displaystyle \min_{\boldsymbol{\pi} \geqslant \boldsymbol{0}} \quad & \displaystyle \sum_{i=1}^{n} \sum_{j=1}^{m} d( \mathbb{Q}_i , \mathbb{P}_j ) \pi_{ij} \\ \text{subject to} \quad & \displaystyle \sum_{j=1}^{m} \pi_{ij} = \mathbb{Q}_i \quad (i = 1,2,\dots, n), \\ &\displaystyle \sum_{i=1}^{n} \pi_{ij} = \mathbb{P}_j \quad (j = 1,2,\dots, m). \end{array} \right. \end{split} \tag{*}\label{LP} \end{align}

ここで，$d(\cdot, \cdot)$ は 2 点間の距離を表します．また，上記の線形計画問題 (\ref{LP}) の最適解 $\boldsymbol{\pi}^\star \in \mathbb{R}^{n\times m}_{\geqslant 0}$ を輸送プランと呼ぶこととします．

本エントリでは特に用いることはありませんが，線形計画問題 (\ref{LP}) に対する双対問題は以下のように定式化されます：

\begin{align} \max_{\boldsymbol{\mu},\boldsymbol{\nu}} \quad & \sum_{i=1}^{n} \mathbb{Q}_i \mu_i + \sum_{j=1}^{m} \mathbb{P}_j \nu_j \\ \text{subject to} \quad & \mu_i + \nu_j \leqslant d( \mathbb{Q}_i , \mathbb{P}_j ) \quad (i =1,2,\dots,n,\, j = 1,2,\dots, m). \end{align}

ここで，$\boldsymbol{\mu}\in\mathbb{R}^{n}$ および $\boldsymbol{\nu}\in\mathbb{R}^{m}$ は等式制約に対するラグランジュ乗数を表します．双対問題の導出方法は例えば以前に書いた記事 (双対問題の作り方―LP の例―) を参照ください．

ここでは例として，最適輸送距離を求めるライブラリ POT で取り上げられている例をモデリング言語 CVXPY を用いて線形計画問題としてモデリングして，付属の Cbc ソルバーで解いてみましょう．

最適輸送距離の Python による求解

点群データを正規分布を用いて生成します．2 点間の距離はユークリッド距離を用いることとします．

n = 5 # number of source samples
m = 10 # number of target samples

np.random.seed(19830119)

mu_s = np.array([0, 0])
cov_s = np.array([[1, 0], [0, 1]])

mu_t = np.array([4, 4])
cov_t = np.array([[1, -.8], [-.8, 1]])

xs = np.random.multivariate_normal(mu_s, cov_s, n)
xt = np.random.multivariate_normal(mu_t, cov_t, m)

Q, P = np.ones((n,)) / n, np.ones((m,)) / m  # uniform distribution on samples

C = distance.cdist(xs, xt, metric='sqeuclidean')

生成した点群データを描画するとした図のようになります．左図が点群のデータ，右図は距離に基づくヒートマップを表します：

さて，青点 (＋) を赤点 (×) に “輸送” することを考えます．ここで，各点の重みは同一である，すなわち， $$ \mathbb{Q}=\left(\frac{1}{n},\dots,\frac{1}{n}\right)^\top,\,\mathbb{P}=\left(\frac{1}{m},\dots,\frac{1}{m}\right)^\top $$ とします．最適輸送距離 $\mathcal{W}_d(\mathbb{Q},\mathbb{P})$ およびそのときの輸送プラン $\boldsymbol{\pi}^{\star}$ を求める実装は例えば以下のように記述できます．

def OT(Q, P, C):
    Q = np.array(Q)
    P = np.array(P)
    n = Q.size
    m = P.size
    X = cp.Variable((n, m))
    obj = cp.sum(cp.multiply(C, X))
    con = [X >= 0]
    for i in range(n):
        con.append(sum([X[i][j] for j in range(m)]) == Q[i])
    for j in range(m):
        con.append(sum([X[i][j] for i in range(n)]) == P[j])
    prob = cp.Problem(cp.Minimize(obj), con)
    prob.solve()
    return (prob.value, X.value)

この OT メソッドを用いて，分布 $\mathbb{Q}$ と分布 $\mathbb{P}$ の最適輸送距離を求めてみましょう：

res_fun, res_x = OT(Q, P, C)
res_x_rounded = np.around(res_x, decimals=2) # round to 2
print('optimal value    = {}'.format(res_fun))
print('optimal solution = \n{}'.format(res_x_rounded))

# optimal value    = 41.150790627650416
# optimal solution = 
# [[0.  0.  0.1 0.  0.  0.  0.1 0.  0.  0. ]
#  [0.1 0.  0.  0.  0.  0.  0.  0.  0.1 0. ]
#  [0.  0.  0.  0.  0.1 0.  0.  0.1 0.  0. ]
#  [0.  0.1 0.  0.  0.  0.1 0.  0.  0.  0. ]
#  [0.  0.  0.  0.1 0.  0.  0.  0.  0.  0.1]]

最適輸送距離はおよそ 41.15 であることが分かりました．また，そのときの輸送プランが得られました．

POT ライブラリを用いる場合，最適輸送距離とそのときの輸送プランは以下のコマンドで得ることができます．

print(ot.emd(Q, P, C))
print(ot.emd2(Q, P, C))
# 41.15079054201769
# [[0.  0.  0.1 0.  0.  0.  0.1 0.  0.  0. ]
#  [0.1 0.  0.  0.  0.  0.  0.  0.  0.1 0. ]
#  [0.  0.  0.  0.  0.1 0.  0.  0.1 0.  0. ]
#  [0.  0.1 0.  0.  0.  0.1 0.  0.  0.  0. ]
#  [0.  0.  0.  0.1 0.  0.  0.  0.  0.  0.1]]

同じように，最適輸送距離はおよそ 41.15 であることが確認できます．

また，POT ライブラリを用いて輸送結果を可視化すると以下のようになります．

最適輸送距離に基づく分布的ロバスト最適化

分布間の距離を測る最適輸送距離が導入できたので，これを用いて分布に関する不確実性集合を定義します．いま，手元にあるデータ集合 $\{(\hat{\boldsymbol{x}}_n, \hat{y}_n)\}_{n=1}^{N} \subset \mathcal{X} \times \mathcal{Y}$ の経験分布を $\hat{\mathbb{P}}_N$ とすると，不確実性集合 $\mathbb{B}_{\rho}(\hat{\mathbb{P}}_N)$ は経験分布からの距離が最適輸送距離の意味で $\rho$ 以下である分布全体の分布として次のように定義されます：

$$ \mathbb{B}_{\rho}(\hat{\mathbb{P}}_N) = \left\{ \mathbb{Q} \left|\, \mathbb{Q}(\Xi) = 1, \, \mathcal{W}_d(\mathbb{Q}, \hat{\mathbb{P}}_N) \leqslant \rho \right\}\right.. $$

ここで，$\Xi$ は $\mathbb{Q}$ の台集合を表します．

この不確実性集合を用いて分布的ロバスト最適化問題は次のように定式化されます：

\begin{align} \min_{\boldsymbol{w}\in\mathbb{R}^d} \max_{\mathbb{Q}\in\mathbb{B}_{\rho}(\hat{\mathbb{P}}_N)} \mathbb{E}_{(\boldsymbol{x}, y)\sim\mathbb{Q}} [\ell (\left<\boldsymbol{w},\boldsymbol{x}\right>, y)]. \tag{1} \label{1} \end{align}

和訳すると，「経験分布 $\hat{\mathbb{P}}_N$ からの距離が最適輸送距離の意味で $\rho$ 以下のすべての分布の中で経験誤差 $\mathbb{E}_{(\boldsymbol{x}, y)\sim\mathbb{Q}} [\ell (\left<\boldsymbol{w},\boldsymbol{x}\right>, y)]$ を最も悪化させるような最悪ケースの分布 $\mathbb{Q}^\star$ に関して，そのときの経験誤差 $\mathbb{E}_{(\boldsymbol{x}, y)\sim\mathbb{Q}^\star} [\ell (\left<\boldsymbol{w},\boldsymbol{x}\right>, y)]$ を最小にするような線形関数の係数ベクトル $\boldsymbol{w}^\star \in \mathbb{R}^d$ を求める」問題です．

この最適化問題 (\ref{1}) は min-max 型の最適化問題であり一般に解くことが難しく，商用の汎用ソルバーでも扱うことができません．典型的な解決策としては，内側の最大化問題と外側の最小化問題を繰り返し解く方法がありますが，そのようにして得られる解は必ずしも最適性の条件を満たすとは限りません．

min-max 問題の凸計画問題への帰着

分布的ロバスト最適化問題を素直に定式化すると min-max 型の最適化問題となることを見ました．このままでは解くことが困難なのですが，この最適化問題を求解可能な形に再定式化する研究が近年進んできています．ここでは，$\mathcal{Y}=\{-1,+1\}$ であるような二値分類問題についての結果を紹介したいと思います．なお，回帰タスク (i.e., $\mathcal{Y}=\mathbb{R}$) の場合やより詳細については [8] を参照ください．誤差関数を $L$ としたとき，解くべき分布的ロバスト最適化問題は以下のように記述されます：

\begin{align} \min_{\boldsymbol{w}\in\mathbb{R}^d} \max_{\mathbb{Q}\in\mathbb{B}_{\rho}(\hat{\mathbb{P}}_N)} \mathbb{E}_{(\boldsymbol{x}, y)\sim\mathbb{Q}} [L (y\left<\boldsymbol{w},\boldsymbol{x}\right>)]. \tag{2} \label{2} \end{align}

また，$\mathcal{X}\times\mathcal{Y}$ 上の点同士の距離を次のように定めます：

\begin{align} d ( (\boldsymbol{x}, y), (\boldsymbol{x}', y') ) = \|\boldsymbol{x} - \boldsymbol{x}' \| + \kappa \boldsymbol{1}_\left\{y \neq y'\right\}. \tag{3} \label{3} \end{align}

ここで，$\kappa>0$ は正の定数であり，ラベルをフリップするために必要なコストを表します．

以上の条件の下で，解くべき分布的ロバスト最適化問題 (\ref{2}) は次と等価となります [8, Theorem 14 (ii)]：

\begin{align} \begin{split} \min_{\boldsymbol{w}, \lambda, \boldsymbol{s}} \quad &\lambda\rho + \frac{1}{N} \sum_{n=1}^{N} s_i \\ \text{subject to} \quad & L(\hat{y}_i \left<\boldsymbol{w},\hat{\boldsymbol{x}}_i\right>) \leqslant s_i \quad &(i = 1,2,\dots,N), \\ & L(-\hat{y}_i \left<\boldsymbol{w},\hat{\boldsymbol{x}}_i\right>) - \kappa\lambda \leqslant s_i \quad &(i = 1,2,\dots,N), \\ & \mathrm{lip}(L) \|\boldsymbol{w}\|_{\ast} \leqslant \lambda. \end{split} \tag{4} \label{4} \end{align}

ここで，$\mathrm{lip}(f)$ は関数 $f$ のリプシッツ定数を表します．この最適化問題は，目的関数は決定変数 $\lambda, \boldsymbol{s}$ に関して線形であり，それぞれの制約関数は決定変数に関して凸のため凸最適化問題となります．

例えば，一変数損失関数 $L$ がヒンジロス関数で距離の定義 (\ref{3}) におけるノルムが無限大ノルム ($\|\cdot\| = \|\cdot\|_{\infty}$) であるとき，問題 (\ref{4}) は以下のようになります：

\begin{align} \begin{split} \min_{\boldsymbol{w}, \lambda, \boldsymbol{s}\geqslant \boldsymbol{0}} \quad &\lambda\rho + \frac{1}{N} \sum_{n=1}^{N} s_i \\ \text{subject to} \quad & 1 - \hat{y}_i \left<\boldsymbol{w},\hat{\boldsymbol{x}}_i\right> \leqslant s_i \quad &(i = 1,2,\dots,N), \\ & 1 + \hat{y}_i \left<\boldsymbol{w},\hat{\boldsymbol{x}}_i\right> - \kappa\lambda \leqslant s_i \quad &(i = 1,2,\dots,N), \\ & \|\boldsymbol{w}\|_{1} \leqslant \lambda. \end{split} \tag{5} \label{5} \end{align}

ここで，無限大ノルムの双対ノルムは $\ell_1$ ノルムであることを用いました．

分布的ロバスト最適化問題 (\ref{2}) と凸最適化問題 (\ref{4}) が等価であることの証明は理論的に興味深いのですが，やや煩雑な議論を要するので概略を述べるだけに留めます．興味のある方だけ目を通してください．

等価性の証明の概略

以下では，分布的ロバスト最適化問題 (\ref{2}) の内側の最大化問題が強双対性により最小化問題に帰着できることを見ます．より一般的な議論や証明の詳細については [8, Appendix] を参照してください．

簡単のため，$\boldsymbol{\xi} := (\boldsymbol{x}, y)$ とし，損失関数を $\mathcal{L}(\boldsymbol{\xi})$，距離 (\ref{3}) を $d(\boldsymbol{\xi}_i, \boldsymbol{\xi}_j) = \| \boldsymbol{\xi}_i - \boldsymbol{\xi}_j \|$ で定義します．また，経験分布は訓練データ集合上で一様，すなわち，$\hat{\mathbb{P}}_N = \frac{1}{N}\sum_{n=1}^{N} \delta_{\boldsymbol{\xi}_n}$ であるとします．

問題 (\ref{2}) の内側の最大化問題を考えます．最適輸送距離に関する制約条件を明示的に記述すると，解くべき分布的ロバスト最適化問題 (\ref{2}) は次のようになります：

\begin{align} \begin{split} \max_{\mathbb{Q} \geqslant 0} \quad &\mathbb{E}_{\boldsymbol{\xi}\sim\mathbb{Q}}[\mathcal{L}(\boldsymbol{\xi})] \\ \text{subject to} \quad &\mathcal{W}_d (\mathbb{Q}, \hat{\mathbb{P}}_N) \leqslant \rho. \end{split} \tag{6} \label{6} \end{align}

この最適化問題の制約条件は，決定変数である確率分布 $\mathbb{Q} = (\underbrace{\mathbb{Q}_1,\dots,\mathbb{Q}_K}_{K})^\top$ と経験分布 $\hat{\mathbb{P}}_N = (\underbrace{\frac{1}{N},\dots,\frac{1}{N}}_{N})^\top$ の最適輸送距離が $\rho$ 以下であるというものです．準備で述べたように，最適輸送距離は線形計画問題 (\ref{LP}) の最適値として与えられるのでした．この最適輸送距離に関する最適化問題の決定変数 $\boldsymbol{\pi}$ をいま検討している分布的ロバスト最適化問題 (\ref{6}) の決定変数に押し込んだ以下の問題を考えます：

\begin{align} \begin{split} \max_{\mathbb{Q} \geqslant 0,\, \boldsymbol{\pi} \geqslant 0} \quad &\mathbb{E}_{\boldsymbol{\xi}\sim\mathbb{Q}}[\mathcal{L}(\boldsymbol{\xi})] = \sum_{k=1}^{K} \mathcal{L}(\boldsymbol{\xi}_k) \mathbb{Q}_k \\ \text{subject to} \quad & \sum_{k=1}^{K} \sum_{n=1}^{N} \| \boldsymbol{\xi}_k - \hat{\boldsymbol{\xi}}_n \| \pi_{kn} \leqslant \rho,\\ &\sum_{n=1}^{N} \pi_{kn} = \mathbb{Q}_k \quad (k = 1,2,\dots, K), \\ &\sum_{k=1}^{K} \pi_{kn} = (\hat{\mathbb{P}}_N)_{n} \quad (n = 1,2,\dots, N). \end{split} \tag{7} \label{7} \end{align}

さらに，$\mathbb{Q}$ に関する等式制約条件から決定変数 $\mathbb{Q}$ を消去することができます．すなわち，問題 (\ref{7}) は次のように再定式化できます：

\begin{align} \begin{split} \max_{\boldsymbol{\pi} \geqslant 0} \quad &\sum_{k=1}^{K} \sum_{n=1}^{N} \mathcal{L}(\boldsymbol{\xi}_k) \pi_{kn} \\ \text{subject to} \quad & \sum_{k=1}^{K} \sum_{n=1}^{N} \| \boldsymbol{\xi}_k - \hat{\boldsymbol{\xi}}_n \| \pi_{kn} \leqslant \rho,\\ &\sum_{k=1}^{K} \pi_{kn} = \frac{1}{N} \quad (n = 1,2,\dots, N). \end{split} \tag{8} \label{8} \end{align}

この問題 (\ref{8}) は決定変数を $\boldsymbol{\pi}$ とする線形計画問題であり，双対問題との間に強双対性が成り立ちます．1 つ目の制約条件に対するラグランジュ乗数を $\lambda \geqslant 0$，2 つ目の制約条件に対するラグランジュ乗数を $\boldsymbol{s}=(s_1,\dots,s_N)^\top$ とすると，問題 (\ref{8}) に対する双対問題は以下で与えられます：

\begin{align} \begin{split} \min_{\lambda\geqslant 0 , \boldsymbol{s}} \quad &\lambda\rho + \frac{1}{N}\sum_{n=1}^{N}s_i \\ \text{subject to} \quad &s_i + \lambda \| \boldsymbol{\xi}_k - \hat{\boldsymbol{\xi}}_n \| \geqslant \mathcal{L}(\boldsymbol{\xi}_k) \quad (k=1,2,\dots,K,\,n=1,2,\dots,N). \end{split} \tag{9} \label{9} \end{align}

以上より，分布的ロバスト最適化問題 (\ref{2}) の内側の最大化問題は最小化問題 (\ref{9}) と (最適値の意味で) 等価であることが分かりました．後は，問題 (\ref{2}) の外側の最小化問題と統合して解いてあげればよいという寸法です．

他の最適化手法との関係性

数値実験の紹介に移る前に，分布的ロバスト最適化問題 (\ref{2}) の特別な場合が (通常の) ロバスト最適化問題や正則化学習問題と一致することを見ます．個人的な感想になりますが，よく知られたロバスト最適化や正則化手法が分布的ロバスト最適化を通して特徴付けられる点は面白いと思いました．

ここでも，$\mathcal{Y}=\{=1,+1\}$ である二値分類問題を考えます．また，距離の定義 (\ref{3}) において $\kappa = +\infty$ としたときについて考えます．

ロバスト 最適化問題との関係性

分布的ロバスト最適化問題 (\ref{2}) の内側の最大化問題 $\max_{\mathbb{Q}\in\mathbb{B}_{\rho}(\hat{\mathbb{P}}_N)} \mathbb{E}_{(\boldsymbol{x}, y)\sim\mathbb{Q}} [L (y\left<\boldsymbol{w},\boldsymbol{x}\right>)]$ は (損失関数 $L$ に関する仮定の下で) 以下の最適化問題と等価となります [8, Theorem 23]：

\begin{align} \max_{\Delta\boldsymbol{x}_i} \quad & \frac{1}{N} \sum_{n=1}^{N} L(\hat{y}_n\left<\boldsymbol{w},\hat{\boldsymbol{x}}_i+\Delta\boldsymbol{x}_i\right>) \\ \text{subject to} \quad & \frac{1}{N}\sum_{n=1}^{N} \|\Delta\boldsymbol{x}_i\| \leqslant \rho. \end{align}

この最適化問題は決定変数に関する摂動を考慮した問題であり，通常のロバスト最適化問題と一致します．制約条件は摂動の平均値が $\rho>0$ で上からバウンドされることを意味します．

正則化手法との関係性

分布的ロバスト最適化問題 (\ref{2}) は (損失関数 $L$ に関する仮定の下で) 以下の最適化問題と等価となります [8, Remark 18]：

\begin{align} \min_{\boldsymbol{w}\in\mathbb{R}^d} \quad & \frac{1}{N} \sum_{n=1}^{N} L(\hat{y}_n\left<\boldsymbol{w},\hat{\boldsymbol{x}}_i \right>) + \rho \|\boldsymbol{w}\|_{\ast}. \end{align}

この最適化問題は通常の経験誤差最小化に正則化項を加えた最適化問題と見ることができ，正則化学習手法は分布的ロバスト最適化問題 (\ref{2}) の特別な場合とみなすことができます．

例えば，損失関数 $L$ がヒンジ関数で $\|\cdot\| = \|\cdot\|_2$ の場合 (i.e., $\|\cdot\|_{\ast} = \|\cdot\|_2$)，ソフトマージン SVM と等価となります．

数値実験

scikit-learn で提供されている breast cancer wisconsin dataset のデータセットを使って二値分類問題を解いてみましょう．

理想的には設定する確実性集合 $\mathbb{B}_{\rho}(\hat{\mathbb{P}}_{N})$ の中にデータを生成する真の分布が含まれるように $\rho > 0$ を設定することが望ましいですが，真の分布は未知であるため適切な $\rho$ の値を事前に知ることはできません．そこで，ここでの数値実験では，様々な $\rho$ に対して分布的ロバスト最適化問題 (\ref{5}) を解き，そのときの訓練データに対する正解率とテストデータに対する正解率を見てみたいと思います．

まず，データをロードします．本実験では，全体の 3 割 (170/569) を訓練データとし，残りの 7 割 (399/569) をテストデータとして用いることにします．

X, y = load_breast_cancer(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.7, random_state=42
)
y_train[y_train == 0] = -1
y_test[y_test == 0] = -1
print(X_train.shape)
print(X_test.shape)
# (170, 30)
# (399, 30)

次に，分布的ロバスト最適化問題 (\ref{5}) を解くための実装を以下に示します．

def fit(X, y, rho, kappa):
    n_sample = X.shape[0] # dimension of decision variable s
    n_dim = X.shape[1] # dimension of decision variable w
    # define decision variabel
    w = cp.Variable(n_dim)
    lam = cp.Variable(1)
    s = cp.Variable(n_sample)
    ones = np.ones(n_sample) / n_sample
    # define objective function
    objective = cp.Minimize(lam * rho + ones.T @ s)
    # define constraints
    constraints = [s >= 0]
    for i in range(n_sample):
        constraints.append(1 - y[i] * X[i,:] @ w <= s[i])
        constraints.append(1 + y[i] * X[i,:] @ w - kappa * lam <= s[i])
    constraints.append(cp.norm(w, 1) <= lam)
    # define optimization problem to be solved
    prob = cp.Problem(objective, constraints)
    prob.solve() # solve
    return w.value
    
def predict(X, w):
    return np.sign(X.dot(w))

def score(X, y, w):
    return accuracy_score(y, predict(X, w))

ここでの実験では，不確実性集合の大きさを表すパラメータ $\rho > 0$ とラベルをフリップするために必要なコストを表すパラメータ $\kappa > 0$ を次のように定めます： \begin{align} \rho &= \{10^{-4},\, 2\times 10^{-4},\, 5 \times 10^{-4},\,10^{-3},\, 2\times 10^{-3},\, 5 \times 10^{-3},\, 10^{-2},\, 2\times 10^{-2},\, 5 \times 10^{-2},\, 10^{-1} \},\\ \kappa &= \{0.1,\, 0.5,\, 1.0,\, 2.0\}. \end{align}

例えば，$\rho = 10^{-4},\, \kappa = 0.1$ のときは以下のように分類器 $h^\star\in\mathbb{H}$ を学習させることができます．

w = fit(X_train, y_train, rho=1e-4, kappa=0.1)
print(score(X_train, y_train, w))
print(score(X_test, y_test, w))
# 1.0
# 0.9298245614035088

各 $(\rho, \kappa)$ に対する分類器を学習し，訓練データに対する正解率 (左図) とテストデータに対する正解率 (右図) をプロットしたものを下図に示します．

左図を見ると，不確実性集合の大きさを表すパラメータが小さいときは訓練データに対して過適合する傾向が強く，逆に $\rho$ の値が大きいときは過適合の傾向が減少することが確認されます．また，右図を見ると，一概にどの $(\rho,\kappa)$ の組合せが良いということは明言できませんが，極端に $\rho$ の値を小さくしたり大きくすると性能は良くないことが伺えます．

おわりに

本エントリでは分布に対する不確実性を考慮した最適化モデリング手法である分布的ロバスト最適化法について扱いました．

本エントリでは仮説集合 $\mathbb{H}$ として最も単純である線形関数全体の集合に焦点を絞りましたが，より一般化した結果として仮説集合 $\mathbb{H}$ を正定値カーネルを用いて非線形化した研究 [8, §3.3] やニューラルネットワークを用いた場合の研究があります [8, §3.4]．

また，数値実験において不確実性集合 $\mathbb{B}_{\rho}(\hat{\mathbb{P}}_{N})$ の大きさを表すパラメータ $\rho>0$ を変えて実験をしましたが，どれくらいのデータ数に対してどれくらいの大きさに設定したら真の分布を含むのかに関して集中不等式を用いた理論的な結果などもあります [8, §4]．

今回はさらっと流しましたが，確率計画法／ロバスト最適化に関する詳細については例えば [9] や [10] を参照してください．分布的ロバスト最適化に関する包括的なレビューは例えば [11] や [12] を参照してください．

冒頭にも述べましたが，本エントリで扱った数値計算については，こちらの Google Colaboratory に集約してあります．参考になれば幸いです．

参考文献

[1] Finlayson, S. G., Bowers, J. D., Ito, J., Zittrain, J. L., Beam, A. L., & Kohane, I. S. (2019). Adversarial attacks on medical machine learning. Science, 363(6433), 1287-1289.
[2] Thys, S., Van Ranst, W., & Goedemé, T. (2019). Fooling automated surveillance cameras: adversarial patches to attack person detection. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops.
[3] Model Hacking ADAS to Pave Safer Roads for Autonomous Vehicles. (link, 2022/12/24 閲覧)
[4] Carlini, N., & Wagner, D. (2018, May). Audio adversarial examples: Targeted attacks on speech-to-text. In 2018 IEEE security and privacy workshops (SPW) (pp. 1-7). IEEE.
[5] 最適輸送入門 (link, 2022/12/24 閲覧)
[6] 最適輸送と自然言語処理 (link, 2022/12/24 閲覧)
[7] 最適輸送理論とその周辺 (link, 2022/12/24 閲覧)
[8] Shafieezadeh-Abadeh, S., Kuhn, D., & Esfahani, P. M. (2019). Regularization via mass transportation. Journal of Machine Learning Research, 20(103), 1-68.
[9] Birge, J. R., & Louveaux, F. (2011). Introduction to stochastic programming. Springer Science & Business Media.
[10] Ben-Tal, A., El Ghaoui, L., & Nemirovski, A. (2009). Robust optimization (Vol. 28). Princeton university press.
[11] Rahimian, H., & Mehrotra, S. (2019). Distributionally robust optimization: A review. arXiv preprint arXiv:1908.05659.
[12] Wiesemann, W., Kuhn, D., & Sim, M. (2014). Distributionally robust convex optimization. Operations Research, 62(6), 1358-1376.
[13] Peyré, G., & Cuturi, M. (2019). Computational optimal transport: With applications to data science. Foundations and Trends® in Machine Learning, 11(5-6), 355-607.

補遺

"What are some of the important points to keep in mind when using machine learning models and deep learning models for various real-world tasks?" に対する ChatGPT の回答：

There are several important points to keep in mind when using machine learning (ML) models and deep learning (DL) models for real-world tasks:

Data quality and quantity: It is important to have a large and diverse dataset to train your ML/DL model on. The data should also be clean and free of errors, as this can significantly affect the model's performance.

Feature engineering: Good feature engineering is crucial for ML/DL models to work well. This involves selecting the most relevant features from the data and transforming them in a way that the model can learn from them effectively.

Model selection and hyperparameter tuning: Choosing the right model and adjusting its hyperparameters can have a big impact on the model's performance. It is important to experiment with different models and hyperparameter settings to find the combination that works best for your task.

Overfitting and underfitting: Overfitting occurs when a model is too complex and has learned the noise in the data rather than the underlying patterns, which leads to poor generalization to new data. Underfitting occurs when the model is too simple and is unable to capture the complexity of the data, leading to poor performance. It is important to avoid both overfitting and underfitting by choosing an appropriate model complexity and using techniques such as regularization and early stopping.

Evaluation and monitoring: It is important to evaluate the model's performance on a separate test set and to monitor the model's performance over time to ensure that it is still performing well on new data.

Ethical considerations: ML/DL models can be biased if the training data is biased, and this can lead to unfair outcomes. It is important to consider ethical implications when building and deploying ML/DL models.

冷めたコーヒー

Weniger, aber besser