サポートベクトルマシン（SVM）を簡単に、わかりやすく説明したい

機械学習統計 matlab

この記事は、MATLAB アドベントカレンダー（その2）の12日目の記事として書かれています。

1 はじめに

本記事では、サポートベクトルマシン（Support Vector Machine: SVM）の仕組みについてまとめます。

サポートベクトルマシンは複数の次元を有するデータを超平面で分離する手法として有名です。色々な場面で利用でき、機械学習の代表的な手法の1つ言えます。以下の図は、3つの変数をもつデータを2つのクラスに分類しているときの図です。各軸が、それぞれの変数に対応しています。また、赤色の局面（超平面）によって、2つのクラスに分類されます。この図は3変数の場合を示していますが、これがより多くの変数の数（次元）になっても問題ありません。このように2つのクラスを超平面（曲面）でうまく分類できると非常に便利そうで、かつおもしろいです。本記事では、このサポートベクトルマシンがどのように計算されるのかを紹介したいと思います。

このサポートベクトルマシンの分離平面を可視化するためのコード (Python、MATLAB）は以下のブログやfile exchangeにて公開しています。もしよろしければ、こちらもご覧ください。

過去の記事：サポートベクターマシン（SVM）の分離平面の可視化 (Python, MATLAB)

kentapt.hatenablog.com

File Exchange: Visualizing a hyper-plane in SVM classifier（SVMの分離境界面の可視化）

jp.mathworks.com

また、本記事を執筆するにあたり、私の理解に基づき、2クラス分類のサポートベクトルマシンを1から実装し（ただしかなりおおばっぱな実装です）、以下のようなそれらしい結果を得ることができました。赤丸と青丸がそれぞれ異なるクラスに属するサンプルで、赤線と青線がそれぞれを分類（分離）する境界です。赤色および青色のサンプルはそれぞれ、同じ色の境界で囲まれていて、うまく分類ができていそうです。

下の図は、本記事で説明する、正則化項の値を値を変えながら学習した時の、境界です。

また、下の図は、ガウシアンカーネルのパラメータを変化させながら学習させたときの境界になります。

本記事の執筆にあたっては、

赤穂先生の「カーネル多変量解析」を非常に参考にさせていただきました。

honto.jp

非常にわかりやすい書籍であり、大変おすすめです。

注意1：説明する流れについて

サポートベクトルマシンについて説明する方法として、

1) 「マージンの最大化」

2) 「損失関数＋正則化項」

といった複数の説明の方法があります。ネットの記事では1で説明するものも多いですが、本記事では、2）のアプローチで説明することを試みます。

注意2: サポートベクトルマシンのという用語について

本記事では、サポートベクトルマシンを2つのクラスに分類するアルゴリズムとして説明しています。多クラスに分類する際の方法も7章にて述べます。

2 サポートベクトルマシンの導入について

まずは、簡単にするために、データの分布と、サポートベクトルマシンの境界を2次元で見てみます。

以下の図を見てください。赤色と青色がそれぞれ異なるクラスに属するサンプルで、x, yがそれぞれのサンプルの持つ変数の値を示します。Boundaryという境界も直線で示されています。この境界の右上が「クラス1」、左下が「クラス-1」を示します。

このように、直線や曲線で境界が定義できれば、簡単にクラスの分類ができそうです。ただ、その境界はどのようにして作ればいいのでしょうか。

以下に示す、私の投稿した過去の記事では、線形回帰（単回帰・重回帰分析）の式の重みを解析的に求めていました。

線形回帰（単回帰分析）を1から実装して理解を深めてみよう

kentapt.hatenablog.com

線形回帰（単回帰分析）を1から実装して理解を深めてみよう: つづき

kentapt.hatenablog.com

重回帰分析の勉強＆実装による検算をしてみよう

kentapt.hatenablog.com

そこでは、以下のような式を定義して、最小二乗法によって、その係数（最小二乗推定量）を求めました。

$y~=~a~+~b_1 x_{1~} +~b_2 x_2 ~+~\cdots +~b_k x_k$ 　...(1)

線形回帰の場合の損失関数は以下のようなものでした。この値が小さいほど、予測の誤差が小さいことを示しています。

$E={\left(y-\omega X\right)}^2$ ...(2)

前回の記事では、正解ラベルyと予測値の値の差分を二乗し、その和を最小にする係数を解析的に計算しました。

今回の、サポートベクトルマシンでは、各サンプルをカテゴリーに分ける（ラベルは1 や -1）ことを目指します。

線形回帰の場合は、誤差は4である、とか、10である、といった値が返され、計算が簡単にできました。

しかし、カテゴリー分けの場合は、その分類が間違ってる or 正解の2択になるため、その2乗誤差による評価は直接的に利用することができません。もちろん、ニューラルネットワークのように2乗誤差による評価を利用可能にすることも考えられますが、ひとまず今は以下に述べる方法によって、サポートベクトルマシンの分類の性能を評価します。

sign (sng)を符号関数を利用します。

符号関数については、以下のページがわかりやすかったです。

mathlandscape.com

...(3)

符号関数（以下、sgnと表す）を用いて、サポートベクトルマシンの演算結果が0以上の場合は、1を、0より小さい場合は-1になるようにすることができます。なお、ここでは、演算結果が0の場合も、1を返すようにしています。なお、 $f(x)$ はサポートベクトルマシンにサンプルを入力した時に利用する関数で、ざっくりいうと、サポートベクトルマシンの中身、といったふうにも言えると思います。

すると、損失関数は以下のように定義することができます。

$E={\left(\frac{y-sgn\left\lbrack f\left(x\right)\right\rbrack }{2}\right)}^2$

...(4)

$sgn\left\lbrack f\left(x\right)\right\rbrack$ の値は、先ほど述べたように、1か-1でした。yは正解ラベルで同様に1か-1を示します。クラスの名前を与えるとすると、クラスAやクラスBというふうに言い換えることもできます。これにより、サポートベクトルマシンの演算結果が正解の場合、損失は0、そうでない（間違えている）場合は1になることがわかります。具体的な例は後で挙げることにして、ここでは式の整理を進めます。

$=~\frac{1-y\times sgn\left\lbrack f\left(x\right)\right\rbrack }{2}$

...(5)

$=~\frac{1-sgn\left\lbrack y\times f\left(x\right)\right\rbrack }{2}$

...(6)

(6)式を得たところで、yのパターンと、そのサンプルの推論結果が正しい/間違い場合の、損失関数の値を確かめてみます。yが1/-1の場合と、それぞれで正解/不正解の場合があるため、4通りの場合が存在します。それぞれについて確認してみます。

