冷めたコーヒー

Weniger, aber besser

研究関連でお世話になった学術書 (2023 年)

tl;dr

本エントリは,2023 年 (主に下半期) にお世話になった学術書を紹介するコーナーです. 昨年は,不確実な環境下における逐次的な意思決定を扱うような業務が多く,それに関連した書籍を主に参照しました. 以下が本エントリで紹介する書籍のリストです:

あたらしい数理最適化

業務最適化でそこそこ参照しました.数理最適化や凸解析の理論的な知識があっても,業務で必要とされるような現実の諸問題をモデル化したり,効率の良い求解方法を考えたりするのは思ったよりも難しいです. 本書が発行されたのは 2012 年ですが,現代でも活かされるようなマニアック (?) なモデリングのコツ (e.g., Big-M における M はなるべく小さく設定する) などについても触れられていて参考になります. また,様々な典型的な問題に対する数理計画問題としての定式化が網羅されていて,実際の現場で必要とされる問題を考える上での出発点として参考になります.

関連書籍

関連ページ

オンライン予測

オンライン最適化(不確実な環境下における逐次的な意思決定をを扱う数理モデル)について知るための一冊目として読みました. オンライン最適化は従来の数理最適化とは異なり,意思決定をするときに最適化すべき目的関数が未知であるような状況を扱います. この書籍も発行されたのは 2016 年と,近年発展が目覚ましい当該分野においては少々古いですが,一冊目として概要を抑えるという観点では良い入門書であると思います. 後に述べますが,より専門的な内容を知りたい場合は,Elad Hazan の Introduction to Online Convex Optimization を参照すると良いです.

バンディット問題の理論とアルゴリズム

オンライン最適化に近い問題設定としてバンディット問題があります.オンライン最適化との主な違いとして,オンライン最適化では意思決定をする度に得られる報酬 (損失) に関する関数形が与えられるが,バンディット問題では関数形が与えられず,実際に得られた報酬(被った損失) のみ与えられるような状況を考えます.バンディット問題に対する有名なアルゴリズムとして,UCB アルゴリズムやトンプソンサンプリングに基づく方策がありますが,それらの紹介にとどまらずリグレットに関する理論評価も本書で扱われています.また,発展として,best arm identification やベイズ最適化に関する内容も盛り込まれており,この一冊でバンディット問題にまつわる内容を幅広く見通すことができます.

Optunaによるブラックボックス最適化

ブラックボックス最適化についても扱う必要があったので,手始めに Optuna について知っておこうと思い読みました . Optuna はご存じの通り Preferred Networks 社が開発しているブラックボックス最適化のためのツールです. 本書では,ブラックボックス最適化に関する最低限の内容を述べた後に,Optuna の使用例が紹介されています.

関連ページ

ベイズ最適化: 適応的実験計画の基礎と実践

ベイズ最適化に関する新しい (2023 年 11 月発行) 書籍です.上記の『Optunaによるブラックボックス最適化』はあくまで “与えられたブラックボックス関数を Optuna を用いてどう最適化するのか” に重点が置かれていましたが,本書ではもう少し理論的な背景が紹介されています.また,Optuna に関する記述もあるだけでなく,BoTorch (Meta 社が開発しているベイズ最適化のためのツール) をバックエンドとした BoTorchSampler による実装例も紹介されています. 理論と応用のバランスが取れた良い書籍ですので,これからベイズ最適化に入門したいと思っている方には是非手に取って頂きたい一冊です.

関連ページ

番外編

オンライン最適化に関しては,その道の大家である Elad Hazan による Introduction to Online Convex Optimization を最終的によく参照しました.ありがたいことに,arXiv 版が こちら から参照できます. また,バンディット問題に関しては,こちらもその道の大家である Tor Lattimore と Csaba Szepesvári による Bandit Algorithms を頻繁に参照しました.こちらに関しても,著者の HP で全文公開されています.

