2020-02-02

20年1月の振り返り

不労所得

正月休みのボーナス期間でWeb収入は好調。が、月末にかけて急激に尻すぼみ状態で、2月はこのブログに記録を残して以来の最低額を更新するかもしれない。

投信のほうはいろいろ事情があって、いったんNISAを除いて全て処分することに。1月上旬に処分したあと、相場は急落。全く予見していたわけではないのですが、結果的にはベストなタイミングで処分できた。あくまでも短期的にはですけど。まあいずれ再開するつもりです。でもいったん処分しちゃうと、再開するなら大きな下落相場を迎えてからがいいなあと思ってしまいます。

20年1月の実績

Web収入	投信前月比評価損益	計
880,615円	117,835円	998,450円

2020-01-22

サポートベクターマシン（ハードマージン）（１）

サポートベクターマシン

サポートベクターマシンにはハードマージン/ソフトマージンという分類があります。訓練データが線形分離可能であるという前提をおいたものをハードマージン、そうでないものをソフトマージンと呼びます。まずはハードマージンを見ていきます。下図はサポートベクターマシンで分類問題を解いた結果例です。緑色で強調した点がサポートベクトルと呼ばれる点です。分類境界を作っているのはこのサポートベクトルのみであり、その他の点は境界に一切影響しないのが特徴です。

f:id:opabinia2:20200122224207p:plain

分類境界が $y(\mathbf{x})=\mathbf{w}^{T}\phi(\mathbf{x})+b$ で与えられるとすれば、訓練データ点 $\mathbf{x}_n$ と境界までの距離 $d_n$ は、直線と点の距離より、

$\displaystyle | d_n | = \frac{|y(\mathbf{x}_n)|}{\|\mathbf{w}\|} = \frac{|\mathbf{w}^{T}\phi(\mathbf{x}_n)+b|}{\|\mathbf{w}\|}\tag{1}$

です。

問題の前提をハードマージン、つまり訓練データが線形分離可能としていますから、データのラベルを $t_n \in \{-1,1\}$ とすれば常に $t_n y(\mathbf{x}_n) \gt 0$ となります。また、 $t_n$ をかけても大きさは変わりませんから、

$\displaystyle d_n = \frac{t_n(\mathbf{w}^{T}\phi(\mathbf{x}_n)+b)}{\|\mathbf{w}\|} \tag{2}$

のように絶対値を外すことができます。

最大マージン分類器で書いた通り、サポートベクターマシンは、分類境界と最も近い点との距離（マージン）を最大にするような分類境界を求めることが目標です。マージンが最大となる $\mathbf{w},b$ は、

$\newcommand{\argmax}{\mathop{\rm arg~max}\limits} \newcommand{\min}{\mathop{\rm min}\limits} \displaystyle \argmax_{\mathbf{w},b} \left( \min_{n} \frac{t_n(\mathbf{w}^{T}\phi(\mathbf{x}_n)+b)}{\|\mathbf{w}\|} \right) \tag{3}$

です。ややこしいですが、内側の $\min$ は、分類境界と最も近い $n$ 番目のベクトルを指し、このようなベクトルが境界から最も遠くなる $\mathbf{w},b$ を選ぶというのが式(3)の意味です。まあ要するにマージンを最大化するということです。

ここで、 $\mathbf{w} \rightarrow \kappa\mathbf{w}$ 、 $b \rightarrow \kappa b$ とすると、式(2)の右辺は

$\begin{eqnarray*} \displaystyle \frac{t_n(\kappa \mathbf{w}^{T}\phi(\mathbf{x}_n)+\kappa b)}{\kappa \|\mathbf{w}\|} &=&\frac{t_n\kappa (\mathbf{w}^{T}\phi(\mathbf{x}_n)+ b)}{\kappa \|\mathbf{w}\|} \tag{4}\\ &=& \frac{t_n(\mathbf{w}^{T}\phi(\mathbf{x}_n)+b)}{\|\mathbf{w}\|} \tag{5}\\ &=&d_n \tag{6} \end{eqnarray*}$