$y$ が1で $sgn\left\lbrack f\left(x\right)\right\rbrack$ が1 (正解）　⇒　 $E$ = 0

$y$ が1で $sgn\left\lbrack f\left(x\right)\right\rbrack$ が-1 (不正解）　⇒　 $E$ = 1

$y$ が-1で $sgn\left\lbrack f\left(x\right)\right\rbrack$ が-1 (正解）　⇒　 $E$ = 0

$y$ が-1で $sgn\left\lbrack f\left(x\right)\right\rbrack$ が1 (不正解）　⇒　 $E$ = 1

...(7)

以上のように、それぞれの場合において、正解の時は損失が0、不正解の時は損失が1になっていることがわかります。

このルールで、適切に誤差を評価できそうです。この誤差が0（ゼロ）であれば、そのサンプル（訓練データ）中のミスもゼロとなります。そのため、この誤差をゼロにするような関数 $f(x)$ を求めればよいことがわかります。

しかし、(6)式では、凸関数になっていないという問題があります。以下の図を見てください。

凸関数では、「コブ」が1つだけ存在します。一方、凸関数でない関数は、「コブ」が複数存在します。

凸関数の場合は、最適解を求めやすいです。直感的に言うと、右と左を見て、その点が下方向にあれば、その点が最小値を取っていることがわかります。一方、凸関数でない場合は、関数全体の値を見渡し、その中での最小値を探す必要があります。

画像出展：凸関数と凸でない関数

www.i.u-tokyo.ac.jp

(6)式と似た形で、かつ凸関数の損失関数を用意するために、ヒンジ関数を導入します。

すると、(6)式は以下のように書き換えることができます。

$\gamma_{\textrm{hinge}} \left(y\;f\left(x\right)\right)=\max \left\lbrace 0,1-y\;f\left(x\right)\right\rbrace$

...(8)

(8)式をより具体的にイメージするために、以下の図を見るとわかりやすいと思います。

(8)は b = 0 の場合と考えてください。0と $1-y\;f\left(x\right)$ の2つの値を比べます。0というのは、サポートベクトルマシンの予測結果が正しい場合です。一方、 $1-y\;f\left(x\right)$ の値が採用される場合は、その推論が正しくない場合です。例えば、y = 1 のサンプルに対して、負の値（正しくない）である -10 を返した場合を考えます。この場合、 $1-y\;f\left(x\right)$ の値は11となり、 $\max \left\lbrace 0,1-y\;f\left(x\right)\right\rbrace$ では、11という値が採用されることになります。

画像出展：ヒンジ関数の意味、損失関数として使えることの大雑把な説明 https://mathwords.net/hinge

この $\gamma_{\textrm{hinge}} \left(y\;f\left(x\right)\right)$ を利用すると、(7)は以下のように書き換えることができます。

$y$ が1で $f\left(x\right)$ が0以上 (正解）⇒ $\gamma_{\textrm{hinge}} \left(f\left(x\right)\right)=0$

$y$ が1で $f\left(x\right)$ が0 より小さい(不正解）⇒ $\gamma_{\textrm{hinge}} \left(f\left(x\right)\right)=1-y\;f\left(x\right)$ 　

$y$ が-1で $f\left(x\right)$ が0 より小さい (正解）⇒ $\gamma_{\textrm{hinge}} \left(f\left(x\right)\right)=0$

$y$ が-1で $f\left(x\right)$ が0以上(不正解）⇒ $\gamma_{\textrm{hinge}} \left(f\left(x\right)\right)=1-y\;f\left(x\right)$

...(9)

このように、予測が正解の場合は $\gamma_{\textrm{hinge}} \left(f\left(x\right)\right)=0$ 、そして、不正解のときは正の誤差値を返すことがわかります。

3. カーネル関数の利用

前節では、 $\gamma_{\textrm{hinge}} \left(y\;f\left(x\right)\right)=\max \left\lbrace 0,1-y\;f\left(x\right)\right\rbrace$ というヒンジ関数を用いた損失関数を定義しました。

多次元の入力 $x$ に対して、ラベルが1の場合は、0以上、ラベルが-1の場合は、0より小さい値を返すとよい、という設定でした。

非常にシンプルに考えると、以下のような線形回帰（単回帰）のような形でもよさそうです。

$f\left(x\right)=\textrm{ax}$

...(10)

ここで、以前のブログ記事で、カーネル法による回帰を紹介しました。

kentapt.hatenablog.com

カーネル法を利用することで、以下の図のようにグネグネとした曲線を描くことができ、

正則化などを利用しない場合は、サンプル（訓練データ）に完全にフィットする曲線を得ることができるのでした。

サポートベクトルマシンの場合について考えます。以下の図の左の図が直線で青と赤のクラスを分類しています。

右の図のような場合はどうでしょうか。左の例のように直線ではきれいに分けることができず、図のように曲がった形の分離境界が必要であることがわかります。サポートベクトルマシンの場合でも、このような2次元や多次元のサンプルに対しては、このようなより複雑な分離境界を利用するほうがより高精度な分類が可能になると考えられます。

カーネル法の便利な点として、何らかのモデル式（例： $\left.y=ax^3 \right)$ を考える必要はなく、

適切なカーネル関数やグラム行列Kを利用することで、複雑な曲線（曲面）を描くことができるのでした。

グラム行列の内容も以前の記事にて記述しています。

https://kentapt.hatenablog.com/entry/entry/matlabAdventDay11

そこで、 $x$ を $K\left(x^{\left(i\right)} ,x\right)$ に読み替え、 $f\left(x\right)$ を以下のように定義します。

$f\left(x\right)=\sum_{i=1}^n \alpha_i \;K\left(x^{\left(i\right)} ,x\right)$

...(11)

また、損失関数は、正則化項を加え、以下のようになります。

$min\sum_{i=1}^n \left(r_{hinge} \left(y^i f\left(x^i \right)\right)+\frac{\lambda }{2}w^T Kw\right)$

...(12)

この正則化項に関しても、以前のブログ記事に詳細の記述があります。

まとめると、この(12)式を最小化する、重み $\alpha_i$ が求めるとよいことがわかります。

そして、予測の際には、

$\sum_{i=1}^n \left(w_i X_i \right)$ を計算する ⇒ 正（負）ならラベル1（-1)と分類します。

このように、損失関数や、最適化したい値を決めることができました。次の章では具体的な計算のステップについて考えていきます。

4. サポートベクトルマシンの学習のための具体的な計算

4.1. スラック変数の導入

さきほどの(12)式を解く方法をより簡単にすることを考えます。(12)式中の $y\;f\left(x\right)$ の値の値に応じて場合分けが生じるため、

最適化の計算が面倒になります。そこで、より式を解きやすくするために、凸二次計画問題 (QP=(convex) quadratic programming) に言い換えることを考えます。線形関数で表された制約条件をみたし，2次の目的関数を最小化する解を求める問題を2次計画問題といい，さらに目的関数が凸関数である場合に凸2次計画問題といいます。

http://www.me.titech.ac.jp/~mizu_lab/lib/pdf/kougisiryou/suuti/handout/11/suuti_kougi11-6.pdf

凸二次計画問題の場合は、既存のライブラリが色々と存在し、比較的簡単に解けるようです。

www.letsopt.com

それでは、凸二次計画問題に言い換えてみます。

$\gamma_{\textrm{hinge}} \left(y\;f\left(x\right)\right)=\max \left\lbrace 0,1-y\;f\left(x\right)\right\rbrace$

の値を $\xi_i$ （スラック変数）とします。すると、求めるものは、

$\xi_i \ge 0$ ...(13)

および

$\xi_i \ge {1-y}^i f\left(x^i \right)=1-y^{(i)} \sum_{j=1}^n \alpha_j K_{ij}$ ...(14)

の両方を満たす $\xi_i$ の中で、最も小さいものである。

「を満たすもの中での最小値」というのは変に聞こえるかもしれません。

例えば、 $\xi_i$ の値が2以上であれば、上の条件を満たす場合、この「を満たすもの中での最小値」という条件がないと、

$\xi_i$ は2でも3でも4でも、何であっても満たしてしまうことになります。

このような言い換えをすることで、式(12)の最適化の問題を、凸2次計画問題に持ち込むことができました。

$\xi_i$ や以下のようなグラム行列 $K$ を利用して、

...(15)

(12)式を書き換えると、最適化する2次関数の式は、以下のようになります。

$min\sum_{i=1}^n \left(\xi_i +\frac{\lambda }{2}\alpha^T K\alpha \right)$

...(16)

4.2. ラグランジュの未定乗数法を用いた凸二次計画問題の求解

本記事では、先ほど求めた、凸二次計画問題を双対問題に帰着させて解いていきます。

式(14)のような制限のもと、式(16)の最小値を求めるには、ラグランジュの未定乗数を利用します。

ラグランジュの未定乗数法は、以下の主成分分析の記事で扱いました。本記事では省略しますが、もし興味がございましたら以下の記事もご覧ください。

kentapt.hatenablog.com

ラグランジュ関数は、以下のようになります。

$\left.~L(\xi ,\alpha ,\beta ,\gamma \right)=~\sum_{i=1}^n \xi_i +\frac{\lambda }{2}\alpha^T K\alpha -\sum_{i=1}^n \beta_i \xi_i -\sum_{i=1}^n \gamma_i \left(\xi_i -1+y^i \sum_{j=1}^n \alpha_j K_{ij} \right)$

...(17)

ここで、ラグランジュ乗数 $\beta_i$ および $\gamma_i$ を導入しました。

$\alpha$ について偏微分をします。この値が0となるので、以下の式を得ます。

$\frac{\partial L}{\partial \alpha }=\sum_{j=1}^n \left(\lambda K_{ij} \alpha_j -\gamma_j y^{\left(i\right)} K_{ji} \right)=0$

...(18)

また、グラム行列は対象行列であるため、

$K_{ij} =K_{ji}$

...(19)

です。そのため、式(18)は

$\lambda \textrm{Ka}-K\hat{r} =0$

...(20)

ただし、 $\hat{\gamma} ={\left({\hat{\gamma} }_i ,{\hat{\gamma} }_2 ,\ldotp \ldotp \ldotp ,{\hat{\gamma} }_n \right)}^{\mathrm{T}}$

...(21)

${\hat{r} }_i =\gamma_i y^{\left(i\right)}$

...(22)

としました。

Kに逆行列が存在する場合、

$\alpha_i =\frac{1}{\lambda }\gamma_i y^{\left(i\right)}$

...(23)

また、 $\xi_i$ について偏微分を行うと、

$1-\beta_i -\gamma_i =0$

...(24)

が得られます。

よって、式(23)と(24)を利用すると、式(17)は以下のように変形できます。

$\left.~L(\xi ,w,\beta ,\gamma \right)=~\sum_{i=1}^n \xi_i +\frac{\lambda }{2}\alpha^T K\alpha -\sum_{i=1}^n \beta_i \xi_i -\sum_{i=1}^n \gamma_i \left(\xi_i -1+y^i \sum_{j=1}^n \alpha_j K_{ij} \right)$

...(17を再掲)

$\gamma_i$ を含む項をひとまとめにし、さらに $\alpha_i =\frac{1}{\lambda }\gamma_i y^{\left(i\right)}$ を代入して、

$=\sum_{i=1}^n \left(1-\beta_i -\gamma_i \right)\xi_i +\frac{\lambda }{2}\frac{1}{\lambda^2 }\sum_{i=1}^n \sum_{j=1}^n y^i y^j \gamma_i \gamma_j K_{ij} +\sum_{i=1}^n \gamma_i -\frac{1}{\lambda }\sum_{i=1}^n \left(y^i \gamma_i \sum_{j=1}^n \gamma_j y^j K_{ij} \right)$

...(25)

最後の項を整理して、

$=\sum_{i=1}^n \left(1-\beta_i -\gamma_i \right)\xi_i +\frac{1}{2\lambda }\sum_{i=1}^n \sum_{j=1}^n y^i y^j \gamma_i \gamma_j K_{ij} +\sum_{i=1}^n \gamma_i -\frac{1}{\lambda }\sum_{i=1}^n \sum_{j=1}^n y^i y^j \gamma_i \gamma_j K_{ij}$

...(26)

$1-\beta_i -\gamma_i =0$ であり、さらに2項目と4項目をまとめることができるため、

$=\sum_{i=1}^n \gamma_i -\frac{1}{2\lambda }\sum_{i=1}^n \sum_{j=1}^n y^{(i)} y^{(j)} \gamma_i \gamma_j K_{ij}$

...(27)

となります。

また、 $\beta_i >0$ および、 $\gamma_i >0$ と $1-\beta_i -\gamma_i =0$ という条件を合わせて、

${0<\gamma }_i <1$

...(28)

を得ることができます。そのため、 ${0<\gamma }_i <1$ の条件下で、

$\sum_{i=1}^n \gamma_i -\frac{1}{2\lambda }\sum_{i=1}^n \sum_{j=1}^n y^{(i)} y^{(j)} \gamma_i \gamma_j K_{ij}$

を最大化する $\gamma_i$ を求めればよいことがわかります。グラム行列Kは与えられたサンプル（訓練データ）で計算できます。また、 $y$ も正解ラベルのため、既知です。そのため、 $\gamma_i$ のみが未知です。

この凸二次計画問題は、SMO(sequential minimal optimization)と呼ばれるアルゴリズムなどを用いて

解くことができます。本記事では、SMOについては触れず、これらの手順で、 $\gamma_i$ を求めたらよいということで本章は終わりにしたいと思います。

4.3. サポートベクトルマシンの計算手順の簡単なまとめ

比較的シンプルな式(27)を得ることができました。これにより、サポートベクトルマシンの計算において、どのような手順を踏めばよいかわかりました。

1) （訓練データとして）与えられたサンプル $x^{\left(1\right)}$ から $x^{\left(n\right)}$ までの値を用いて、カーネル法で用いるグラム行列 $K_{\textrm{ij}}$ 　を計算する

2) 式(27)を ${0<\gamma }_i <1$ という条件下で最適化するパラメータ $\gamma_i$ を求める

3) 式(23)の $\alpha_i =\frac{1}{\lambda }\gamma_i y^{\left(i\right)}$ に $\gamma_i$ の値を代入し、 $\alpha_i$ の値を求める

4) 式(11)の $f\left(x\right)=\sum_{i=1}^n \alpha_i \;K\left(x^{\left(i\right)} ,x\right)$ という値を求める（ここでは、テストしたいデータを推論している）

5) この値が0以上であれば、クラスA、そうでなければクラスBというふうに、分類を行うことができる

5. サポートベクトルについて

4章にて、具体的なサポートベクトルマシンの学習に必要な計算について述べました。それでは、サポートベクトルマシンの名前の中の、「サポートベクトル」とは何のことでしょうか。さらに、式(27)を見ると、グラム行列Kを計算していますが、それはサンプル数×サンプル数の形をしているため、訓練データが増えれば増えるほど、非常に大きなデータになってしまいます。それでは、運用上大変になりそうです。

(17)式の $\gamma_i$ に関して偏微分を行います。この値は0になるので、

$\gamma_i \not= 0$ のとき、

$\xi_i =1-y^{\left(i\right)} f{\left(x\right)}^{(i)}$

...(29)

が成り立ちます。

厳密には、カルーシューキューンータッカー (Karush-Kuhn-Tucker) 定理について触れる必要がありそうですが、ひとまずここでは、ざっくりとした議論を進めます。

式(23)の $\alpha_i =\frac{1}{\lambda }\gamma_i y^{\left(i\right)}$ から、

$\gamma_i =0$ のときは $\alpha_i =0$ となることがわかります。

つまり、 $\alpha_i \not= 0$ となるのは、 $\xi_i =1-y^{\left(i\right)} f{\left(x\right)}^{(i)}$ の場合に限られます。

よって、サポートベクトルマシンの計算においては、以下のヒンジ関数の図の、

折れ線の折れている点か、その左側の点のみを利用すればよいことになります。

また、右側の点は、サポートベクトルマシンの計算では登場しません。

そのため、サポートベクトルマシンの計算では、全体のサンプルが使われるわけではなく、

一部の点のみを保持すればよいことになります（スパース性）。一般には、サポートベクターは、もとの訓練サンプル数に比べてかなり少ないと言われています（栗田先生サポートベクターマシン入門）。

6. マージン最大化について

冒頭でも述べた通り、サポートベクトルマシンを説明する場合、「マージン最大化」が利用される場合も多いです。

赤と青のサンプルを「スパっと」わける超平面を想定し、その超平面と赤と青のサンプル（のうち最も近いもの）の距離を

最大化しようという考え方です。

図出展：最初に学ぶ 2分類SVMモデルの基本 www.acceluniverse.com

そして、その文脈において、以下の図のように、直線や曲面で分けられない時は、

カーネル法を利用しながら、別の特徴空間にサンプルの写像をして、その特徴空間上で、線形分離をする。そうすると、以下のようなサンプルに対してもうまく分類を行うことができる、といったような説明がよくあると思います（図の出展のブログでは別の例を使ってわかりやすく説明されていました）。個人的には、サポートベクトルマシンをはじめて知ったときは、この説明はわかりやすいものの、実際の計算が想像しにくいという印象でした。しかし、本記事のような、「損失関数と正則化項の最適化」というアプローチを知ったうえでもう一度見てみると、より具体的に想像ができる気がしました。

図出展：最初に学ぶ 2分類SVMモデルの基本

www.acceluniverse.com

そこで、今回の説明では採用しなかった、マージン最大化という観点で考えてみたいと思います。

簡単に、マージン最大化をもとに、サポートベクトルマシンを説明してみます。

以下は、冒頭で掲載した図を少し変更したものです。四角の点がそれぞれサポートベクターを示しています。

超平面の式を $f\left(x\right)=\omega^{\mathrm{T}} \phi \left(x\right)$ とし、各サンプルの値（特徴ベクトル）を $\phi \left(x^{\left(i\right)} \right)$ とすると、

超平面と、そのサンプルの距離は

$d_i =\frac{\left|\omega^{\mathrm{T}} \phi \left(x^{\left(i\right)} \right)\right|}{\|w\|}$

...(30)

で表すことができます。

さて、冒頭の

に合わせて、

$f\left(x\right)=\omega^{\mathrm{T}} \phi \left(x\right)$ の値が1より大きい場合は、ラベル1、

小さい場合はラベル-1を割り当てるとします。（こうなるようにωの値を等倍してうまく選ぶようにします）

このとき、式(30)で定義される、分離境界（超平面）と訓練データのサンプルの距離の最小値は、

$f\left(x\right)=\omega^{\mathrm{T}} \phi \left(x\right)=1$ のときの、

$\frac{1}{\left\|w\right\|}$

...(31)

となることがわかります。

このマージンを最大化する、を、別の言い方で表現すると、この式の逆数の二乗である、 ${\left\|w\right\|}^2$ を最小化することになります。

この状態では、すべての点をうまく分類するための、超平面を作成することを念頭に置いています。

しかし、カーネル法による回帰の記事にもあったように、外れ値に対しても敏感に反応する超平面を作成すると、

他のデータに対して対応できない（過学習）ものになってしまいます。

そのため、誤分類を許さない今回の想定（ハードマージン）ではなく、いくつかの誤分類は許容し（ソフトマージン）、より滑らかな超平面を作成するために、正則化項に相当するものを追加します。

式(14)のスラック変数と同様に、

${y^{\left(i\right)} \omega }^{\mathrm{T}} \phi \left(x^{\left(i\right)} \right)>1-\xi_i$

...(32)

という制約を設けて、

$\min_{\omega ,\xi } {\left\|\omega \right\|}^2 +\frac{1}{\lambda }\sum_{i=1}^n \xi_i$

...(33)

を満たす、ωを見つければよいことになります。λは、誤分類が発生した時に加えるペナルティ項の係数（正則化項の係数）です。

ここで、

$\omega =\sum_{i=1}^n \alpha \varphi \left(x^{\left(i\right)} \right)$

...(34)

と置くと、

${\left\|\omega \right\|}^2 =\alpha^{\mathrm{T}} K\alpha$

...(35)

となります。

4.1で述べた、最小化する損失関数は以下のようになっていました。

$min\sum_{i=1}^n \left(\xi_i +\frac{\lambda }{2}\alpha^T K\alpha \right)$

...(16の再掲)

式(33)を式(35)を利用して変形すると、

$\min_{\omega ,\xi } {\left\|\omega \right\|}^2 +\frac{1}{\lambda }\sum_{i=1}^n \xi_i$

$\underset{\omega ,\xi }{=\min } {\;\alpha }^{\mathrm{T}} K\alpha +\frac{1}{\lambda }\sum_{i=1}^n \xi_i$

$=\frac{1}{\lambda }\min \left(\sum_{i=1}^n \xi_i +{\;\lambda \;\alpha }^{\mathrm{T}} K\alpha \right)$

...(36)

となり、 $\lambda =\frac{1}{2}\lambda^{\prime }$ とすれば

4章で示した方法と、等価な式を得ることができました。

このように、「損失関数＋正則化項」を最小化する枠組みで説明した場合でも、「マージンの最大化」を考えた場合も、同じ計算をしていることになります。おもしろいですね。

7. 複数クラスの分類について

これまで述べてきたように、サポートベクトルマシンは2つのクラスの分類に利用可能です。しかし、3つ以上のクラスに分類したい場合も多く存在します。

多クラスの分類として知られている方法の例を挙げます。

one vs all (1対他方式)
one vs one (1対1方式)
誤り訂正出力符号

図出展：カーネル多変量解析 p.119

7.1. one vs all (1対他方式)について

上図の(a)にあるように、各クラスとその他のクラス、という分類器をそれぞれ用意します。

利点

シンプルでわかりやすい。

欠点

全ての分類器で、all側になる場合や、複数のパターンで、one側になる（1つの分類器だけ、そのクラスを支持してくれると嬉しいがそうならない場合が多い。分離境界までの距離を基準にすることも考えられるが、サンプル数がクラス間で異なっていたりすると、フェアに判断することが難しくなる。）
クラス間のサンプル数が不均衡である場合、学習が難しくなる。

7.2. one vs vs one (1対1方式)について

上図の(b)にあるように、1クラスどうしを比較するように、学習器をたくさん用意する。

利点

one vs allのようなサンプル数が不均衡である場合でも対応しやすくなる。

欠点

クラスの数が増えるに従い、用意すべき学習器の数が膨大になる。

7.3. 誤り訂正出力符号を用いた方法

詳しくは以下のサイトをご覧ください。

朱鷺の杜Wiki　誤り訂正出力符号

https://ibisforest.org/index.php?%E8%AA%A4%E3%82%8A%E8%A8%82%E6%AD%A3%E5%87%BA%E5%8A%9B%E7%AC%A6%E5%8F%B7

サポートベクトルマシンの多クラス分類への適用

https://deus-ex-machina-ism.com/?p=23045

8. MATLABでのパラメータ設定を見てみよう

サポートベクトルマシンの仕組みや、そこで必要な設定値をまとめてきました。実際に、サポートベクトルマシンを実行するときに、以上のパラメータ調整が可能であることを確認していきます。

今回は、ドキュメントが充実している、MATLABを例に考えます。以下に、主要なものを挙げます。

以下のドキュメントを参照します。

fitcsvm関数

2クラスの分類、または、1クラス (one-class SVM。本記事では、紹介しません)の分類に利用可能なサポートベクトルマシンの関数。

https://jp.mathworks.com/help/stats/fitcsvm.html

fitcecoc関数

「サポートベクターマシンまたはその他の分類器向けのマルチクラスモデルの近似」

https://jp.mathworks.com/help/stats/fitcecoc.html

8.1. ボックス制約

以下のような説明があります。

ボックス制約は、マージンに違反している観測値に課せられる最大ペナルティを制御するパラメーターであり、過適合の防止 (正則化) に役立ちます。ボックス制約の値を大きくすると、SVM 分類器が割り当てるサポートベクターは少なくなります。ただし、ボックス制約の値を大きくすると、学習時間が長くなる場合があります。

式(12)で導入した、正則化項の係数であると思われます。これを大きくすることで、過適合（overfitting)を防ぐことができますが、同時に、分離境界の超平面が滑らかになりすぎて、認識能力が下がる可能性もあるのでうまくバランスをとる必要があります。

