ノンパラメトリックモデルのためのオンライン学習（理論編）

## 初めに

オンライン学習（online learning）に関する多くの文献は，パラメトリックモデル（parametric model）に焦点を当てています．そこで，理論編，実装編，実験編の三編に渡って，ノンパラメトリックモデル（nonparametric model）のためのオンライン学習について紹介します．取り上げる論文は，次の三つです．

- Kivinen, J., Smola, A. J., and Williamson, R. C. [Online learning with kernels.](https://papers.nips.cc/paper/2128-online-learning-with-kernels) In Advances in Neural Information Processing Systems (NIPS) 14, pp. 785-792, 2001.
- Shalev-Shwartz, S., Singer, Y., and Srebro, N. [Pegasos: Primal Estimated sub-GrAdient SOlver for SVM.](http://doi.acm.org/10.1145/1273496.1273598) In Proceedings of the 24th International Conference on Machine Learning (ICML), pp. 807-814, 2007.
- Wang, Z., Crammer, K. and Vucetic, S. [Breaking the curse of kernelization: Budgeted stochastic gradient descent for large-scale svm training.](http://www.jmlr.org/papers/v13/wang12b.html) Journal of Machine Learning Research (JMLR), Vol. 13, pp. 3103-3131, 2012.

## 問題設定

実数の集合を $\mathbb{R}$ とします．また，太字の小文字で記した $\boldsymbol{x} \in \mathbb{R}^{d}$ を $d$ 次元列ベクトル，その第 $i$ 成分を $x_{i}$ と表します．そして，上付きの $\top$ をベクトルの転置とします．したがって，$\boldsymbol{x}^{\top}$ は，行ベクトルとなります．

$\mathcal{X}$ を入力空間（input space），$\mathcal{Y}$ を出力集合，$\mathcal{H}$ を仮説集合（hypothesis set）とします．教師あり学習は，入力（input） $\boldsymbol{x} \in \mathcal{X}$ と出力（output） $y \in \mathcal{Y}$ を関係付ける仮説（hypothesis） $f: \mathcal{X} \to \mathcal{Y}$ を訓練事例集合（training instances） $\mathcal{S} \subseteq \mathcal{X} \times \mathcal{Y}$ から求める問題で，回帰（regression）と判別（classification）に大きく分けられます．

```math
\begin{align}
\mathcal{S} = \{ (\boldsymbol{x}_{t}, y_{t}) \}_{t = 1}^{n}
\end{align}
```

以下，説明を簡潔にするため，出力 $y$ が連続値を取る回帰について議論します．

```math
\begin{align}
\mathcal{X} &= \mathbb{R}^{d} \\
\mathcal{Y} &= \mathbb{R}
\end{align}
```

はじめに，仮説集合 $\mathcal{H}$ として線形モデル（linear model）を仮定します．

```math
\begin{align}
\mathcal{H} = \{ f (\boldsymbol{x}) = \boldsymbol{w}^{\top} \phi (\boldsymbol{x}) + b \mid \boldsymbol{w} \in \mathbb{R}^{p}, b \in \mathbb{R} \}
\end{align}
```

ここで，$\phi:\mathcal{X} \to \mathbb{R}^{p}$ は，任意の非線形関数です．

背後にある確率分布 $P$ は，未知であるため，経験損失（empirical loss） $R_{\mathrm{emp}} (\boldsymbol{w}, b)$ を最小化することで，最適解 $(\boldsymbol{w}^{\star}, b^{\star})$ を得ます．ただし，訓練事例数 $n$ に対して次元 $p$ が大きいと，汎化性能が劣化する過剰適合（overfitting）を起こすことが知られています．そこで，経験損失 $R_{\mathrm{emp}} (\boldsymbol{w}, b)$ に罰則（penalty） $\Phi (\boldsymbol{w})$ を加えた正則化損失（regularized loss） $ R_{\mathrm{reg}} (\boldsymbol{w}, b)$ を最小化することにします．最も単純な罰則 $\Phi (\boldsymbol{w})$ は，パラメータ $\boldsymbol{w}$ の L2 ノルムの二乗です．

```math
\begin{align}
(\boldsymbol{w}^{\star}, b^{\star}) &= \mathrm{argmin} R_{\mathrm{reg}} (\boldsymbol{w}, b) \\
R_{\mathrm{reg}} (\boldsymbol{w}, b) &= R_{\mathrm{emp}} (\boldsymbol{w}, b) + \lambda \Phi (\boldsymbol{w}) \\
R_{\mathrm{emp}} (\boldsymbol{w}, b) &= \frac{1}{n} \sum_{t = 1}^{n} l (f (\boldsymbol{x}_{t}), y_{t}) \\
\Phi (\boldsymbol{w}) &= \frac{1}{2} \| \boldsymbol{w} \|^{2}
\end{align}
```

ここで， $l: \mathcal{Y} \times \mathcal{Y} \to \mathbb{R}$ は，任意の損失関数（loss function），$\lambda \geq 0$ は，正則化の強弱を調節する正則化係数（regularization parameter）です．

### squared loss function

二乗損失関数（squared loss function） $l_{\mathrm{squared}}$ は，最小二乗法（least squares method）で使用される損失関数です．

```math
\begin{align}
l_{\mathrm{squared}} (f (\boldsymbol{x}), y) = \frac{1}{2} (f (\boldsymbol{x}) - y)^{2}
\end{align}
```

### epsilon-insensitive loss function

$\varepsilon$-許容損失関数（epsilon-insensitive loss function） $l_{\varepsilon\mathrm{-insensitive}}$ は，サポートベクトルマシン (Support Vector Machine) で使用される損失関数です．

```math
\begin{align}
l_{\varepsilon\mathrm{-insensitive}} (f (\boldsymbol{x}), y) &= \max (0, \left | f (\boldsymbol{x}) - y \right | - \varepsilon )
\end{align}
```

### Huber loss function

Huber 損失関数（Huber loss function） $l_{\mathrm{huber}}$ は，外れ値に頑健な損失関数です．

```math
\begin{align}
l_{\mathrm{huber}} (f (\boldsymbol{x}), y) &= \begin{cases}
\frac{1}{2} (f (\boldsymbol{x}) - y)^{2} & \textrm{if} \quad \left |f (\boldsymbol{x}) - y \right | \leq \sigma \\
\sigma \left | f (\boldsymbol{x}) - y \right | - \frac{1}{2} \sigma^{2} & \textrm{otherwise}
\end{cases}
\end{align}
```

## 最適化手法

オンライン学習の場合，時刻 $t$ に事例（instance） $(\boldsymbol{x}\_{t}, y\_{t}) \in \mathcal{X} \times \mathcal{Y}$ を観測する度，パラメータ $(\boldsymbol{w}, b)$ を逐次更新します．そこで，時刻 $t$ に観測した事例 $(\boldsymbol{x}\_{t}, y\_{t})$ に対する損失 $R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}\_{t}, y\_{t})$ を次のように定義します．

```math
\begin{align}
R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}_{t}, y_{t}) = l (f (\boldsymbol{x}_{t}), y_{t}) + \lambda \Phi (\boldsymbol{w})
\end{align}
```

### Stochastic Gradient Descent

損失関数 $l$ が凸関数である場合，時刻 $t$ に観測した事例 $(\boldsymbol{x}\_{t}, y\_{t})$ に対する損失 $R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}\_{t}, y\_{t})$ の最急降下方向へパラメータ $(\boldsymbol{w}, b)$ を逐次更新すれば，いずれ最適解 $(\boldsymbol{w}^{\star}, b^{\star})$ に収束します．この最適化手法は，Stochastic Gradient Descent (SGD) と呼ばれています．