となり、結局距離 $d_n$ は変わりません。 $\kappa$ によって任意に $\mathbf{w},b$ をスケーリングできるわけですから、境界に最も近い点に対して、

$t_n(\mathbf{w}^{T}\phi(\mathbf{x}_n)+b) = 1 \tag{7}$

となるように $\mathbf{w},b$ を選ぶことができます。これを式(3)に代入すれば、

$\newcommand{\argmax}{\mathop{\rm arg~max}\limits} \displaystyle \argmax_{\mathbf{w},b} \frac{1}{\|\mathbf{w}\|} \tag{8}$

となります。式(8)を最大にする $\mathbf{w},b$ は、 $\|\mathbf{w}\|^{2}$ を最小にする $\mathbf{w},b$ と同じです。したがって求める $\mathbf{w},b$ は

$\newcommand{\argmin}{\mathop{\rm arg~min}\limits} \displaystyle \argmin_{\mathbf{w},b} \frac{1}{2}\|\mathbf{w}\|^{2} \tag{9}$

と書けます。 $\frac{1}{2}$ をかけているのは、微分したときにきれいな形にするためのもので、あってもなくても求められる解は同じです。また、わざわざ2乗しているのは、後で書くように2次計画問題にもっていきたいからです*1。なお $b$ は式(7)の制約がありますから、消えているわけではありません。

式(7)は、分類境界から最も近い点までの距離が1であるということですから、その他の点までの距離は1より大きいといえます。

以上より、

$\displaystyle \frac{1}{2}\|\mathbf{w}\|^{2} \tag{10}$

を、

$t_n(\mathbf{w}^{T}\phi(\mathbf{x}_n)+b) \ge 1 \quad\quad (n=1,\cdots,N) \tag{11}$

という不等式制約条件下での2次関数の最小問題を解けば良いことになります。制約条件が線形で、目的関数が凸の2次形式であるこのような問題を2次計画と呼ぶようです。

続き：サポートベクターマシン（ハードマージン）（２）

*1:詳しく調べていないのですが、そうしないと解を求めるのが大変、、、なのだと思います。

2020-01-14

直線と点の距離

線形代数・ベクトル

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

直線 $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}$

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

2020-01-12

最大マージン分類器

サポートベクターマシン

2クラスの分類問題を

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

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

f:id:opabinia2:20200112121236p:plain

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

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

*1:参考：線形識別の概要

2020-01-04

カーネル法の目次

カーネル法

カーネル法の概要

カーネル回帰分析

ガウス過程

必要な数学知識

ガウス過程による分類（７）実験結果

カーネル法

ガウス過程による分類（６）の続きです。

ガウス過程による分類（１）～（６）までの長い道のりを経て、ようやく以下の式が得られました。

$\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 結果の図だけ見るとカーネル回帰分析（２）実験結果と似たような感じですね。

訓練データを２点にした結果を見ると、何が行われているのかよくわかります。点を中心としてガウスカーネルが配置されているようなイメージです。中心が少しずれているのは互いの点の影響を受けているためと思います。 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()

2020-01-03

ガウス過程による分類（６）

カーネル法

ガウス過程による分類（５）の続きです。

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

$\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}$

を求めるため、右辺後半の項を計算している途中です。ラプラス近似やニュートン法などを駆使し、ようやくガウス過程による分類（２）の式(10)、ガウス過程による分類（５）の式(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}$

です。ここでガウス過程による分類（５）の式(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}$

です。ここでガウス過程による分類（５）より、 $\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)は、

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

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

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

20年1月の振り返り

20年1月の実績

サポートベクターマシン（ハードマージン）（１）

直線と点の距離

最大マージン分類器

カーネル法の目次

カーネル法の概要

カーネル回帰分析

ガウス過程

ガウス過程による回帰

ガウス過程による分類

必要な数学知識

関連

ガウス過程による分類（７）実験結果

ガウス過程による分類（６）