洋食の日記

「だ・である」調ではなく「です・ます」調で書きはじめれば良かったなと後悔してる人のブログです

SVMKitにPipeline機能を追加した

はじめに

SVMKitで、一通りベーシックな機械学習アルゴリズムの実装を終えたので、しばらく便利機能の追加を予定している。バージョン0.7.2ではPipelineを実装した。Pipelineを使うことで、正規化して主成分分析してSVMで分類といった連結処理を定義できる。

svmkit | RubyGems.org | your community gem host

使い方

gemコマンドでSVMKitをインストールする。

$ gem install svmkit

例で使うデータセットLIBSVM Dataからpendigitsデータセットを取得する。

$ wget https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multiclass/pendigits

RBFカーネル近似を行い、線形SVM分類器で分類することをPipelineを使って実装すると次のようになる。

require 'svmkit'

# pendigitsデータを読み込む。
samples, labels = SVMKit::Dataset.load_libsvm_file('pendigits')
samples = Numo::DFloat.cast(samples)

# カーネル近似と線形SVMによる分類器を構成する。
# Pipelineの各ステップはHashで定義する。
rbf = SVMKit::KernelApproximation::RBF.new(gamma: 0.0001, n_components: 800, random_seed: 1)
svc = SVMKit::LinearModel::SVC.new(reg_param: 0.0001, max_iter: 1000, random_seed: 1)
pipeline = SVMKit::Pipeline::Pipeline.new(steps: { hoge: rbf, fuga: svc })

# 5-交差検定を実施する。
kf = SVMKit::ModelSelection::StratifiedKFold.new(n_splits: 5, shuffle: true, random_seed: 1)
cv = SVMKit::ModelSelection::CrossValidation.new(estimator: pipeline, splitter: kf)
report = cv.perform(samples, labels)

# 結果を出力する。
mean_accuracy = report[:test_score].inject(:+) / kf.n_splits
puts("5-CV mean accuracy: %.1f %%" % (mean_accuracy * 100.0))

これを実行すると次のようになる。

$ ruby pipeline.rb
5-CV mean accuracy: 99.2 %

Pipelineは、他の機械学習アルゴリズム同様にMarshal.dumpとloadで、モデルの書き出しと読み込みが行える。

おわりに

SVMKitも多機能になり、名前と内容が一致していないため、改名を考えている。0.7系でユーティリティ的なのを実装したら変えるつもりです。

Red DatasetsとSVMKitを使ってIrisデータセットでの線形SVMの分類精度を確認する

はじめに

Red Datasetsは、IrisやMNISTといった公開されているデータセットを、Rubyで簡単に扱えるようにするプロジェクトである(Pythonでいえば、scikit-learnのsklearn.datasetsや、Kerasのkeras.datasetsに近い)。本記事では、Red DatasetsでIrisデータセットを読み込み、SVMKitで線形SVMによる分類精度の交差検定を行う。SVMKitは、データをNumo::NArrayで扱うので、そこの変換が必要になる。

インストール

Red DatasetsとSVMkitともに、gemで簡単にインストールできる。

$ gem install red-datasets svmkit

Red Datasetsの簡単な使い方

使い方は、Red DatasetsのUsageがわかりやすい。 Red Datasetsは、データセットを、関係データベースの様な複数レコードからなるテーブルで表現する。

require 'datasets'

# Irisデータセットをnewする。
iris = Datasets::Iris.new

# Irisデータセットは、花のアヤメの種類を表すラベルと、ガク片と花びらの長さ・幅による特徴量からなる。
# eachでそれらを1つずつ見ていく。
iris.each do |r| 
  puts "#{r.label}, #{r.sepal_length}, #{r.sepal_width}, #{r.petal_length}, #{r.petal_width}"
end

これを実行すると以下のようになる。便利!!

$ ruby red_test.rb
Iris-setosa, 5.1, 3.5, 1.4, 0.2
Iris-setosa, 4.9, 3.0, 1.4, 0.2
Iris-setosa, 4.7, 3.2, 1.3, 0.2

... (省略)

コード

Red Datasetsで読み込んだIrisデータセットで、線形SVMの分類精度の交差検定を行う。テーブルによるデータの取り出しを試してみた。

require 'datasets'
require 'svmkit'
require 'numo/narray'

# Irisデータセットを読み込む。
iris = Datasets::Iris.new

# テーブルを取得する。
iris_table = iris.to_table

# ラベルと特徴量に分けてとりだす。
iris_labels = iris_table[:label]
iris_attrs = iris_table.fetch_values(
  :sepal_length, :sepal_width, :petal_length, :petal_width).transpose

