洋食の日記

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

RandSVDで乱数のシードを固定できるように修正した

RandSVDをマイナー（タイニー？）バージョンアップした。RandSVD（というよりも乱択アルゴリズムの特異値分解）では、ランダム行列との積により行列を小さくする処理を含んでいる。なので、乱数のシードを固定できるようにしたほうが親切だと考えた。 ※ NMatrixに、numpy.random.seedとかみたいにグローバルに固定できる仕組みがあると、勝手に思い込んでいたのがそもそもの間違いでして…

require 'randsvd'

mat = NMatrix.rand [5, 5]
nb_singular_values = 2
nb_power_iterations = 3
random_seed = 1

u, s, vt = RandSVD.gesvd(mat, nb_singular_values, nb_power_iterations, random_seed)

単純にメソッド引数を後ろに増やしました。最初からキーワード引数とかなにかで、柔軟に引数増やせるように設計すれば良かったと反省…

randsvd | RubyGems.org | your community gem host

RandSVDを使った主成分分析によるMNISTデータセットの可視化

Ruby NMatrix 機械学習

RandSVDを作ったので、特異値分解（固有値分解）をベースにした機械学習アルゴリズムを実装していこうと思っている。 まずは、主成分分析によるMNISTデータセットの可視化を行った。 ※特異値分解と主成分分析の関係は、そのまま「主成分分析 特異値分解」でググると沢山でてくるので、ネットの海を汚さないためにも、ここでは割愛させて頂きました。

MNISTデータセットは、LIBSVMのページからダウンロードした。

$ wget https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multiclass/mnist.scale.t.bz2
$ bunzip2 mnist.scale.t.bz2

LibSVMLoaderで、MNISTデータセットを読み込み、中心化したデータセットを、RandSVDで特異値分解する。 得られた射影行列で２次元に射影し、GnuplotRBでプロットする。

require 'libsvmloader'
require 'randsvd'
require 'gnuplotrb'

# データを読み込む。
samples, labels = LibSVMLoader.load_libsvm_file('mnist.scale.t', stype: :dense)

# 中心化する。
nb_samples, = samples.shape
mean_vec = samples.mean(0)
centered = samples - mean_vec.repeat(nb_samples, 0)

# 特異値分解する。
nb_subspace_dimensions = 2
nb_power_iterations = 3
_u, _s, vt = RandSVD.gesvd(centered, 
                           nb_subspace_dimensions, nb_power_iterations)
projection_matrix = vt.transpose
#
# ※もちろん共分散行列を作っても良い。
# covariance_matrix = centered.transpose.dot(centered) / (nb_samples - 1)
# u, _s, _vt = RandSVD.gesvd(covariance_matrix, 
#                            nb_subspace_dimensions, nb_power_iterations)
# projection_matrix = u
#

# 低次元空間に射影する。
projected = centered.dot(projection_matrix)

# データをラベル別に色分けしてプロットする。
datasets = labels.uniq.map do |lid|
  x = projected.col(0).map.with_index { |e, n| e if labels[n] == lid }
  y = projected.col(1).map.with_index { |e, n| e if labels[n] == lid }
  GnuplotRB::Dataset.new([x.to_a, y.to_a], title: lid.to_s)
end

plot = GnuplotRB::Plot.new(*datasets, with: 'points', title: 'MNIST')
plot.to_png('mnist.png', size: [640, 480])

結果はこんな感じで上手くいってる。

f:id:yoshoku:20170908003948p:plain

自分がRubyで実現したいもので、共通化できる処理をGemにすると、自宅だけでなく会社とか他所のPCで使えて便利だし、もしかしたら同じ問題で悩んでいる遠い友人の手助けにもなるかもしれない。

乱択アルゴリズムな打ち切り特異値分解のgemを公開した

Ruby NMatrix 機械学習

※Truncated SVD (Singular Value Decomposition) って和訳ありますか？勝手に打ち切り特異値分解としました。

はじめに

