冷めたコーヒー

Weniger, aber besser

歯科検診へ行こう

これは何?

歯科検診を定期的に受診することを薦めるエントリです.わたしは 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

SciPy==1.9.0 で使用可能になった混合整数計画問題のソルバー scipy.optimize.milp に関する備忘録

はじめに

この記事の目的:

  • SciPy==1.9.0 で導入された scipy.optimize.milp の簡単な説明
  • scipy.optimize.milp の使い方の解説

この記事で扱わないこと

  • scipy.optimize.milp で採用されている最適化手法
  • 他のソルバーとの比較

[追記]

  • 公式ドキュメントにナップザック問題を解くチュートリアルが掲載されたので参考 [3] にリンクを追記しました (2022/09/08)

SciPy==1.9.0 のリリース

Atsushi さんによる以下の記事を参照してください.

myenigma.hatenablog.com

先行記事

書籍 [2] で挙げられている問題を scipy.optimize.milp を用いて解く方法について解説されています.

github.com

scipy.optimize.milp で解ける最適化問題

scipy.optimize.milp では以下の形式の混合整数計画問題を解くことができる:

\begin{align} \min_{\mathbf{x} \in \mathbb{R}^n} \quad &\mathbf{c}^\top\mathbf{x} \\ \text{subject to} \quad&\mathbf{b}_l \leqslant\mathbf{A}\mathbf{x} \leqslant \mathbf{b}_u,\\ &\mathbf{l} \leqslant \mathbf{x} \leqslant \mathbf{u},\\ &x_i \in \mathbb{Z},\, i \in \mathcal{I} \subseteq [n]. \end{align}

ここで,$\mathcal{I}$ は添字集合 $[n] := \{1,2,\dots,n\}$ の部分集合 (e.g., $\mathcal{I} = \{2, 3, 5 \}$) を表す.

残念ながら,scipy.optimize.milp を用いて上記の形式でない混合整数計画問題を解くことは今のところできません.頑張って上記の形式に書き直すか,他の最適化モデリング言語などを用いましょう.

scipy.optimize.milp の使い方

導入

以下で scipy==1.9.0 のインストールをしましょう:

python -m pip install scipy==1.9.0

以下では,scipy.optimize に含まれるパッケージを使用します.ついでに numpy も使用するので import しておきましょう.

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

例 1

まず,簡単な問題として Mathworksintlinprog のドキュメントで扱われている以下のような混合整数計画問題を考えます:

\begin{align} \min_{x_1,x_2} \quad &8x_1 + x_2 \\ \text{subject to} \quad&x_1 + 2x_2 \geqslant -14,\\ &-4x_1 - x_2 \leqslant -33,\\ &2x_1 + x_2 \leqslant 20,\\ &x_2 \in \mathbb{Z}. \end{align}

この問題はベクトルと行列を用いて次のように表すことができます:

\begin{align} \min_{x_1,x_2} \quad &\begin{bmatrix}8 \\ 1\end{bmatrix}^\top \begin{bmatrix}x_1 \\ x_2\end{bmatrix} \\ \text{subject to} \quad&\begin{bmatrix}-1 & -2 \\ -4 & -1 \\2 & 1\end{bmatrix}\begin{bmatrix}x_1 \\ x_2\end{bmatrix} \leqslant \begin{bmatrix}14 \\ -33 \\ 20\end{bmatrix},\\ &x_2 \in \mathbb{Z}. \end{align}

この問題を scipy.minimize.milp を用いて解けるように係数ベクトルや行列を以下のように定めます.なお,この問題は $(x_1,x_2)^\top$ に関する上下限制約はないので特別指定する必要はありません.また,$\mathbf{Ax}$ に関する下限制約もないので,以下のように b_l の各要素に -np.inf を設定します.整数制約は $x_2$ に対してのみ課されるため,以下のように integrality = np.array([0, 1]) と設定します.なお,0 は連続変数,1 は整数制約を表します.他に 23 を設定することもできますが,それらの説明は 公式ドキュメント を参照してください.また,integrality を指定しなかった場合,デフォルトでは 0 が指定されます.

