読者です 読者をやめる 読者になる 読者になる

nykergoto’s blog

機械学習とpythonをメインに、雑多な内容をとりとめなく扱うブログです。

pythonによる固有値固有ベクトル計算

python固有値固有ベクトルを計算しようと思うと numpy.linalg (線形代数の便利関数のモジュール)以下にある numpy.linalg.eigs を使うことになります。

numpy.linalg.eig — NumPy v1.11 Manual

この挙動が個人的な感覚とちょっとずれていて、そこでかなり消耗してしまったのでメモとして残しておきます。

nump.linalg.eigs でハマる

numpy にはnumpy.linalg.eigsという関数が用意されています。
返り値として第一変数に固有値(eigenvalues)、第二変数に正規化された(すなわち \( \|u_i\|=1 \) )固有ベクトル(eigenvector)を返します。

X = # 適当な(n,n)のnp.array
eig_vals,eig_vec = np.linalg.eigs(X)

ここで1点注意点があります。

http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eig.html
Returns:
w : (..., M) array The eigenvalues, each repeated according to its multiplicity. (以下略
v : (..., M, M) array The normalized (unit “length”) eigenvectors, such that the column v[:,i] is the eigenvector corresponding to the eigenvalue w[i].

重要なのは、index を i とした時に固有値 eig_vals[i] に対応する固有ベクトルeig_vec[i]ではないということ!!

正しくはeig_vec[:,i]が対応します。

この問題のやっかいな点がデバックすることが困難である点です。

なぜなら、インデックスを入れる次元を間違えて eig_val[i] としても、その配列の長さは正しい場合と同じ M になり、ベクトルの長さでは判定できません。また「固有ベクトル同士は直行する」という条件を使ってデバッグをしようとしても、間違えたベクトルに対しても成り立ってしまいます。すなわちi=jの時だけ

eig_vec[i].dot(eif_vec[j])=1

が成り立ちます。ですから間違えたindexを与えた際にも固有ベクトルの特徴である「お互いが直交する」という性質も満たしてしまっています。

したがって、計算した結果がなんだかおかしい、ということ以外でのデバッグが不可能になります。僕はこれで3時間ぐらい持って行かれました。おとなしく最初にドキュメントを読めばいいのですが。ちょっと慣れてくると適当にやってもそれっぽく出来てしまうのはダメだなと改めて感じました。ドキュメント大事。

までも普通に考えて、i番目の固有ベクトルだからeig_vec[i]でいいと思うしソッチのほうが自然な定義な気がしてならないんですが、なぜそうなっていないのか……

勾配ブースティングを実装してみた

勾配ブースティングとは

機械学習アルゴリズムはいろいろなものがありますが、その中でも木構造を使ったアンサンブル学習の一つとして勾配ブースティング法と呼ばれるものがあります。

勾配ブースティング法は学習器を一つづつ追加していく加法的モデルと呼ばれるモデルに分類されます。ざっくりとしたアルゴリズムとしては

  1. 今ある学習器の結果と目的の値との差を縮めるような学習器を作り
  2. それを今ある学習器に追加する

ということを何回も繰り返して、段々と目的の値を得られるような学習器を生成します。

勾配ブースティング法は、ブースティングの学習の時に目的関数の勾配およびヘシアン情報を使うもので

  • 学習ステップを増やしても過学習を起こしにくい
  • 非常に柔軟な予測モデルを生成できる
  • 木構造ゆえ変数の正規化が必要ない

という特徴から人気者です。実際自分で使っていても大概の問題でうまく行ってる印象があります。

高速な実装としては、xgboostなどが有名です。(内部でC++で実装されているようです

普通に使う文にはxgboost使いましょうという話なんですが思い立って実装してみました。

github.com

使い方

基本的にはreadme.mdに書いてあるのでそちらを見てほしいのですが、

の二値分類であればgbdtreeフォルダを実行ファイルと同じ階層に置いて

import gbdtree as gb

clf = gb.GradientBoostedDT()
x_train,t_train = ~~ # 適当なトレーニングデータ
clf.fit(x=x_train,t=t_train)

pred = clf.predict(x_test) # x_test の予測を出力

というような感じで学習ができます。

やってみた

連続変数に対する回帰問題をやってみます

import pandas as pd
import numpy as np
import gbdtree as gb
import matplotlib.pyplot as plt

sample_size = 100
x_train = np.linspace(-np.pi,np.pi,sample_size)
t_train = np.sin(x_train) + np.random.normal(scale=.5,size=sample_size)

# xの形は(N,feature_dim)の形に直す
x_train = x_train.reshape(sample_size,1)

regobj = gb.LeastSquare()
loss = gb.leastsquare

clf = gb.GradientBoostedDT(regobj,loss)
clf.fit(x=x_train,t=t_train,num_iter=20)

pred = clf.predict(x_train)
plt.plot(x_train,t_train,'o')
plt.plot(x_train,np.sin(x_train))
plt.plot(x_train,pred)
plt.show()

f:id:dette:20160624075124p:plain

いい感じですね。適当にclf.fitに与えるパラメータを変えて挙動がどういうふうに変わるかを見るのも面白いかもしれません。

二値分類に関しては先程のgithubの中に2つ

  • sample.pyに適当なガウス分布に対する分類問題
  • mnist.pyにMNISTの手書き文字(3と8)の二値分類問題

のサンプルが入っているのでよければどうぞ。

gitのコミットメッセージの書き方

最近ようやくgitをまともに使うようになってきました。

コミットメッセージの書き方がよくわかってなかったのでまとめます。

 

基本的なコミットメッセージのスタイル

  • 英語で記述すること!!
  • ピリオドなし
  • 現在時制
  • 文頭の単語は大文字から開始

基本フォーマット

  • 一行目:変更内容の要約(タイトルと概要)
  • 二行目:空白
  • 三行目:変更の理由、詳細

大事な一行目

三行目は適当に英作文すれば良いのでやはり大事なのは一行目。コミットの種類と要約が人目でわかるように記述します。

フォーマット

[コミット種別]+"要約"

ex). [fix] Fix some Bugs in hogehog method

コミット種別

  • fix:バグ修正
  • hotfix:クリティカルなバグ修正
  • add:新規(ファイル)機能追加
  • update:機能修正(バグではない)
  • change:仕様変更
  • clean:整理(リファクタリング等)
  • disable:無効化(コメントアウト等)
  • remove:削除(ファイル)
  • upgrade:バージョンアップ
  • revert:変更取り消し

普通に書いていると hotfix とか revert とか書かないので覚えとこう

 参考url

git commit時のコメントを英語で書くための最初の一歩 | hiro345

PRML 第3章のEvidence近似

PRML第3章では、線形回帰モデルを扱っています。
その後半で、これまで固定だとしていた、重みwの分布{N(w|0,\alpha^{-1} I)}の精度パラメータ{\alpha}と、実際に観測される値tのモデルである{N(t|w^{T}\phi(x),\beta^{-1})}に現れるノイズの精度を表す{\beta}の値もいい感じの数値にしちゃいましょうという話が出てきます。

これを真面目にやろうとすると、α,βに関する事前分布を定義してやって、事後分布である
{p(\alpha,\beta|t) \propto p(t|\alpha,\beta)p(\alpha,\beta)}
を最大化する{\alpha,\beta}を探すことになります。

じゃあαとβの分布はどうするの?ということになります。
ここで、仮定を追加して、われわれは、αやβについて何も情報を持っていない、とします。

これはどういうことかというと、どのαやβに対しても、だいたい同じぐらいの確率があるよね、という風に考えるということです。
なので、(α,β)=(0.2,0.1)だと思っているモデルM1も、(3.0,5.0)だと思っているモデルM2も、同じぐらいの確率で選ばれるよね、とするわけです。

こうなると、p(α,β)はどのαやβに対しても同じ値となるので、先ほどの最大化すべき関数が
{p(t|\alpha,\beta)}
となり、これは単にwを積分消去すれば計算できて、
{p(t|\alpha,\beta)=\int{p(t|w,\beta)p(w|\alpha)dw}}
となり、特に{p(t|w,\beta),p(w|\alpha)}ガウス分布に従うと仮定しているので、これは解析的に解くことが出来て、とても嬉しいねという流れになります。これをEvidence近似とか経験ベイズとか第二種最尤推定とか言ったりします。第二種最尤推定と呼ばれるのは、普通の最尤推定がwの分布が存在していることまで仮定しているのに対して、今度はαやβも固定されているのではなくて、(平らという条件はついていますが)確率的な分布があると考えているからです。(これをもっと発展させて、αやβも分布をもたせることも可能で、それはPRML下巻の変分推論のあたりに書いてあります。)

さてでは今度は実際に解析的に解けたものを最大化していこうと思いますが、実はこれにも二つのアプローチがあります。一つは、今回とる解析的に解いて微分をする方法。もう一つがEMアルゴリズムを用いる方法です。EMアルゴリズムは、最大化したい関数が隠れ変数の積分で表されている時に使えるので、今回だとwは積分によって消えているので、EMアルゴリズムを適用することが出来ます。
でも今回は、せっかく解析的に解くことができるので、前者のアプローチを取ることにします。(計算すると、EMアルゴリズムと解析的に解く方法は一致することも確認できます。PRML下巻を参照)

どうやって解くのかを式変形を書いても良いのですが、それは何処にでも載っているので、ここではその気持について書きたいと思います。(自分が後で見て復習(あの時の考え方ってどうだったんだろうとか)に使える為)
まず大事なのが有効パラメータ数(well-determined parameter)の事です。有効パラメータとは、ざっくり言ってしまうと、特徴量ベクトル{phi(x)}によって計算できる{\beta*\sum_i{\phi(x)\phi(x)^T}}固有値{\lambda_i}のわりと大きい次元です。

この固有値が小さい方向というのは、実際のxから射影された先の特徴量空間において、あまり動きのない方向であると言えます。
だから、その方向のデータはあまり意味を成さないというか、{phi(x)}の動きがそちら方向に乏しいので、あてにできないのです。なので、大きな固有値を持つ固有ベクトルの方向に大きく動いて、それ以外はあまり動かない(確証がないから)という作戦を取りたくなります。

その度合を表すのが、有効パラメータ数{\gamma = \sum_i{\frac{\lambda_i}{\alpha+\lambda_i}}}になります。
これはwの事前分布の精度αと、先ほどのデータの固有値(重要度)λとの比になっていて、αに比べて大きければ各∑の中身は1になり、小さければ0になります。なので、すべての和を取ると、今あるデータ点は、αとくらべて重要度が高いものが何個ぐらいあるか、という指標になります。

今回は、ガウス基底と多項式基底を用いた時に、

・普通にMAP推定をやった時の予測平均値
・Evidence近似によって尤度を最大化した時に得られる予測平均値
・その時の有効パラメータγの変化

をプロットするやーつを作ってみました。

まずはガウス基底から。基底関数には、ガウス基底を50次元用意しています。
訓練データは、y=sin(x/2)+ノイズ(精度パラメータ5.)で生成していて、グラフ上には青い点として表示しました。青の実線はノイズの無い目的の値y=sin(x/2)です。
f:id:dette:20150922144805p:plain

・データ点が多くなると、有効パラメータγも大きくなること
・MAP推定だけだと、選んだαとβに依存するので、今回初期の値としてα=0.1としているために、データ点にめちゃくちゃ引っ張られていること
・データ点が多くなると、両者は一致してくること
とかがわかると思います。

データ点が多くなればML推定に全部近づいていくので問題ないのですが、次元に比べてデータ点が少ない時(今回は基底に50個のガウス基底を用いていますから、10,20,30の時とかが顕著だと思います)に、Evidence近似では、データにひっぱられることなくうまく平滑化ができています。

つぎは多項式基底。10次式まで用意しました。こちらはy=sin(x)がターゲットです。
f:id:dette:20150922150110p:plain

データ点が少ない時は繰り返し回数が足りてませんね。雰囲気サンプル数が20以下の時は、γは最終的に0になりそうです。
(データが足りないからこれからは何も言えないので全部零にしちゃえということでしょうか)