(2024/01/07: 追記) ガウス過程について習熟するために『ガウス過程と機械学習』も頻繁に参照しました(が,2023 年以前から参照していたので上記のリストには入れませんでした).ガウス過程の原理について初歩から分かりやすく書かれていて一冊目として強くお勧めできる書籍です.また,ガウス過程については,Williams, C. K., & Rasmussen, C. E. の "Gaussian Processes for Machine Learning" が こちらのリンク から入手できます.同様に,ガウス過程のバンディット設定の下での解析については [Srinivas et al. (ICML '09)] の "Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design" (arXiv) をよく参照しました.

最後に

最近,本エントリで紹介した周辺分野で論文を書いたり投稿したりしていますが,なかなか難しさを感じています.議論してくださる方がいたらぜひ議論させてください.(@mirucaaura)

【凸解析】平滑性に関する諸性質

これは何?

機械学習や数理最適化に関する論文を読んでいると,アルゴリズムの収束性の議論をするときなどに関数の平滑性 (Smoothness) を仮定することがあります.

先日,友人氏 (@IamTakala) から平滑性に関する性質をまとめた記事 [1] において,[1] -> [0] は成立しないのか,という旨の連絡を受け,凸性を仮定すれば成立することを確認したので,その内容をまとめておこうと思います.記事 [1] の [1] -> [0] がどのような主張かというと,「関数 $g(x) = \frac{L}{2} \|x\|^2 - f(x)$ が凸関数であるとき,関数 $f$ は $L$-平滑である」というものです.

なお,本記事では,ベクトルであっても太字を用いないこととし,$\|\cdot\|$ は適当なノルムを表すものとします.

また,本記事は 数理最適化のカレンダー | Advent Calendar 2023 - Qiita の第 11 日目の代理として書かれたものです(何か不都合があれば,@mirucaaura までご連絡頂ければと思います).

準備

まず,関数 $f$ に関する平滑性は次のように定義されるのであった:

定義 1 (平滑性) 連続的微分可能な関数 $f:\mathbb{R}^n\to\mathbb{R}$ が $L$-平滑 ($L$-smooth) であるとは,任意の $x, y\in\mathbb{R}^n$ に対して, $$ \begin{align} \|\nabla f(x) - f(y) \| \leq L \|x - y \| \end{align} $$ が成り立つことを言う.

上記の平滑性の定義において,関数 $f$ に対する凸性は必要ないことに注意されたい.また,上記の $f$ の平滑性は,$\nabla f(x)$ に関する Lipschitz 連続性と等価であり,$L$ は Lipschitz 定数に他ならない.

次に,関数 $f$ が凸関数であるとき,平滑であることと等価な条件を確認する.証明は [2] などを参照されたい.

命題 1 (平滑性の必要十分条件) [2, 命題 2.3.5] 関数 $f:\mathbb{R}^n\to\mathbb{R}$ を連続的微分可能な凸関数とし,$L>0$ とする.このとき,$f$ が $L$-平滑であるための必要十分条件は次の (1)-(3) である:
(1) $f(y) \leq f(x) + \left<\nabla f(x), y-x\right> + \frac{L}{2} \|x-y\|^2$
(2) $f(y) \geq f(x) + \left<\nabla f(x), y-x\right> + \frac{1}{2L} \|x-y\|^2$
(3) $\left<\nabla f(x) - \nabla f(y), x-y\right> \geq \frac{1}{L} \|\nabla f(x) - \nabla f(y)\|^2$

本題

さて,懸案事項は,凸関数 $f$ が $L$-平滑であるとき,以下のような関係が成立するかどうかである (実際に成立するので命題という形で書いている):

命題 関数 $g(x) = \frac{L}{2} \|x\|^2 - f(x)$ は凸関数であるとする.このとき,$f$ は $L$-平滑,すなわち, 任意の $x, y\in\mathbb{R}^n$ に対して, $$ \begin{align} \|\nabla f(x) - f(y) \| \leq L \|x - y \| \end{align} $$ が成り立つ.

証明の方針:凸関数 $f$ が命題の条件を満たすとき,$f$ に対する $L$-平滑性と等価な条件式 (命題 1-(2)) が成立することを確かめる.

証明.まず,関数 $\phi_x:\mathbb{R}^n\to\mathbb{R}$ を次のように定める:

\begin{align} \phi_x(z) := f(z) - \left<\nabla f(x), z\right>. \end{align}

この関数 $\phi_x$ の勾配ベクトルは,

\begin{align} \nabla \phi_x(z) = \nabla f(z) - \nabla f(x) \end{align}

であるから,$z=x$ において最小となる ($\because \nabla \phi_x(x) = \nabla f(x) - \nabla f(x) = 0$).

いま,命題の条件より,関数 $h(z) := \frac{L}{2} \|z\|^2 - \phi_x(z)$ も凸である.実際,この関数 $h:\mathbb{R}^n\to\mathbb{R}$ は

\begin{align} h(z) = \underbrace{\frac{L}{2} \|z\|^2 - f(z)}_{=g(z)} + \left<\nabla f(x), z\right> \end{align}

であり,凸関数 $g$ に線形項を加えた関数であり,したがって凸である.

さて,関数 $h$ の凸性より,任意の $z, y \in \mathbb{R}^n$ に対して次の不等式が成り立つ:

\begin{align} h(z) \geq h(y) + \left<\nabla h(y), z-y\right>. \end{align}

つまり,任意の $y, z\in\mathbb{R}^n$ に対して,

\begin{align} &\frac{L}{2}\|z\|^2 - \phi_x(z) \geq \frac{L}{2}\|y\|^2 - \phi_x(y) + \left< L y - \nabla \phi_x(y), z-y \right> \\ \iff\quad& \phi_x(z) \leq \phi_x(y) + \left<\nabla \phi_x(y), z-y\right> + \underbrace{\frac{L}{2}\|z\|^2 - \frac{L}{2}\|y\|^2 - \left< Ly, z-y\right>}_{=\frac{L}{2}\|z-y\|^2}. \end{align}

が成り立つ.上式において,左辺の $\phi_x(z)$ は $z=x$ において最小となる.また,右辺全体を $\varphi(z)$ とおくと,$\nabla \varphi(z) = \nabla \phi_x(y) + L(z-y)$ であるから,右辺 $\varphi(z)$ は $z=y-\frac{1}{L}\nabla\phi_x(y)$ で最小となる.したがって,

\begin{align} f(x) - \left< \nabla f(x), x \right> &\leq \phi_x(y) + \left< \nabla \phi_x(y) , -\frac{1}{L} \nabla \phi_x(y) \right> + \frac{L}{2} \left\| -\frac{1}{L} \nabla \phi_x(y) \right\|^2 \\ &= \phi_x(y) - \frac{1}{L} \| \nabla \phi_x(y) \|^2 + \frac{1}{2L} \| \nabla \phi_x(y) \|^2 \\ &= f(y) - \left< \nabla f(x), y\right> - \frac{1}{2L} \| \nabla f(y) - \nabla f(x)\|^2. \end{align}

整理すると,所望の不等式 (命題 2-(2)):

\begin{align} f(y) - f(x) - \left<\nabla f(x), y-x \right> \geq \frac{1}{2L} \| \nabla f(y) - \nabla f(x) \|^2 \end{align}

が得られる.したがって,関数 $g(x) = \frac{L}{2} \|x\|^2 - f(x)$ が凸であるとき,関数 $f$ は $L$-平滑である.(証明終

参考文献

[1] Lipschitz continuous gradient · Xingyu Zhou's blog (2023/12/26 閲覧)
[2] 飯塚秀明 (2023)『連続最適化アルゴリズム』, オーム社

歯科検診へ行こう

これは何?

歯科検診を定期的に受診することを薦めるエントリです.わたしは 2-3 ヵ月に 1 度のペースで受診しています.周りの人と話していると,年に一度も歯医者を受診しない話をしばしば耳にします.わたし自身,歯科医師でも歯科衛生士でも歯科技工士でも何でもないのですが,定期的に自分の歯の状態を専門家に診てもらうことは有意義であり,コスパの良い自己投資であると思います.

統計情報

厚生労働省「平成 28 年国民健康・栄養調査」[1] によると「過去 1 年間に歯科検診を受けた者の割合」は 52.9% ($N=13,104$) であり,年齢階級別にみると,男女ともに年齢が高い層で増加の割合が高い結果となっています.例えば,20-29 歳では 43.3 % ($N=774$),60-69 歳では 58.1 % ($N=3,185$) となっています.

歯科検診で何をするのか

これはあくまで自分が受診している歯科医院の事例ですが,30-40 分くらいかけて歯のクリーニング,フッ素塗布,歯石の除去,歯周病・虫歯の有無の検査などを行います (ただし,初診ではレントゲンを撮って虫歯の有無や歯周病の進行状態などを確認することが一般的であると思います).また,普段の歯磨きにおいてどのあたりに磨き残しが多くあるのかを教えてもらったりします.わたしは普段 Philips の DiamondClean 9000 シリーズ (e.g., HX9911/67) を使っていますが,いくら電動歯ブラシといえども磨き残しが発生するということを痛感しています.

継続的に通院するコツ

継続的に通院するためのコツとして,なるべく徒歩圏内の歯医者を受診するべきであると思います.また,受診した当日に次回の診療日の予約をしてしまうと良いと思います.気持ちの問題として,「歯科医に何か言われるのが怖い」というのがありますが,とにかく一度近くの歯科医院で評判の良さそうなところに行ってみると良いと思います.

さいごに

歯科検診は圧倒的にコスパの良い自己投資だと思います.即効性のある投資ではないですが,積み重ねによる影響が大きい類の投資だと思います.本エントリが読者の皆様の一助となれば幸いです.

参考

[1] 厚生労働省「平成28年国民健康・栄養調査」

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

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

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

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

はじめに

様々な現実のタスクに対して機械学習/深層学習モデルが人間に匹敵,あるいは,凌駕するような性能を発揮することがごく当たり前の世の中になっています. 本エントリを執筆時点 (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:

  1. 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.

  2. 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.

  3. 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.

  4. 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.

  5. 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.

  6. 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.

混合整数計画問題のソルバー scipy.optimize.milp と Python-MIP でのデフォルトのソルバー CBC の簡単な比較

前回の記事で,SciPy==1.9.0 で新たに導入された scipy.optimize.milp に関する簡単な紹介をしました.

mirucacule.hatenablog.com

本記事では,上記の記事で扱った問題と同様の問題を Pyhton-MIP で求解する方法を紹介します.Python-MIP のより詳細な使用方法については公式ドキュメント [1] を参照してください.

解くべき混合整数計画問題(再掲)

次の混合整数計画問題を考えます:

\begin{align} \min_{\mathbf{x}} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{A}\mathbf{x} \leqslant \mathbf{b},\\ &\mathbf{x} \in \mathbb{Z}^n,\,\mathbf{x} \geqslant \mathbf{0}. \end{align}

ここで,$\mathbf{c}\in\mathbb{R}^n, \, \mathbf{A}\in\mathbb{R}^{m\times n}, \, \mathbf{b}\in\mathbb{R}^{m}$ とします.

scipy.optimize.milp による求解(再掲)

必要なパッケージを import しておきます.

import numpy as np
from scipy.optimize import milp
from scipy.optimize import LinearConstraint
from scipy.optimize import Bounds

解くべき問題の係数行列などを定義します.

np.random.seed(42)
n, m = 40, 80
c = np.random.rand(n)
x = np.random.randint(0, 5, n)
A = np.random.randn(m, n)
s = 10 * np.ones_like(m)
b = A @ x + s

制約条件に関わる部分を定義します.

l = np.zeros_like(c)
u = np.full_like(c, np.inf)
b_u = b
b_l = np.full_like(b_u, -np.inf)

bounds = Bounds(l, u)
constraints = LinearConstraint(A, b_l, b_u)
integrality = np.ones_like(c)

求解します.

%%time
res = milp(
    c=c,
    bounds=bounds,
    constraints=constraints,
    integrality=integrality,
)
print(res)

結果は以下の通りです.

            fun: 7.441844504063867
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 7.44111388697805
        mip_gap: 9.8176881473139e-05
 mip_node_count: 49538
         status: 0
        success: True
              x: array([1., 0., 1., 0., 0., 1., 1., 0., 0., 0., 1., 0., 0., 4., 0., 0., 0.,
       0., 0., 0., 2., 1., 0., 0., 0., 0., 3., 1., 1., 0., 1., 3., 0., 0.,
       1., 0., 0., 1., 0., 0.])
CPU times: user 17.8 s, sys: 665 ms, total: 18.4 s
Wall time: 17.6 s

Python-MIP の CBC による求解

必要なパッケージを import しておきます.

import numpy as np
import mip

同様の問題を扱うため,係数行列などの定義は同様です.

np.random.seed(42)
n, m = 40, 80
c = np.random.rand(n)
x = np.random.randint(0, 5, n)
A = np.random.randn(m, n)
s = 10 * np.ones_like(m)
b = A @ x + s

次に,決定変数,目的関数,制約条件をそれぞれ定義して求解を行います.

model = mip.Model()
# define decision variable
x = {}
for i in range(n):
    x[i] = model.add_var(lb=0, ub=mip.INF, var_type=mip.INTEGER)
# define objective function
model.objective = mip.minimize(mip.xsum(c[i] * x[i] for i in range(n)))
# define constraints
for i in range(m):
    model.add_constr(
        mip.xsum(A[i,:][j] * x[j] for j in range(n)) <= b[i]
    )
# solve
model.optimize()

model.optimize() を実行すると以下のような出力が得られます.実際には,最適化計算の各反復における情報が出力されますが,ここでは省略します.

Cbc0032I Strong branching done 29094 times (142226 iterations), fathomed 305 nodes and fixed 796 variables
Cbc0041I Maximum depth 33, 3997 variables fixed on reduced cost (complete fathoming 1136 times, 157339 nodes taking 667811 iterations)
Total time (CPU seconds):       8.99   (Wallclock seconds):       9.63

<OptimizationStatus.OPTIMAL: 0>

得られた最適値と最適解を確認してみましょう.

print('Objective value: {model.objective_value:.15}'.format(**locals()))
print('Solution: ', end='')
for v in model.vars:
    if v.x > 1e-5:
        print('{v.name} = {v.x}'.format(**locals()))
        print('          ', end='')

実行すると以下のような結果が得られます.

Objective value: 7.44184450406387
Solution: var(0) = 1.0
          var(2) = 1.0
          var(5) = 1.0
          var(6) = 1.0
          var(10) = 1.0
          var(13) = 4.0
          var(20) = 2.0
          var(21) = 1.0
          var(26) = 3.0
          var(27) = 1.0
          var(28) = 1.0
          var(30) = 1.0
          var(31) = 3.0
          var(34) = 1.0
          var(37) = 1.0

最適値,最適解ともに scipy.optimize.milp で計算した結果と整合していることが分かります.なお,求解にかかった時間は CPU 時間だけで比較すると,milp は 17.8 秒,Python-MIP (cbc) は 8.99 秒でした.

参考

[1] Python MIP Documentation