c = np.array([8, 1])
A = np.array([[-1, -2], [-4, -1], [2, 1]])
b_u = np.array([14, -33, 20])
b_l = np.full_like(b_u, -np.inf)
constraints = LinearConstraint(A, b_l, b_u) $b_l <= Ax <= b_u
integrality = np.array([0, 1]) # x_2 is integer

求解は以下のように行います.

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

結果は以下のようになります.

            fun: 59.0
        message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'
 mip_dual_bound: 59.0
        mip_gap: 0.0
 mip_node_count: 1
         status: 0
        success: True
              x: array([6.5, 7. ])
CPU times: user 5.49 ms, sys: 0 ns, total: 5.49 ms
Wall time: 3.24 ms

Mathworksドキュメント に掲載されている結果と同じ結果が得られていることを確認できます.

例 2

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

\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}$ とします.この問題の実行可能領域 $\mathcal{F}=\left\{ \mathbf{x} \in \mathbb{Z}^n \mid \mathbf{A}\mathbf{x} \leqslant \mathbf{b},\, \mathbf{x} \geqslant \mathbf{0} \right\}$ に実行可能解が存在することを保証するために,$\overline{\mathbf{x}} \in \mathbf{Z}^n$ をランダムに生成し,その $\overline{\mathbf{x}} \in \mathbf{Z}^n$ が実行可能領域 $\mathcal{F}$ に含まれるように $\mathbf{A}$ と $\mathbf{b}$ を設定します.

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

決定変数 $\mathbf{x}$ に対する下限制約が存在するため,scipy.optimize.Bounds を使用します.決定変数 $\mathbf{x}$ に対する上限制約は存在しないので,以下のように u の各要素に np.inf を設定します.また,今回の問題では決定変数のすべての要素は整数であることを課しているため,integrality のすべての要素を 1 とします.

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

決定変数の次元 $n=40$ の問題を解いてみました.およそ 20 秒弱の時間で最適解を計算することができました.

おわりに

この記事では SciPy==1.9.0 で新たに導入された scipy.optimize.milp を用いて混合整数計画問題を解く方法について解説をしました.

参考

[1] scipy.optimize.milp
[2] 『Pythonではじめる数理最適化: ケーススタディモデリングのスキルを身につけよう』

[3] Knapsack problem example

リーマン多様体上の最適化―特異値分解の例を通して―

はじめに

以前(2019 年 11 月)に「リーマン多様体上の最適化の初歩と Pymanopt による数値実験」 という記事を投稿した.

mirucacule.hatenablog.com

記事内で用いた最適化 Toolbox である Pymanopt のバージョンアップに伴い,実行方法に変更が生じたため改めて書き直そうと思ったのが本記事を執筆するに至った経緯である.また,リーマン多様体上での共役勾配法を簡単に勉強したので,それを用いて行列の特異値分解を行う問題について扱いたい.

なお,本記事を執筆するに際して参考にした文献は 参考文献 にまとめてある.特に,本記事の内容は [1] の Chapter 4 に書かれている内容の重要と思われる箇所を抜粋してまとめたものであり,より深く知りたい方は [1] を参照されたい.本記事で扱っていない Newton 法や,共役勾配法と Newton 法を組み合わせたアルゴリズムに関する記述もあり楽しめると思う.

また,Pymanopt の使い方についてのみ知りたい読者は Pymanopt による求解 の節を確認されたい.より詳細な使い方や例については以下を参照されたい.

特異値分解

一般に,行列 $A \in \mathbb{R}^{m \times n},\ (m \ge n)$ は直交行列 $U \in \mathbb{R}^{m \times p},\, V \in \mathbb{R}^{n \times p}$ を用いて $$ A = U D V^\top \tag{1} \label{eq:1} $$

