機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

直線と点の距離

下図に示すような点\mathbf{a}から、直線y(\mathbf{x})=\mathbf{w}^{T}\mathbf{x}+b=0におろした垂線との交点\mathbf{h}との距離を考えます。

f:id:opabinia2:20200114225421p:plain


直線y(\mathbf{x})の法線ベクトルは\nabla y=\mathbf{w}(参考:法線ベクトルと勾配)ですから、


\displaystyle \mathbf{a} - \mathbf{h} = d \frac{\mathbf{w}}{\|  \mathbf{w} \|} \tag{1}

となる実数dが存在します。式(1)において\displaystyle \frac{\mathbf{w}}{\|  \mathbf{w} \|}\mathbf{w}方向への単位ベクトルです。したがって|d|は点\mathbf{a}から\mathbf{h}までの距離を表します。使っている文字が1次関数y=ax+bと似ててややこしいですが、y(\mathbf{x})=w_1 x_1+w_2 x_2+b=0を想定しています。

式(1)より


\displaystyle \mathbf{h} =\mathbf{a}- d \frac{\mathbf{w}}{\|  \mathbf{w} \|} \tag{2}

です。ここで\mathbf{h}は直線y(\mathbf{x})=0上の点ですから、


\begin{eqnarray*}
\displaystyle \mathbf{w}^{T} \left( \mathbf{a} - d \frac{\mathbf{w}}{\|\mathbf{w}\|} \right) + b &=& 0 \tag{2} \\
\displaystyle \mathbf{w}^{T}\mathbf{a} - d \frac{\mathbf{w}^{T}\mathbf{w}}{\|\mathbf{w}\|} + b &=& 0 \tag{3}
\end{eqnarray*}

です。ベクトルのノルムの式(5)より \mathbf{w}^{T}\mathbf{w} = \| \mathbf{w} \|^{2}ですから、


\begin{eqnarray*}
\mathbf{w}^{T}\mathbf{a} - d \| \mathbf{w}\| + b &=& 0 \tag{4} \\
\displaystyle d &=& \frac{\mathbf{w}^{T}\mathbf{a}+b}{ \|\mathbf{w}\| } \tag{5} \\
&=& \frac{y(\mathbf{a})}{\|\mathbf{w}\|} \tag{6}
\end{eqnarray*}

です。したがって、点\mathbf{a}から直線y(\mathbf{x})におろした垂線の距離dは、


\displaystyle |d| = \frac{|y(\mathbf{a})|}{\|\mathbf{w}\|} \tag{6}

です。直線でなくとも、曲線でも平面でも同じ式で距離が求められます。

最大マージン分類器

2クラスの分類問題を


y(\mathbf{x}) = \mathbf{w}^{T}\phi(\mathbf{x})+b \tag{1}

で識別することを考えます。*1 このとき、分離できる解は複数存在します。例えばパーセプトロンの場合、係数の初期値によって下図のように複数の解が求まります。

f:id:opabinia2:20200112121236p:plain

最大マージン分類器では、分離境界から最も近くの訓練データまでの距離をマージンと定義し、このマージンが最大となるような解を求めます。そして、分離境界からマージンだけ離れた点をサポートベクトルと呼びます。下図のようなイメージです。

f:id:opabinia2:20200112122721p:plain


このように、マージンを最大化するという考えのもとに分類する手法をサポートベクターマシンと呼ぶようです。サポートベクトルマシンとも呼ぶようですが、ベクターのほうがかっこいい。次回以降、解の求め方を確認していきますが、識別モデルはカーネル関数を使って表されます。つまりサポートベクターマシンはカーネル法と同様、無限次元を扱うことができるモデルです。(参考:ガウスカーネル) そして、カーネル法では新たな入力に対する結果を求めるのに、全ての訓練データを用いる必要がありましたが、サポートベクターマシンでは、計算に使われるのはサポートベクトルのみです。メモリ上に多くのデータを持つ必要がないなどのメリットがあるのかな。この特性を、疎な解を持つ、あるいはスパースな解、と呼ぶようです。

カーネル法の目次

カーネル法の概要

カーネル回帰分析

ガウス過程


必要な数学知識

関連

ガウス過程による分類(7)実験結果

ガウス過程による分類(6)の続きです。

ガウス過程による分類(1)~(6)までの長い道のりを経て、ようやく以下の式が得られました。


\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )\simeq \sigma \left(   \frac{\mathbf{k}^{T}\mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N})}{1+\displaystyle \frac{\pi}{8}\left\{  c-\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k}  \right\} } \right) \tag{1}

今回はこれを使って予測分布を求めてみたいと思います。おさらいですが、各文字は以下のようなものでした。

\mathbf{C}_{N}は各要素が式(2)の行列。


C(x_{n},x_{m}) = k(x_{n},x_{m}) + \nu \delta_{nm} \tag{2}


