2019-01-06

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）などでは、最大の特異値・特異ベクトルだけが欲しいというアルゴリズムもある。そういったものには、べき乗法が適していると思う。

2018-12-17

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

ポエム

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

f:id:yoshoku:20181217122510p:plain

2018-08-26

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

Ruby 機械学習データ分析

久しぶりに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

2018-06-30

SVMKitにクラスタリングと行列分解を実装した

Ruby SVMKit データ分析機械学習

svmkit | RubyGems.org | your community gem host

クラスタリングはK-MeansとDBSCAN、行列分解は主成分分析（Principal Component Analysis, PCA）と非負値行列因子分解（Non-negative Matrix Factorization, NMF）を実装した。K-Meansは、K-Means++による初期値設定を実装した。PCAはベーシックなアルゴリズムでは固有値分解が必要になるが、Numo::NArrayには固有値分解はない（Numo::Linalgにある）ので、今回は不動点法（Fixed-point algorithm）によるものを実装した。DBSCANとNMFはベーシックな手法を実装した。

これ以前に、0.4系では、RMSPropやNADAMといったOptimizerを実装した。Optimizerは深層学習で使われるイメージが強いが、最適化法が確率的勾配降下法（Stochastic Gradient Descent, SGD）であれば、SVMやロジスティック回帰でも有効である。動機としては、SGDなFactorization Machine（FM）が安定せず、その解決策に導入した。FMは、勾配に二乗が含まれるため、学習率のスケジュールがシンプルだと、一時的にでも大きな値が入ると、NaNが出まくるということがあった。また学習率の初期値が適切ではないと、ものすごく低い分類精度になったりした（これら最初は実装誤りを疑ったがlibfmなどでも同じ現象が見られた）。RMSProp系のOptimizerは、勾配を正規化するような形の学習率を求めるので、FMにおいて多少勾配が大きくなっても、NaNが出まくるような暴走はしなくなった。

これで、SVMKitでは、教師あり学習（分類・回帰）と教師なし学習（クラスタリング・次元圧縮）の基本的なことができるようになった。これからはコード整備とパフォーマンスの向上をメインに据えつつ、ちょこちょことアルゴリズムを増やしていこうと思う。

2018-05-26

SVMKitでの回帰の実装をだいたい終えました