Frank-Wolfeアルゴリズム

#はじめに
過学習を防ぐためには正則化を行います．例えばリッジ回帰ではパラメータの $L_2$ ノルムを抑えますし，LASSO ではこれが $L_1$ ノルムになります．いずれも，二乗誤差のような何かの損失関数 $f(\boldsymbol{x})$ を最小にする $\boldsymbol{x} ; {\rm s.t.} ; ||\boldsymbol{x}||_p \le R$ を探そうという問題に定式化されます．実際はノルムを罰則項としてくっ付けた $f(\boldsymbol{x})+\lambda ||\boldsymbol{x}||_p$ の最小化として定式化されることも多いですが，実はある $R$ に対応する $\lambda$ が存在し2つの定式化は等価です．

また行列 $X$ を変数とする推薦システムのような問題も，損失関数 $f(X)$ を最小にする $X$ を正則化の下で求める問題に定式化されます．例えば低ランク性が期待できる $M\times N$ 行列 $X^\star$ の一部の要素 $\Omega$ が観測できたとき，行列 $X^\star$ の近似を次式で求めます．
$$\arg\min_X f(X) := \sum_{(i,j) \in \Omega} (X_{i,j}-X^\star_{i,j})^2 \quad {\rm s.t.} \quad {\rm rank} X \le K$$
$X$ の低ランク性を担保し過学習を防ぐために ${\rm rank}X$ を制約しています．実際にこれを解く時には $X=UV^\top$ と分解した形で最適化問題を書き直し，$U$ と $V$ について交互に最適化を回す方法があり，変数の数が $(M+N)K$ で済む利点があります．しかし，こうするとランク制約のせいで非凸最適化問題になってしまい，局所解や初期値依存の問題があります．そこでランク制約の凸緩和として $X$ のトレースノルム $||X||_{\rm Tr}$（＝ $X$ の特異値の和）を ${\rm rank}X$ の代わりに制約することがあります．この方法では変数の数は $MN$ に増えますが，凸最適化問題になるので局所解と初期値依存が解消でき，得られる解も特異値の分布が疎＝低ランク性が期待できます．これは $L_0$ ノルム（非ゼロ要素の個数）の代わりに $L_1$ ノルムを抑えるような話です．

実行可能領域（例えば $L_p$ ノルム球）を $\mathcal{K}$ とすると，上述の問題はどれも制約付き凸最適化問題 $\arg\min_{\boldsymbol{x}\in\mathcal{K}} f( \boldsymbol{x} )$ の形で書けます．これを解くとき，罰則項をくっつけた形に変形して Proximal Gradient Descent を用いると Proximal operator を計算して射影する必要があります．$L_1$ 正則化や Group LASSO 等では射影は簡単に計算できますが，変数 $\boldsymbol{x}$ が行列 $X$ で，トレースノルム正則化する時は射影に特異値分解（SVD）が必要になり重い！ということがあります（後ほど詳しく述べます）．

しかし，$\nabla f$ が計算できて $\mathcal{K}$ 上での線形最適化が解けるなら， Frank-Wolfe アルゴリズムが使えるかもしれません．このアルゴリズムを使う機会があったため，本稿で調べたことをまとめました．

#定義
凸最適化問題とは，凸集合 $\mathcal{K}$ と凸関数 $f : \mathcal{K}\to \mathbb{R}$ について， $\arg\min_{\boldsymbol{x} \in \mathcal{K}} f(\boldsymbol{x})$ を求める問題です．凸性の定義は以下の通りです．

次式の条件が満たされるとき $\mathcal{K}$ は凸集合と呼ばれます．
$$\forall \boldsymbol{x}, \boldsymbol{y} \in \mathcal{K}, \forall \alpha \in [0,1]
\quad
\alpha \boldsymbol{x} + (1-\alpha) \boldsymbol{y} \in \mathcal{K}$$
直感的には，$\mathcal{K}$ 内の任意の2点を結んだ線分が $\mathcal{K}$ からは}み出ないと言っています．
また，次式の条件を満たすとき $f$ は凸関数です．
$$\forall \boldsymbol{x}, \boldsymbol{y} \in \mathcal{K}, \forall \alpha \in [0,1]
\quad
f(\alpha \boldsymbol{x} + (1-\alpha) \boldsymbol{y}) \le
\alpha f(\boldsymbol{x}) + (1-\alpha) f(\boldsymbol{y}) $$
こちらは，$f$ のグラフ上の任意の2点を結んだ線分が $f$ のグラフの上にあるということを言っています．