機械学習アルゴリズムで特異値分解（もしくは固有値分解）をするときは、だいたい、大きい（もしくは小さい）方からk個だけ特異値/特異ベクトルを取り出すことが多い。NMatrixに実装されている特異値分解は、LAPACKのラッパーで、全ての特異値を求めるスタイル（Full SVD）なので、分解したい行列が大きい場合、実行時間が大きなものとなる。そこで、Randomized AlgorithmなTruncated SVDを実装した。

randsvd | RubyGems.org | your community gem host

NMatrixがARPACKとかに対応すると、存在意義を失うかも。

インストール

行列の扱いにはSciRubyのNMatrixを使用する。また、アルゴリズム中で特異値分解とQR分解を行うため、NMatirx-Lapackeに依存する。

$ gem install randsvd

使い方

内部の特異値分解で、NMatrix::LAPACK.gesvdを使うものと、NMatrix::LAPACK.gesddを使うもので二種類のメソッドを用意した。わかりやすさのため、メソッド名はNMatrixのものと同じにした。gesddの方が、分割統治法で速いらしいっすよ。

require 'randsvd'

# 適当に、分解する行列として、m x nの大きさの行列を用意する。
m = 1000
n = 100
input_matrix = NMatrix.rand([m, n])

# 欲しい特異値の数を定義する。
k = 10

# 特異値分解を実行する。分解したい行列と、欲しい特異値の数を指定する。
# 特異値が大きいほうからk個の特異値と特異ベクトルが得られる。
# ここが u, s, vt = RandSVD.gesdd(input_matrix, k) でも良い。
u, s, vt = RandSVD.gesvd(input_matrix, k)

# 上の例では、
#   u は m x kの大きさの左特異ベクトルが並ぶ行列、
#   s は k個の特異値が並ぶベクトル、
#   vt は k x nの大きさの右特異ベクトルが並ぶ行列となる。
# これら特異値と特異ベクトルで元の行列を復元する場合、以下の様になる。
reconstructed_matrix = u.dot(NMatrix.diag(s).dot(vt))

ベンチマーク

Core m3 1.1GHzで8GBメモリのMacBookで、ベンチマークとってみた。 10,000 x 1,000の大きさでランクが10な行列の特異値分解を行う。 欲しいのは特異値の大きい方から10個となる。 NMatrix::LAPACKの特異値分解は、実際には、必要な特異値を取り出す処理（例えば u=u[0..m - 1, 0..k - 1] みたいなの）が必要だけど、 速度に大した影響ないだろうと思って省略した。

require 'benchmark'
require 'randsvd'

m = 10000
n = 1000
k = 10

mat = NMatrix.rand([m,k]).dot(NMatrix.rand([k,n]))

Benchmark.bm 10 do |r|
  r.report 'NMatrix' do
    u, s, vt = mat.gesvd
  end

  r.report 'RandSVD' do
    u, s, vt = RandSVD.gesvd(mat, k)
  end
end

実行結果は以下のとおり。totalを見ると、RandSVDの方が約10倍くらい速い。

$ ruby bench_svd.rb
                 user     system      total        real
NMatrix     64.910000   1.300000  66.210000 ( 31.076849)
RandSVD      6.790000   0.860000   7.650000 (  5.485305)

アルゴリズム

乱択アルゴリズムな特異値分解は、色んなバリエーションが提案されている。 NIPSのワークショップの論文に、良い感じに要約されたアルゴリズムが掲載されていたので、主にこれを参考にして実装した。

P.-G. Martinsson et. al, “Normalized power iterations for the computation of SVD,” In Proc. of NIPS Workshop on Low-Rank Methods for Large-Scale Machine Learning, 2011.

乱択アルゴリズムな特異値分解では、ランダム行列由来の直交行列との積で、元の行列のサイズを縮めるということをやっている。 上記の論文では、（近似精度の高い?）直交行列を得るのに normalized power iterations というのを行う。 が、これをしなくても、そんなに精度的には変わらなかった（けど実行速度は落ちる）ので、RandSVDでは、デフォルトでは切っている。 normalized power iterationsを試したければ、gesvd/gesddメソッドで、繰り返し回数を指定するとよい（デフォルトでは0回になっている）。

