pythonで「異常検知と変化検知」を読む5 サポートベクトルデータ記述法による異常検知

「異常検知と変化検知」は異常検知の本です。アルゴリズム部分をpythonで実装していきたいと思います。たくさんの方が同じ内容をより良い記事でとして投稿しています。

個人的な勉強のメモ書きとなります。
私には難しい内容であり、コードがかなり汚いです。申し訳ありません。
まとめることが目的ではないので本文について参考書とほぼ同じ表現となっています。問題があればお教えください。
興味を持ったり、導出や詳細などを知りたい方は「異常検知と変化検知」を読んでいただければと思います。

参考

前回

混合分布モデルによる逐次更新型異常検知

サポートベクトルデータ記述法による異常検知

データを囲む最小の球

$D={\boldsymbol{x}^{(1)},\cdots,\boldsymbol{x}^{(N)}}$のような、ラベルが与えられていないデータを考える。
$D$の中には異常標本が含まれていないか、含まれていたとしても圧倒的少数であるとする。
「標本集合のほぼ全体を囲む球を作り、その球に入りきらなかったものを異常と判定する」という考え方は自然である。

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(1)

# 参考書内の図を作成

def make_circle(x, y, r, theta):
    """円の座標を出力
    """
    return np.concatenate([[x + r * np.cos(theta)], [y + r * np.sin(theta)]], axis=0).T

def find_near_k(x1, x2, k):
    """近傍点と距離の計算(近傍法の使いまわし)
    """
    dist = np.sum((x1 - x2)**2, axis=1)
    near = x1[np.argsort(dist)[:k]]
    dist = dist[np.argsort(dist)[k-1]]
    return near, dist

# データの準備
x1 = np.random.rand(20, 2)
x2 = np.mean(x1, axis=0)

plt.scatter(x1[:,0], x1[:,1]);
plt.scatter(x2[0], x2[1], color=['red']);

# 近傍点とそれを囲む円の描画
near, dist = find_near_k(x1, x2=x2, k=20)
theta = np.linspace(0, 2*np.pi, 100)
c = make_circle(x2[0], x2[1], r=np.sqrt(dist), theta=theta)
plt.plot(c[:,0], c[:,1], '--', color='red');

データを含むという条件のもとで、できる限り小さい球を求める。
これは次のような最適化問題で表せる。
$$
\min_{R^2,\boldsymbol{b}}{R^2}\ ただし条件\ |\boldsymbol{x}^{(n)}-\boldsymbol{b}|^2\leq R^2
$$
ここで、$n=1,\cdots,N$で$|\boldsymbol{a}|^2=\boldsymbol{a}^T\boldsymbol{a}$であり、$R$は球の半径、$\boldsymbol{b}$は中心の位置ベクトルを表す。
$\boldsymbol{x}^{(n)}$に対して、半径2乗に$u^{(n)}$分だけの「遊び」を許したうえで次の問題を解くことにする。

\min_{R^2,\boldsymbol{b},\boldsymbol{u}}\biggl\{R^2+C\sum_{n=1}^Nu^{(n)} \biggr\}\ subject\ to\ |\boldsymbol{x}^{(n)}-\boldsymbol{b}|^2\leq R^2+u^{(n)}

$C$は正則化定数と呼ばれる。
また、$u^{(n)}\geq 0$でなければならない。
解が${R^{2*},b^*,u^*}$と求まっているとする。
直感的に異常度を「球からはみ出した長さ」として