のように分解することができる.ここで,$D \in \mathbb{R}^{n\times n}$ は対角行列であり,その対角成分を特異値と呼ぶ.また,分解 (\ref{eq:1}) を行列 $A$ の特異値分解と呼ぶ.具体的な例としては [6] などを参照されたい.

特異値分解最適化問題

特異値分解は以下の最適化問題と密接な関係がある.

$$ \begin{align} (\mathrm{P}_1)\qquad \min_{U,V}\quad& -\mathrm{tr}(U^\top AVN)\\ \text{subject to}\quad& U\in \mathbb{R}^{m\times p},\, V \in \mathbb{R}^{n\times p},\, U^\top U = V^\top V = I_p. \end{align} $$

ここで,$A \in \mathbb{R}^{m\times n}$ は所与の行列,$N \in \mathbb{R}^{p\times p}$ は所与の対角行列をそれぞれ表す. 特に,$N$ の対角成分は $\mu_1 > \mu_2 > \dots > \mu_p > 0$ を満たす.また,$1 \le p \le n$ である.このとき,以下の命題が成り立つ.

命題 [4, Proposition 2.1] $m \ge n$ とする.最適化問題 ($\mathrm{P}_1$) の最適解を $(U_\star, V_\star)$ とする.このとき,$U_\star$ と $V_\star$ の列ベクトルたちは,それぞれ $A$ の小さい方から $p$ 個の特異値に属する左特異ベクトル,右特異ベクトルとなる.

この命題より,$U_\star = (\mathbf{u}_1,\dots,\mathbf{u}_p),\, V_\star = (\mathbf{v}_1,\dots,\mathbf{v}_p)$ としたとき,行列 $A$ の特異値分解

$$ A = \sum_{i=1}^{p} d_i \mathbf{u}_i \mathbf{v}_i^\top $$

と表される.

リーマン多様体上での特異値分解

ユークリッド空間上の最適化問題 ($\mathrm{P}_1$) は次のリーマン多様体上の最適化問題と等価である.

$$ \begin{align} (\mathrm{P})\qquad \min_{U,V}\quad& F(U, V) = -\mathrm{tr}(U^\top AVN)\\ \text{subject to}\quad& (U,V) \in \mathrm{St}(p,m) \times \mathrm{St}(p,n). \end{align} $$

ここで,制約条件に現れる $\mathrm{St}(p,n)$ は Stiefel 多様体といい,

$$\mathrm{St}(p,n) := \left\{Y \in \mathbb{R}^{n \times p} \mid Y^\top Y = I_p \right\}$$

で定義される.また,$\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ は $\mathrm{St}(p,m)$ と $\mathrm{St}(p,n)$ の積多様体を表す.

$\mathrm{St}(p,m)\times\mathrm{St}(p,n)$ 上の最適化

リーマン多様体 $\mathcal{M}$ 上での最適化を考える上で必要となるレトラクション,勾配について述べる.一般的な定義や導出方法については [1] や [2] を参照されたい.ここでは単に Stiefel 多様体の積多様体におけるレトラクションと勾配について,その結果についてのみ提示するに留める.以下の図は,リーマン多様体上の最適化手法で生成される点列のイメージである.リーマン多様体 $\mathcal{M}$ 上の点 $\mathbf{x}_k$ を探索方向 $\mathbf{\xi}_k \in T_{\mathbf{x}_k} \mathcal{M}$ 方向に $t_k>0$ でスケールして動かし,レトラクション $R_{\mathbf{x}_k}$ という操作によってリーマン多様体 $\mathcal{M}$ 上に写すことを繰り返し行いながら探索を行う.

リーマン多様体上のレトラクション

接空間

多様体 $\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ 上の点 $(U, V)$ における接空間 $T_{(U,V)}(\mathrm{St}(p,m) \times \mathrm{St}(p,n))$ は次で与えられる.