8.2. カーネル関数

式(11)で導入した、カーネル関数の中身のことです。再掲になりますが、以前のカーネル関数の記事もご覧ください。i番目とj番目のサンプルを何らかのルールで計算させてグラム行列を作るのでした。私のブログ記事では、ガウシアンカーネルを利用していました。線形カーネル、多項式カーネル、カスタマイズしたカーネルなど、色々と利用が可能です。

kentapt.hatenablog.com

8.3. カーネル関数のスケール

上の過去のブログ記事でも、ガウシアンカーネルのパラメータについて記述していました。

上の記事では、式3-1に相当する内容です。

例えば、ガウシアンカーネルを以下のように設定したとします。

その場合は、以下のように設定するとよいそうです。

mdlSVM = fitcsvm(..., 'KernelScale', 1/sqrt(gamma));

出展：MATLAB answers "Gaussian kernel scale for RBF SVM"

https://jp.mathworks.com/matlabcentral/answers/328796-gaussian-kernel-scale-for-rbf-svm

8.4. 最適化ルーチン

式(28)周辺で、 ${0<\gamma }_i <1$ の条件下で、

$\sum_{i=1}^n \gamma_i -\frac{1}{2\lambda }\sum_{i=1}^n \sum_{j=1}^n y^{(i)} y^{(j)} \gamma_i \gamma_j K_{ij}$

を最大化する $\gamma_i$ を求めるという、凸二次計画問題を解けばよい、ということを示しました。

本文中で説明した、SMOに加え、他のアルゴリズムも利用可能です。

8.5. キャッシュサイズ

カーネル法の際に用いる、グラム行列Kは、サンプルの数が多くなればなるほど、2乗のスケールで大きくなり、計算が大変そうだという、話をしていました。サポートベクトルマシンの計算をする際に、確保するメモリサイズに関するパラメータだと思います。

8.6. 多クラス分類をする際の方式

7章で述べたような、本来2クラスの分類に利用可能なサポートベクトルマシンを用いて、多クラスの分類をする際の方式です。上で述べた、one vs oneやone vs allに加えて、色々な方式が利用可能です。

8.7. Standarize

入力データを標準化する/しないを設定するもの。

8.8. 誤分類のコスト

例えば、AというものをBと間違えてしまった時のコスト（ペナルティ）の値を制御します。例えば、犬と猫を見分ける場合を考えます。訓練データ内に犬が90匹、猫が10匹存在するとします。その場合、学習せずとも、「全て犬である」と答えると精度は90%になってしまいます。このようなクラス間でデータ数が不均衡である場合などに、このコストを調整するといいと思われます。例えば、猫なのに、犬と訓練すると、非常に大きなペナルティ（損失）を与えるようにします。簡単にするために、クラス数の逆比を用いて、猫を犬と判断した場合に、9倍のペナルティを与えるようにすると、このような「全て犬という」事態はさけられそうです。

9. まとめ

かなり長くなりましたが、ここまで読んでいただき誠にありがとうございました。本記事では、

サポートベクトルマシンの概要と、その学習時の計算について説明を行いました。
損失関数と正則化項の最適化に焦点を当てて、サポートベクトルマシンの計算方法について述べました
別の切り口である、「マージン最大化」で説明を行っても、等価な数式に帰着することを述べました
MATLABを例に、ソフトウェアで実行するときのパラメータについて述べました
この記事が、サポートベクトルマシンを理解するうえでのきっかけとなれば幸いです。

注意点としては、

学習の基本的なところは、単回帰や、重回帰に関する過去の記事でカバーしていますのでそちらもご覧ください。
カーネル法に関しても、過去のカーネル法による回帰の記事でカバーしていますのでそちらもご覧ください。
筆者も細部については、勉強中のところが多く、わかりにくかったり、修正すべき点があるかもしれません。その場合は教えていただけますと幸いです。

参考文献

「カーネル多変量解析」赤穂昭太郎先生

honto.jp

「サポートベクターマシン」赤穂昭太郎先生 https://www.ism.ac.jp/~fukumizu/ISM_lecture_2006/svm-ism.pdf

「サポートベクターマシン入門」栗田多喜夫先生 https://home.hiroshima-u.ac.jp/tkurita/lecture/svm.pdf

不等式条件下におけるラグランジュの未定乗数法

laid-back-scientist.com

SVMの2次計画問題に関する解法の考察

http://numaf.net/Z9/Z9a/html/THESIS/H15/abst_toda.pdf

サポートベクターマシンとは［ソフトマージンサポートベクターマシン］

sudillap.hatenablog.com

最初に学ぶ 2分類SVMモデルの基本

www.acceluniverse.com

本文中に参考文献としての記載のあるものは、ここでは一部省略しています。

2022-12-30

重回帰分析の勉強＆実装による検算をしてみよう

統計機械学習 matlab

この記事は、MATLAB アドベントカレンダー（その2）の6日目の記事として書かれています。

qiita.com

1 はじめに

前回の記事では、1つの変数のみで、線形回帰（単回帰）を行うときの、数式の導出や自身の実装による検算を行いました。

線形回帰（単回帰分析）を1から実装して理解を深めてみよう

kentapt.hatenablog.com

線形回帰（単回帰分析）を1から実装して理解を深めてみよう: つづき

kentapt.hatenablog.com

そこでは、y = ax + bという非常にシンプルな場合を扱いました。

今回は、

$y~=~a~+~b_1 x_{1~} +~b_2 x_2 ~+~\cdots +~b_k x_k$ 　

といったような、変数の数を増やした場合について考えたいと思います。

そこで、本記事では、変数の数を増やして線形回帰（重回帰分析）をしたときの、計算方法について簡単にまとめます。さらに、MATLABを用いてコーディングをし、ここでの理解による実装の結果と、MATLABのregress関数の結果が一致することを確かめます。本記事は、タイプミスや間違いなどがあるかもしれません。その場合教えていただけますと幸いです。また、本記事は、前回の記事をもとに、執筆をしており、説明が比較的少なくなっています。もしよろしければ、先述した前回の記事も見ていただけますと幸いです。

また、本記事の原稿や用いたコードは以下のページにアップロードされています。

github.com

本記事の執筆においては、以下の資料が非常に参考になりました。

回帰分析（重回帰） http://fs1.law.keio.ac.jp/~aso/ecnm/pp/reg2.pdf

2 重回帰分析について

重回帰モデル (multiple regression model) とは、単回帰とは異なり、説明変数が複数ある回帰モデルのことです。

モデル式は、線形回帰に変数をさらに加えて以下のようになります。

$y~=~a~+~b_1 x_{1~} +~b_2 x_2 ~+~\cdots +~b_k x_k$

...(1)

単回帰のときと同様に、誤差を2乗して、足したもの（残差平方和）を求めます。

...(2)

前回の記事と同様に、各変数で偏微分を行います。

$X~=~a~+~b_1 x_{1~} +~b_2 x_2 ~+~\cdots +~b_k x_k$

のようにおいて、合成関数の微分を行っています。

aおよび $b_i$ に対して、偏微分を行った時の式を以下に記述します。

$-2\sum_{i=1}^n \left(y_i -a-b_1 x_1 -b_2 x_2 \cdots -b_k x_k \right)=0$

$-2\sum_{i=1}^n \left(y_i -a-b_1 x_1 -b_2 x_2 \cdots -b_k x_k \right)x_1 =0$

$-2\sum_{i=1}^n \left(y_i -a-b_1 x_1 -b_2 x_2 \cdots -b_k x_k \right)x_2 =0$

$\cdots$

$-2\sum_{i=1}^n \left(y_i -a-b_1 x_1 -b_2 x_2 \cdots -b_k x_k \right)x_j =0$

$\cdots$

$-2\sum_{i=1}^n \left(y_i -a-b_1 x_1 -b_2 x_2 \cdots -b_k x_k \right)x_k =0$

...(3)

yを左辺に移動させて、両辺を-2で割り、整理すると以下のようになります。

この式を正規方程式（normal equation）と呼びます。

$\sum_{i=1}^n y_i =\;\sum_{i=1}^n \left(a+b_1 x_1 +b_2 x_2 \cdots +b_k x_k \right)$

$\sum_{i=1}^n x_1 y_i =\;\sum_{i=1}^n \left(ax_1 +b_1 x_1 \times x_1 +b_2 x_2 \times x_1 \cdots +b_k x_k \times x_1 \right)$

$\cdots$

$\sum_{i=1}^n x_j y_i =\;\sum_{i=1}^n \left(ax_j +b_1 x_1 \times x_j +b_2 x_2 \times x_j \cdots +b_k x_k \times x_j \right)$

$\cdots$

$\sum_{i=1}^n x_j y_i =\;\sum_{i=1}^n \left(ax_j +b_1 x_1 \times x_k +b_2 x_2 \times x_k \cdots +b_k x_k \times x_k \right)$

...(4)

なお，変数の数が1つの単回帰の場合については，(4)の最初の方程式を他の方程式に代入することで、連立方程式を解くことができました。しかし、変数が多くなると、このように解析に解くことは非常に困難です。そこで行列を用いて解いていきます。

(4)の左辺を行列式で表すと以下のようになります。

$\left\lbrack \begin{array}{c} y_1 +y_2 +\cdots +y_n \\ x_1 y_1 +x_1 y_2 +\cdots +x_1 y_n \\ \cdots \\ x_k y_1 +x_k y_2 +\cdots +x_k y_n \end{array}\right\rbrack =\left\lbrack \begin{array}{cccc} 1 & 1 & \cdots & 1\\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & \vdots & \vdots \\ x_{k1} & x_{\textrm{k2}} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack \;\left\lbrack \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right\rbrack ={\left\lbrack \begin{array}{cccc} 1 & x_{11} & \cdots & x_{k1} \\ 1 & x_{12} & \cdots & x_{\textrm{k2}} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1n} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack \;\;}^{\prime } \;\left\lbrack \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array}\right\rbrack =X^{\prime } y$

...(5)

なお、 $y_i$ はi番目のサンプルのyの値を示し、 $x_{\textrm{kl}}$ は、k番目の変数のl番目の値を示しています。

また、

$X=\left\lbrack \begin{array}{cccc} 1 & x_{11} & \cdots & x_{k1} \\ 1 & x_{12} & \cdots & x_{\textrm{k2}} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1n} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack$

と置いています。

さらに、(4)の右辺を行列式で表します。

$\left\lbrack \begin{array}{cccc} n & \sum_{i=1}^n x_{1i} & \cdots & \sum_{i=1}^n x_{\textrm{ki}} \\ \sum_{i=1}^n x_{1i} & \sum_{i=1}^n {x_{1i} }^2 & \cdots & \sum_{i=1}^n x_{1i} x_{\textrm{ki}} \\ \vdots & \vdots & \cdots & \vdots \\ \sum_{i=1}^n x_{\textrm{ki}} & \sum_{i=1}^n x_{\textrm{ki}} x_{1i} & \cdots & \sum_{i=1}^n {x_{\textrm{ki}} }^2 \end{array}\right\rbrack \;\left\lbrack \begin{array}{c} a\\ b_1 \\ \vdots \\ b_k \end{array}\right\rbrack =\left\lbrack \begin{array}{cccc} 1 & 1 & \cdots & 1\\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & \vdots & \vdots \\ x_{k1} & x_{\textrm{k2}} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack \;\left\lbrack \begin{array}{cccc} 1 & x_{11} & \cdots & x_{k1} \\ 1 & x_{12} & \cdots & x_{\textrm{k2}} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1n} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack \;\left\lbrack \begin{array}{c} a\\ b_1 \\ \vdots \\ b_k \end{array}\right\rbrack$

...(6)

(6)のようにきれいな形で表すことができました。

さらに、(6)の式をXを用いて表すと、以下のようになります。

$={\left\lbrack \begin{array}{cccc} 1 & x_{11} & \cdots & x_{k1} \\ 1 & x_{12} & \cdots & x_{\textrm{k2}} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1n} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack }^{\prime } \;\left\lbrack \begin{array}{cccc} 1 & x_{11} & \cdots & x_{k1} \\ 1 & x_{12} & \cdots & x_{\textrm{k2}} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1n} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack \;\left\lbrack \begin{array}{c} a\\ b_1 \\ \vdots \\ b_k \end{array}\right\rbrack =X^{\prime } X\;b$

...(7)

(4)の左辺は(5)で、(4)の右辺は(6)のようになりました。

そこで、それらを合わせて、

$X^{\prime } y=X^{\prime } X\;b$

...(8)

となります。

$X^{\prime } X$ の逆行列が存在すれば、

$b={\left(X^{\prime } X\right)}^{-1} X^{\prime } y$

..(9)

となり、係数の値（最小二乗推定量）を得ることができました。

3 回帰直線がサンプルの平均を通ることの確認

単回帰分析の記事で、回帰直線がサンプルの平均を通ることを示しました。

kentapt.hatenablog.com

重回帰分析でも同様です。

(4)の1つ目の式で以下のことが示されていました。

$\sum_{i=1}^n y_i =\;\sum_{i=1}^n \left(a+b_1 x_1 +b_2 x_2 \cdots +b_k x_k \right)$

この式をサンプル数nで割ると、以下の式を得ることができます。

${\hat{y} =a+b}_1 {\bar{x} }_1 +b_2 {\bar{x} }_2 +\cdots +b_k {\bar{x} }_k$

...(10)

なお、yの平均を $\hat{y}$ 、xの平均を $\bar{x}$ としています。

この式から、回帰直線は、必ずサンプルの平均（重心）点を通ることがわかります。

また、単回帰分析の場合と同様に、残差の和は０であることや、残差ベクトルとが説明変数xは直交することも導くことができます。