本稿では内積を $\cdot$ で書くことにします．行列同士の内積もベクトルの内積と同様に $X \cdot Y = \sum_{i,j} X_{i,j}Y_{i,j} = {\rm Tr}(XY^\top)$ で定義されます．

#Frank-Wolfe Algorithm
##アルゴリズムの説明
$f$ の勾配 $\nabla f(\boldsymbol{x})$ が計算でき，$\mathcal{K}$ 上での線形最適化問題 $\arg\min_{\boldsymbol{x} \in \mathcal{K}} ; \boldsymbol{x}\cdot \boldsymbol{g}$ は解けると仮定します．Frank-Wolfeアルゴリズムは，適当な初期値 $\boldsymbol{x}_0 \in \mathcal{K}$ から始めて，収束するまで以下を繰り返します．

\begin{align}
&{\rm initialize} \; \boldsymbol{x}_0 \in \mathcal{K} \\
&{\rm for} \; k=0,1,2,\dots \\
&\quad \boldsymbol{s}_{k} = \arg\min_{\boldsymbol{s} \in \mathcal{K}}
 \;\boldsymbol{s} \cdot \nabla f(\boldsymbol{x}_k) \\
&\quad \gamma_k = \frac{2}{2+k} \quad {\rm or} \quad \gamma_k = \arg\min_{\gamma \in [0,1]} \; f(\boldsymbol{x}_k  + \gamma(\boldsymbol{s}_k -\boldsymbol{x}_k ) )\\
&\quad \boldsymbol{x}_{k+1} = (1-\gamma_k) \boldsymbol{x}_k + \gamma_k \boldsymbol{s}_k \\
&{\rm end \; for}
\end{align}

$\boldsymbol{s}_k$ は $\nabla f(\boldsymbol{x}_k)$ との内積を最小にするように決めるので，$\nabla f(\boldsymbol{x}_k)$ とは逆=勾配降下方向で $\mathcal{K}$ の端にいます（[1]のp.1の図がわかりやすい）．$k$ 回目のループで現在地 $\boldsymbol{x}_k$ から勾配降下方向に動くとき， $\mathcal{K}$ に入っているハズの $\boldsymbol{x}_k$ と $\boldsymbol{s}_k$ の内分点を取るように更新することで，$x_{k+1}$ が $\mathcal{K}$ の外に出て行っちゃうこと（とそれに伴う射影）を回避しています．

##ステップ幅の決め方
ステップ幅 $\gamma_k$ は $2/(2+k)$ で決め打つ方法と line search して決める方法があります．Line search の式は一見ごついですが，簡単のために

\boldsymbol{x}_k+\gamma(\boldsymbol{s}-\boldsymbol{x}_k)=\boldsymbol{x}_\gamma

と書いて合成関数の微分をすると，

\begin{align}
\frac{\partial}{\partial \gamma} f(\boldsymbol{x}_\gamma)
=
\frac{\partial}{\partial \boldsymbol{x}_\gamma} f(\boldsymbol{x}_\gamma) \cdot \frac{\partial (\boldsymbol{x}_k  + \gamma(\boldsymbol{s} -\boldsymbol{x}_k ))}{\partial \gamma}
=
\nabla f(\boldsymbol{x}_\gamma) \cdot (\boldsymbol{s} -\boldsymbol{x}_k ) = 0
\end{align}

を満たす $\gamma$ を見つけることになり，見た目がちょっとすっきりします（実際に解けるかどうかは $\nabla f$ 次第ですし， $\gamma \in [0,1]$ も考慮しないといません）．次節で細かく説明しますが，わざわざ line search せずに $\gamma_k=2/(2+k)$ の決め打ちでも収束に関する理論保証はあり，いずれの場合も $O(1/k)$ です．ただ，[2]によると決め打ちは収束が遅く，可能なら line search が良さそうです．また line search だと $f(\boldsymbol{x}_k)$ が $k$ について単調減少するようになります．

#収束性
Frank-Wolfe アルゴリズムの収束について， $f$ の curvature constant $C_f$ というものが効いてきます． $C_f$ の定義は次式の通りです．

C_f =
\sup_{\boldsymbol{x},\boldsymbol{s} \in \mathcal{K}, \gamma \in [0,1], \boldsymbol{y} = \boldsymbol{x} + \gamma (\boldsymbol{s}-\boldsymbol{x})}
\frac{2}{\gamma^2}
(f(\boldsymbol{y}) - f(\boldsymbol{x}) - (\boldsymbol{y}-\boldsymbol{x})\cdot \nabla f(\boldsymbol{x}))