\begin{align} T_{(U,V)}(\mathrm{St}(p,m) \times \mathrm{St}(p,n)) &\simeq T_{U}(\mathrm{St}(p,m)) \times T_{V}(\mathrm{St}(p,n)) \\ &= \left\{ (\xi, \zeta) \in \mathbb{R}^{m \times p} \times \mathbb{R}^{n \times p} \mid \xi^\top U + U^\top\xi = \zeta^\top V + V^\top\zeta = O_p \right\}. \end{align}

ここで,$O_p$ は $p$ 次零行列を表す.

レトラクション $R_{(U,V)}$

レトラクションは接空間 $T_{(U, V)} (\mathrm{St}(p,m) \times \mathrm{St}(p,n))$ からリーマン多様体 $\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ への写像である.図の緑色の破線がレトラクションに対応する.点 $(U,V)\in\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ における $\mathrm{St}(p,m) \times \mathrm{St}(p,n)$ のレトラクションは次のように定められる.

$$ R_{(U, V)} (\xi, \zeta) = \left( \mathrm{qf}( U + \xi ), \mathrm{qf}( V + \zeta ) \right), \quad (\xi, \zeta) \in T_{(U, V)} (\mathrm{St}(p,m) \times \mathrm{St}(p,n)). $$

ここで,$\mathrm{qf}:\mathbb{R}^{m \times n} \to \mathrm{St}(p, n)$ であり,具体的にはフルランク行列 $X \in \mathbb{R}^{n \times p}$ に対して QR 分解を行ったときの $Q \in \mathrm{St}(p, n)$ を返す写像として定義される.ここで,$R$ は $p$ 次上三角行列であり,その対角成分はすべて正とする.

勾配 $\mathrm{grad} F(U,V)$

目的関数 $F(U,V) = -\mathrm{tr}(U^\top AVN)$ の点 $(U, V) \in \mathrm{St}(p,m)\times\mathrm{St}(p,n)$ における勾配は以下で定義される.

$$ \mathrm{grad} F(U, V) = \left( U \mathrm{sym}(U^\top AVN) - AVN, V \mathrm{sym}(V^\top A^\top UN) - A^\top UN \right). $$

ここで,$\mathrm{sym}(X)$ は正方行列 $X \in \mathbb{R}^{k\times k}$ に対して $\mathrm{sym}(X) = \frac{X + X^\top}{2}$ で定義される.

リーマン多様体上での共役勾配法

ユークリッド空間上での最適化アルゴリズムとして共役勾配法(Conjugate gradient descent method; CG 法)がある.この手法は,暫定解における探索方向を上手く選ぶことによって最急降下法のようなジグザグとした挙動を抑えようというアルゴリズムである.もともとは線形方程式系 $Ax = b$ を解く手法として考案されたアルゴリズムだったが,関数を最小化する問題に対しても適用できるということで様々な拡張が考えられてきた.この辺の詳しい話は,例えば [3] の Chapter 5 に詳しく書かれているので興味がある方は参照されたい. また,リーマン多様体上での共役勾配法アルゴリズムに関する和文記事は例えば [5] を参照されたい.

Pymanopt による求解

以下では,最適化問題 (P) を Pymanopt で解く手順について説明する.Pymanopt はリーマン多様体上での最適化を行うフレームワークである.本記事における各モジュールのバージョンは以下の通りである.

  • pymanopt == 2.0.0
  • autograd == 1.4
  • numpy == 1.23.0

モジュールのインポート

まずは必要なモジュールをインポートする.今回,リーマン多様体として Stiefel 多様体と Stiefel 多様体の積多様体を考えるため,pymanopt.manifolds から StiefelProduct をインポートする.また,今回は最適化手法として共役勾配法を用いるため,pymanopt.optimizers から ConjugateGradient をインポートする.