4 重回帰分析を自分で実装して検算してみよう

3章まででまとめた重回帰分析に対する理解が正しいことを確認するために、4章では、自分で重回帰分析を実装します。まずは、MATLABを用いて計算し、さらに自分の手で実装します。そして、それらの結果が一致することを確認します。

4.1. MATLABのfitlm関数による計算

はじめに、MATLABの関数を用いて、重回帰分析を計算します。carsmallデータをロードし、さらに、fitlm関数を利用します。

load carsmall
X = [Weight,Horsepower,Acceleration];
% fitlm を使用して、線形回帰モデルを当てはめます。
mdl = fitlm(X,MPG)

mdl = 
線形回帰モデル: 
    y ~ 1 + x1 + x2 + x3

推定された係数: 
                    Estimate        SE          tStat        pValue  
                   __________    _________    _________    __________

    (Intercept)        47.977       3.8785        12.37    4.8957e-21
    x1             -0.0065416    0.0011274      -5.8023    9.8742e-08
    x2              -0.042943     0.024313      -1.7663       0.08078
    x3              -0.011583      0.19333    -0.059913       0.95236

観測数: 93、誤差の自由度: 89
平方根平均二乗誤差: 4.09
決定係数: 0.752、自由度調整済み決定係数: 0.744
F 統計量 - 定数モデルとの比較: 90、p 値は 7.38e-27 です

このように、切片47.977、および、各変数の係数を得ることができました。

4.2. 自分の手で実装してみる

Xは以下のように定義されました。

$X=\;\left\lbrack \begin{array}{cccc} 1 & x_{11} & \cdots & x_{k1} \\ 1 & x_{12} & \cdots & x_{\textrm{k2}} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{1n} & \cdots & x_{\textrm{kn}} \end{array}\right\rbrack$

縦に、サンプル、横に変数が並んでいて、データセットと同じ形をしています。そのため、

一番左の列に1を足すだけでよく、以下のようにするとXを定義できます。

validIdx = ~isnan(MPG)&~isnan(X(:,2));
X = X(validIdx,:);
MPG = MPG(validIdx);
myX = [ones(size(X,1),1),X];

さらに、重回帰分析の式は以下のように解くことができることがわかりました。

$b={\left(X^{\prime } X\right)}^{-1} X^{\prime } y$

さきほどのXと合わせて、以下のように計算できます。

b = inv(myX'*myX)*myX'*MPG

すると、4.1と同じ結果を得ることができました。

4.3. 重回帰分析を3次元プロットする

4.2では3変数を利用しましたが、次は、2変数のみを利用し、その結果を3次元プロットします。

以下のページを参考にしています。

https://jp.mathworks.com/help/stats/regress.html

load carsmall
x1 = Weight;
x2 = Horsepower;    % Contains NaN data
y = MPG;
% 
X = [ones(size(x1)) x1 x2];
b = regress(y,X)

scatter3(x1,x2,y,'filled')
hold on
x1fit = min(x1):100:max(x1);
x2fit = min(x2):10:max(x2);
[X1FIT,X2FIT] = meshgrid(x1fit,x2fit);
YFIT = b(1) + b(2)*X1FIT + b(3)*X2FIT;
mesh(X1FIT,X2FIT,YFIT)
xlabel('Weight')
ylabel('Horsepower')
zlabel('MPG')
view(50,10)
hold off

上の図のように、説明変数をx、y軸とし、目的変数をz軸としてプロットすることができました。

5 まとめ

本記事では、重回帰分析について簡単にまとめました。行列式を用いることで、単回帰と同じ流れで簡単に導出することができました。さらに、MATLABでも実装し、自身の理解が正しいことを確認できました。

今回は、具体的な利用については述べませんでした。重回帰分析は様々な変数を入れることができるため、便利です。しかし多重共線性などについても考慮しながら、うまく運用することが重要です。今後も引き続き勉強し、うまくこのようなモデルを扱えるようになりたいです。

多重共線性への対応として、主成分分析を利用することも考えられます。主成分分析により、各変数を無相関化することができます。主成分分析については、以下の記事をご覧ください。

kentapt.hatenablog.com

sorb.co.jp

参考文献（本文中に記載のないもの）

Albertさまデータ分析基礎知識サイト

www.albert2005.co.jp

丹羽時彦さま放課後の数学入門

www.kwansei.ac.jp

MATLAB ドキュメンテーション fitlm

jp.mathworks.com

2022-12-22

iPhone LiDARで取得した3次元点群の境界を書き出してみよう

3次元点群 GIS matlab

この記事は、MATLAB アドベントカレンダー22日目の記事として書かれています。

qiita.com

この記事の前編として以下の記事も投稿しています。ぜひご覧ください。

kentapt.hatenablog.com

この記事では、iPhone LiDARで取得した3次元点群をMATLABにて読み込み、さらに、その領域の境界のデータをKMLファイルとして出力する方法を紹介します。言語はMATLABを使用します。

3次元点群を取得するためのアプリは3D Scanner Appを用いました。

詳しい使い方は、以下の記事がわかりやすかったです。

note.com

また、以下の記事でもこちらのアプリの使い方を紹介しています。

note.locusblue.com

ここでは、すでにiPad LiDARにて計測し、PCに転送した状態から始めます。

なお、ここで用いたコードや、点群ファイルは以下のページにアップロードされています。

github.com

点群ファイルの読み込み

LiDAR toolboxの機能を用いて、点群データの読み込みを行います。以下のように、lasFileReaderやreadPointCloud関数を用いて、簡単に点群ファイルを読み込むことができます。まずは、GPS情報を有した点群（LAS）ファイルをそのまま読み込み、表示してみます。なお、上のgithubのページでは、ファイルサイズを小さくするため、サンプルファイルのファイル形式をLAZ形式に圧縮してアップロードしています。（3D Scanner Appでは、LAS形式で点群がエキスポートされます。）MATLABに読み込む方法などは変わりません。

clear;clc;close all
% las/lazファイルを読み込むためのlasFileReader objectを作成する
lasReader = lasFileReader('sample.laz');
% las/lazファイルの、xyzおよび色情報を読み込み
ptCloud = readPointCloud(lasReader);
% 可視化
figure;pcshow(ptCloud)

しかし、単に読み込むだけでは、このように、うまく表示できません。2022年12月現在、3D Scanner Appから、Geo-referenceされた（位置情報を持っている）LASファイルをエキスポートすると、その座標は経度・緯度の値を持っています。経度緯度の値は、小数点以下の数が非常に多く、単にビューワーに読み込むだけではうまく行かない場合も多いです。

座標系の変換

geocrsオブジェクトの作成

MATLABのpcshow関数で点群を閲覧するために、経度緯度のxy座標を投影座標に変換します。

% geojcrsオブジェクトを作成する
geocrs4326 = geocrs(4326);

ここで、作成した、geocrsオブジェクトの中身について確認します。

% 確認のため、geocrsを表示
disp(geocrs4326)

  geocrs のプロパティ:

             Name: "WGS 84"
            Datum: "World Geodetic System 1984"
         Spheroid: [1x1 referenceEllipsoid]
    PrimeMeridian: 0
        AngleUnit: "degree"

測地系（Datum）の情報や、Spheroid（楕円体）の情報格納されています。

PrimeMeridianは、本初子午線のことで、ここの数字は、緯度の基準が旧グリニッジ天文台からどれくらいの差分があるかを示しています。今回は0なので特に気にする必要はありません。

また、Spheroidをクリックすると、その中身を確認することができます。

一部のみを紹介します。

Code: 割り当てられている識別コード EPSG:7030
SemimajorAxis: 赤道半径 6,378,137 m
InverseFlattening: 扁平率の逆数 298.257223563

扁平率とは、（今回は）地球を表す楕円体が、完全な球と比べた時にどれくらい、扁平であるか（つぶれているか）を示します。全くつぶれていない（＝球と同じ）場合は、0となり、つぶれているほど、1に近づきます。

楕円体の長半径をa、短半径をbとすると、扁平率は、

$1-\frac{b}{a}\;=\frac{a-b}{a}\;$

で表すことができます。

実際にその逆数を計算してみると、

f = 1-geocrs4326.Spheroid.SemiminorAxis/geocrs4326.Spheroid.SemimajorAxis;
1/f

ans = 298.2572

となり、InverseFlatteningの値と一致することが確認できます。

また、geocrsオブジェクトを用いて、そのWKTを簡単に作成することもできます。

WKT（Well-Known Text）は、ジオメトリの情報をテキスト形式で記述するものです。

str = wktstring(geocrs4326)

str = "GEOGCRS["WGS 84",DATUM["World Geodetic System 1984",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]"

ここでも、さきほどの、geocrsオブジェクトのような情報が確認できます。

projfwd関数を用いた、座標の変換

それでは、経度緯度の座標を投影座標に変換したいと思います。

このサンプルデータは、東京都で取得されたため、EPSGコードは6676を利用します。

変換自体は、projfwd関数を用いて簡単に行うことができます。

lon = ptCloud.Location(:,1);
lat = ptCloud.Location(:,2);
proj = projcrs(6676);
[x,y] = projfwd(proj,lat,lon);

pcshow関数を用いた、点群の閲覧

pointCloudProj = pointCloud([x,y,ptCloud.Location(:,3)],"Color",ptCloud.Color);
figure;pcshow(pointCloudProj)

以下の動画のように、MATLAB上でiPhone LiDARにて取得した点群を可視化することができました。手軽にこのようなきれいな3D情報が取れてすごいですね。

3次元マップでの点群の可視化

次に、投影したxy座標をもとに、iPhone LiDARで取得した点群を、3次元のマップ上で可視化してみます。

このセクションの内容は、以前の投稿と同様です。詳細は以下の記事をご覧ください。

kentapt.hatenablog.com

% 25点のうち1点のみを表示する。重くなってしまうため。
mskip = 1:50:length(lat);
% 3Dで可視化するための図を用意する
uif = uifigure;
g = geoglobe(uif);
% プロット
geoplot3(g,lat,lon,ptCloud.Location(:,3)-mean(ptCloud.Location(:,3)),'r','MarkerSize',2,'Marker','o','HeightReference','terrain','LineStyle','none','MarkerIndices',mskip)

境界の表示の書き込み

さきほどは、点群ファイル全体を3Dのマップ上に表示させました。次に、この点群ファイルの境界を計算し、その境界のみを2Dの地図上に可視化してみようと思います。

境界の計算は、boundary関数で簡単に行うことができます。

ここで、計算を軽くするために、さきほど用いた間引きの間隔を用いて、50点ごとに境界を計算する点を抽出します。

% 点の間引き
latSub = lat(mskip);
lonSub = lon(mskip);
% 境界の計算
k = boundary(latSub,lonSub);
% 可視化は、geoplot関数が便利です
figure;geoplot(latSub(k),lonSub(k))
% ベースマップの種類を指定できます
geobasemap streets

kmlファイルへの書き込み

境界をkmlファイルとして書き込みを行います。

filename = 'myBoundary.kml';
kmlwriteline(filename, latSub(k), lonSub(k), 'Color','black','LineWidth', 3);

Google Earth Proでの表示

Windows環境で、かつ、Google Earth Proをインストールしている場合、以下のコマンドで、さきほど保存したkmlファイルをGoogle Earth上で確認することができます。Google Earth Proを別途開くことなく、MATLABからコマンド1行で起動できるのは非常に便利ですね。

Google Earth Pro（パソコン用）のダウンロードページ

www.google.com

winopen(filename)

以下の動画（GIF）では、winopen(filename)とコマンドウィンドウで入力しています。さきほど書き出した点群の境界をGoogle Earth上で可視化できていることがわかります。

まとめ

この記事では、iPhone LiDARで取得した点群をMATLABを用いて、地図上で可視化しました
投影座標系から経度・緯度への変換も簡単に行うことができました
3Dのベースマップ上で、点群を可視化することができました
表示する点数が多くなると、非常にビューワーが重くなるため、間引くことが必要です

参考にしたページ（本文中に記載がないもの）

とらりもん：測地系 (datum)

pen.envr.tsukuba.ac.jp

Wikipedia: 扁平率

ja.wikipedia.org

tohka383様：WKT にみる座標系について

qiita.com

MATLAB ドキュメンテーション：Transform Coordinates to a Different Projected CRS

jp.mathworks.com

MATLAB ドキュメンテーション：kmlwriteline

jp.mathworks.com

2022-12-18

深層学習を用いた物体検出の手法（YOLOv2）についてのまとめ

深層学習画像解析 matlab

この記事は、MATLAB アドベントカレンダー（その2）の21日目の記事として書かれています。

qiita.com

はじめに

この記事では、物体検出の有名な手法である、YOLOv2について説明を行います。物体検出の手法を用いて、以下のように画像中から対象を自動的に検出することができます。

YOLOにより物体を検出しているときの動画の例も挙げます。以下の動画をご覧ください。

YOLOモデルは物体検出において非常に有名で、様々なバージョンが公開されています。この中でも、ネットワークがシンプルで、比較的高速かつ高精度に物体検出ができる、YOLOv2についてまとめたいと思います。YOLOなどの物体検出ネットワークの仕組みを説明した記事は多く存在します。本記事では、もとの論文の和訳をして、作者による説明を把握しつつ、適宜、補足説明を行います。これにより、論文をベースにした情報に触れながらも、（まずはぼんやりと）内容の理解の進む記事になればよいなと思いました。

　そこで、本記事では、以下のリンクにある Redmon and Farhadi (2017) の論文の、Introductionと、Betterのセクションの和訳を行い、さらに、適宜、筆者の理解による補足を行います。特にBetterセクションでは、図を追加して、YOLOv2の原理について説明を行います。これにより、元の論文の流れをベースとして、YOLOv2の概要をつかみやすくすることを目指します。本記事を通して、物体検出の方法の理解度が上がれば嬉しく思います。なお、本論文では、以下の論文から図を引用しております。

[論文情報]

Redmon, J., & Farhadi, A. (2017). YOLO9000: better, faster, stronger. In Proceedings of the IEEE conference on computer vision and pattern recognition (pp. 7263-7271).

[論文のリンク]

arxiv.org

なお、本記事で用いたコードや本文の原稿は以下のページにアップロードされています。

github.com

論文の和訳

それでは、さっそく論文のIntroductionから和訳を行っていきます。また、和訳に関しては、筆者の勝手な理解に基づいた意訳を多く含みます。ご注意ください。（）内の内容は筆者による補足説明です。

Introduction

物体検出は、高速かつ正確で、さまざまな物体を認識できる必要がある。ニューラルネットワークの利用が広まり、物体検出のネットワークは、高速かつ高精度になってきている。しかし、ほとんどの物体検出手法は、まだ少数のオブジェクトに限定されており、現在の物体検出のデータセットは、分類やタグ付けといった他のタスクのデータセットに比べ、（データ数が）限られている。最も一般的な物体検出用データセットは、数千から数十万の画像のタグが含まれている[3] [10] [2]。一方で、分類データセットには、数百万の画像と数万から数十万のカテゴリが含まれる[20] [2] 。（このように、物体検出用のデータセットの量が少ないことがわかる。）

　本論文では、検出を物体分類のレベルにまで拡張することを考える。しかし、物体検出のためのラベリングは、分類やタグ付けのためのラベリングよりもはるかに時間とコストがかかる。そのため、分類と同規模の物体検出用データセットが近い将来に登場する可能性は低い。我々は、既に存在する大量の分類データを利用し、現在の検出システムの適用範囲を拡大する新しい手法を提案する。本手法は、異なるデータセットを一緒に組み合わせるために、物体の分類ラベルに対して階層的な見かたをする（例：食べもの⇒お肉系⇒焼肉）。また、物体検出と画像分類の両方のデータセットを用いて、物体検出のためのネットワークを学習することが可能なアルゴリズムも提案する。本手法は、（分類用の）ラベル付けをされたデータセットを利用して物体の正確な位置情報の推定方法を学習する。それに加え、分類のためのデータセットを利用して分類できる画像の種類を増やしたり、頑健性を向上させる。この方法を用いて、9000以上の異なる物体のカテゴリを検出でき、さらにリアルタイムで推論のできるネットワーク（YOLO9000とよぶ）を学習させる。まず、ベースとなるYOLOとよばれるネットワークを改良し、YOLOv2の開発を行った。次に、ImageNetの9000クラス以上の画像分類のデータや物体検出のデータを用いてモデルの学習を行う。

Better

YOLOは、その他の物体検出手法と比較して、様々な欠点を抱えている。YOLOとFast R-CNNを比較した場合、YOLOはより多くの位置推定の誤差を有することがわかった（物体を検出したとしても、その位置に誤差が含まれる）。さらに、YOLOはregion proposal-basedな手法と比較して、再現率が低い（検出漏れが多い）ことがわかった。したがって、YOLOv2では分類精度を維持しつつ、主に再現率と位置推定の精度を向上させることを目指す。（以降、再現率をrecallとする）

　近年は深層学習を用いた物体検出の精度を向上させるために、ネットワークをより大きく、より深くする傾向にある[6] [18] [17]。より良い性能を得るためには、複数のモデルを構築し、それらの結果の平均などを取ることもある（アンサンブル）。しかし、アンサンブルは計算コストも多くなるため、YOLOv2では、より正確かつ高速なネットワークを実現する。ネットワークのサイズを大きくする代わりに、ネットワークを単純化する。今回の提案では、過去の研究から得たアイデアと独自のコンセプトを融合させた。それらの結果の要約を表2に示す（この表の見方としては、チェックがあるところが、YOLOと比べて新たに追加した技術だと思います。例えば、左から2つ目の列だと、Batch Normalizationを追加した時はmAPが65.8で、そうでない場合は63.4である、といった形です）。

Batch Normalization

バッチ正規化により、他の正則化手法を使うことなく、学習を大幅に収束させやすくなる[7]。YOLOに含まれる全てのの畳み込み層の後にバッチ正規化層を加えることで、mAPが2%以上改善された。また、バッチ正規化はネットワークの正則化に役立つ。バッチ正規化により、ネットワークの過学習を防ぎ、ネットワークからドロップアウト層を除去することもできる。

バッチ正規化については、以下の記事がわかりやすかったです。

cvml-expertguide.net

High Resolution Classifier

最先端の検出手法は、すべてImageNet [16]で事前に学習させた分類器を用いている。AlexNetを始めとして、ほとんどの分類器は、入力の画像のサイズが、256×256より小さいことが必要である [8]。オリジナルのYOLOでは，224×224の入力サイズで分類器ネットワークを学習し、検出のために画像サイズがを448×448になるように調整している。つまり、物体検出においては、画像分類から物体検出という異なるタスクのために、新たに学習が必要になるだけでなく、新しい入力解像度に適応する必要がある。そこで、YOLOv2では、まずImageNet上で10エポックの間、448×448の画像サイズで分類器ネットワークを微調整する。これにより、ネットワークはより画像サイズの大きい場合にも学習できるように各種重みなどをを調整するしやすくなる。この操作により、4%近いmAPの増加が得られた。

Convolutional With Anchor Boxes (その1)

YOLOは、畳み込み層の後に全結合層を置き、バウンディングボックスの座標を直接予測する。一方、Faster R-CNNでは座標を直接予測するのではなく、手動で選んだ事前分布を使用してバウンディングボックスを予測する[15]。そして、領域提案ネットワーク（RPN）を用いて、そのボックスのオフセットと信頼度を予測する。全結合層を用いて、直接バウンディングボックスの座標を予測するものとは異なり、畳み込み層の各グリッドが各ボックスのオフセットや信頼度を計算している。

==========論文の和訳をひとやすみ==========

補足事項1: YOLOv2のおおまかな流れについて

ちょうど上の、Convolutional With Anchor Boxesという章の前半で、YOLOv2の処理の内容が少し記述されていました。

そこで、論文の和訳は中断し、ここで、YOLOv2の大まかな流れを以下の図をもとに説明します。そして、このおおまかな流れの説明の後で、論文に戻り、どのような経緯でこのような方法に行きついたかを見ていきたいと思います。

1) 各画像をS×Sの領域に分割する (a)