```math
\begin{align}
\boldsymbol{w}_{t + 1} &= \boldsymbol{w}_{t} - \eta_{t} \left . \frac{\partial R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}_{t}, y_{t})}{\partial \boldsymbol{w}}  \right　|_{\boldsymbol{w} = \boldsymbol{w}_{t}} \\
b_{t + 1} &= b_{t} - \eta_{t} \left . \frac{\partial R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}_{t}, y_{t})}{\partial b}  \right　|_{b = b_{t}}
\end{align}
```

ここで，$0 < \eta_{t} < \frac{1}{\lambda}$ は，勾配降下の歩幅を表した学習率（learning rate）です．

時刻 $t$ に観測した事例 $(\boldsymbol{x}\_{t}, y\_{t})$ に対する損失 $R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}\_{t}, y\_{t})$ の勾配は，次のように計算することができます．

```math
\begin{align}
\frac{\partial R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}_{t}, y_{t})}{\partial \boldsymbol{w}} &= \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial \boldsymbol{w}} + \lambda \frac{\partial \Phi (\boldsymbol{w})}{\partial \boldsymbol{w}} \\
\frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial \boldsymbol{w}} &= \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f (\boldsymbol{x}_{t})} \frac{\partial f (\boldsymbol{x}_{t})}{\partial \boldsymbol{w}} \\
&= \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f (\boldsymbol{x}_{t})} \phi (\boldsymbol{x}_{t}) \\
\frac{\partial \Phi (\boldsymbol{w})}{\partial \boldsymbol{w}} &= \boldsymbol{w} \\
\frac{\partial R_{\mathrm{inst}} (\boldsymbol{w}, b, \boldsymbol{x}_{t}, y_{t})}{\partial b} &= \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial b} + \lambda \frac{\partial \Phi (\boldsymbol{w})}{\partial b} \\
\frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial b} &= \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f (\boldsymbol{x}_{t})} \frac{\partial f (\boldsymbol{x}_{t})}{\partial b} \\
&= \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f (\boldsymbol{x}_{t})} \\
\frac{\partial \Phi (\boldsymbol{w})}{\partial b} &= 0
\end{align}
```

以上より，パラメータ $(\boldsymbol{w}, b)$ に関する更新式を変形することができます．

```math
\begin{align}
\boldsymbol{w}_{t + 1} &= (1 - \eta_{t} \lambda) \boldsymbol{w}_{t} - \eta_{t} \left . \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f (\boldsymbol{x}_{t})} \right|_{f = f_{t}} \phi (\boldsymbol{x}_{t}) \\
b_{t + 1} &= b_{t} - \eta_{t} \left . \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f (\boldsymbol{x}_{t})} \right|_{f = f_{t}}
\end{align}
```