\mathbf{k} = ( k(\mathbf{x_{1}},\mathbf{x_{N}}) , \cdots ,  k(\mathbf{x_{N}},\mathbf{x_{N}}))^{T} \tag{3}


c = k(\mathbf{x_{N+1}},\mathbf{x_{N+1}}) + \nu \tag{4}


\boldsymbol {\sigma}_{N} = (\sigma(a_{1}), \cdots, \sigma(a_{N})) \tag{5}


\mathbf{w}_{N} =  \left(
    \begin{array}{cccc}
      \sigma(a_{1})(1-\sigma(a_{1})) &  & & 0 \\
       &  & \ddots &  \\
       0 &  &  & \sigma(a_{N})(1-\sigma(a_{N}))
    \end{array}
  \right) \tag{6}

式(5)、(6)は、ニュートン法により求める \nabla \Psi(\mathbf{a}_{N}^{\prime})=\mathbf{0}となる点 \mathbf{a}_{N}^{\prime} を使います。

では結果です。 f:id:opabinia2:20200104134317p:plain ガウスカーネルの\sigmaを小さくしてみると次のようになりました。 f:id:opabinia2:20200104134324p:plain 結果の図だけ見るとカーネル回帰分析(2)実験結果と似たような感じですね。

訓練データを2点にした結果を見ると、何が行われているのかよくわかります。点を中心としてガウスカーネルが配置されているようなイメージです。中心が少しずれているのは互いの点の影響を受けているためと思います。 f:id:opabinia2:20200104134334p:plain

入力を非線形変換したロジスティック回帰でもガウス基底を使って識別をしていましたが、これは基底を自分で適切に設定しなければなりませんでした。一方カーネル法では無限次元の特徴を持つガウスカーネルを使えば、基底をどう設定するか、ということに頭を使わなくてよくなってしまうんですね。

今回のコードです。ニュートン法のイテレーションは適当に10回固定としています。

# ガウス過程による分類
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


# 教師データ
def create_data(x,y):
    if y > -2*x +1 and y > 2*x -1 :
        return 1
    else:
        return 0


# カーネル関数
def kernel(x0, x1):
    sigma = 0.3
    return np.exp(- np.linalg.norm(x0-x1) ** 2 / (2 * sigma ** 2))


# グラム行列
def gram_matrix(x):
    gram = np.empty([N, N])
    for i in range(N):
        for j in range(N):
            gram[i, j] = kernel(x[:,i], x[:,j])
    return gram


def sigmoid(a):
    return 1/(1+np.exp(-a))


# 予測分布
def p_tt(x2, y2):
    # kを求める
    k = np.zeros(N)
    for i in range(N):
        k[i] = kernel(x[:, i], np.array([x2, y2]))

    mu = k @ (t-sigma_N)
    v = kernel(np.array([x2, y2]), np.array([x2, y2])) - k @ np.linalg.inv(np.linalg.inv(WN)+CN) @ k
    return sigmoid(mu/np.sqrt(1+(np.pi*v)/8))


# ランダムシードを固定
np.random.seed(0)

# 行列Cに加えるノイズ
NU = 0.1

# 訓練データ点数
N = 100

# 2クラス分のデータを作成
x = np.random.uniform(0, 1, [2, N])
t = np.empty([N])
t = np.vectorize(create_data)(x[0, :], x[1, :])

# a' をニュートン法で求める

# a の初期化
a = np.random.uniform(-1,1,N)

# CNを求める
CN = gram_matrix(x) + NU * np.eye(N)

# ニュートン法の繰り返し
for i in range(10):
    # sigma_N
    sigma_N = sigmoid(a)
    # WNを求める
    WN = np.diag(sigmoid(a)*(1 - sigmoid(a)))
    # aの更新
    a = CN @ np.linalg.inv((np.eye(N) + WN @ CN)) @ (t - sigma_N + WN @ a)


#新たな入力
x2, y2 = np.meshgrid(np.linspace(0, 1, 100), np.linspace(0, 1, 100))

# 予測値を計算
result = np.vectorize(p_tt)(x2, y2)


plt.xlim(0, 1)
plt.ylim(0, 1)
plt.contourf(x2, y2, result,30, cmap=cm.coolwarm)
plt.scatter(x[:,np.where(t==1)][0], x[:,np.where(t==1)][1], color="red", alpha=0.5)
plt.scatter(x[:,np.where(t==0)][0], x[:,np.where(t==0)][1], color="blue", alpha=0.5)
plt.show()

ガウス過程による分類(6)

ガウス過程による分類(5)の続きです。

何をやっている途中かと言うと、最終目標である


\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )= \int p(t_{N+1}=1 |a_{N+1} )p(a_{N+1}|\mathbf{t}_{N}) d a_{N+1} \tag{1}