上図の例では、(a)の画像が、4×4に分割されています
実際は入力の画像をこのように直接的に分割するのではなくて畳み込み層を通して、結果的に入力が小さくなります。その畳み込み層などの結果も、RGBの3チャンネルではなくて、およそ数千ものチャンネル数があります。

2) 1)で作成した各グリッドに物体があった場合、それはどのクラスか（クラス確率）を計算する (b)

(b)を見てください。各グリッドがカラフルに色付けされています。例えば、黄色が樹木、緑が植生、青が車、赤が地表面付近の物体だとします。このように各グリッドのクラスを計算します。厳密には、各クラスの数値を計算して、その中で最も確率の高いクラスをそのグリッドのクラスとします。
ここでは、物体があった場合の、条件付確率を求めています。つまり、物体があったと仮定したら、それぞれのクラスの確率はいくらであるかということです。
「そのグリッドに物体がある」とは、その物体の中心がそのグリッドに存在する場合を指します。

（引用：YOLOの論文）

If the center of an object falls into a grid cell, that grid cell is responsible for detecting that object.

3) 各グリッドにおいて、物体がありそうなバウンディングボックスの位置とその信頼度を計算

(c)を見てください。検出するための箱（バウンディングボックス）が複数出力されていることがわかります
ここで、(バウンディングボックスの情報) = (中心のx, y座標＋横幅、縦幅）で表すことができます。なぜなら、それぞれのバウンディングボックスの中心（や左上の隅）の座標を基準として、さらに横幅と縦幅がわかれば、バウンディングボックスは一つに定まるからです。
前方の樹木の幹や後方の樹木、自動車を検出したいのですが、(c)では、それらと関係のないものも複数出ており、この時点では正確な検出ができていないことがわかります。
推論時は、信頼度スコアqのみを予想。損失を計算して、重みを更新していく際は、右辺の値をもとに誤差が求まる。物体がない場合は0で、ある場合はIoUの値と等価

4) . (b), (c)より、信頼度×クラス確率を計算（[d]）

⇒各クラスごとの信頼度スコア（どれだけ推論が正確そうか）がわかる

（4’. 上の信頼度スコアやボックスの大きさが閾値以下なら消去）

5) 複数のバウンディングボックスが出力されないようにNMS を実行

NMSについて補足

1.Scoreの高いボックスから選択