from pymanopt.manifolds import Stiefel, Product
from pymanopt.optimizers import ConjugateGradient
from pymanopt import Problem
import pymanopt
import autograd.numpy as np

なお,リーマン多様体に関するその他の選択肢は以下の通りである.

また,最適化手法に関するその他の選択肢は以下の通りである.

解くべき最適化問題の定義

問題を定義する方法について述べる.今回は行列 $A$ の左特異ベクトル,右特異ベクトルが既知となるように直交行列 $U_r \in \mathbb{R}^{m\times n}, \, V_r \in \mathbb{R}^{n\times n}$ をランダムに生成し,特異値を対角成分に持つ対角行列 $S$ は $S = \mathrm{diag} [n, n-1, \dots, 1]$ となるように定義する.これによって,アルゴリズムによって得られた解を最適解と比較することができる.

np.random.seed(2022)

# setting problem
m, n, p = 5, 4, 3
N = np.diag(np.arange(p, 0, -1))
Ur, _ = np.linalg.qr(np.random.normal(size=(m, n)))
Vr, _ = np.linalg.qr(np.random.normal(size=(n, n)))
S = np.diag(np.arange(n, 0, -1))
A = Ur @ S @ Vr.T

次に,解くべき最適化問題を Pymanopt の記法に従って記述する.まず,多様体は次のように定める.今回は,Stiefel 多様体の積多様体として記述されるリーマン多様体上での最適化問題を扱うため,次のように書くことができる.

# setting manifold
stiefel_1 = Stiefel(m, p)
stiefel_2 = Stiefel(n, p)
manifold = Product([stiefel_1, stiefel_2])

目的関数は次のように記述することができる.以前のバージョンと異なり,@pymanopt.function.autograd(manifold) といったデコレータが必要である.このデコレータを付けないと ValueError: Function 'cost' must be decorated with a backend decorator というエラーが返ってくる.

# define cost function
@pymanopt.function.autograd(manifold)
def cost(U, V):
    return -np.trace(U.T @ A @ V @ N)

以上の下で解くべき最適化問題Problem クラスのインスタンスとして次のように定める.引数にリーマン多様体と目的関数を与えればよい.

# define problem
problem = Problem(manifold=manifold, cost=cost)

最適化手法の定義

最適化問題を解くための最適化手法を定める.最適化アルゴリズム共役勾配法を扱うため ConjugateGradient クラスのインスタンスとして生成する.引数は以下の通りで,verbosity は出力されるログレベルを表し,0 が出力なしを表し,2 が最大限出力することを課す.詳しくは コード を確認されたい.

# define solver
solver = ConjugateGradient(
    max_time=100,
    max_iterations=100,
    min_gradient_norm=1e-6,
    min_step_size=1e-10,
    max_cost_evaluations=5000,
    verbosity=2,
    log_verbosity=0,
    beta_rule="PolakRibiere",
)

初期値は以下のように与える.初期値を自前で定義して与えることもできるが,pymanopt.manifolds から import したクラス(e.g., Stiefel, Product)は random_point というメソッドを持っており,これを使うと対象とするリーマン多様体上の変数をランダムに生成できる.

# initial guess
U0, _ = np.linalg.qr(np.random.rand(m, p))
V0, _ = np.linalg.qr(np.random.rand(n, p))
x0 = np.array([U0, V0])
# x0 = manifold.random_point()

最後に実行する.引数には,Problem クラスのインスタンスと初期値を与えればよい.

# run
res = solver.run(problem, initial_point=x0)
print(res)

出力内容

出力内容は以下の通り.反復回数 60 回で終了条件を満たし停止していることが確認できる.

Optimizing...
Iteration    Cost                       Gradient norm
---------    -----------------------    --------------
  1          +4.7770270621793562e+00    9.86333970e+00
  2          -4.1494305543650158e+00    9.87110012e+00
  3          -1.1198367922895978e+01    6.04524074e+00
  4          -1.3925445689443663e+01    3.76110809e+00
  5          -1.5989745840394253e+01    2.93949033e+00