$f$ の点 $\boldsymbol{x}$ での一次近似は，点 $\boldsymbol{y}$ では $f$ とどのくらいズレるかな？という量（を $\gamma$ が小さいとこのズレも小さくなるので $1/\gamma^2$ したもの）です．値が大きいほど「$f$ は曲がっている」と言えます．

ちょっと実例で計算してみましょう．例えば $L_2$ ノルム制約

\mathcal{K} = \{ \boldsymbol{x} \in \mathbb{R}^N : ||\boldsymbol{x}||_2 \le D \}

の上で，
ベクトル $\boldsymbol{g}$ と対称行列 $H$ で定義される関数 $f(\boldsymbol{x}) = \boldsymbol{g} \cdot \boldsymbol{x} + \frac{1}{2} \boldsymbol{x} \cdot H \boldsymbol{x}$ の場合，$\nabla f(\boldsymbol{x}) = \boldsymbol{g} + H\boldsymbol{x}$ なので，

\begin{align}
f(\boldsymbol{y}) - f(\boldsymbol{x}) - (\boldsymbol{y}-\boldsymbol{x})\cdot \nabla f(\boldsymbol{x})
&=
\boldsymbol{g} \cdot (\boldsymbol{y} - \boldsymbol{x}) + \frac{1}{2} (\boldsymbol{y} \cdot H \boldsymbol{y} - \boldsymbol{x} \cdot H \boldsymbol{x})
 - (\boldsymbol{y}-\boldsymbol{x})\cdot(\boldsymbol{g} + H\boldsymbol{x}) \\
&=
  \frac{1}{2} (\boldsymbol{y} \cdot H \boldsymbol{y}
- 2\boldsymbol{y} \cdot H \boldsymbol{x}
 + \boldsymbol{x} \cdot H \boldsymbol{x} )  \\
&=
 \frac{1}{2} (\boldsymbol{y}-\boldsymbol{x})
 \cdot H (\boldsymbol{y}-\boldsymbol{x}) \\
 &=
  \frac{\gamma^2}{2} (\boldsymbol{s}-\boldsymbol{x})
  \cdot H (\boldsymbol{s}-\boldsymbol{x})
\le \gamma^2 D^2 \lambda_{\rm max} (H)
\end{align}

となり， $C_f \le D^2 \lambda_{\rm max} (H)$ です（$\lambda_{\rm max}(H)$ は $H$ の最大固有値）．

さてこの $C_f$ を使って，Frank-Wolfeアルゴリズムの収束率は，$\gamma_k=2/(2+k)$ とする時と line search する時，いずれの場合も

f(\boldsymbol{x}_k) - f(\boldsymbol{x}^\star) \le \frac{2C_f}{2+k}

となり，繰り返し回数について $O(1/k)$ となることが知られています[1]．Accelerate なしの Proximal Gradient Descent とオーダー的には一緒です．

#使いどころ
[1] でも書かれていますが，低ランク制約の凸緩和としてよく用いられるトレースノルム制約下での凸最小化問題

\arg\min_X f(X) \quad \rm{s.t.} \quad X \in
\mathcal{K} := \{X \in \mathbb{R}^{M\times N} \: : \:||X||_{\rm Tr} \le \tau \}

が Frank-Wolfe の効く一例です．冒頭で書いた行列の穴埋め問題のように，$f$ が観測できた要素との2乗誤差で定義されていれば Soft-Impute [3] で楽に解けるのですが，一般の凸関数のとき Proximal gradient を適用するとなかなか込み入ったことになり大変です．この場合は制約ではなく罰則項としてノルムを追加した $\arg\min_X f(X)+ \lambda ||X||_{\rm Tr}$ を考え（冒頭でも述べた通り，ある $\tau$ に対応する $\lambda$ があり，2つの最適化問題は等価です），更新式は以下の通りになります．

\begin{align}
X_{k+1} &= {\rm prox}_{\eta \lambda||\cdot||_{\rm Tr}}(X_k - \eta \nabla f(X_k)) \\
{\rm prox}_{\lambda ||\cdot||_{\rm Tr}}(Y) &=
{\rm argmin}_X \left(\lambda ||X||_{\rm Tr} + \frac{1}{2} || X - Y ||_{\rm Fr}^2 \right)
\end{align}

${\rm prox}$ を計算すると $X$ の小さい特異値を 0 にする計算である Singular Value Thresholding [4] が出てきますが，これには各ステップで SVD が必要になり，特に $X$ の疎性等を仮定できないとき計算量は $M,N$ に関して $O(NM^2)$ です（$N\le M$ と仮定）．一方，Frank-Wolfe だと $-\nabla f(X_k)$ の最大特異値に対応する特異ベクトル $\boldsymbol{u}_1,\boldsymbol{v}_1$ さえ求めることが出来ればよく，