2.そのボックスとIoUの高いもの (e.g. 0.5）を消去⇒繰り返し

==========論文の和訳を再開==========

上の章で、YOLOv2の大まかな流れを紹介しました。これにより、だいたいの雰囲気がつかめました。そこで、論文の和訳を再開し、このような手法い至った経緯を確認したいと思います。

Convolutional With Anchor Boxes (その2)

YOLOでは、全連結層を用いて、直接バウンディングボックスの位置を推定するのではなく、アンカーボックスを使って、各グリッドに存在する、バウンディングボックスの位置を予測する。ここではまず、プーリング層を1層削除し、ネットワークの畳み込み層の出力の解像度を高める（≒特徴マップの縦と横の大きさを大きくする）

　またYOLOに入力する画像のサイズを小さくする。従来の448×448ピクセルから416×416ピクセルに縮小した。この理由は、畳み込みなどの結果で得られる特徴マップの縦横のピクセル数を奇数にしたいためである。提案するYOLOの畳み込み層により、入力画像を32分の1にダウンサンプリングする。つまり、416×416ピクセルの入力画像をYOLOにより処理すると、13×13の出力特徴マップが得られる。

　アンカーボックスを用いた物体検出においては、各物体のクラスの予測と、各物体の位置の予測は別々に行われる。アンカーボックスごとにクラスとそのセルの物体らしさをYOLOでは予測したが、本手法もそれと同様のプロセスを踏む。物体らしさの予測の方法として、YOLOと同様に、グランドトゥルースとバウンディングボックスのIoUを予測する。また、その物体のクラスの予測においては、「そのバウンディングボックス内に物体が存在する場合における、その物体のクラスの確率」という、条件付き確率を予測する。

　しかし、アンカーボックスを用いると、精度が少し低下してしまうという課題も存在する。YOLOは1画像あたり98個のボックスを予測するが、今回のようなYOLOv2にてアンカーボックスを用いた検出を行うことで、バウンディングボックスは1000個以上生成されるする。本発表の過程で試作した、アンカーボックスを用いない場合のモデルでは、mAPが、69.5 mAP、recallは81%であった。一方、今回のようなアンカーボックスを用いた場合、mAPは69.2 mAP、recallは88%であり、mAPは若干低下していることがわかる。しかし、recallの値は高く、mAPが減少しているとはいえ、本モデルにも改良の余地があると言える。

Dimension of Clusters

　YOLOでアンカーボックスを使用する場合、2つの課題が生ずる。1つ目は、バウンディングボックスの寸法を手動にて選択する必要があるということだ（バウンディングボックスの寸法は自動では決められないということだ）。先述したように、ネットワークはバウンディングボックスの大きさ（の差分）を出力すること（や学習すること）が可能ではあるが、事前により確からしいバウンディングボックスのサイズを入力できれば、より精度の良い検出ができると考えられる。

　完全に手動で（勘や経験に近い形で）事前分布を選択するのではなく、学習セットにてラベリング情報にあるバウンディングボックスの寸法に対してk-meansクラスタリングを実行することを考える。これにより、よりよいバウンディングボックスの大きさの設定が可能となる。標準的な、ユークリッド距離を用いたk-means法を使用した場合、大きなバウンディングボックスは小さなバウンディングボックスよりも、より大きな誤差を生ずる（単に寸法が大きいため、例えば同じ10%の誤差でも絶対量としては大きくなる）。ここでの目標は、より高いIoUを得るためのバウンディングボックスの大きさを決定することである。そこで、k-meansにて用いる「距離」の定義として、以下のものを用いる。

以下の図2に示すように、様々なクラスタ数（kの値）を用いて、k-meansを実行した。そして各クラスタにて、

1 最もそのクラスタの中心に近いバウンディングボックスの寸法を求める

2 そのクラスタの全てのバウンディングボックスと1で特定したバウンディングボックスのIoUを求める

3 全てのバウンディングボックスに対する、そのIoUの値を平均する

4 そのkの値のIoUの平均（3の値）をプロットする

ということを行う。

本論文では、k = 5にてClusterの数とIoUの値のバランスがとれていると判断し、クラスタの数は5とした。その結果、図2の右側のパネル中にあるように、様々な種類のバウンディングボックスを得ることができた。

==========論文の和訳を終了==========

YOLOv2の論文の一部しか和訳ができていませんが、これまでの内容で、おおよその仕組みや流れが紹介できたと思います。最後に、テストデータにて推論を行い、検出を行ったときに、その検出の良しあしを評価する方法について述べます。

物体検出の評価方法について補足

物体検出の評価指標として、平均適合率（mean Average Precision: mAP）が有名です。

具体的な計算方法については、以下の過去の記事で紹介しています。ぜひご覧ください。

kentapt.hatenablog.com

アンカーボックスボックスの寸法を求めてみよう

後半は、YOLOv2の学習を行う上で、必要なステップについて、コードとともに確認していきたいと思います。

和訳の最後の部分で、アンカーボックスの寸法を、kmeans法により求めるとありました。論文中の図を見るだけでは、少し難しかったので、ここで例を用いて確認をしたいと思います。

なお、2022年12月現在では、estimateAnchorBoxes関数で、1行でアンカーボックスの計算を行うことができます。

jp.mathworks.com

データセットのダウンロード

この例では、295 枚の車両の画像から構成されるデータセットを使用します。Caltech の Cars 1999 データセットおよび Cars 2001 データセットからのものです (Caltech Computational Vision の Web サイトで入手可能)であり、Pietro Perona 氏によって作成されました。

以下のコードにて、ダウンロードすることができます。

clear;clc;close all
pretrainedURL = 'https://www.mathworks.com/supportfiles/vision/data/yolov2ResNet50VehicleExample_19b.mat';
websave('yolov2ResNet50VehicleExample_19b.mat',pretrainedURL);
unzip vehicleDatasetImages.zip
data = load('vehicleDatasetGroundTruth.mat');
vehicleDataset = data.vehicleDataset;

バウンディングボックスの分布を可視化

アンカーボックスの寸法を求めるときに、学習データのバウンディングボックスの情報を用います。本データの様子を眺めるため、面積とアスペクト比を用いてプロットします。

% バウンディングボックスの情報を抽出
allBoxes = cell2mat(vehicleDataset.vehicle);
% バウンディングボックスの面積とアスペクト比のプロット
% アスペクト比の計算
aspectRatio = allBoxes(:,3) ./ allBoxes(:,4);
% 面積の計算
area = prod(allBoxes(:,3:4),2);
% 図示
figure
scatter(area,aspectRatio)
xlabel("バウンディングボックスの面積")
ylabel("アスペクト比 (width/height)");
title("面積 vs. アスペクト比")

バウンディングボックスをクラスタリングし、アンカーボックスを決定

例えば、クラスタの数を6に設定し、論文と似た方法で、バウンディングボックスの寸法を決定したいと思います。

論文では、kmeans法が紹介されていました。しかし、ここでは、kmedoidsを用います。kmedoids関数では、デフォルトの設定として、kmeans++法を実行します。

まずは、例として、クラスタの数を6としてみます。また、距離は論文にあったとおり、ユークリッド距離ではなく、IoUに関する値（1-IoU）でした。距離の定義をカスタムに設定するため、'Distance'の引数を設定します。また、具体的な距離の計算の仕方は、iouDistanceMetric.mというコードに記述します。

% 6つのグループにクラスタリング
numAnchors = 6;

% K-Medoidsを使ってクラスタリング
[clusterAssignments, anchorBoxes, sumd] = kmedoids(allBoxes(:,3:4),numAnchors,'Distance',@iouDistanceMetric);

% クラスタリングした平均アンカーボックスのサイズ
disp(anchorBoxes);

% クラスタリング結果のプロット
figure
gscatter(area,aspectRatio,clusterAssignments);
title("K-Mediodsで"+numAnchors+"個クラスタリング")
xlabel("面積")
ylabel("アスペクト比 (width/height)");
grid

% 累積加算
counts = accumarray(clusterAssignments, ones(length(clusterAssignments),1),[],@(x)sum(x)-1);

% 平均IoUを計算
meanIoU = mean(1 - sumd./(counts));
disp("平均IoU : " + meanIoU);

平均IoU : 0.86311

このように、kmedoids法とバウンディングボックス情報から、クラスタ数が6のときの平均IoUを求めることができました。最後に、論文にあったように、kの値を変えながら、平均IoUとアンカーボックスの数のバランスを可視化しようと思います。

アンカーボックスの数と平均IoUの関係の確認

上のセクションのコードをfor文で回します。

アンカーボックスを増やすと平均IoUは改善するが計算量は増大します。

このデータセットでは以下のような関係を得ることができました。論文での記述と合わせると、だいたいk=5くらいが良いのではないかと思われます。

% 1から15までアンカーボックスを増やしたときに平均IoUの改善がどうなるかを計算する
maxNumAnchors = 15;
for k = 1:maxNumAnchors
    fprintf('number of anchor box: %d\n',k)
    % Estimate anchors using clustering.
    [clusterAssignments, ~, sumd] = kmedoids(allBoxes(:,3:4),k,'Distance',@iouDistanceMetric);
    
    % 平均IoUを計算
    counts = accumarray(clusterAssignments, ones(length(clusterAssignments),1),[],@(x)sum(x)-1);
    meanIoU(k) = mean(1 - sumd./(counts));
end

number of anchor box: 1
number of anchor box: 2
number of anchor box: 3
number of anchor box: 4
number of anchor box: 5
number of anchor box: 6
number of anchor box: 7
number of anchor box: 8
number of anchor box: 9
number of anchor box: 10
number of anchor box: 11
number of anchor box: 12
number of anchor box: 13
number of anchor box: 14
number of anchor box: 15

figure
plot(1:maxNumAnchors, meanIoU,'-o')
ylabel("平均IoU")
xlabel("アンカー数")
title("アンカー数 vs. 平均IoU")

まとめ

この記事では、YOLOv2の仕組みを論文の和訳とともに紹介しました。また、後半では、アンカーボックスの寸法の決定方法を、小規模なデータセットを用いて確認しました。 YOLOをはじめとする、深層学習を用いた物体検出の手法は非常に多くの分野で応用が進められています。農学分野での応用として、以下のような利用を発表させていただきました。こちらもご覧いただけますと幸いです。

www.youtube.com

2022-12-17

iPadで数式を書いて、MATLAB livescriptでブログ記事を書いてみよう

tips matlab

この記事は、MATLAB アドベントカレンダー（その2）の19日目の記事として書かれています。

qiita.com

はじめに

私が数式を含むブログ記事を投稿するときの手順を紹介したいと思います。あくまで私が気に入っている方法であり、他にもよい方法があると思います。参考になれば幸いです。

例えば、以下のような数式が多い記事をブログ投稿したい場合を考えます。

kentapt.hatenablog.com

そのような場合は、ここで紹介する方法が効果的でした。

大まかな手順として、以下の流れを想定しています。

1 内容をまとめ、ブログに掲載する数式を用意する

2 記事を執筆する

3 ブログ投稿がしやすい形で記事や画像をエキスポートする（マークダウン形式）

4 記事を投稿する

ここでは以下のものが必要です。

ペンシルで、すらすらメモがかけるタブレット端末（≒iPadとappleペンシル）
MATLAB

手順

手順1 内容をまとめ、ブログに掲載する数式を用意する

はじめに、ブログに書くような数式やその流れを考える必要があります。ここではその作業が完了していることを想定します。私は、以下のアプリを利用しています。もし、iPadをお持ちの方はぜひダウンロードしてみてください。

Nebo: メモ & PDF 注釈

MyScript
仕事効率化
無料

apps.apple.com

使い方は簡単です。以下の赤枠で示されているように、数式を入れるモードがあるので、そこで、以下のようにどんどん数式を書いていきます。

ある程度数式がかけたら、右上の、「三つの点・・・」のボタンを選択して、数式に変換してください。

すると、以下の動画のように、高精度で数式を認識して、自動変換してくれます。これよりも複雑なものでも問題ありませんでした。

このような数式の書き込みを重ねて、記事に必要な数式が書けた場合を考えます。

次のステップとしては、このiPad上の数式をブログ投稿できる形でエキスポートします。以下のボタンをタッチしてください。

出力形式はWordを選択します。

このように、簡単なタッチ操作で、iPadで書いた数式をWordにエキスポートすることができます。（もしかしたらこのあたりで少し課金が必要かもしれません）

この作業の後、Google Driveなどを経由して、Wordファイルとしてエキスポートが可能になります。

こちらの、2値化の記事の数式の例です。このように、さきほどよりも、ごちゃごちゃした式でもうまく数式に変換してくれます。

kentapt.hatenablog.com

手順2 記事を執筆する

次に、MATLABを開き、livescriptにて記事を書いていきます。下の「テキスト」を押すと、文字を打つことができます。ここで、どんどん原稿を書いていくとよいです。

1でエキスポートした数式のコピペも簡単です。Wordにエキスポートしたファイルを開いてください。

右下の▼を押すと、行形式に変換することができます。その変換した数式をコピーしてください。

以下の動画は、Wordにて、さきほどのエキスポートしたファイルを開いたときの様子です。

MATLABに数式をペーストするために、挿入タブ⇒式⇒LaTeX式を選択してください。

さきほどコピーした数式を張り付けると数式になります。このように数式を挿入しながら記事を書き進めます。

大変な執筆作業が終了したとします。手順3に移ります。

手順3 ブログ投稿がしやすい形で記事や画像をエキスポートする

私が普段利用している、はてなブログでは、マークダウン記法による記事の執筆が可能です。そのため、マークダウン形式でlive scriptをエキスポートすることができれば、ブログ投稿をスムーズに行うことができます。また、私は利用したことがないのですが、QiitaやZennといった他のサービスでも、マークダウン形式にて記事が執筆できるそうです。

live scriptをマークダウン形式でエキスポートするには、追加のパッケージをインストール必要があります。ホームタブ=>アドオン をクリックしてください。

以下のLive Script To Markdown Converterを検索し、インストールしてください。ボタンを押せば簡単にインストールできます。

また、githubからコードをダウンロードして、コードを実行することも可能です。

github.com

上記のアドオンがインストールできたら、コマンドウィンドウから、live scriptをマークダウン形式のファイルに変換します。コードは以下の1行で完了です。

すると、以下のようにマークダウン形式の原稿（.md）と、画像ファイルが別々でエキスポートされます。

手順4 記事を投稿する

以下の図は、さきほどのマークダウン形式のファイルをVisual Studio Code で開いたときの様子です。ここで用いるソフトウェアは他のものでも構いません。

azure.microsoft.com

ここでの文章を再確認しながら、はてなブログに貼り付けていきます。

プレビューモードにすると数式もうまく変換されていることがわかります。画像などを貼り付けて、その他の設定が完了すれば、投稿準備完了です。

今回の方法で嬉しいこと

この章では、今回の方法でブログ記事を書いた場合に嬉しいことを述べます。あくまで個人的な意見です。また、ほかにもよりよい方法があるかもしれません。

1 初めからマークダウンで書くよりも、簡単にかける

マークダウンで書く場合は、太字にする場合や、下線を引く場合、画像を添付する場合に、それぞれの書き方があります。例えば、太字にしたい場合は、太字とし、画像を挿入したい場合は、

![``image_0.png``](``README_images/image_0.png``)

などとします。

個人的には、ひとつひとつこのように書いていくのは少し時間がかかってしまいます。一方、MATLABのlive scriptでは、画像は、Ctrl + V 貼り付けができ、非常に簡単です。

このような色付けも、Ctrl + M で簡単にできます。

さらに、MATLABのlive scriptでは、Wordと同様に、太字は、Ctrl + B、斜体はCtrl + IのようにWordと同じ方法で文字の調整ができます。

2 コードを書く＆実行する、と同時にブログ執筆ができる

技術関連のブログ投稿は、プログラミングコードも一緒に記述する場合が多いと思います。

私が普段行っているブログ投稿のように、技術的なことを説明しつつも、時折、コードを実行しながら理解を深める場合は、live scriptが非常に便利です。そして、ここで紹介した方法では、図（figure; plot(x, y))も一緒にエキスポートすることができます。

3 Github（Git）と連携しながら執筆できる

2と関連して、コードを書きながら記事の執筆を進める場合、間違いや、コードの修正は避けられません。一度書いたコードや原稿を、修正したあとに、「やっぱりこれがよかったな」ということはありうると思います。そのようなバックアップ用などにGit連携の機能は有用です。また、他のPCで執筆を継続する場合でも、git pullをすると、最新の状態から執筆を継続できます。

MATLAB上での操作も簡単です。

git連携については以下のページがわかりやすそうでした。

ただ、私の場合は、他の言語などを使う場合に、git bashを利用しており、

MATLABの場合も、リポジトリを作成した直後の作業では、git bashを併用しています。

qiita.com

github.com

注意

iPad (nebo) =>word=>matlabでうまく数式が反映されない場合があります
iPadにて、うまく数式が認識されない場合があります
私の環境では、MATLAB 2022bでは、日本語入力がやりにくいです。左上に入力した文字が表示されてしまいます。そのため、22aを利用しています

まとめ

本記事では、iPadで数式を書いて、Wordでエキスポートし、MATLABのlivescriptでブログ記事の準備をする方法について紹介しました。記事を書くのも大変な作業ですが、便利なツールを使い分けながら進めると楽しみながらできるのではないかと思います。ブログを投稿してみたい方がいらっしゃれば、その参考になれば幸いです。

先日の講演会でも、ブログ投稿のためのMATLAB利用についても質疑応答セッションで、無理やり紹介したかったのですが、時間（と関係ないことをぶっこむ勇気）が足りませんでした。アドベントカレンダーの記事として紹介できてよかったです。

note.locusblue.com

2022-12-11

サポートベクトルマシンを理解するうえで重要なカーネル法に関して勉強してみよう

機械学習統計 matlab

この記事は、MATLAB アドベントカレンダー（その2）の11日目の記事として書かれています。

qiita.com

本記事の執筆にあたっては、

赤穂先生の「カーネル多変量解析」を非常に参考にさせていただきました。

カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学) | 赤穂昭太郎 |本 | 通販 | Amazon

非常にわかりやすい書籍であり、大変おすすめです。

なお、本記事では、機械学習の中でも非常に有名な、サポートベクトルマシン（SVM）を理解するうえで非常に重要な手法（カーネル法）を扱います。機械学習を学習している方の参考になれば幸いです。本記事は自分の勉強のまとめであり、厳密でない記述や間違いなどあるかもしれません。その場合は教えていただけますと幸いです。

この記事の原稿や、コードは以下のページに格納されています。勉強の助けになれば幸いです。

github.com

1. はじめに

これまでの記事で、線形回帰について、詳しく述べてきました。そこで用いた図の例としては以下のようなものがあります。しかし、この例では、決定係数も約0.3とあまりうまくフィッティングすることができていません。また、次数を上げて、3次式などでフィッティングしても、ある程度の残差出てしまうことが予想されます。

上の図のような、y = ax+bで表されるような、線形的な回帰ではなく、非線形なフィッティングを考えたいと思います。

ここで、本記事にて、勉強した内容をもとに、カーネル法によるフィッティングを筆者が1から実装してみました。その結果の様子をご覧ください。

サンプルの点（青）に対して、うまくフィッティングできている（赤）ことがわかります。

さきほどの線形回帰よりも、ぐにゃぐにゃとした線でフィッティングされていることがわかります。

このような複雑な形状でフィッティングするにはどのような方法を取ればよいのでしょうか。

2章や3章にて、このフィッティングをするための式について考えていきます。

2. 線形回帰（複数の変数）を行列で表現する

2章では、以前のブログ記事で行った、線形回帰の式の導出を、行列を用いて行います。また、この場合は、変数の数が多くなっても、うまく対応することができます。

線形回帰のときは y=ax+bという式を用意し、以下のように、各yの値と、予測した値の差の2乗の値が最も小さくなるような、aとbの値を求めました。

$E=\;\frac{1}{2}\sum_{i\;=\;1}^N {\left(y_i -\left({\textrm{ax}}_i +b\right)\right)}^2$

...(2-1)

そして、これらを解くことで、以下のような、aを求める式を得ました。

$a\;=\;\frac{\sum_{i\;=\;1}^N x_i y_i \;-\frac{1}{N}\sum_{i\;=\;1}^N x_i \;\sum_{i\;=\;1}^N y_i }{\sum_{i\;=\;1}^N {x_i }^2 -\frac{1}{N}{\left(\sum_{i\;=\;1}^N x_i \;\right)}^2 }$

...(2-2)

線形回帰については、こちらの記事をご覧ください。

kentapt.hatenablog.com

この場合は、変数が1つの場合でしたが、今回は、変数がもっと多くなっても対応できるようにしたいと思います。

そのために、さきほどの線形回帰の傾きaや切片bを求めるための式を行列で書いてみることにします。

xやyは以下のようになります。これは、サンプル数が4つ（データが4件）あることを意味します。xの場合、データの値は、2、6、4、2です。

$x=\left\lbrack \begin{array}{c} 2\\ 6\\ 4\\ 2 \end{array}\right\rbrack$

...(2-3)

yについても同様です。

$y=\left\lbrack \begin{array}{c} 4\\ 4\\ 5\\ 1 \end{array}\right\rbrack$

...(2-4)

傾きのaについては、どのデータに対しても同じ値が用いられるため、単純に、以下のようにあらわせます。

$a\mathrm{\ }=\ 3$

...(2-5)

次に、変数の数が多くなった場合を考えます。y = ax1+bx2+cx3+dx4+....のような形です。

以下の式(6)の縦方向がサンプルの数、横方向が変数の次元数であるとします。

...(2-6)

例えば、以下の式を見てください。先ほどと同様にデータが4つあるとします。1つ目の変数のデータの値は、2、6、4、2です。横方向は3列ありますが、変数が3つあることを表しています。

$x=\left\lbrack \begin{array}{ccc} 2 & 3 & 5\\ 6 & 1 & 7\\ 4 & 6 & 8\\ 2 & 3 & 12 \end{array}\right\rbrack$