a(\boldsymbol{x}')=|\boldsymbol{x}'-\boldsymbol{b}^*|＾2-R^{2*}

のように定義できる。
2乗を展開して、

a(\boldsymbol{x}')=K(\boldsymbol{x}',\boldsymbol{x}')-2K(\boldsymbol{b}^*,\boldsymbol{x}')+K(\boldsymbol{b}^*,\boldsymbol{b}^*)-R^{2*}

と表しておく。
ただし、$K(・,・)$は引数同士の内積を表す。(あとで拡張した定義を与える。)
最適化問題を解くと、異常と正常を分ける球が、ごく少数の訓練標本で表されるという結果になる。
それが、サポートベクトルデータ記述法(支持ベクトルデータ記述法) と呼ばれる由来である。
この手法は、1クラスサポートベクトルマシンとも呼ばれる。

双対問題への変換とカーネルトリック

非線形制約を持つ最適化問題を「双対問題」に変換することで、非線形制約を消去する手法がある。
それを使えないか検討する。
ラグランジュ関数$L$を次のように定義する。

\begin{align}
L(R^2, \boldsymbol{b}, \boldsymbol{u}, \boldsymbol{\alpha}, \boldsymbol{\beta})&=R^2+C\sum_{n=1}^Nu^{(n)}-\sum_{n=1}^N\beta_nu^{(n)}-\sum_{n=1}^N\alpha_n\bigl\{R^2+u^{(n)}-|\boldsymbol{x}^{(n)}-\boldsymbol{b}|^2 \bigr\}
\end{align}

もとの変数$R^2$、$\boldsymbol{b}$、$\boldsymbol{u}$について最大化するため以下の計算を行う。

0=\frac{\partial L}{\partial R^2}=1-\sum_{n=1}^N\alpha_n\\
\boldsymbol{0}=\frac{\partial L}{\partial \boldsymbol{b}}=2\sum_{n=1}^N\alpha_n\boldsymbol{b}-2\sum_{n=1}^N\alpha_n\boldsymbol{x}^{(n)}\\
0=\frac{\partial L}{\partial u^{(n)}}=C-\beta_n-\alpha_n

これを使って$R^2$、$\boldsymbol{b}$、$\boldsymbol{u}$を消去する。

\begin{align}
l(\boldsymbol{\alpha},\boldsymbol{\beta})&=\min_{R^2,\boldsymbol{b},\boldsymbol{u}}L(R^2, \boldsymbol{b}, \boldsymbol{u}, \boldsymbol{\alpha}, \boldsymbol{\beta})\\
&=\sum_{n=1}^N\alpha_nK_{n,n}-\sum_{n=1}^N\alpha_n\alpha_{n'}K_{n,n'}
\end{align}

ただし、$K_{n,n'}=K(\boldsymbol{x}^{(n)},\boldsymbol{x}^{(n')})$と略記している。
ラグランジュ乗数$\boldsymbol{\alpha},\boldsymbol{\beta}$については、KKT条件から、非負条件$\alpha_n\geq0,\beta_n\geq0$が付される。
明らかに
$$
0\leq\beta_n=C-\alpha_n
$$
が成り立つので、$\beta_n$は$\alpha_n$が求まれば決定できる。
$\beta_n$に対する非負条件から$\alpha_n\leq C$という条件を付す。
以上まとめると、解くべき問題は、