S_k = \arg\min_{S\in\mathcal{K}} \; S \cdot \nabla f(X_k) = \tau \boldsymbol{u}_1\boldsymbol{v}_1^\top

と計算でき SVD よりは早く各ステップで $M,N$ に関して $O(NM)$ の計算時間で済みます[5]．色んな論文でサクッと「$S_k$ はこうなるよ！」と書いてあり，確かに何だかそれっぽいけどそんなに自明か？と思ったので証明してみます．$-\nabla f(X_k)$ の特異値分解を $-\nabla f(X_k)=U\Sigma V^\top = \sum_{i}^{N} \sigma_i \boldsymbol{u}_i \boldsymbol{v}_i^\top$ とします．

\begin{align}
S_k &= \arg\min_{S\in\mathcal{K}} \; S \cdot \nabla f(X_k) \\
&= \tau \arg\max_{S \in \mathcal{K}} \frac{1}{\tau} S \cdot (U\Sigma V^\top) \\
&= \tau \arg\max_{\{\bar S \;{\rm s.t.}\; ||\bar S||_{\rm Tr}\le 1\}}
 \left(
   \sigma_1 \boldsymbol{u}_1 \cdot \bar S \boldsymbol{v}_1 +
   \sigma_2 \boldsymbol{u}_2 \cdot \bar S \boldsymbol{v}_2 + \dots +
   \sigma_N \boldsymbol{u}_N \cdot \bar S\boldsymbol{v}_N
 \right) \\
 &= \tau \arg\max_{\{\bar S \;{\rm s.t.}\; ||\bar S||_{\rm Tr}\le 1 \}}
   (\sigma_1,\sigma_2,\dots,\sigma_N) \cdot
   (\boldsymbol{u}_1 \cdot \bar S \boldsymbol{v}_1 ,
    \boldsymbol{u}_2 \cdot \bar S \boldsymbol{v}_2 , \dots ,
    \boldsymbol{u}_N \cdot \bar S\boldsymbol{v}_N )
  \end{align}

$\bar S = S/\tau$ としました．こうするとベクトルの内積として書き下すことが出来ます．ここで，制約 $||\bar S||_{\rm Tr} \le 1$ の下でこの内積を最大にするには $\sigma_1$ に掛かっている $\boldsymbol{u}_1 \cdot \bar S \boldsymbol{v}_1$ を最大にするように $\bar S$ を決めればよく， $\bar S = \boldsymbol{u}_1\boldsymbol{v}_1^\top$ が導かれます．

もっと効果があるのが変数 $X$ がテンソルの時です．そもそもテンソル相手にランクやノルムをどう定義するのかという疑問があるのですが，[6] では Scaled latent trace norm なるややこしいノルムを正則化に使っています．これは，ざっくりと言うと $X$ を展開（unfolding）して得られる行列のトレースノルムを次元のサイズに応じてスケーリングしたものです．ノルムが最小になるように unfolding するので，このノルムの定義そのものがすでに最小化問題になっており，なにこれ怖い近寄らんとこ......という気持ちになりますが，実はこの制約下でうまい具合に線形最適化が解けて Frank-Wolfe が適用できるよという話 [7] があります．これはそれなりに重たい話になる（たぶん）のでまたの機会に書ければと思います．

#参考文献
[1] Martin Jaggi, Revisiting Frank–Wolfe: Projection-Free Sparse Convex Optimization, 2013
[2]Makoto Yamada, Wenzhao Lian, Amit Goyal, Jianhui Chen, Kishan Wimalawarne, Suleiman A Khan, Samuel Kaski, Hiroshi Mamitsuka and Yi Chang, Convex Factorization Machine for Regression, 2015
[3] Rahul Mazumder, Trevor Hastie, and Robert Tibshirani, Spectral Regularization Algorithms for Learning Large Incomplete Matrices, 2010
[4] Jian-Feng Cai, Emmanuel J. Candes and Zuowei Shen, A Singular Value Thresholding Algorithm for Matrix Completion, 2008
[5] Zeyuan Allen-Zhu, Elad Hazan, Wei Hu and Yuanzhi Li, Linear Convergence of a Frank-Wolfe Type Algorithm over Trace-Norm Balls
[6] Kishan Wimalawarne, Masashi Sugiyama and Ryota Tomioka, Multitask learning meets tensor factorization: task imputation via convex optimization, 2014
[7] Xiawei Guo, Quanming Yao and James T. Kwok, Efficient sparse low-rank tensor completion using the Frank-Wolfe algorithm,2017