この記事は「数理最適化 Advent Calendar 2020」の24日目の記事です．
Advent Calendar 初参戦です．今回は「数理最適化」がテーマということで，凸計画問題に対する近接勾配法について書こうと思います*1．連続最適化や機械学習に馴染みのある方にとっては新鮮味のない話題かもしれませんが，どうぞよろしくお願いします．
はじめに
数理最適化において，凸集合*2上で凸関数*3を最小化する*4問題は「凸計画問題」 あるいは「凸最適化問題」などと呼ばれる．線型計画問題，二次計画問題，半正定値計画問題，二次錐計画問題などは凸計画問題の代表的な例であり，馴染み深いと感じる方も多いと思う．凸計画問題は非線形最適化問題に分類されるが，局所的最適解が大域的最適解に一致するなどの性質上，理論・応用の両面から古くから研究されてきた最も重要な最適化問題のクラスの一つである．
凸計画問題において，制約条件が凸集合に限定されず空間全体*5であるとき，問題は次のように記述される：
\begin{align}
\min_{x\in\mathbb{R}^n} f(x)\tag{$\diamond$}
\end{align}
ここに，$f\colon\mathbb{R}^n\to\mathbb{R}$ は凸関数であり，一般に目的関数と呼ばれる．もし関数 $f$ が微分可能であるならば，後述する勾配降下法などのアルゴリズムによって最適解に収束する点列を構成することができる．
一方，目的関数が微分できない点を含む場合，関数の勾配の情報が使えないため，そのようなアルゴリズムを適用することができないことは想像に容易いだろう．例えば，機械学習などで頻出するLASSO（$\ell_1$ 正則化）における正則化項や，行列の補完問題（Matrix completion problem）における核ノルム（nuclear / trace norm）は微分できない点を含むことでよく知られている．
なお，本稿は数理最適化がテーマなので，$n$ 次元ベクトルが頻出する．多くの文献ではベクトルを $\mathbf{x}$ あるいは $\boldsymbol{x}$ と表記されることが多い*6が，本稿では太字にせず $x$ と表すこととする．また，本稿を通して，$\|\cdot\|_2$ はユークリッドノルム，$\|\cdot\|_1$ は $\ell_1$ ノルムを表す*7．行列 $X\in\mathbb{R}^{m\times n}$ に対して，フロベニウスノルム $\|X\|_{F}$ および核ノルム $\|X\|_{\mathrm{tr}}$ はそれぞれ
\begin{align}
\|X\|_{F} := \sqrt{\sum_{i=1}^{\min \{m,n\}} \left\{ \theta_i (X) \right\}^2 },\quad
\|X\|_{\mathrm{tr}} := \sum_{i=1}^{\min \{m,n\}} \theta_i (X)
\end{align}
で定義される．ここで，$\theta_i(X)$ は $X$ の特異値を表す*8．
最後に，指示関数（indicator function）を定義しておく．閉凸集合 $\mathcal{C}$ に対して，
\begin{align}
1_\mathcal{C}(x) := \left\{\begin{array}{ll}
0, &\quad \mathrm{if} \; x \in \mathcal{C}, \\
\infty, &\quad \mathrm{otherwise}
\end{array}\right.
\end{align}
と定義される関数 $1_\mathcal{C} \colon \mathbb{R}^n \to \mathbb{R} \cup \{\infty\}$ を集合 $\mathcal{C}$ の指示関数という*9．すなわち，入力ベクトル $x$ が集合 $\mathcal{C}$ に含まれていたら $0$ を，含まれていなかったら $\infty$ を返す関数である．この指示関数を用いると，制約条件を含んだ最適化問題 $\min_{x\in\mathcal{C}} f(x)$ が次のように等価に書き換えられる*10．
\begin{align}
\min_{x\in\mathbb{R}^n} f(x) + 1_{\mathcal{C}}(x)
\end{align}
問題設定
本稿で取り上げる最適化問題は以下のように定式化される問題である：
\begin{align}
\min_{x \in \mathbb{R}^n} f(x) := g(x) + h(x) \tag{$\ast$}
\end{align}
ここに，$g\colon\mathbb{R}^n\to\mathbb{R}$ は微分可能な凸関数，$h\colon\mathbb{R}^n\to\mathbb{R}$ は必ずしも微分可能とは限らない凸関数とする．
先に例に挙げたLASSOを例にとると，次のような対応関係になる．
\begin{align}
\min_{\beta\in\mathbb{R}^{n}}f(\beta) = \underbrace{\frac{1}{2} \| y - X \beta \|_2^2}_{g(\beta)} + \underbrace{\lambda \| \beta \|_1}_{h(\beta)}
\end{align}
また，行列の補完問題を例にとると，次のような対応関係になる．
\begin{align}
\min_{B\in\mathbb{R}^{m\times n}}f(B) = \underbrace{\frac{1}{2} \sum_{(i,j)\in\Omega} \left(Y_{ij} - B_{ij}\right)^2 }_{g(B)} + \underbrace{\lambda \| B \|_{\mathrm{tr}}}_{h(B)}
\end{align}
収束解析をする上では，関数 $g$ や $h$ にもう少し条件を付け加える必要があったりするのだが，今回は概説ということなので詳しく触れないこととする．
勾配降下法
問題($\diamond$)の目的関数が微分可能な凸関数である場合，次のような手続きによって，最適解に収束するような点列を生成することができる．
\begin{align}
x^{k+1} = x^{k} + t_k \cdot d_k, \quad k=0,1,2,\dots
\end{align}
ここで，$t_k$ は「ステップ幅」と呼ばれる正の実数であり，$d_k$ は降下方向を表すベクトルである．ベクトル $d_k$ として勾配ベクトル $-\nabla f(x^k)$ を採用したアルゴリズムが「最急降下法」に相当する．
少し違った見方をしてみよう．現在の点を $\overline{x}$ として，関数 $f$ の点 $\overline{x}$ の周りで Taylor 展開して得られる関数：
\begin{align}
f(x) \simeq f(\overline{x}) + \nabla(\overline{x})^\top (x-\overline{x}) + \frac{1}{2} (x-\overline{x})^\top \nabla^2 f(\overline{x}) (x-\overline{x})
\end{align}
は点 $\overline{x}$ の周りにおける関数 $f$ の近似した関数になっている．上式においてヘッセ行列 $\nabla^2 f(\overline{x})$ を $tI_n$ とおいた式：
\begin{align}
\tilde{f}(x) = f(\overline{x}) + \nabla f(\overline{x})^\top (x-\overline{x}) + \frac{1}{2t} \|x-\overline{x}\|_2^2
\end{align}
は $x$ に関する二次式であり，最小となる点が存在する．そこで，その最小となる点を求めて次の点として採用することを考えてみる．
具体的には，以下のように式変形をすればよい．
\begin{align}
x^{\star}
&= \mathop{\mathrm{arg~min}}_{x} \left\{ \tilde{f}(x) \right\} \\
&= \mathop{\mathrm{arg~min}}_{x} \left\{ f(\overline{x}) + \nabla f(\overline{x})^\top (x-\overline{x}) + \frac{1}{2t} \|x-\overline{x}\|_2^2 \right\} \\
&= \mathop{\mathrm{arg~min}}_{x} \left\{ 2t \cdot \nabla f(\overline{x})^\top x + \|x \|^2_2 - 2\overline{x}^\top x \right\} \\
&= \mathop{\mathrm{arg~min}}_{x} \left\{ \|x \|^2_2 - 2 (\overline{x} - t \cdot \nabla f(\overline{x}))^\top x \right\} \\
&= \mathop{\mathrm{arg~min}}_{x} \left\{ \|x - (\overline{x} - t \cdot \nabla f(\overline{x})) \|^2_2 \right\}
\end{align}
以上より，関数 $\tilde{f}$ を最小とする $x$ は以下のようになる：
\begin{align}
x^{\star} = \overline{x} - t \cdot \nabla f (\overline{x})
\end{align}
これは勾配法において勾配ベクトル $d$ を $-\nabla f(x)$ とした最急降下法に相当する．
近接勾配法
問題($\ast$)に対する近接勾配法は以下の更新式によって点列を更新するアルゴリズムである．
\begin{align}
x_{k+1} &= \textrm{prox}_{t_k,h} ( x_k - t_k \nabla g(x_k) ),\\
\textrm{prox}_{t,h}(y) &= \mathop{\mathrm{arg~min}}_{x} \left\{ h(x) + \frac{1}{2t} \|x-y\|_2^2 \right\}.
\end{align}
上式における $\textrm{prox}$ は近接作用素（proximal operator / mapping）と呼ばれる作用素*11である．意味合いとしては，現在地 $y$ の近くで関数を最小化する問題である．
本稿では上記の更新式がどのような経緯を辿って導出されたものなのか見ていきたい．
解くべき問題について，関数 $h$ は微分不可能な点を含むが，関数 $g$ については微分可能なので $g$ に関する勾配の情報のみを用いることを考える．先のアナロジーにより，現在の点 $x^k$ における関数 $f$ の二次式による近似を考えると，次のような手続きを考えたらよさそうである：
\begin{align}
x^{k+1} &= \mathop{\mathrm{arg~min}}_{x} \left\{ g(x^{k}) + \nabla g(x){^\top} (x^{k}-x) + \frac{1}{2t} \|x^{k}-x\|_2^2 + h(x) \right\}
\end{align}
この式を変形していくと次のようになる：
\begin{align}
x^{k+1} &= \mathop{\mathrm{arg~min}}_{x} \left\{ g(x^{k}) + \nabla g(x){^\top} (x-x^{k}) + \frac{1}{2t_k} \|x-x^{k}\|_2^2 + h(x) \right\} \\
&= \mathop{\mathrm{arg~min}}_{x} \left\{ \frac{1}{2t_k} \left( \| x - x^k \|_2^2 + 2t_k \nabla g (x^k)^\top x \right)+ h(x) \right\} \\
&= \mathop{\mathrm{arg~min}}_{x} \left\{ \frac{1}{2t_k} \left\| x - (x^k - t_k \nabla g (x^k)) \right\|_2^2 + h(x) \right\}
\end{align}
ここで，近接作用素の定義を用いると，
\begin{align}
x^{k+1} = \mathrm{prox}_{t_k,h} (x^k - t_k \nabla g (x^k))
\end{align}
となり，所望の更新式が得られた．
近接勾配法の更新式において，近接作用素 $\mathrm{prox}$ を逐次計算する必要がある．
この計算が重たいと効率的に解の更新を行うことができない．
本稿では応用上重要な以下の 3 つの関数に対する近接作用素が具体的にどのように構成されるのかについて見ていく．
まず，$\ell_1$ ノルムについて見ていく．関数 $h$ が $h(x) = \lambda\|x\|_1$ で与えられるとき，近接作用素は次のように表される．
\begin{align}
\mathrm{prox}_{t,h} (x) = \mathop{\mathrm{arg~min}}_{y\in\mathbb{R}^n} \left\{ \frac{1}{2t} \left\| x - y \right\|_2^2 + \lambda \|y\|_1 \right\}
\end{align}
簡単な変形により，上記の右辺を最小にするようなベクトル $y^\star$ は解析的に求めることができ，一般にソフト閾値作用素（soft-thresholding operator）と呼ばれる以下の関数と等価となる．
\begin{align}
\left[ \mathrm{prox}_{t,h} (x) \right]_i = \mathrm{sgn} (x_i) \max \{ |x_i| - \lambda , 0 \}
\end{align}
ここに，$\mathrm{sgn}(\cdot)$ は符号関数である．
次に，核ノルムについてである．関数 $h$ が $h(X) = \|X\|_{\mathrm{tr}}$ で与えられるとき，近接作用素は次のように表される．
\begin{align}
\mathrm{prox}_{t,h} (X) = \mathop{\mathrm{arg~min}}_{Y\in\mathbb{R}^{m\times n}} \left\{ \frac{1}{2t} \left\| X - Y \right\|_{F}^2 + \lambda \|Y\|_{\mathrm{tr}} \right\}
\end{align}
上式は行列 $X$ の特異値に対するソフト閾値作用素に帰着される．すなわち，行列 $X$ を $X = U\Theta V^\top$ と特異値分解したとき，近接作用素は
\begin{align}
\mathrm{prox}_{t,h} (X) = U \tilde{\Theta}_{\lambda} V^\top
\end{align}
で与えられる．ここで，$\tilde{\Theta}_{\lambda}$ は対角成分が $\left(\tilde{\Theta}_{\lambda} \right)_{ii} = \max \{ \theta_i (X) - \lambda , 0 \} , \; (i=1,2,\dots, \min \{m,n\} )$ で与えられる対角行列を表す．
最後に，指示関数について見ていこう．関数 $h$ が $h(x) = 1_{\mathcal{C}} (x)$ で与えられるとき，近接作用素は次のように表される．
\begin{align}
\mathrm{prox}_{t,h} (X) &= \mathop{\mathrm{arg~min}}_{y\in\mathcal{C}} \left\{ \frac{1}{2t} \left\| x - y \right\|_{2}^2 + 1_{\mathcal{C}}(y) \right\} \\
&= \mathop{\mathrm{arg~min}}_{y\in\mathcal{C}} \left\| x - y \right\|_{2}^2 =: P_{\mathcal{C}}(x)
\end{align}
要するに，入力ベクトル $x$ とのユークリッド距離が最小となるような $\mathcal{C}$ 内の点を見つける問題に相当し，一般に凸射影と呼ばれる写像 $P_{\mathcal{C}}$ に他ならない*12．
おわりに
本稿では，凸計画問題に対する近接勾配法に関する初歩的な内容について解説した．途中で断ったように，今回はアルゴリズムの収束解析について触れていない．そのため，目的関数の条件について凸性と微分可能性についてしか言及していない．様々な文献をみていくと，収束解析をするために，関数 $g$ にリプシッツ連続性や平滑性を仮定として置いている場合がある．その上で，リプシッツ定数などの情報を用いてステップ幅などを決定するアルゴリズムを構成し，収束解析がなされているようである．
Refference
本稿を記述する上で参考にした資料を以下に挙げる．より詳しい内容を知りたい場合は以下の文献を参照されたい．
更新ログ
- （12/24 8:18）指示関数の記述を整理しました。集合論などで出てくる指示関数と最適化の文脈で出てくる指示関数を混同してしまっていました。今回の場合は、入力ベクトルが集合に含まれていたら $0$ を返し、含まれていなかったら $\infty$ を返す関数の方が正しいです。指摘していただいた @weary_samurai さん、 @paruma184 さん、@jirohhana2 さん、ありがとうございました。