通常，訓練事例集合を $\mathcal{S}$ を一巡しても，パラメータ $(\boldsymbol{w}, b)$ は，滅多に収束しないため，何巡も学習する必要があります．

### Naive Online R Minimization Algorithm (Kivinen et al., 2001)

SGD は，パラメトリックモデルのための最適化手法でした．続いて，SGD をノンパラメトリックモデルのための最適化手法に拡張した Naive Online R Minimization Algorithm (Norma) を紹介します．

表現定理（representer theorem）は，「最適解 $\boldsymbol{w}^{\star}$ は，$\phi (\boldsymbol{x})$ の線形結合の形で書くことができる」という定理です．したがって，時刻 $t$ におけるパラメータ $\boldsymbol{w}_{t}$ を次のように表します．表現定理の証明は，[赤穂昭太郎著「カーネル多変量解析」](https://www.amazon.co.jp/%E3%82%AB%E3%83%BC%E3%83%8D%E3%83%AB%E5%A4%9A%E5%A4%89%E9%87%8F%E8%A7%A3%E6%9E%90%E2%80%95%E9%9D%9E%E7%B7%9A%E5%BD%A2%E3%83%87%E3%83%BC%E3%82%BF%E8%A7%A3%E6%9E%90%E3%81%AE%E6%96%B0%E3%81%97%E3%81%84%E5%B1%95%E9%96%8B-%E3%82%B7%E3%83%AA%E3%83%BC%E3%82%BA%E7%A2%BA%E7%8E%87%E3%81%A8%E6%83%85%E5%A0%B1%E3%81%AE%E7%A7%91%E5%AD%A6-%E8%B5%A4%E7%A9%82-%E6%98%AD%E5%A4%AA%E9%83%8E/dp/4000069713)に譲ります．

```math
\begin{align}
\boldsymbol{w}_{t} = \sum_{i = 1}^{t - 1} \alpha_{i}^{(t)} \phi (\boldsymbol{x}_{i})
\end{align}
```

すると，時刻 $t$ における仮説 $f_{t}$ は，次のようになります．

```math
\begin{align}
f_{t} (\boldsymbol{x}) &= \sum_{i = 1}^{t - 1} \alpha_{i}^{(t)} \phi (\boldsymbol{x}_{i})^{\top} \phi (\boldsymbol{x}) + b_{t} \\
&= \sum_{i = 1}^{t - 1} \alpha_{i}^{(t)} k (\boldsymbol{x}_{i}, \boldsymbol{x}) + b_{t}
\end{align}

```

ここで，カーネル法（kernel method）を用いました．$\phi (\boldsymbol{x})$ の内積を置き換えることで，非線形関数 $\phi$ に関する計算を明示的に経ることなく計算量を抑えることができます．$k:\mathcal{X} \times \mathcal{X} \to \mathbb{R}$ は，任意のカーネル関数（loss function）です．

以上より，パラメータ $\boldsymbol{w}$ に関する更新式を係数 $\alpha_{i}$ に関する更新式に変形することができます．

```math
\begin{align}
\alpha_{i}^{(t + 1)} = \begin{cases}
(1 - \eta_{t} \lambda) \alpha_{i}^{(t)} & i < t \\
- \eta_{t} \left . \frac{\partial l (f (\boldsymbol{x}_{t}), y_{t})}{\partial f(\boldsymbol{x}_{t})} \right|_{f = f_{t}} & i = t
\end{cases}
\end{align}
```

学習率は，次のように減衰させます．

```math
\begin{align}
\eta_{t} = \frac{\eta}{\sqrt{t}}
\end{align}
```

特に，係数 $\alpha_{i}$ が非零となる入力 $\boldsymbol{x}_{i}$ をサポートベクトル（support vector）と言います．

### Primal Estimated sub-GrAdient SOlver for SVM (Shalev-Shwartz et al., 2007)

Primal Estimated sub-GrAdient SOlver for SVM (Pegasos) は，Norma と同じく，SGD の派生手法です．

次のように学習率を減衰させることで，Norma より速く収束することができます．

```math
\begin{align}
\eta_{t} = \frac{1}{\lambda t}
\end{align}
```

### Budgeted SGD (Crammer and Vucetic, 2012)

Norma や Pegasos は，学習が進むにつれてサポートベクトルの数が増え，計算が困難になるという問題を抱えています．この問題は，カーネル化の呪い（curse of kernelization）と呼ばれています．

そこで，予算（budget） $B > 0$ を定め，学習中にサポートベクトルの数を調整します．最も単純な調整方法は，古い，または，係数の絶対値が小さいサポートベクトルから順に破棄（removal）することです．その他の調整方法として，射影（projection），併合（merging）があります．この枠組みは，Budgeted SGD (BSGD) と呼ばれています．

## 終わりに

次回は，BSGD を Python で実装するノンパラメトリックモデルのためのオンライン学習（実装編）をお届けします．