...(2-7)

yについては、さきほどと変わりません。

$y=\left\lbrack \begin{array}{c} 4\\ 4\\ 5\\ 1 \end{array}\right\rbrack$

...(2-8)

傾きの値は、重みとして、wで表すことにします。変数が3つあるのでそれに対応する形で、3つの要素があります。

$w=\left\lbrack \begin{array}{c} 1\\ 3\\ 4 \end{array}\right\rbrack$

...(2-9)

切片に相当するものはここでは考えないとすると、(2-1)式を(2-7)や(2-9)式を用いて書くと以下のようになります。

$E\left(w\right)=\sum_{i=1}^n {\left(y^{\left(i\right)} -x^{\left(i\right)} w\right)}^2$

...(2-10)

i番目のデータに対して、それぞれの対応する重みを掛けて、それとi番目のyの値との差を取ります。さらにその値を2乗します。

そして、それをデータの数だけ繰り返し、足していくことで、予測値と実際の値の差の2乗の和、つまり偏差平方和を求めることができます。

これを行列で書くと以下のようになります。

$E(w)=(y-x{w)}^T (y-xw)$

...(2-11)

${\left(y-xw\right)}^T$ という所が少し分かりにくいと思います。以下に例を上げさせてください。

予測の誤差を2乗して、足すという計算が、行列の転置を使うと、簡単に表現することができます。

そして、式(2-11)の値が最小になる、重みwの値を求めたいので、wで微分をし、その値が0になる値を求めたいと思います。

$\frac{\partial E\left(w\right)}{\partial w}={\left(y-xw\right)}^T (y-xw)=0$

...(2-12)

ここで、転置を含むときの微分の式は以下のようになります。

$\frac{\partial }{\partial w}(Aw+c)^T \left(Bw+d\right)=2A^T Bw+A^T d+B^T c$

...(2-13)

式(2-12)と(2-13)を見比べると、

$A=B=-x$

...(2-14)

$c=d=y~$

...(2-15)

$\frac{\partial E(w)}{\partial w}=2{\left(-x\right)}^T \left(-x\right)w-x^T y-x^T y$

...(2-16)

これが0になるとき、

$2x^T xw-2x^T y=0$

...(2-17)

となって、wについて解くと、

$w={\left(x^T x\right)}^{-1} x^T y$

...(2-18)

となります。

このように、wについて、変数が複数ある場合でも、解析的に、誤差を最小にする重みwを求めることができることがわかりました。

3. カーネル法の利用

3.1. カーネル法を用いた時の誤差の最小化

このように、線形回帰を行うときの、重みwを求める方法はわかりましたが、この方法では、うまくフィッティングできるデータにも限りがあります。より柔軟なフィッティングができるように、カーネル法というものを導入したいと思います。

以下の例の、①や②では、足し算と引き算を行っています。別な言い方をすると、足し算や引き算をするときのルールのもと、計算しています。③でも同様に、カーネルaという筆者がここで定めた方法で計算して、その結果0という値を得ました。ちなみに、ここでは、カーネルaという計算のルールでは、2つのベクトルの内積を計算する、というふうにしました。

よく用いられるカーネルとして、ガウシアンカーネルというものがあります。その式を以下に示します。

$K\left(a,b\right)=\exp \left(-\beta {\left\|a-b\right\|}^2 \right)$

...(3-1)

βは任意に設定できるパラメータです。 $\|\;\|$ の計算については以下の例をご覧ください。

要素の2乗をして、その和を取り、その平方根を求めています。つまり、

$\sqrt{3^2 +1^2 }=\sqrt{10}$

...(3-2)

としています。

この式をグラフに表すと以下のようになります。横軸は、 $a-b$ の値です。

さきほどの、y=ax+bのような形式は以下のように計算できました。つまり、i番目の変数に対する重みと、i番目の変数をかけ合わせ、そしてそれらの和を計算していました。

カーネル法を用いた場合も同様です。ただ、この場合は、カーネル法の計算のKの箇所が挟まります。

$y=\sum_{i=1}^n w_i {K(x}_{i,} x)$

...(3-3)

カーネル法では、インプットしたxに対して、他の各サンプルxiとのカーネルの値に変換し、それを、重みwと掛け合わすことをしています。

例えば、サンプル数3で、変数の数（次元）が2のデータがあったとします。このデータをもとに以下のデータをテストデータとして計算する時を考えます。

...(3-4)

テストデータは1つなので、サンプル数は1、次元は同様に2です。

そして、式3-2の

${K(x}_{i,} x)$

の部分では、

(2,3)と(2,8)を使って計算をして、

(9,1)と(2,8)を使って計算をして、

(-3,1)と(2,8)を使って計算をして、

それらを最終的に以下のように計算するイメージです。

数式で表すと式3-3のようになります。

そして1章で示した、線形回帰の誤差を計算するための式を再掲します。

$E=\;\frac{1}{2}\sum_{i\;=\;1}^N {\left(y_i -\left({\textrm{ax}}_i +b\right)\right)}^2$

...(1-1)

線形回帰の場合は上のような式でしたが、

同様に、カーネル法を用いた場合でも以下のように、偏差平方和を定義できます。

$E=(y-\sum_{i=1}^n w_i {K(x}_{i,} x))^2$

...(3-5)

この式を最小化する、wを求めたいと思います。

まず、その下準備を行います。

...(3-6)

は、xのi番目と、j番目の値を用いて、カーネルの値を計算するとします。

これをもとに、行列の (i, j) 成分をその値とする行列Kを準備します。つまり、

を定義します。これをグラム行列（Gram matrix）と呼びます。

グラム行列の定義や性質については以下の記事がわかりやすかったです。

mathlandscape.com

ここでも同様に、

の結果を重み付けするために、

$w=\left\lbrack \begin{array}{c} w_1 \\ w_2 \\ \vdots \\ w_d \end{array}\right\rbrack$

...(3-8)

を用意します。

すると、

式3-5は、

$=(y-Kw)^T (y-Kw)$

...(3-9)

と表すことができます。

ここで、線形回帰の場合は、以下の式を解きました。

$\frac{\partial E\left(w\right)}{\partial w}={\left(y-xw\right)}^T (y-xw)=0$

...(3-12)

つまり、xがKになっただけで、同じ形になっていることがわかります。

そこで、2-18のxとKを入れ替えて、

この3-9の式が最小となる、wの値は、

$w={\left(K^T K\right)}^{-1} K^T y$

...(3-10)

であることがわかります。

3.2. グラム行列の特性を生かして解き進める

ここで、2-1のように、今回は以下のようなカーネル関数を用いていました。

ここでは、aとbの値を引いて、2乗の計算をするため、aとbを入れ替えても同じ値になります。

つまり、

${K(a}_, b)={K(b}_, a)$ 　

...(3-11)

です。

そのため、2-7からわかるように、グラム行列は

$K^T =K$

...(3-12)

を満たします。これを利用して、3-10を解き進めると、

$w={\left(KK\right)}^{-1} Ky=K^{-1} y$

...(3-13)

となることがわかります。

3.3. カーネル法を利用した時の誤差の和（残差平方和）について

線形回帰の場合は、誤差の二乗の和が最小となるwを選んでも、完全にそのデータに適合する直線は得られない場合がほとんどでした。そして、そのときの直線の当てはまり具合を決定係数などで評価していました。今回のようなカーネル法を利用した場合の、残差平方和はどうなるでしょうか。3-9に、 $w=K^{-1} y$ を代入します。すると、

$=(y-Kw)^T (y-Kw)$

$=0$

となります。つまり、

$E(w)=0$ 　で、カーネル法を用いてフィッティングしたときのトレーニングデータに対する

誤差はゼロとなります。

つまり、以下の図のように、青色のサンプルに対して、赤線のようにフィッティングを行うと、

与えられたサンプル全てを漏れなく通る、式を得ることができます。

3.4. 正則化項について

これまでは誤差が全くない（ゼロ）のフィッティングを行うことができることを勉強してきました。

しかし、サンプルの数（次元）が多くなるほど、そのトレーニングデータに過適合してしまいます。例えば、トレーニングデータ中にノイズの点が含まれていると、その点に引っ張られてしまいます。

つまり、トレーニングしたデータに対してはうまく行くが、それがうまくいきすぎて、逆に他の新たなデータに対しては、うまく行かないという問題が考えられます。つまり、低い汎化性が問題になります。

そこで正則化項という、このフィッティングの具合を下げるためのものを考えます。

以下のように、誤差の式に $\lambda ~w^T Kw~$ を足します。

$E(w)=(y-Kw)^T (y-Kw)+~\lambda ~w^T Kw$

...(3-14)

$\lambda ~w^T Kw~$ は2次形式であるため、微分の結果もシンプルになり、

$2\lambda ~Kw$ となります。これにより、重みの計算もやりやすくなりそうです。

manabitimes.jp

また、この微分に関しては、さきほどの、1-13式に、c=0, d=0, A=λ, B=Kを代入しても同じ結果が得られます。

$\frac{\partial }{\partial w}(Aw+c)^T \left(Bw+d\right)=2A^T Bw+A^T d+B^T c$

...(2-13)

$2\lambda ~Kw$ という項が入ると、さきほどのように、トレーニングデータに対して、完全にフィッティングすることはできません。しかし、それにより、トレーニングデータに過剰に適合することを抑えることができます。

3-14を微分して、

$\frac{\partial R\left(w\right)}{\partial w}=-2K\left(y-Kw\right)+2\lambda ~Kw$

...(3-15)

この値が0になるときは、以下の式を満たす。

$\left(y-Kw\right)-\lambda w=0$

...(3-16)

$\left(K+\lambda I_n \right)w=y$

...(3-17)

$w={\left(K+\lambda I_n \right)}^{-1} y$

...(3-18)

このような重みwを選ぶことで、トレーニングデータに過適合しないようなフィッティング結果が得られます。しかし、この正則化の程度をコントロールする、係数λを適切に選ぶ必要があります。

3.5. カーネルについてもう少し考えてみる

先ほどは、カーネルの関数として、ガウシアンカーネルを利用しました。

（これがカーネルとして利用できるかは置いておいて）、以下のような関数をカーネルとして選んだとします。

$K\left(x,~y\right)={\left(1+x^T y\right)}^2$

...(3-19)

すると、以下のように6次元のベクトルの内積として書き換えることができます。

なお、この節での説明は以下のページを参考にさせていただきました。

How to intuitively explain what a kernel is?

stats.stackexchange.com

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

www.iwanttobeacat.com

カーネルとは直感的に説明するとなんなのか？

qiita.com

次に、ガウシアンカーネルについて、同様に、n次元のベクトルの内積に書き換えられないかを考えてみます。

ガウシアンカーネルは以下のように定義されていました。

$K\left(a,b\right)=\exp \left(-\beta {\left\|a-b\right\|}^2 \right)$

ここで、

$a-b=X$

とおきます。これをテーラー展開 (Taylor expansion)すると以下のようになります。

$=1-\beta ~x+\frac{\beta^2 X^2 }{2}-\frac{\beta^3 X^3 }{6}+\frac{\beta^4 X^4 }{24}+\cdots +{\left(-1\right)}^n \frac{\beta^n X^n }{n!}$

...(3-20)

ガウシアンカーネルも、無限次元のベクトルの内積に書き換えられる気がしてきました。

ここから、ガウシアンカーネルを無次元特徴量として表すための、式展開は、以下のようになります。

なお、この式展開については、

佐久間淳先生　機械学習(6)カーネル/確率的識別モデル1の講義資料をそのまま掲載させていただいております。

非常にわかりやすい資料を公開していただき、ありがとうございました。

https://ocw.tsukuba.ac.jp/wp/wp-content/uploads/2019/10/e279a56fc6a000e67d6370309f9374c1.pdf

このように、ガウシアンカーネルを用いた時は、無限次元の特徴量を計算していることと等しいです。

次元が無限だけあるので、これを全て1からコツコツ計算していくことは、（真正面から計算する場合は）できません。

しかし、ガウシアンカーネルを用いた上の計算をすると、この無限次元の計算をした場合と等価になります。

このようにカーネルとして利用できる関数をもってくれば，無限次元の特徴ベクトルの内積を非常に簡単に計算できます。これは「カーネルトリック」と呼ばれています。

4. カーネル法を用いたフィッティングを実装してみよう

冒頭で述べた通り、4章で用いたコードは以下のページにアップロードされています。

github.com

4.1. 正則化項を使用しない場合

3章までは、カーネル法を用いたフィッティングを行うことで、線形回帰よりもはるかに柔軟性のあるフィッティングを行うことができることを示してきました。4章では、これまでの理解を確認および定着させるために、実際にそれらをコーデイングして、うまくいくか確認してみたいと思います。

言語はMATLABを用いました。3章までの数式をもとに実装しました。いったん数式を解いてしまえば、比較的短い行で書くことができました。正則化項を使用しない場合は、上で説明したとおり、サンプルデータに完全にフィットする結果を得ることができました。

clear;clc;close all
%% パラメータ設定
% ガウシアンカーネル中のベータの値
beta = 1;
% 正則化のパラメータ
L2regu = 0;
%% サンプルデータの作成と計算の準備
rng('default') % For reproducibility
x_observed = linspace(0,10,21)';
y_observed1 = x_observed.*sin(x_observed);
y_observed2 = y_observed1 + 0.5*randn(size(x_observed));
numSample = numel(x_observed);
K = zeros(numSample, numSample);
%% グラム行列の計算
% わかりやすくfor文を2つ用いて実装します
for i =1:numSample
    for j = 1:numSample
        diff = (x_observed(i)-x_observed(j))^2;
        K(i,j) = exp(-beta*diff);
    end
end
%% 重みwの計算
w = inv((K+L2regu*eye(numSample)))*y_observed2;
%% 結果の可視化
% サンプルデータのプロット
figure;scatter(x_observed,y_observed2);hold on
% フィッティングした結果のプロット
numSample_test = 100;
x_test = linspace(0,10,numSample_test)';
kOut = zeros(numSample_test,1);
for i =1:numSample_test
    x_i = x_test(i);
    diff_test = (x_observed-x_i).^2;
    kOut(i) = w'*exp(-diff_test.*beta);
end

plot(x_test,kOut)
title(sprintf('L2正則化の係数=%d',L2regu));hold off

4.2. 正則化項を変えながら実行する

4.1では、正則化項を利用しませんでした。以下の図は、正則化項のパラメータを変えながら実行した時の様子です。この図のコードもgithubにアップロードしています。正則化項の係数が小さくなる（影響を小さくする）と、真ん中下の、急に値が小さくなっている点にも適合しようと、曲線全体がより複雑になっていることがわかります。

5. まとめ