(中略)
 56          -1.9999999999970512e+01    8.87200818e-06
 57          -1.9999999999987086e+01    7.52846533e-06
 58          -1.9999999999997979e+01    3.67079919e-06
 59          -1.9999999999999311e+01    1.58361854e-06
 60          -1.9999999999999545e+01    9.85344173e-07
Terminated - min grad norm reached after 60 iterations, 0.09 seconds.

OptimizerResult(point=[array([[-3.27751202e-04, -2.00772692e-01,  6.05821509e-02], 
       [ 1.75156668e-01,  7.33257570e-01, -5.71602349e-04],
       [ 2.31897948e-01,  5.11838715e-01,  4.82436752e-01],
       [ 7.09519651e-01, -3.96506778e-01,  4.94409705e-01],
       [ 6.41969813e-01,  5.31711091e-02, -7.20515996e-01]]), array([[-0.13814535, 
-0.50517496,  0.67497762],
       [-0.88146048, -0.32414775, -0.21981553],
       [ 0.14490487, -0.1186366 ,  0.56149162],
       [-0.4277217 ,  0.79098157,  0.42522171]])], cost=-19.999999999999545, iterations=60, stopping_criterion='Terminated - min grad norm reached after 60 iterations, 0.09 seconds.', time=0.09300017356872559, cost_evaluations=60, step_size=2.3126464564821429e-07, gradient_norm=9.853441734587333e-07, log={'optimizer': 'ConjugateGradient', 'stopping_criteria': {'max_time': 100, 'max_iterations': 100, 'min_gradient_norm': 1e-06, 'min_step_size': 1e-10, 'max_cost_evaluations': 5000}, 'optimizer_parameters': {'beta_rule': 'PolakRibiere', 'orth_value': inf, 'line_searcher': <pymanopt.optimizers.line_search.AdaptiveLineSearcher object at 0x0000022E9CA75C70>}, 'iterations': None})

今回は,行列 $A$ の構成の仕方から Pymanopt で得られた(近似)最適解の質を評価することができる.

def F_opt(S, N):
    p = N.shape[0]
    S = S[:p, :p]
    return -np.trace(S @ N)

U_star, V_star = res.point[0], res.point[1]

print('Optimal value obtained by conjugate descent: {}'.format(cost(U_star, V_star)))
print('Optimal value: {}'.format(F_opt(S, N)))

結果は以下の通り.

Optimal value obtained by conjugate descent: -19.999999999999545
Optimal value: -20

おわりに

本記事では,特異値分解を Stiefel 多様体の積集合上で最適化する問題について扱った.また,最適化問題 $(\mathrm{P})$ をリーマン多様体上の最適化を行うパッケージである Pymanopt を用いて求解する方法について述べた.

参考文献

[1] Sato, H. (2013). Riemannian Optimization Algorithms and Their Applications to Numerical Linear Algebra.
[2] Absil, P. A., Mahony, R., & Sepulchre, R. (2009). Optimization algorithms on matrix manifolds. In Optimization Algorithms on Matrix Manifolds. Princeton University Press.
[3] Nocedal, J., & Wright, S. J. (Eds.). (1999). Numerical optimization. New York, NY: Springer New York.
[4] Sato, H., & Iwai, T. (2013). A Riemannian optimization approach to the matrix singular value decomposition. SIAM Journal on Optimization, 23(1), 188-212.
[5] 【Python】リーマン多様体上の共役勾配法
[6] https://manabitimes.jp/math/1280
[7] https://pymanopt.org/docs/stable/

おまけ

最適化問題 $(\mathrm{P})$ を最急降下法共役勾配法で解いたときの勾配ノルムの推移を描いた結果を載せておく.

勾配ノルムの推移