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]
でいいと思うしソッチのほうが自然な定義な気がしてならないんですが、なぜそうなっていないのか……
勾配ブースティングを実装してみた
勾配ブースティングとは
機械学習のアルゴリズムはいろいろなものがありますが、その中でも木構造を使ったアンサンブル学習の一つとして勾配ブースティング法と呼ばれるものがあります。
勾配ブースティング法は学習器を一つづつ追加していく加法的モデルと呼ばれるモデルに分類されます。ざっくりとしたアルゴリズムとしては
- 今ある学習器の結果と目的の値との差を縮めるような学習器を作り
- それを今ある学習器に追加する
ということを何回も繰り返して、段々と目的の値を得られるような学習器を生成します。
勾配ブースティング法は、ブースティングの学習の時に目的関数の勾配およびヘシアン情報を使うもので
という特徴から人気者です。実際自分で使っていても大概の問題でうまく行ってる印象があります。
高速な実装としては、xgboostなどが有名です。(内部でC++で実装されているようです
普通に使う文にはxgboost使いましょうという話なんですが思い立って実装してみました。
使い方
基本的には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()
いい感じですね。適当に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
PRML 第3章のEvidence近似
PRML第3章では、線形回帰モデルを扱っています。
その後半で、これまで固定だとしていた、重みwの分布の精度パラメータと、実際に観測される値tのモデルであるに現れるノイズの精度を表すの値もいい感じの数値にしちゃいましょうという話が出てきます。
これを真面目にやろうとすると、α,βに関する事前分布を定義してやって、事後分布である
を最大化するを探すことになります。
じゃあαとβの分布はどうするの?ということになります。
ここで、仮定を追加して、われわれは、αやβについて何も情報を持っていない、とします。
これはどういうことかというと、どのαやβに対しても、だいたい同じぐらいの確率があるよね、という風に考えるということです。
なので、(α,β)=(0.2,0.1)だと思っているモデルM1も、(3.0,5.0)だと思っているモデルM2も、同じぐらいの確率で選ばれるよね、とするわけです。
こうなると、p(α,β)はどのαやβに対しても同じ値となるので、先ほどの最大化すべき関数が
となり、これは単にwを積分消去すれば計算できて、
となり、特にがガウス分布に従うと仮定しているので、これは解析的に解くことが出来て、とても嬉しいねという流れになります。これをEvidence近似とか経験ベイズとか第二種最尤推定とか言ったりします。第二種最尤推定と呼ばれるのは、普通の最尤推定がwの分布が存在していることまで仮定しているのに対して、今度はαやβも固定されているのではなくて、(平らという条件はついていますが)確率的な分布があると考えているからです。(これをもっと発展させて、αやβも分布をもたせることも可能で、それはPRML下巻の変分推論のあたりに書いてあります。)
さてでは今度は実際に解析的に解けたものを最大化していこうと思いますが、実はこれにも二つのアプローチがあります。一つは、今回とる解析的に解いて微分をする方法。もう一つがEMアルゴリズムを用いる方法です。EMアルゴリズムは、最大化したい関数が隠れ変数の積分で表されている時に使えるので、今回だとwは積分によって消えているので、EMアルゴリズムを適用することが出来ます。
でも今回は、せっかく解析的に解くことができるので、前者のアプローチを取ることにします。(計算すると、EMアルゴリズムと解析的に解く方法は一致することも確認できます。PRML下巻を参照)
どうやって解くのかを式変形を書いても良いのですが、それは何処にでも載っているので、ここではその気持について書きたいと思います。(自分が後で見て復習(あの時の考え方ってどうだったんだろうとか)に使える為)
まず大事なのが有効パラメータ数(well-determined parameter)の事です。有効パラメータとはカジュアルに言えば特徴量ベクトルによって計算できる行列 の $n$ 個の固有値 のうち事前分布 $\alpha$ と比べて大きい値を持つものの数といえます。
この行列の固有値が小さい固有ベクトルというのは特徴量の空間において、あまり動きのない方向であると言えます。この固有値はデータ数が増えると単調に増加していくということを考慮すれば、まだデータが十分に無いためにその方向の情報量がない(別の固有ベクトルで十分に説明がつく) という状況を表しています。すなわち各データ がそちら方向の変動に乏しい、ということです。
一方で固有値が大きいというのは上記の逆で、その方向への変動が大きかったということを表します。これらを加味すると、行列の固有値のうち比較的数値が大きい物の数、すなわちデータを表すのに必要な次元の数を計算して、その数の次元程度でデータを予測するのが良さそうです。
これを表すのが有効パラメータ数になります。
これはwの事前分布の精度αと、先ほどのデータの固有値(重要度)λとの比になっていて、αに比べて大きければ各∑の中身は1になり、小さければ0になります。なので、すべての和を取ると、今あるデータ点は、αとくらべて重要度が高いものが何個ぐらいあるか、という指標になります。
・普通にMAP推定をやった時の予測平均値
・Evidence近似によって尤度を最大化した時に得られる予測平均値
・その時の有効パラメータγの変化
をプロットするやーつを作ってみました。
まずはガウス基底から。基底関数には、ガウス基底を50次元用意しています。
訓練データは、y=sin(x/2)+ノイズ(精度パラメータ5.)で生成していて、グラフ上には青い点として表示しました。青の実線はノイズの無い目的の値y=sin(x/2)です。
・データ点が多くなると、有効パラメータγも大きくなること
・MAP推定だけだと、選んだαとβに依存するので、今回初期の値としてα=0.1としているために、データ点にめちゃくちゃ引っ張られていること
・データ点が多くなると、両者は一致してくること
とかがわかると思います。
データ点が多くなればML推定に全部近づいていくので問題ないのですが、次元に比べてデータ点が少ない時(今回は基底に50個のガウス基底を用いていますから、10,20,30の時とかが顕著だと思います)に、Evidence近似では、データにひっぱられることなくうまく平滑化ができています。
つぎは多項式基底。10次式まで用意しました。こちらはy=sin(x)がターゲットです。
データ点が少ない時は繰り返し回数が足りてませんね。雰囲気サンプル数が20以下の時は、γは最終的に0になりそうです。
(データが足りないからこれからは何も言えないので全部零にしちゃえということでしょうか)