はじめに
一般的な線形最小二乗法
$$
\text{argmin}_\mathbf{x} \bigl\vert \mathbf{y} - \mathrm{A} \mathbf{x} \bigr\vert^2
$$
の解 $\mathbf{x}$ を求めるというようなシーンはいたるところで出てきます。なおこの解は
$$
\mathbf{x}^\mathrm{opt} = \bigl(\mathrm{A}^\top \mathrm{A}\bigr)^{-1} \mathrm{A}^\top \mathbf{y}
$$
となります。実装する際は、既存のライブラリを用いることが多いと思います。例えば python の場合、
np.linalg.solve(A.T @ A, A.T @ y)
とすればよいでしょう。
（ここでは本題でないので触れませんが、逆行列 $(\mathrm{A^\top A})^{-1}$ を計算してから $\mathbf{A^\top y}$ との積を求めるより、直接解くほうが安定性・スピード両面で有利です。）
$\mathbf{x}$ を 次元 $n$ のベクトルだとすると、上記の方程式を解くには $\mathcal{O}(n^3)$ の計算量・メモリが必要です。そのため、次元 $n$ が大きくなると、この問題を解くための計算コストは非常に大きくなります。
このような線形問題を解く際のもう1つの選択肢は、最急降下法のような反復法を用いるものです。
後に述べるように、反復法は行列和と行列積のみで行うものがほとんどですが、行列和・行列積の計算量は $\mathcal{O}(n^2)$ なので、充分近い初期値から始めることができると計算量とメモリを節約することができます。
python の場合、 scipy の線形代数モジュール scipy.linalg に様々な反復アルゴリズムが実装されています。一般には各自でアルゴリズムを実装するよりもこういったパッケージを使う方が速度・安定性の点でよいことが多いですが、特殊な場合には自分で実装するほうがよいこともあるでしょう。
さらにこの方法は、Lasso
$$
\mathrm{argmin}_\mathbf{x}\Bigl\{
\bigl\vert \mathbf{y} - \mathrm{A} \mathbf{x} \bigr\vert ^2
+ \lambda \bigl\vert \mathbf{x} \bigr\vert_1
\Bigr\}
$$
のような、解析解のない問題にも応用することができるので、知っておいても良いと思います。
ここでは、最も簡単な最急降下法について述べたいと思います。特に、収束が保証される学習率をどのようにして決めるかなどを少し理論的に述べたいと思います。
線形問題を最急降下法で解く
ロス関数を
$$
\mathcal{L}(\mathbf{x}) = \frac{1}{2} \bigl\vert \mathbf{y} - \mathrm{A} \mathbf{x} \bigr\vert^2
$$
と定義すると、最急降下法は、以下の更新を繰り返すことになります。
$$
\mathbf{x} \gets \mathbf{x} - \alpha \frac{d\mathcal{L}(\mathbf{x})}{d\mathbf{x}}
$$
ここで、
$$
\frac{d\mathcal{L}(\mathbf{x})}{d\mathbf{x}} =
-(\mathbf{y} - \mathrm{A}\mathbf{x})^\top\mathrm{A}
$$
から、python での擬似コードは以下のようになるでしょう。
for iter in range(maxiter):
res = y - A @ x
x += alpha * (A.T @ res)
$\alpha$ をどうやって決めるかや、収束を早くするためにはどのような工夫があるかが問題となります。
最急降下法を導出する
直感的には最急降下法は最適解に収束しそうな気がしますが、それをもう少し掘り下げて考えてみましょう。
まず、初期値 $\mathbf{x}^{(t)}$ が手元にあるとします。
ロス関数を $\mathbf{x}^{(t)}$ の周りでTaylor展開すると、以下のようになります。
$$
\mathcal{L}(\mathbf{x}) = \mathcal{L}(\mathbf{x}^{(t)})
+ \frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}(\mathbf{x} - \mathbf{x}^{(t)})
+ \frac{1}{2}(\mathbf{x} - \mathbf{x}^{(t)})^\top
\mathrm{A^\top A} (\mathbf{x} - \mathbf{x}^{(t)})
$$
$\mathcal{L}(\mathbf{x})$ に よく似た 関数を最適化することで、新しい $\mathbf{x}^{(t+1)}$ を求めることにしましょう。
具体的には、以下の性質を満たす関数 $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ を考えます。
$\mathcal{L}(\mathbf{x}^{(t)}) = \mathcal{G}(\mathbf{x}^{(t)}; \mathbf{x}^{(t)})$ ：
$\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ は $\mathbf{x} = \mathbf{x}^{(t)}$ で $\mathcal{L}(\mathbf{x})$ と交わる$\mathcal{L}(\mathbf{x}) \le \mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ ：
その他の点では、$\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ は $\mathcal{L}(\mathbf{x})$ より上にある
図にすると以下のような状態です。
このような関数 $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ を最小化する $\mathbf{x}$ は、
必ず元のロス関数も小さくすることがわかります。
そのため、
- $\mathcal{L}(\mathbf{x})$ を上から抑える $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ を求める
- $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ を最小化する $\mathbf{x}$ を求めて $\mathbf{x}^{(t+1)}$ とする
を繰り返すことができれば、最終的には元のロス関数を最小化することができるでしょう。
補助関数で上から抑える
この性質を満たす $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ のうち、以下のような形式のものを探すことにしましょう。
$$
\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)}) =
\mathcal{L}(\mathbf{x}^{(t)})
+
\frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}(\mathbf{x} - \mathbf{x}^{(t)})
+
\frac{\rho}{2}\bigl\vert\mathbf{x}-\mathbf{x}^{(t)}\bigr\vert^2
$$
この形式では、
$\mathcal{L}(\mathbf{x}^{(t)}) = \mathcal{G}(\mathbf{x}^{(t)}; \mathbf{x}^{(t)})$
は自明です。
では、どのような $\rho$ が $\mathcal{L}(\mathbf{x}) \le \mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ を満たすでしょうか。
しかし $\rho$ の条件を求める前に、まずは $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ の最適値を求めましょう。
上で示した $\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)})$ は、簡単に平方完成できる形式になっています。
$$
\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)}) =
\frac{\rho}{2} \Bigl\vert
\mathbf{x} - \bigl(
\mathbf{x}^{(t)} - \frac{1}{\rho} \frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}
\bigr)\Bigr\vert^2
$$
そのため、この補助関数を用いた更新式は以下のようになります。
$$
\mathbf{x}^{(t+1)} \gets
\mathbf{x}^{(t)} - \frac{1}{\rho} \frac{d\mathcal{L}(\mathbf{x}^{(t)})}{d\mathbf{x}}
$$
この形は、 $\alpha = \frac{1}{\rho}$ を学習率とする再急降下法の更新式そのものであることがわかります。
学習率の条件を求める
ここで、2つの関数の差を見てみましょう。
$$
\mathcal{G}(\mathbf{x}; \mathbf{x}^{(t)}) - \mathcal{L}(\mathbf{x}) =
\frac{1}{2}
(\mathbf{x} - \mathbf{x}^{(t)})^\top (\alpha \mathrm{I} - \mathrm{A^\top A})
(\mathbf{x} - \mathbf{x}^{(t)})
$$
この値が全ての $\mathbf{x}$ で0以上となるためには、行列 $(\alpha \mathrm{I} - \mathrm{A^\top A})$ が(半)正定値（全ての固有値がゼロ以上）であればよいことがわかります。つまり、 $\mathrm{A^\top A}$ の最大の固有値が $\rho$ として取りうる最小の値になります。
ただし、固有値分解も基本的には $\mathcal{O}(n^3)$ の計算量が必要なので、もっと簡便に最大の固有値の大きさを見積もることが必要です。その上限を簡単に見積もる方法としてゲルシュゴリンの定理(Gershgorin circle theorem)があります。
この定理から、以下の関係がわかります。
$$
\lambda \le \max_j \sum_i\bigl\vert{\mathrm{A^\top A}}_{ij}\bigr\vert
$$
$\mathrm{A^\top A}$ の要素の絶対値を列方向に総和をとったもののうちの最も大きいものより、固有値は必ず小さくなります。python の擬似コードは以下のようになります。
def supremum_eigen(A):
return np.max(np.sum(np.abs(A), axis=0))
そのため、学習率 $\alpha$ に $\frac{1}{\max_j\sum_i\bigl\vert{\mathrm{A^\top A}}_{ij}\bigr\vert}$ を採用することが1つの目安になるでしょう。
まとめ
この記事では、線形問題における最急降下法の導出と、確実に解が収束するための学習率の条件を求めました。
なお今回紹介した方法は、確実に解に収束することは保証されていますが、収束速度は一般に早いとは言えないかもしれません。$\mathrm{A}$ が単位行列に近ければ収束は速いですが、単位行列から異なる行列になればなるほど、つまり固有値のばらつきが大きければ大きいほど、収束は遅くなります。
これらを簡便に改善する方法に、スケーリングと加速法があります。これらについては次回まとめたいと思います。