# 例えば5回おこなう場合:
t = 5
RandSVD.gesvd(mat, k, t)

おわりに

つまらないものですが、よろしくお願い致します。

LIBSVM/SVMLight形式のデータの読み込み/書き込みをするgemを公開した

Ruby NMatrix 機械学習

はじめに

LIBSVM形式のデータセットファイルを読み込み、データとラベルをNMatrixで返すライブラリが欲しかったので、作成してgemパッケージを公開した。

libsvmloader | RubyGems.org | your community gem host

Pythonのsvmlight-loaderにお世話になっているので、リスペクトのつもりで「LibSVMLoader」としたけど、Wekaに同名のクラスがあった。Loaderと名付けたけど、svmlight-loaderと同じく書き込みもできる。

インストール

読み込んだデータやラベルをNMatrixで返すのでNMatrixに依存する。

$ gem install libsvmloader

使い方

読み込みにはload_libsvm_fileメソッドを、書き込みにはdump_libsvm_fileメソッドを用いる。ドキュメントはこちら。

require 'libsvmloader'

# 分類タスクを読み込む場合:
# samplesは、サンプル数x次元数のfloat64のNMatrix
# labelsは、サンプル数x1のint32のNMatrix
samples, labels = LibSVMLoader.load_libsvm_file('hoge.t')

# 回帰タスクを読み込む場合:
# ※ラベルに相当するNMatrixのdtypeを指定できるのでfloat64として実数とする
# samplesは、サンプル数x次元数のfloat64のNMatrix
# target_valuesは、サンプル数x1のfloat64のNMatrix
samples, target_values = LibSVMLoader.load_libsvm_file('hoge.t', label_dtype: :float64)

# libsvmファイルの特徴ベクトルの添字が「0」から始まる場合:
samples, labels = LibSVMLoader.load_libsvm_file('hoge.t', zero_based: true)

# 書き込み：
LibSVMLoader.dump_libsvm_file(samples, labels, 'foo.t')

# 書き込み（特徴ベクトルの添字を「0」から始めたい場合）:
LibSVMLoader.dump_libsvm_file(samples, labels, 'foo.t', zero_based: true)

線形SVMによる分類のサンプル

例として、libsvmloaderでデータを読み込んで、liblinear-rubyで線形SVMの訓練・テストをするのは以下の様になる。

require 'libsvmloader'
require 'liblinear'

# 訓練データを読み込む．
tr_samples, tr_labels = LibSVMLoader.load_libsvm_file('pendigits')

# 線形SVMによる分類器を訓練する．
Liblinear.quiet_mode
model = Liblinear.train(
  { solver_type: Liblinear::L2R_L2LOSS_SVC },
  tr_labels.to_flat_a,
  tr_samples.to_a
)

# テストデータを読み込む．
ts_samples, ts_labels = LibSVMLoader.load_libsvm_file('pendigits.t')

# テストデータのラベルを推定する．
pr_labels = ts_samples.to_a.map { |fv| Liblinear.predict(model, fv).to_i }

# 正解数をカウントする．
n_hits = pr_labels.map.with_index { |l, n| 1 if l == ts_labels[n] }.compact.sum

# Accuracyを出力する．
n_ts_samples = ts_labels.size
accuracy = 100.0 * (n_hits / n_ts_samples.to_f)
puts(format('Accuracy = %.4f%% (%d/%d)', accuracy, n_hits, n_ts_samples))

おわりに

つまらないものですが、よろしくお願い致します。

Liblinear-Rubyによる線形SVMを試すついでにカーネル近似も試した

Ruby NMatrix 機械学習

はじめに