\max_{\alpha}\biggl\{\sum_{n=1}^N\alpha_nK_{n,n}-\sum_{n,n'=1}^N\alpha_n\alpha_{n'}K_{n,n'} \biggr\}\\
subject\ to\ 0\leq\alpha_n\leq C(n=1,\cdots,N)

この双対問題は、通常の二値分類問題に使われる支持ベクトル分類器の双対問題と本質的に同じである。
通常、LIBSVMなどの既存のプログラムを用いて数値計算を行う。
具体的な最適化手法としては、SMO法と双対座標降下法が代表的である。

内積$K$を、たとえばRBF(動径基底関数)カーネルを使って、
$$
K(\boldsymbol{x},\boldsymbol{x'})←\exp{-\sigma|\boldsymbol{x}-\boldsymbol{x'}|^2}
$$
のように置き換えたと考える。
これは、元の座標を、内積が上記の関数で与えるような座標に非線形変換したとに対応している。
内積を与えることで非線形変換をしたことにする考え方をカーネルトリックと呼ぶ。
内積として与える関数$K(\boldsymbol{x},\boldsymbol{x'})$のことをカーネル関数と呼ぶ。

解の性質と分類

最適解$\alpha^*$が求まったら、KKT条件を逆に使い、元の解を求める。
まず、

\boldsymbol{b}^*=\sum_{n=1}^N\alpha_n^*\boldsymbol{x}^{(n)}

が得られる。
また、KKT条件の式に対応して、次の2種類の式が成り立つ。

\begin{align}
\alpha_n\bigl\{R^2+u^{(n)}-|\boldsymbol{x}^{(n)}-\boldsymbol{b}|^2 \bigr\}&=0\\
(C-\alpha_n^*)u^{(n)}&=0
\end{align}

この条件式から、以下の3つの場合があることが分かる。

1. $\alpha^*$が$0$でも$C$でもない場合、$u^{(n)}=0$となり、$R^{2*}=|\boldsymbol{x}^{(n)}-\boldsymbol{b}^*|^2$が成り立つ。これは球面の上に乗っていることを意味する。支持ベクトルやサポートベクトルと呼ばれる。
2. $\alpha^*=0$となる場合、$u^{(n)}=0$となり、$R^{2*}-|\boldsymbol{x}^{(n)}-\boldsymbol{b}^*|^2<0$が成り立つ。球の内部にあることを意味する。
3. $C<1$に対して$\alpha^*_n=C$となる場合、$u^{(n)}>0$となり、球の外にあることを意味する。

$C$を1未満の正値にすることにより、球の外に出る標本の数を制御できる。

$0<\alpha_n^*<C$を満たす$n$を1つ選びそれを$n'$とすると、$R^{2*}$は$|\boldsymbol{x}^{(n')}-\boldsymbol{b}^*|^2$とるが、ない席を展開すると、
$$
R^{2*}=K_{n',n'}+\sum_{n_1,n_2=1}^N\alpha_{n_1}^*\alpha_{n_2}^*K_{n_1,n_2}-2\sum_{n=1}^N\alpha_n^*K_{n,n'}
$$
と表せる。
1クラスサポートベクトル分類器における異常度が次のように表せる。

a(\boldsymbol{x}')=K(\boldsymbol{x}',\boldsymbol{x}')+\sum_{n_1,n_2=1}^N\alpha^{*}_{n_1}\alpha^*_{n_2}K_{n_1,n_2}-2\sum_{n=1}^N\alpha^*_nK(\boldsymbol{x}',\boldsymbol{x}^{(n)})-R^{2*}

データを準備する。
ここでは参考書に倣って、2つの2次元ガウス分布からサンプリングする。

# データの準備
X1 = np.random.multivariate_normal(mean=[0,0], cov=np.array([[1,0],[0,1]]), size=60)
X2 = np.random.multivariate_normal(mean=[3,3], cov=np.array([[1,0],[0,1]]), size=60)
X = np.concatenate([X1, X2])

X = (X - X.mean(axis=0))/(X.std(axis=0))
plt.plot(X[:,0], X[:,1], 'o')

one-class SVMを実装する。
ここでは、最適化問題を解くためにcvxoptというライブラリを使用する。

import cvxopt as cv

def rbf(x1, x2, beta):
    """RBFカーネル
    """
    return np.exp(-beta*np.sum((x1 - x2.reshape(-1,1,2))**2, axis=2))

# パラメータ
C = 0.2
beta = .5
N = len(X)

# 各行列の計算
gram_matrix = rbf(X, X, beta=beta)

# 以下のような最適化問題を解く関数を使用
# 上記の問題と比べて値をきめていく
#    minimize    (1/2)*x'*P*x + q'*x
#    subject to  G*x <= h
#                A*x = b.

# 目的関数
P = cv.matrix(gram_matrix*2)
q = cv.matrix(-np.diag(gram_matrix))

# 不等式の条件
G = cv.matrix(np.r_[np.identity(N), -np.identity(N)])
h = cv.matrix(np.r_[C*np.ones(N).T, np.zeros(N).T])

# 等式の条件
A = cv.matrix(np.array([np.ones(N)]))
b = cv.matrix(1.0)

# 最適化問題を解かせる
sol = cv.solvers.qp(P, q, G=G, h=h, A=A, b=b, kktsolver='ldl2') # ldl, qr, chol

# 解の取得
alpha = np.array(sol['x'])

# バイアスの計算
b = np.sum(alpha*X, axis=0)

# サポートベクトル
index_list = np.where((alpha>=1e-6)&(alpha<C))[0]
support_vector = X[index_list]
print(index_list)

異常度を計算し、可視化する。
赤い点がサポートベクトルを表している。

xx, yy = np.meshgrid(np.arange(-3,3, 0.1), np.arange(-3,3, 0.1))
xy = np.concatenate([xx.ravel().reshape(-1,1), yy.ravel().reshape(-1,1)], axis=1)

# カーネル行列の計算
gram_matrix_d = rbf(X, xy, beta=beta)
gram_matrix_dd = rbf(xy, xy, beta=beta)

# 異常度の計算
R2 = gram_matrix[index_list[0],index_list[0]]+np.sum(alpha@alpha.T * gram_matrix) - 2*np.sum(alpha.T*gram_matrix[index_list[0],:], axis=1)
a = np.diag(gram_matrix_dd)+np.sum(alpha@alpha.T * gram_matrix) - 2*np.sum(alpha.T*gram_matrix_d, axis=1) - R2
a_r = np.array([1 if _a>=0 else -1 for _a in a])

# 可視化
plt.contourf(xx,yy,a.reshape(len(xx), len(yy)), cmap='coolwarm');
plt.contour(xx,yy,a_r.reshape(len(xx), len(yy)), levels=[0,1], colors='red');
plt.scatter(X[:,0], X[:,1]);
plt.plot(support_vector[:,0], support_vector[:,1], 'o', color='red')

one-class SVMはschikit-learnでも実装されている。
ほぼ同様の結果が得られた。

from sklearn.svm import OneClassSVM

osvm = OneClassSVM(nu=1/(N*C), gamma=.5)
osvm.fit(X)
pred_osvm = osvm.predict(X)
support_vector_osvm = osvm.support_vectors_

pred_xy = osvm.decision_function(xy)
pred_osvm_r = np.array([1 if p>=0 else -1 for p in pred_xy])
plt.contourf(xx,yy,pred_xy.reshape(len(xx), len(yy)), cmap='coolwarm');
plt.contour(xx,yy,pred_osvm_r.reshape(len(xx), len(yy)), levels=[0,1], colors='red');
plt.scatter(X[:,0], X[:,1]);
plt.scatter(support_vector_osvm[:,0], support_vector_osvm[:,1], color='red');

以上となります。
世の中のデータサイエンティストと呼ばれる人たちにとっては基本的なことなんでしょうが、やはり難しいです。
なんとなく数式が読めても、自分で導けるようになれるとは思いません。

次回

方向データの異常検知