# Irisデータセットの文字列によるラベルを整数値のラベル (Numo::Int32) に変換する。
encoder = SVMKit::Preprocessing::LabelEncoder.new
labels = encoder.fit_transform(iris_labels)

# Irisデータセットの特徴量をNumo::DFloatに変換する。
samples = Numo::DFloat[*iris_attrs]

# 線形SVMの5-fold分割による交差検定を定義する。
svc = SVMKit::LinearModel::SVC.new(
  reg_param: 0.0001, fit_bias: true, max_iter: 3000, random_seed: 1)
kf = SVMKit::ModelSelection::StratifiedKFold.new(n_splits: 5, random_seed: 1)
cv = SVMKit::ModelSelection::CrossValidation.new(estimator: svc, splitter: kf)

# 交差検定を実行する。
report = cv.perform(samples, labels)

# 平均Accuracyを出力する。
mean_accuracy = report[:test_score].inject(:+) / kf.n_splits
puts("Mean Accuracy: %.1f%%" % (100.0 * mean_accuracy))

これを実行すると以下のとおり。Irisデータセットにおける線形SVMの分類精度を、5交差検定で確認すると94.7%であるとわかる。

Mean Accuracy: 94.7%

おわりに

公開されているデータセットは、ファイル形式が独自のバイナリだったりして、まず扱えるようにするまでが大変だったりするものも多い。Red Datasetsのように、データセットを統一された使い方で扱えるのはとてもありがたい。何かデータセットを追加して貢献できたらな〜。

Numo::Linalgで特異値の数を指定できる特異値分解を実装する

はじめに

Numo::Linalgは、Rubyで行列のノルムの計算や固有値分解などの線形代数計算を行うライブラリである。Pythonのnumpy.linalgやscipy.linalgに相当する。Numo::Linalgのバージョン0.1.3から、固有値分解で求める固有値固有ベクトルの範囲を指定できる(主成分分析など機械学習アルゴリズムでよくある、固有値の小さい/大きい方からk個...ができる)。

一方で特異値分解は範囲を指定できない。これは、numpy.linalgやscipy.linalgも同様であり、ライブラリが叩いているLAPACKのgesvd関数が、特異値の範囲を指定できないことに起因する(ちなみにscikit-learnのTruncatedSVDはARPACKを利用した特異値の個数指定に対応している)。

特異値分解は、固有値分解で書き下すことができる。そこで、Numo::Linalgの固有値分解を行うeighメソッドを利用して、特異値の大きい方からk個の特異値・特異ベクトルを求める、打ち切り特異値分解(Truncated Singular Value Decomposition)を実装する。

Numo::Linalgのインストールと固有値分解

以降は、macOS Mojave 10.14.2上で行った。まず、バックグラウンドライブラリのBLAS/LPACKのためにOpenBLASをbrewでインストールする。

$ brew install openblas --with-openmp

そして、Numo::Linalgをインストールする。依存関係でベースとなるベクトル・行列ライブリのNumo::NArrayも一緒にインストールされる。これで準備完了となる。

$ gem install numo-linalg

例として、Numo::Linalgのeighメソッドを用いて、3×3の大きさの対称行列の小さい方から2つの固有値と対応する固有ベクトルを求めた。irb上で行った。

> require 'numo/linalg/autoloader'
=> true

> a = Numo::DFloat.new(3, 4).rand
=> Numo::DFloat#shape=[3,4]
[[0.0617545, 0.373067, 0.794815, 0.201042],
 [0.116041, 0.344032, 0.539948, 0.737815],
 [0.165089, 0.0508827, 0.108065, 0.0687079]]

> b = a.dot(a.transpose)
=> Numo::DFloat#shape=[3,3]
[[0.815142, 0.713004, 0.128882],
 [0.713004, 0.967739, 0.145705],
 [0.128882, 0.145705, 0.046242]]

> evals, evecs = Numo::Linalg.eigh(b, vals_range: 0...2)
=> [Numo::DFloat#shape=[2]
[0.022371, 0.174381], Numo::DFloat#shape=[3,2]
[[-0.0739081, 0.744638],
 [-0.0972296, -0.667395],
 [0.992514, -0.00993004]]]

固有値分解と特異値分解の関係

固有値分解は、正方行列Aを、固有値による対角行列Dと固有ベクトルからなる行列Eに分解する。また、固有ベクトルは直交するので、行列Eは直交行列となる。

 A=EDE^\top, \hspace{1em} E^\top E=I