Liblinearは、線形SVMの実装として有名なライブラリ/ツールである。 これをRubyから叩くライブラリとして、Liblinear-Rubyがある。 Rubyの配列で表現された、特徴ベクトルとラベルをわたすだけで、線形SVMの訓練・テストが行える。 これで線形SVMによる分類器が簡単にできるので、ついでにカーネル近似も試した。

インストール

Liblinear-Rubyでは、内部的にswigを使っているので、HomebrewでLiblinearだけでなくswigもインストールする。

$ brew install liblinear swig
$ gem install liblinear-ruby

libsvm形式ファイルの読み込み

Rubyでlibsvmフォーマットのファイルを読み込むメソッド書いていたが、どうにも汚くなってしまったので、 PyCallでScikit-Learnのload_svmlight_fileを叩くことにした。 トリッキーだけど、PyCallで得られたnumpy.arrayをNMatrixに変換することを、考えてみたかったのもある。 あんまりいい感じにならなかったけど。

require 'nmatrix/lapacke'

require 'pycall/import'
include PyCall::Import
pyfrom 'sklearn.datasets', import: 'load_svmlight_file'

def load_libsvm_file(filename)
  # PyCallで呼び出したload_svmlight_fileでlibsvm形式のデータを読み込む。
  py_samples, py_labels = load_svmlight_file.(filename, zero_based: false)
  # サンプル数と次元数を得て、サンプル用のNMatrixとラベル用のNVectorを用意する。
  n_samples, n_dimensions = py_samples.shape
  samples = NMatrix.zeros [n_samples, n_dimensions]
  labels = NVector.zeros n_samples
  # NMatrixとNVectorに変換する。
  n_samples.times do |r|
    labels[r] = py_labels[r].to_i
    n_dimensions.times do |c|
      samples[r,c] = py_samples[r,c].to_f
    end
  end
  return samples, labels
end

Liblinear-Rubyでの分類

まず、動作確認のためのサンプルデータとして、libsvmのサイトから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

これらを読み込んで、訓練・テストを行うコードは次の様になる。

# 訓練データを読み込む。
samples, labels = load_libsvm_file 'pendigits'

# Liblinear-Rubyで線形SVMな分類器の訓練を行う。
# データの入力にはRubyの配列を用いるので、
# to_aやto_flat_aでNMatrixを配列に変換する。
Liblinear.quiet_mode
model = Liblinear.train(
  { solver_type: Liblinear::L2R_L2LOSS_SVC },
  labels.to_flat_a,
  samples.to_a)

# テストデータを読み込む。
test_samples, test_labels = load_libsvm_file 'pendigits.t'

# テストデータのラベルを推定する。
# predictメソッドでは、テストデータをまとめて与えるのではなく、
# 特徴ベクトルの一つ一つを与える形になる。
pred_labels = test_samples.to_a.map do |smpl|
  Liblinear.predict(model, smpl).to_i
end

# 正解しているラベルをカウントする。
n_correctly_preds = 0
n_test_samples = test_labels.size
n_test_samples.times do |n|
  n_correctly_preds += 1 if test_labels[n] == pred_labels[n]
end

# Accuracyを計算して出力する。
accuracy = n_correctly_preds / n_test_samples.to_f
puts sprintf("Accuracy = %.4f%% (%d/%d)", 
  accuracy * 100.0, n_correctly_preds, n_test_samples)

これを実行すると、次の様な結果を得られた。 これは、当然ながら、Liblinearのコマンドで訓練・テストした場合のAccuracyと一致した。

$ ruby hoge.rb
Accuracy = 87.8788% (3074/3498)

カーネル近似による非線形分類

Liblinearで得られるのは線形分類器なので、特徴空間上で非線形に分布するデータは、正しく分類できない。 これを解決する方法で有名なのが、カーネル法だが、データ数が多いと重くなったり、下手をすると動かなかったりする。 この問題に対して、カーネルを近似するような変換を考えよう、というアプローチがある。 最も有名なのが、ガウスカーネルを近似した、Random Fourier Features（RFF）である。 RFFでは、D個の正規ランダムベクトルwと特徴ベクトルxとの積による変換で表される。