を求めるため、右辺後半の項を計算している途中です。ラプラス近似やニュートン法などを駆使し、ようやくガウス過程による分類(2)の式(10)、ガウス過程による分類(5)の式(16)により、


\begin{eqnarray*}
p(a_{N+1}|\mathbf{t}_{N}) &=& \int p(a_{N+1}|\mathbf{a}_{N}) p(\mathbf{a}_{N}|\mathbf{t}_{N}) d\mathbf{a}_{N} \tag{2} \\
&\simeq& N(a_{N+1}| \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{a}_{N}, c- \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k})N(\mathbf{a}_{N}|\mathbf{a}_{N}^{\prime}, \mathbf{H}^{-1})  \tag{3}
\end{eqnarray*}

と求めることができました。今回はこれを計算してきます。

計算するために、以下の公式を使います。


\begin{eqnarray*}
p(\mathbf x) &=& N(\mathbf x | \boldsymbol\mu, \boldsymbol \Lambda^{-1}) \tag{4} \\
p(\mathbf y | \mathbf x) &=& N(\mathbf y | \mathbf A \mathbf x + \mathbf b,\mathbf{L}^{-1}) \tag{5}
\end{eqnarray*}

であるとき、


p(\mathbf y) = N(\mathbf y | \mathbf A \mathbf \mu + \mathbf b , \mathbf{L}^{-1} + \mathbf A \boldsymbol\Lambda^{-1} \mathbf A^{T}) \tag{6}

ここで、\boldsymbol\mu = \mathbf{a}_{N}^{\prime}\boldsymbol{\Lambda}^{-1} = \mathbf{H}^{-1}\mathbf{A}=\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{L}^{-1}=c-\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}とすれば、


\mathbf{A}\boldsymbol{\mu} + \mathbf{b} = \mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{a}_{N}^{\prime} \tag{7}

です。ここでガウス過程による分類(5)の式(3)より、


\begin{eqnarray*}
\mathbf{t}_{N}-\boldsymbol {\sigma}_{N}-\mathbf{C}_{N}^{-1}\mathbf{a}_{N}^{\prime} &=& 0 \tag{8}\\
\mathbf{a}_{N}^{\prime} &=& \mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N}) \tag{9}
\end{eqnarray*}

です。なお式(9)の右辺の\boldsymbol {\sigma}_{N}の中にaが入っていますからこれで\mathbf{a}_{N}^{\prime}を求めることはできてません。式(7)に代入すれば、


\mathbf{A}\boldsymbol{\mu} + \mathbf{b} = \mathbf{k}^{T}\mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N}) \tag{10}

となり、これがp(a_{N+1}|\mathbf{t})の平均値です。次に分散を求めますと、


\mathbf{L}^{-1} + \mathbf A \boldsymbol\Lambda^{-1} \mathbf A^{T} = c-\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}+\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{H}^{-1}(\mathbf{k}^{T}\mathbf{C}_{N}^{-1})^{T} \tag{11}

です。ここでガウス過程による分類(5)より、\mathbf{H}^{-1}=\mathbf{W}_{N}+\mathbf{C}_{N}^{-1}です。さらに転置行列の定理転置行列の逆行列\mathbf{C}_{N}は共分散行列だから対称行列であることに注意し、式(11)は


\mathbf{L}^{-1} + \mathbf A \boldsymbol\Lambda^{-1} \mathbf A^{T} = c-\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}+\mathbf{k}^{T}\mathbf{C}_{N}^{-1}(\mathbf{W}_{N}+\mathbf{C}_{N}^{-1})^{-1}\mathbf{C}_{N}^{-1}\mathbf{k} \tag{12}

となります。

さて、Woodburyの公式より、


(\mathbf{A}+\mathbf{B}\mathbf{C}\mathbf{D})^{-1} = \mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{B}(\mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B})^{-1}\mathbf{D}\mathbf{A}^{-1} \tag{13}

でした。式(13)において\mathbf{B}=\mathbf{D}=\mathbf{I}とすれば、


(\mathbf{A}+\mathbf{C})^{-1} = \mathbf{A}^{-1} - \mathbf{A}^{-1}(\mathbf{C}^{-1}+\mathbf{A}^{-1})\mathbf{A}^{-1} \tag{14}

となります。式(14)を(\mathbf{C}_{N}^{-1}+\mathbf{W}_{N})^{-1}に対して適用すれば、式(12)は


\begin{eqnarray*}
式(12) &=& c-\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}+\mathbf{k}^{T}\mathbf{C}_{N}^{-1} \left\{   \mathbf{C}_{N}-\mathbf{C}_{N}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{C}_{N} \right\}\mathbf{C}_{N}^{-1}\mathbf{k} \tag{15} \\
&=& c-\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}+\mathbf{k}^{T}\mathbf{C}_{N}^{-1}\mathbf{k}  -\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k} \tag{16} \\
&=& c-\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k} \tag{17} \\
\end{eqnarray*}