特異値分解は、大きさm×nの行列Bを、m×mの大きさの特異値による対角行列Sと、特異ベクトルからなるm×mの大きさの行列U, n×mの大きさの行列Vに分解する。また、特異ベクトルも直交する。

 B=USV^\top, \hspace{1em} U^\top U=I, V^\top V=I

ここで、行列Bを自身の転置行列との内積による対象行列Cを考えると、次のようになる。

 C=BB^\top =USV^\top VSU^\top=US^{2}U^\top

つまり、行列Cの固有値平方根は行列Bの特異値に、行列Cの固有ベクトルは行列Bの右特異ベクトルに対応する。さらに、特異値と右特異ベクトルがわかれば、左特異ベクトルを以下のように、求めることができる。

 S^{-1}U^\top B=V^\top

Numo::Linalgのeighメソッドは、対称行列の固有値固有ベクトルを範囲を指定して求めることができる。以上の関係を利用することで、打ち切り特異値分解を実装することができる。

実装

上記の式を実装していく。欲しい特異値の数が、行列の大きさを超えた場合のエラー処理などは省いている。

# a: 入力行列
# k: 大きい方から得る特異値・特異ベクトルの数
def svd(a, k)
  # 入力行列の大きさを得る。
  n_rows, = a.shape

  # 対称行列を計算する。
  b = a.dot(a.transpose)

  # 対称行列を固有値の範囲を指定して固有値分解する。
  # Numo::Linalg.eighメソッドでは固有値は昇順に並ぶ。
  # 大きい方から得られるように値の範囲の指定に注意する。
  vals_range = (n_rows - k)...n_rows
  evals, evecs =  Numo::Linalg.eigh(b, vals_range: vals_range)

  # 固有値・固有ベクトルから特異値・左右の特異ベクトルを求める。
  # reverseメソッドで降順にしている。
  s = Numo::NMath.sqrt(evals.reverse.dup)
  u = evecs.reverse(1).dup
  vt = (1.0 / s).diag.dot(u.transpose).dot(a)
  
  [s, u, vt]
end

Numo::Linalgのsvdメソッドと結果を比較することで、簡単に動作を確認する。

require 'numo/linalg/autoloader'

# ここに上記のsvdメソッドがあるとする

a = Numo::DFloat.new(3, 5).rand

pp Numo::Linalg.svd(a, job: 'S')

puts '--'

pp svd(a, 2)

これを実行すると以下のようになる。大きい方から2個分の特異値と対応する特異ベクトルを求めることができている。

$ ruby svd_test.rb
[Numo::DFloat#shape=[3]
[1.64956, 0.495033, 0.164851],
 Numo::DFloat#shape=[3,3]
[[-0.543985, -0.139848, -0.827359],
 [-0.564919, -0.668036, 0.48435],
 [-0.620441, 0.73087, 0.284398]],
 Numo::DFloat#shape=[3,5]
[[-0.17883, -0.333785, -0.85485, -0.302866, -0.184692],
 [-0.322162, -0.732599, 0.114647, 0.427095, 0.404915],
 [0.8873, -0.167401, -0.261489, 0.301801, 0.158795]]]
--
[Numo::DFloat#shape=[2]
[1.64956, 0.495033],
 Numo::DFloat#shape=[3,2]
[[-0.543985, -0.139848],
 [-0.564919, -0.668036],
 [-0.620441, 0.73087]],
 Numo::DFloat#shape=[2,5]
[[-0.17883, -0.333785, -0.85485, -0.302866, -0.184692],
 [-0.322162, -0.732599, 0.114647, 0.427095, 0.404915]]]

パフォーマンス比較実験

打ち切り特異値分解が効果を発揮するのは、入力行列が大きな場合である。例えば、1,000次元の特徴ベクトルが、5,000サンプルあり、その50次元の主成分を求めることを考える。

require 'numo/linalg/autoloader'
require 'benchmark'

# ここに上記のsvdメソッドがあるとする

# 1,000次元の特徴ベクトルが、5,000サンプルあるイメージの行列
# ※本当に、特異値分解を用いて主成分分析を行うなら、事前に平均を減じて中心化する。
a = Numo::DFloat.new(1000, 5000).rand

Benchmark.bm 10 do |r|
  # Numo::Linalg.svdによる特異値分解
  r.report 'Numo::Linalg' do
    s, u, vt = Numo::Linalg.svd(a, job: 'S')
  end

  # svdメソッドによる特異値分解
  r.report 'svd' do
    s, u, vt = svd(a, 50)
  end