\displaystyle \phi(\mathbf{x})=\sqrt{1/D}[ \sin(\mathbf{w}_{1}^{\top}\mathbf{x}),\ldots,sin(\mathbf{w}_{D}^{\top}\mathbf{x}), \cos(\mathbf{w}_{1}^{\top}\mathbf{x}),\ldots,cos(\mathbf{w}_{D}^{\top}\mathbf{x}) ]

これは、正規ランダムベクトルを並べた正規ランダム行列Wを考えると、処理自体は行列-ベクトル積となるので重くない。

\displaystyle W=[\mathbf{w}_{1},\ldots,\mathbf{w}_{D}]^{\top}

このRFFの改良版が、NIPS'16においてGoogle ResearchのF X. Yuらが提案した、Orthogonal Random Features（ORF）である。 ORFでは、ランダム行列Wを、カイ分布からサンプルされた要素からなる対角行列S、一様分布なランダム直交行列Qで表す。 ※Qは正規ランダム行列GをQR分解することで得られるそうな。

\displaystyle W=\frac{1}{\sigma}SQ

また、次の簡単化したモノも提案されている。

\displaystyle W=\frac{\sqrt{d}}{\sigma}Q

直交化することで、カーネル近似の精度が上がるとのこと。 さらに論文では、Structured ORFとしてWalsh-Hadamard行列を用いたものも提案されている。

ORFの実装

今回、ORFを試そうと思ったのは、NMatrixのQR分解を試したかったからである。 NMatrixには、正規ランダム行列を生成するメソッドがないので、まずこれを実装する。よくあるBox-Muller法を用いた。

def randn(shape, mu=0.0, sigma=1.0)
  x = NMatrix.random shape
  y = NMatrix.random shape
  ((x.log * -2.0).sqrt * (y * 2.0 * Math::PI).sin) * sigma + mu
end

さらに、正規ランダム行列をQR分解して直交行列を得るメソッドと、特徴変換を行うメソッドは、次のようになる。

# Random Fourier Featuresのための変換行列を得る。
def rff_transform_mat(samples, sigma=100.0)
  n_samples, n_dimensions = samples.shape
  rand_mat = randn [n_dimensions, n_dimensions]
  transform_mat = rand_mat * (1.0 / sigma)
end

# Orthogonal Random Featuresのための変換行列を得る。
def orf_transform_mat(samples, sigma=100.0)
  n_samples, n_dimensions = samples.shape
  rand_mat = randn [n_dimensions, n_dimensions]
  orth_rand_mat, r = rand_mat.factorize_qr
  transform_mat = orth_rand_mat * (Math.sqrt(n_dimensions) / sigma)
end

# 特徴変換を行う。
def mapping(samples, transform_mat)
  n_samples, n_dimensions = samples.shape
  features = samples.dot transform_mat.transpose
  mapped = (features.sin.hconcat features.cos) * Math.sqrt(1.0 / n_dimensions.to_f)
end

これを、訓練データとテストデータに適応する。

# 訓練データを読み込む。
samples, labels = load_libsvm_file 'pendigits'

# 変換行列を計算し、特徴変換を行う。
trans_mat = orf_transform_mat samples
samples = mapping samples, trans_mat

...

# テストデータを読み込む。
test_samples, test_labels = load_libsvm_file 'pendigits.t'

# 特徴変換を行う。
test_samples = mapping test_samples, trans_mat

...

これで実験してみると、RFFのAccuracyが91.84%±1.41で、ORFのAccuracyが93.94%±0.67となった。 たしかに、ORFの方がわずかに良い。 データセットによりけりだろうけど、 元の特徴ベクトルのAccuracyが87.88%なので、カーネル近似による特徴変換により分類精度が向上することが確認できた。

おわりに

思いのほか長くなってしまったので、Liblinear-Rubyとカーネル近似の話で、記事を分けたほうが良かった感。