となります。

平均と分散が式(10)、式(17)の通り求まり、


p(a_{N+1}|\mathbf{t}_{N}) \simeq N(\mathbf{a}_{N+1} | \mathbf{k}^{T}\mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N}) ,c-\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k}  ) \tag{18} \\

です。

p(t_{N+1}=1|a_{N+1})=\sigma(a_{N+1})でしたから、これと式(18)を求めたかった式(1)に代入すれば、


\displaystyle p(t_{N+1} 1 | \mathbf{t}_{N} )\simeq \int \sigma(a_{N+1})N(\mathbf{a}_{N+1} | \mathbf{k}^{T}\mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N}) ,c-\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k}  ) d a_{N+1} \tag{19}

となります。

ここで以下の近似を使います。


\displaystyle \int \sigma(a)N(a|\mu,\sigma^{2})d\sigma \simeq \sigma\left(\frac{\mu}{\sqrt {1+\displaystyle\frac{\pi\sigma^{2}}{8} }} \right)   \tag{20}

この近似の導出は追っていないのですが、とりあえず既知のものとして使います。

すると式(19)は、


\displaystyle p(t_{N+1}=1 | \mathbf{t}_{N} )\simeq \sigma \left(   \frac{\mathbf{k}^{T}\mathbf{C}_{N}(\mathbf{t}_{N}-\boldsymbol {\sigma}_{N})}{1+\displaystyle \frac{\pi}{8}\left\{  c-\mathbf{k}^{T}(\mathbf{W}_{N}^{-1}+\mathbf{C}_{N})\mathbf{k}  \right\} } \right) \tag{21}

となり、非常に長い道のりでしたが、これで新たな入力に対する予測値を計算することができます。式(21)を用いた実験結果は次回:ガウス過程による分類(7)実験結果

分散共分散行列の半正定値性

共分散行列\mathbf{\Sigma}は半正定値行列であることを確認します。

共分散行列\mathbf{\Sigma}ij成分を\Sigma_{ij}を考えると、共分散行列の定義(参考:多変量正規分布)より、


\begin{eqnarray*}
\Sigma_{ij} &=&\mathrm{Cov}(X,Y) \tag{1}\\
 &=& E[(X_{i}-E[X_{i}])(X_{j}-E[X_{j}])] \tag{2} \\ 
\displaystyle  &=& \frac{1}{N}\sum_{k=1}^{N} (x_{i}^{(k)}-E[X_{i}])(x_{j}^{(k)}-E[X_{j}]) \tag{3} \\ 

\end{eqnarray*}

です。 ここで、y_{i}^{(k)} = x_{i}^{(k)}-E[X_{i}]とし、\mathbf{y}_{i} = (y_{i}^{(1)},\cdots,y_{i}^{(N)})^{T} とすれば、


\displaystyle \Sigma_{ij} = \frac{1}{N} \mathbf{y}_{i}^{T}\mathbf{y}_{j} \tag{4}

と書けます。

したがって共分散行列は、式(4)より、


\begin{eqnarray*}
\displaystyle \mathbf{\Sigma} &=& \frac{1}{N} \left(
    \begin{array}{ccc}
      \mathbf{y}_{1}^{T}\mathbf{y}_{1} & \cdots & \mathbf{y}_{1}^{T}\mathbf{y}_{N}  \\
        \vdots  &   \ddots &  \vdots   \\
       \mathbf{y}_{N}^{T}\mathbf{y}_{1}    & \cdots & \mathbf{y}_{N}^{T}\mathbf{y}_{N} 
    \end{array}
  \right) \tag{5} \\
&=&  \left(
    \begin{array}{c}
      \mathbf{y}_{1}^{T} \\
       \vdots  \\
      \mathbf{y}_{n}^{T}
    \end{array}
  \right)
(\mathbf{y}_{1},\cdots,\mathbf{y}_{n})  
 \tag{6}
\end{eqnarray*}

と書けます。ここで\mathbf{Y}=(\mathbf{y}_{1},\cdots,\mathbf{y}_{n})  とすれば、


\displaystyle \mathbf{\Sigma} =\frac{1}{N} \mathbf{Y}^{T}\mathbf{Y} \tag{7}

です。

ここで任意のベクトル\mathbf{z}に対する二次形式を考え、転置行列の定理を使えば、


\displaystyle \frac{1}{N} \mathbf{z}^{T} \mathbf{Y}^{T}\mathbf{Y}\mathbf{z} = \frac{1}{N} (\mathbf{Y}\mathbf{z})^{T}\mathbf{Y}\mathbf{z} \ge 0 \tag{8}

です。二次形式が\ge 0ですので、定義より半正定値行列であることがわかります。(参考:正定値行列