洋食の日記

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

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