end

これを実行すると以下のようになる。おおよそ10倍ほど速いことがわかる。Numo::Linalg.svdメソッドでは、全ての特異値・特異ベクトルを求めるため、それと比べて速いのは当然のようだが...欲しい特異値・特異ベクトルの数が明確な場合、本記事で実装したsvdメソッドで速くなることが確認できた。

$ ruby svd_bench.rb
                 user     system      total        real
Numo::Linalg 14.219540   1.152062  15.371602 ( 11.201133)
svd          1.530432   0.291214   1.821646 (  0.828056)

おわりに

打ち切り特異値分解の手法としては、べき乗法やQR分解を利用したもの(ML Wikiというものにいい感じにまとまっていた)がある。いずれもNumo::NArrayで実装してみたが、本記事にある固有値分解を用いたものを超えるパフォーマンスを得られなかった。一方で、行列補完(Matrix Completion)などでは、最大の特異値・特異ベクトルだけが欲しいというアルゴリズムもある。そういったものには、べき乗法が適していると思う。

はてなブログProを解約してしまった

筆不精なので「お金を払えば書くようになるかな?」と思ったら、そうでもなかったです。ちなみに、解約してスグに元に戻るのではなくて、払った分はProでいられます。良心的ですね。

f:id:yoshoku:20181217122510p:plain

RubyでLIBSVM形式のファイルを扱うLibSVMLoaderをバージョンアップした

久しぶりにLibSVMLoaderをアップデートした。

$ gem install libsvmloader

これまでNMatrixで特徴ベクトルとラベルを表現していたが、これをRubyのArrayに変えた。Arrayにすることで、NMatrixだけでなくNumo::NArrayでも使用できる。

require 'libsvmloader'
require 'nmatrix/nmatrix'

# LIBSVM形式のファイルを読み込む。
samples, labels = LibSVMLoader.load_libsvm_file('foo.t')

# ArrayをもとにNMatrixによるサンプル行列・ラベルベクトルを生成する。
samples_nm = N[*samples]
labels_nm = N[*labels]
require 'libsvmloader'
require 'numo/narray'

# LIBSVM形式のファイルを読み込む。
samples, labels = LibSVMLoader.load_libsvm_file('foo.t')

# ArrayをもとにNumo::NArrayによるサンプル行列・ラベルベクトルを生成する。
samples_na = Numo::NArray[*samples]
labels_na = Numo::NArray[*labels]

また、liblinear-rubyは、RubyのArrayでデータをやりとりするので、より使いやすくなる。

$ gem install liblinear-ruby

例として、pendigitsデータセットを分類する。

$ wget https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multiclass/pendigits
$ wget https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multiclass/pendigits.t

訓練は次のようになる。

require 'libsvmloader'
require 'liblinear'

# pendigitsの訓練データセットを読み込む。
samples, labels = LibSVMLoader.load_libsvm_file('pendigits')

# Liblinearで線形サポートベクターマシンを学習する。
Liblinear.quiet_mode
model = Liblinear.train(
  { solver_type: Liblinear::L2R_L2LOSS_SVC },
  labels,
  samples
)

# 学習したモデルを保存する。
model.save('pendigits.model')

そしてテストは次のようになる。

require 'libsvmloader'
require 'liblinear'

# pendigitsのテストデータセットを読み込む。
samples, labels = LibSVMLoader.load_libsvm_file('pendigits.t')

# 学習したモデルを読み込む。
model = Liblinear::Model.load('pendigits.model')

# テストデータのラベルを推定する(predictメソッドにはサンプルを1つずつ渡す)。
predicted = samples.map { |s| Liblinear.predict(model, s).to_i }

# 正解ラベル数をカウントする。
n_samples = labels.size
n_corrects = Array.new(n_samples) { |n| labels[n] == predicted[n] ? 1 : 0 }.count(1)

# 分類の正確度を出力する。
puts "Accuracy: %.4f" % (n_corrects.fdiv(n_samples))

これらのスクリプトを実行すると、正確度が出力される。約9割が正解している。

Accuracy: 0.8788

だいたい以上のような感じ。インターフェースはあまり変えなかった。label_dtypeやvalue_dtypeオプションを使用すると、ラベルや特徴ベクトルの型を指定できる。また、LIBSVM形式では、特徴ベクトルの添字が1から始まるが、zero_basedオプションにtrueを渡せば、0から始まるとして読み込みを行う。詳しくはドキュメントで。

https://www.rubydoc.info/gems/libsvmloader/0.2.0/LibSVMLoader