カーネル法を用いたフィッティングの方法についてまとめました
トレーニングデータに対して、誤差のない、複雑な曲線を描くことができます
自身の手でも実装し、勉強したとおりのフィッティングができることを確認しました
機械学習において、非常に有名な、サポートベクトルマシンを理解するうえで非常に重要なトピックを扱いました。サポートベクトルマシンを学習しているかたの参考になれば幸いです。

参考ページ（本文中に記載のないもの）

線形な手法とカーネル法（回帰分析）

https://qiita.com/wsuzume/items/09a59036c8944fd563ff

2022-12-06

t-SNEの勉強メモ

機械学習 matlab

この記事は、MATLAB/Simulink Advent Calendar 2022の6日目の記事として書かれています。

qiita.com

1章はじめに

t-SNEと呼ばれる方法を用いて、高次元データを、2次元平面や3次元空間にプロットすることができます。

例えば、以下の図は、MNISTという0から9の手書き数字の画像の情報を2次元平面にプロットしたときの様子です。

0から9の画像のサンプルデータが、それぞれクラスタを形成しており、うまく可視化できていることがわかります。

PCA (主成分分析)と呼ばれる方法を用いて、次元圧縮を行い、上のような2次元や3次元上でのプロットを得ることもできます。しかし、主成分分析では、高次元空間上で線形的な構造を有しているものに対して、有効であり、非線形な構造を有しているデータに対しては、うまく低次元空間にマッピングできないという欠点があります。主成分分析については、以下の私の過去の記事をご覧ください。

kentapt.hatenablog.com

本記事の執筆においては、以下のページを参考にさせていただきました。ありがとうございました。

t-SNE 解説

qiita.com

Parametric t-SNEの理論とKerasによる実装

qiita.com

t-SNE を用いた次元圧縮方法のご紹介

blog.albert2005.co.jp

この記事では、初めにt-SNEについて簡単にまとめ、その後、パラメータをいくつか変えながらt-SNEについて勉強していこうと思います。

なお、この記事の原稿やコードは以下のページにて公開されています。

github.com

2章 t-SNEについての自分用まとめ

ここからは、t-SNEの処理の流れについておおまかにまとめていきます。より正しい/詳しい解説にて勉強したい方は上の参考ページをご覧ください。

2.1. t-SNEの目標

n次元上に複数の点が存在するといます。その空間上の任意の点のペアを作り、それらを2次元平面（や3次元空間）に写像したときのことを考えます。

(i) 近い点：2次元平面（や3次元空間）に点を写像したときも、それらの点が近いままにする

(ii) 遠い点：2次元平面（や3次元空間）に点を写像したときも、それらの点が遠いままにする

(i) (ii)のように、高次元と低次元空間上で、それらの点の近さの関係をできるだけ保つようにしていきます。まずは、2.2で、この近さについてどのように定義していくかについて述べていきます。

2.2. 近さ（類似度）の定義について

2.2.1. 高い次元の場合

高い次元の空間（高次元空間）での任意の点どうしの距離を、正規分布に従う確率分布でモデル化します。

はじめに、正規分布の確率密度関数の式は以下のように表すことができます。

$f(x)=\frac{1}{\sqrt{2\pi \sigma }}\exp (-\frac{(x-\mu )^2 }{2\sigma^2 })$

...(1)

そもそも、なぜ確率密度関数の式がこのようになるのか、ということについては、以下のページがわかりやすかったです。

ai-trend.jp

ここで、2つの点、 $x_i$ と $x_j$ を考えます。先ほどの、正規分布の確率密度関数の式をもとに、以下のように近さ（類似度）を表してみます。

直感的には、分母で、全てのペアの距離を足して、それで割ることで正規化しています。確率密度関数の式のexpの直前の値は、分母分子で共通しているのでキャンセルされています。

$p_{i|j} =\frac{exp\left(-~\parallel x_i -x_j \parallel^2 /2\sigma_i^2 \right)}{\sum_{k\not= i} exp\left(-~\parallel x_i -x_k \parallel^2 /2\sigma_i^2 \right)}$

...(2)

$\sigma_i^2$ はユーザーの指定する必要があります。

上で定義した近さについては、SNEと呼ばれる、t-SNEのもととなった手法での近さを示しています。t-SNEでは、近さを以下のように定義します。

$p_{ij} =\frac{p_{i|j} +p_{j|i} }{2n}$

$n$ はデータ数を示します。

そのため、t-SNEにおいては、

$p_{ij} =p_{ji}$

であると言うことができます。ひとまず、高次元空間の任意の点の距離を正規分布の確率密度関数の式をもとに、定義しました。

2.2.2. 低い次元の場合

低次元空間上の類似度を以下のStudentのt-分布(自由度1)で表現します

$q_{ij} =\frac{1+{\left({\left\|y_i -y_j \right\|}^2 \right)}^{-1} }{\sum_{k,l\not= k} {\left(1+{\left\|y_k -y_l \right\|}^2 \right)}^{-1} }$

...(3)

ここで、t分布の一般式は以下のとおりです。

$x$ を含まない部分は、分母と分子で共通するためキャンセルされます。残った $n$ には1を代入し、 $x$ に

を代入すると上の式になることがわかります。

高次元空間上では、正規分布を使った一方で、ここでは、t分布を用いています。その理由については、Crowing problemがあるそうです。元の論文にも記述されていますが、以下のブログの解説がわかりやすかったです。

blog.albert2005.co.jp

2.3. 高次元空間上の点を低次元空間上に写像するときの考え方

高次元空間上および写像後の点の近さ（遠さ）を保てるようにするため、KLダイバージェンスを用います。

KLダイバージェンスを用いた損失関数としては、以下のようなものを用います。

$C=\sum_{k,k\not= l} p_{ij} log\frac{p_{ij} }{q_{ij} }$

....(4)

ここでは、点 i と点 jの、高次元空間と低次元空間上の距離が似ている（割り算すると1になる）ほど、logの中身が1に近づきます。それにより、式(4)の値自体も0に近づきます。

勾配法によって、低次元空間に写像する点を少しずつ変化させていき、KLダイバージェンスの値をどんどん小さくしていきます。

2.4. KLダイバージェンスの式の微分について

式(4)で示した、KLダイバージェンスの式を微分していきます。

先ほど示したような、以下の式を微分するために、logの中身を展開していきます。

$C=\sum_{k,k\not= l} p_{ij} log\frac{p_{ij} }{q_{ij} }$

$=\sum_{k,l\not= k} \left(\left.~p_{lk} logp_{lk} -p_{lk} logq_{lk} \right)\right.$

...(5)

ここで、式(3)の通り、

$q_{ij} =\frac{1+{\left({\left\|y_i -y_j \right\|}^2 \right)}^{-1} }{\sum_{k,l\not= k} {\left(1+{\left\|y_k -y_l \right\|}^2 \right)}^{-1} }$

でした。

$E_{i,j} =1+{\parallel y_i -y_j \parallel }^2$ とすると、この式は

$=\frac{{E_{i,j} }^{-1} }{\sum_{k,l\not= k} {E_{i,l} }^{-1} }$

と表すことができます。さらに分母の値を $z$ として

$=\frac{{E_{i,j} }^{-1} }{z}$

...(6)

とします。

これを用いると、式(5)は、

$C=\sum_{k,l\not= k} p_{lk} log~p_{lk} -p_{lk} log~E_{lk}^{-1} +p_{lk} log~z$

...(7)

と書き換えることができます。低次元空間上にマッピングしたときの値である $y_i$ を求めるために、

低次元空間上の $y_i$ について微分を行う。

...(8)

第一項目は、yが含まれないため、0となる。

この式(8)の2項目について考える。

$p_{ij} =p_{ji}$ 　および　 $E_{ij} =E_{ji}$ であることから、

ここで、 $E_{ij} =E_{ji}$

$\sum_{k,l\not= k} -p_{lk} \partial logE_{lk}^{-1} =-2\sum_j p_{ij} \partial log{E_{i,j} }^{-1}$

...(9)

とすることができる。

さらに、

$E_{i,j} =1+{\parallel y_i -y_j \parallel }^2$ なので、

（以下、式どうしに説明があったりして、スキャンした画像のまま掲載させてください）

（スキャンした画像終わり）

(11)および(17)で計算した、(8)の第一項と第二項を足し合わせると、

$\frac{\partial C}{\partial y_i }=4\sum_j {\left(p_{ij} -q_{ij} \right){E_{i,j} }^{-1} \left(y_i -y_j \right)}$

...(18)

を得ることができる。

$=4\sum_j \left(p_{ij} -q_{ij} \right){\left(1+{\left\|y_i -y_j \right\|}^2 \right)}^{-1} (y_i -y_j )$

...(19)

となり、 $y$ で微分すると、このようにきれいな式を得ることができました。

この式をもとに勾配法などで更新していくことで、2次元平面や3次元空間上に写像した $y$ の座標を求めることができます。

3章 t-SNEで遊んでみよう

この章では、2章で簡単にまとめた、t-SNEいろいろ試してみたいと思います。以下の公式ドキュメントを参考にします。

jp.mathworks.com

3.1. データのダウンロード

以下の、THE MNIST DATABASE of handwritten digitsにあるデータをダウンロードします。以下のページにアクセスしてください。手書き数字のデータセットをダウンロードすることができます。以下の左下の赤枠をご覧ください。

yann.lecun.com

以下のようなtrain-からはじまるファイルを2つダウンロードしてください。そして、7zipなどを用いて、ダウンロードした gz ファイルを解凍します。すると、以下のように-ubyteで終わるファイルを2つ得ることができます。

imageFileName = 'train-images.idx3-ubyte';
labelFileName = 'train-labels.idx1-ubyte';

processMNISTdata.mというコードを用いて、t-SNEによる解析を行うための前準備を行います。

[X,L] = processMNISTdata(imageFileName,labelFileName);

Read MNIST image data...
Number of images in the dataset:  60000 ...
Each image is of 28 by 28 pixels...
The image data is read to a matrix of dimensions:  60000 by  784...
End of reading image data.

Read MNIST label data...
Number of labels in the dataset:  60000 ...
The label data is read to a matrix of dimensions:  60000 by  1...
End of reading label data.

計算時間短縮のため、以下のコードでサンプルを減らします。必要なければ、===で囲まれた範囲をコメントアウトしてください

% =============
idx = randperm(numel(L),round(numel(L)*0.1));
X = X (idx,:);
L = L(idx);
% =============

画像がどのようなものか確認します。以下のように手書き数字のデータセットであることが確認できます。

dispIdx = randperm(numel(L),9);
dispCell = cell(9,1);
for i=1:9
    X_i = X(dispIdx(i),:);
    dispCell{i} = reshape(X_i,[28 28]);

end
figure;montage(dispCell)

3.2 t-SNEによる次元圧縮

さっとく、さきほどの手書き数字のデータセットをt-SNEを用いて解析し、2次元平面にデータをプロットしたいと思います。t-SNE関数を用いて簡単に実行することができます。

Barnes-Hutアルゴリズムを用いて、大規模なデータを効率的に処理します
主成分分析（PCA）を用いて、次元を784から50に削減してから、t-SNEを実行します。主成分分析については、以下の私の過去の記事をご覧ください。

kentapt.hatenablog.com

rng default % for reproducibility
% t-SNEの実行時間を計測
tim = tic;
Y = tsne(X,'Algorithm','barneshut','NumPCAComponents',50,'Perplexity',30);
toc(tim)

経過時間は 176.732598 秒です。

numGroups = length(unique(L));
% クラス数に従い、色を定義
clr = hsv(numGroups);
% 可視化
figure;gscatter(Y(:,1),Y(:,2),L,clr);title('Default Figure')

t-SNEを用いて、各数字をクラスタに分けることができました。比較的、ミスしている点も少なく、きれいに可視化できていることがわかります。

3.3. パラメータの変更

3.3.1. Perplexityについて

以下のように、Perplexityの値を変更していきます。パープレキシティが大きくなると、tsne はより多くの点を最近傍として使用します。データセットが大きい場合は、Perplexity の値を大きくします。一般的な Perplexity の値の範囲は、5 ～ 50 だそうです。

rng default % for fair comparison
Y50 = tsne(X,'Algorithm','barneshut','NumPCAComponents',50,'Perplexity',50);
figure
gscatter(Y50(:,1),Y50(:,2),L,clr)
title('Perplexity 50')

rng default % for fair comparison
Y4 = tsne(X,'Algorithm','barneshut','NumPCAComponents',50,'Perplexity',4);
figure
gscatter(Y4(:,1),Y4(:,2),L,clr)
title('Perplexity 4')

perplexityを50に設定した場合は、比較的、はじめの図と似ていると思います。一方、4にすると、各クラスターが込み合ったような図になっています。

3.3.2. Learning Rate（学習率）について

次に、学習率を変えて比較したいと思います。

以下のコマンドでは、学習率を5に設定します。

rng default % for fair comparison
YL5 = tsne(X,'Algorithm','barneshut','NumPCAComponents',50,'LearnRate',5);
figure
gscatter(YL5(:,1),YL5(:,2),L,clr)
title('Learning Rate 5')

次に学習率を2000に設定します。

rng default % for fair comparison
YL2000 = tsne(X,'Algorithm','barneshut','NumPCAComponents',50,'LearnRate',2000);
figure
gscatter(YL2000(:,1),YL2000(:,2),L,clr)
title('Learning Rate 2000')

以上のように、学習率を変えると、また異なる結果を得ることができました。データによるため一概には言えませんが、学習率が小さすぎると、最適な解にたどり着けず、十分分離されていない結果が得られると思います。一方、大きすぎても、最適な場所に解が落ち着かないため、きれいな分離結果を得ることができません。

4章簡単なまとめと注意点

簡単なまとめ

この記事では、t-SNEの仕組みを簡単に説明し、さらにMATLABを用いて、その動作を確認しました。Perplexityや学習率などのパラメータがあり、それらの設定によって結果が異なることがわかりました。

その結果を使用して新しいデータを分類することはできない

多くの場合に t-SNE はデータをクラスターにきれいに分離するので、新しい手元のデータにも試したくなります。しかし、t-SNE は新しいデータ点をクラスターに分類できません。t-SNEによるマッピングは、データに依存する非線形写像で、新たに低次元空間に新しい点を埋め込むためには再度その点を含めてt-SNEを実行しなおす必要があります。一方主成分分析（PCA）による方法では、マッピングするための重みなどが求まるため、新たなデータに対しても適用することができます。

参考ページ（本文中に記載のないもの）

1 Yann LeCun (Courant Institute, NYU) and Corinna Cortes (Google Labs, New York) hold the copyright of MNIST dataset, which is a derivative work from original NIST datasets. MNIST dataset is made available under the terms of the Creative Commons Attribution-Share Alike 3.0 license

creativecommons.org

2 MATLAB ドキュメンテーション：t-SNEとは

jp.mathworks.com

3 Federico Errica: Step-By-Step Derivation of SNE and t-SNE gradients

https://diningphil.github.io//media/sne_tsne.pdf