Rubyで疎行列演算を高速に行う

背景

Rubyの環境で行列演算を高速に行う必要があった。
行列計算ライブラリを利用して、疎行列の演算についての知見が少し得られたのでメモ。

※疎行列: ほとんどの要素が0の行列, スパースな行列

行列計算ライブラリ - NMatrix

Rubyで行列計算を高速に扱えるライブラリを調べるとNArrayとNMatrixの二つが候補に上がってくる。
NArrayの方は現時点で疎行列計算に対応していないらしく、今回はNMatrixを利用する。
NMatrixは3種類のデータ構造が利用でき、そのうち2種類が疎行列に適している。

NMatrixのインストール

$ gem install nmatrix

基本的にはgemでインストールできるのだが、C++11オプションが使えるgccが利用できないとエラーになる模様。
特にMacの場合はgccのバージョンが低い場合がよくある。
gcc4.7以上なら大丈夫らしいので、必要に応じてインストールする。

NMatrixで疎行列を扱う

基本的なNMatrixの使い方はこちら。

Getting started · SciRuby/nmatrix Wiki
https://github.com/SciRuby/nmatrix/wiki/Getting-started

以下では、特に疎行列を扱う際のTipsを挙げる

行列の初期化

NMatrixでは3種類のデータ構造が利用でき、stypeで指定する。
stypeとしてyaleを利用すると、疎行列演算が高速に行える。

# 20000行30000列の行列をyaleで初期化
n = N.new([20000, 30000], stype: :yale)

stype毎の演算速度比較

stypeには、yale・dense・listの3種類があるが
以下のようなコードでそれぞれの演算速度を計測してみた。

require 'nmatrix'
require 'benchmark'

# 5000行, 1000列の行例で計測
rows, cols = 5000, 1000
puts "rows: #{rows}, cols: #{cols}"

Benchmark.bm(10) do |x|
  [:yale, :dense, :list].each do |stype|
    # stypeを指定して初期化
    c = N.new([rows, cols], dtype: :byte, stype: stype)
    # 各行10要素ずつ1を入れる
    c[0..-1, 0..9] = 1
    # 加算
    x.report("#{stype}:add") do
      c + c
    end
    # 乗算
    x.report("#{stype}:multi") do
      c * c
    end
  end
end

yaleが最も高速に演算できる。
listは0を含め全ての要素が演算対象になるので遅い。

rows: 5000, cols: 1000
                 user     system      total        real
yale:add     0.120000   0.000000   0.120000 (  0.128647)
yale:multi   0.120000   0.010000   0.130000 (  0.146060)
dense:add    1.230000   0.060000   1.290000 (  1.384727)
dense:multi  1.150000   0.020000   1.170000 (  1.185778)
list:add    32.600000   0.160000  32.760000 ( 32.908609)
list:multi  30.430000   0.190000  30.620000 ( 31.503494)

yaleタイプの説明

The storage of a matrix - How to create an NMatrix

効率的に非0要素にアクセスする

NMatrixでは、メモリ上に保持されている要素にのみアクセスするメソッドがある。
疎行列で非0の要素全てにアクセスしたい場合は、全要素でループを回すよりもこちらを使った方が速い。
※保持されているからといって値が0でないとは限らない点に注意

• each_stored_with_indices
• each_ordered_stored_with_indices

# 保持されている要素のみでループを回す
n.each_stored_with_indices do |value, row, col|
  puts "[#{row}, #{col}] = #{value}"
end

サンプル

require 'nmatrix'

# 10000行, 10000列の行列を準備
n = N.new([10000, 10000], dtype: :byte, stype: :yale)

# 6要素だけ1を入れる
n[10..11, 50..52] = 1

count = 0
n.each_ordered_stored_with_indices do |value, row, col|
  # ループ数をカウント
  count += 1
  # 初期化した時点で0が保持されている要素もある
  next if value == 0
  # 非0要素を出力
  puts "[#{row}, #{col}] = #{value}"
end
# ループ数(=メモリ上に保持されていた要素の数)
puts "count: #{count}"

以下の出力を見ると、値を設定した6要素全てにアクセスできた。
実際にブロックが実行されたのは、1億(10000*10000)要素中の10006要素だけ。
効率的にループを回せていることがわかる。

[10, 50] = 1
[10, 51] = 1
[10, 52] = 1
[11, 50] = 1
[11, 51] = 1
[11, 52] = 1
count: 10006

参考

SciRuby/nmatrix: Dense and sparse linear algebra library for Ruby via SciRuby
https://github.com/sciruby/nmatrix

class NMatrix - RDoc Documentation
http://sciruby.com/nmatrix/docs/NMatrix.html

Numerical Ruby NArray
https://masa16.github.io/narray/index.ja.html