Ruby拡張でNumo::NArrayのデータをポインタで取得する

はじめに

Ruby拡張で、Numo::NArrayのデータをC言語の配列のように扱えないかな〜と思っていたら、na_get_pointer_for_read（for_writeもある）という関数が用意されていて簡単にできた。

準備

もろもろ自動で用意されて便利なので、Ruby拡張を含むgemを作る形で進める。bundle install前に、gemspecファイル中のTODOを削除したり、URLの部分はコメントアウトしたりする必要がある。

$ bundle gem hoge --ext
Creating gem 'hoge'...
MIT License enabled in config
Code of conduct enabled in config
      create  hoge/Gemfile
      create  hoge/lib/hoge.rb
... （省略）

$ cd hoge
$ vim hoge.gemspec
... 「TODO:」を削除する
... spec.homepageなどURLを書く部分をコメントアウトする
... Numo::NArrayを依存関係に追加する
spec.add_runtime_dependency 'numo-narray'
... （省略）

$ bundle install
...

extconfにNumo::NArray関連の記述を追加する。

$  vim ext/hoge/extconf.rb

require 'mkmf'
require 'numo/narray'

$LOAD_PATH.each do |lp|
  if File.exist?(File.join(lp, 'numo/numo/narray.h'))
    $INCFLAGS = "-I#{lp}/numo #{$INCFLAGS}"
    break
  end
end

create_makefile('hoge/hoge')

ライブラリでもNumo::NArrayを読み込むようにする。

$ vim lib/hoge.rb

require 'numo/narray'
require 'hoge/version'
require 'hoge/hoge'

module Hoge
  class Error < StandardError; end
  # Your code goes here...
end

拡張の作成

関連ファイルのinclude

まずは、Numo::NArray関連のヘッダーファイルをincludeする。

$ vim ext/hoge/hoge.h

#ifndef HOGE_H
#define HOGE_H 1

#include "ruby.h"
#include "numo/narray.h"
#include "numo/template.h"

#endif /* HOGE_H */

Numo::NArrayのデータをポインタで取得する

Numo::DFloatの中身を表示するメソッドをもつモジュールを作る。

$ vim ext/hoge/hoge.c

#include "hoge.h"

VALUE rb_mHoge;

static VALUE
print_mat(VALUE self, VALUE vnary)
{
  narray_t *nary;
  size_t i, j;
  size_t n_rows, n_cols;
  double* data;

  /* VALUEなNArrayからnarray_t構造体を得る */
  GetNArray(vnary, nary);

  /* ホントはNumo::DFloatか確認するコードがあると良い
   * if (CLASS_OF(vnary) != numo_cDFloat) { ... とか
   */

  /* ホントは次元数が2か（行列かどうか）のコードがあると良い
   * if (NA_NDIM(na) != 2) { ... とか
   */

  /* NArrayの行数・列数を得る */
  n_rows = NA_SHAPE(nary)[0];
  n_cols = NA_SHAPE(nary)[1];
  printf("(%zd, %zd)\n", n_rows, n_cols);

  /* Numo::DFloatから配列の先頭ポインタを得る */
  data = (double*)na_get_pointer_for_read(vnary);

  /* 配列の中身を表示する */
  for (i = 0; i < n_rows; i++) {
    printf("[ ");
    for (j = 0; j < n_cols; j++) {
      printf("%.4f ", data[i * n_cols + j]);
    }
    printf("]\n");
  }

  return Qnil;
}

void
Init_hoge(void)
{
  rb_mHoge = rb_define_module("Hoge");
  rb_define_module_function(rb_mHoge, "print_mat", print_mat, 1);
}

動作確認

コンパイルして動作を確認する。渡したNumo::DFloatの中身を取得して表示することができた。

$ rake compile
...

$ bin/console

irb(main):001:0> a = Numo::DFloat.new(3,2).rand
=> Numo::DFloat#shape=[3,2]
[[0.0617545, 0.373067],
 [0.794815, 0.201042],
 [0.116041, 0.344032]]
irb(main):002:0> Hoge.print_mat(a)
(3, 2)
[ 0.0618 0.3731 ]
[ 0.7948 0.2010 ]
[ 0.1160 0.3440 ]
=> nil

おわりに

Ruby拡張を作っていると、外部ライブラリとアレコレするとか、完全にC言語の世界に行きたい場合があって、Arrayがポインターで取れたら良いのにな〜と思ってたら簡単にできた。

Rumaleにガウス混合モデルによるクラスタリングを追加した

はじめに

Rumaleにガウス混合モデル（Gaussain Mixture Model, GMM）によるクラスタリングを追加して、ver. 0.12.2としてリリースした。

rumale | RubyGems.org | your community gem host

GMMは、データ分布をいくつかの正規分布の重み付き線形和で表現しようというもので、このいくつかの正規分布がクラスタとなる。GMMは、与えられたデータから、混合比・平均ベクトル・共分散行列で表されるクラスタを解析する。共分散行列の表現方法でいくつかパターンがあるが、Rumaleでは対角要素のみを用いる方法を実装した。これは、混合比を求めるために、共分散行列の行列式と逆行列を計算する必要があり、共分散行列が対角行列であると、これら計算を簡略化できることから選択した。

使い方

Rumaleはgemコマンドでインストールできる。Numo::NArrayに依存している。

$ gem install rumale

データセットの読み込みでred-datasets、データのplotにNumo::Gnuplotを使いたいので、これらもインストールする。※別途gnuplotをインストールする必要がある

$ brew install gnuplot
$ gem install numo-gnuplot red-datasets-numo-narray

Scikit-LearnのGMMの例を試してみる。アヤメデータをGMMでクラスタリングして、1次元目と2次元目の特徴を使って散布図で可視化する。

require 'numo/linalg/autoloader'
require 'numo/gnuplot'
require 'datasets-numo-narray'
require 'rumale'

# 楕円の大きさを求める.
def ellipse_size(cov11, cov22)
  cov = Numo::DFloat.zeros(2, 2)
  cov[0, 0] = cov11
  cov[1, 1] = cov22
  v, = Numo::Linalg.eigh(cov)
  sz = 2 * Math.sqrt(2) * Numo::NMath.sqrt(v)
  sz.to_a
end

# Numo::NArray形式でIRISデータセットを読み込む.
dataset = Datasets::LIBSVM.new('iris').to_narray
labels = Numo::Int32.cast(dataset[true, 0])
samples = Numo::DFloat.cast(dataset[true, 1..-1])

# GMMでクラスタリングする.
gmm = Rumale::Clustering::GaussianMixture.new(n_clusters: 3, random_seed: 5)
cluster_ids = gmm.fit_predict(samples) # ※クラスタラベルは今回の可視化では使用しない.

# Numo::GnuplotでPNGファイルに書き出す.
## サンプル
x = samples[true, 0]
y = samples[true, 1]
plots = labels.to_a.uniq.sort.map { |l| [x[labels.eq(l)], y[labels.eq(l)], t: l.to_s] }
## 各クラスタの平均ベクトル
plots.push([Numo::DFloat[gmm.means[2, 0]], Numo::DFloat[gmm.means[2, 1]], t: 'center 1'])
plots.push([Numo::DFloat[gmm.means[0, 0]], Numo::DFloat[gmm.means[0, 1]], t: 'center 2'])
plots.push([Numo::DFloat[gmm.means[1, 0]], Numo::DFloat[gmm.means[1, 1]], t: 'center 3'])
## 各クラスタの共分散による楕円とともにプロットする.
Numo.gnuplot do
  set(terminal: 'png')
  set(output: 'iris.png')
  set(:object, 1, 'ellipse',
      center: [gmm.means[2, 0], gmm.means[2, 1]],
      size: ellipse_size(gmm.covariances[2, 0], gmm.covariances[2, 1]),
      angle: 0.0, front: true, fs: true, empty: true, bo: 1)
  set(:object, 2, 'ellipse',
      center: [gmm.means[0, 0], gmm.means[0, 1]],
      size: ellipse_size(gmm.covariances[0, 0], gmm.covariances[0, 1]),
      angle: 0.0, front: true, fs: true, empty: true, bo: 2)
  set(:object, 3, 'ellipse',
      center: [gmm.means[1, 0], gmm.means[1, 1]],
      size: ellipse_size(gmm.covariances[1, 0], gmm.covariances[1, 1]),
      angle: 0.0, front: true, fs: true, empty: true, bo: 3)
  plot(*plots)
end

これを実行すると、以下のような画像を得られる。Scikit-LearnのGMMの例でいえば、covariance_typeをdiagにした図に相当する。対角要素のみで当てはめると、楕円の回転が考慮されない（例えば左上のクラスタの楕円は斜め45度くらい傾いていると分布によりフィットした感じになるがそれがない）。それだけ分布にフィットできていないことになるわけだが、親戚的な手法であるK-Meansがクラスタ中心のみによる表現であったり、あんまりフィットし過ぎても良くない場合もあったりすることから、クラスタリングの実利用としては割とこれで十分だったりもする。

おわりに

Rumaleのver. 0.12.2では、確率モデルによるクラスタリング手法であるガウス混合モデル（Gaussain Mixture Model, GMM）を実装した。また、GMMの他に、オマケでカテゴリ変数を整数値に変換するOrdinalEncoderも追加した。こういったfeature engineeringな前処理クラスを足していっても良いかもな〜とも思ったり。

github.com

Rumaleの乱数生成を無難なものにした

はじめに

機械学習アルゴリズムでは、乱数でベクトルを初期化したり、ランダムサンプリングしたりなど、乱数生成をアルゴリズム中に含むものが多い。Rumaleの多くもそんな感じで、クラスのインスタンス変数にRandomクラスによる乱数生成器を持っている。これの扱いを良い感じにした。

rumale | RubyGems.org | your community gem host

良い感じにしたとは？

Rumaleのクラス内で、乱数生成器は、簡単に書くと以下のようになっていた。

class Hoge
  attr_reader :weight

  def initialize(random_seed: 0)
    @rng = Random.new(random_seed)
  end

  def fit
    # 重みを初期化する
    @weight = @rng.rand
  end
end

インスタンス生成として乱数生成器を使いまわすかたちになっている。例えば、乱数で重みベクトルを初期化するようなアルゴリズムの場合に、学習データが同じであっても、fitメソッドを呼び出すたびに、学習結果が異なる状態にあった。

> h = Hoge.new
> h.fit
> h.weight
=> 0.5488135039273248
> h.fit
> h.weight
=> 0.7151893663724195

同じデータ・パラメータでfitメソッドを何度も呼び出すコトがあるか？は置いておいて、同じデータを入力したら同じ結果が出て欲しいのが人情なので、fitメソッドを以下のように修正した。

def fit
  sub_rng = @rng.dup
  # 重みを初期化する
  @weight = sub_rng.rand
end

initializeで作成した乱数生成器をコピーするようにした。これで、fitメソッドを呼び出すたびに、学習結果（乱数での初期化やランダムサンプリング）が変化することはなくなった。

> h = Hoge.new
> h.fit
> h.weight
=> 0.5488135039273248
> h.fit
> h.weight
=> 0.5488135039273248

他にこの現象の回避方法はいくらでもあるだろうが、すでにMarshal.dumpされた学習済みモデルがあることを考えると、最短距離の対処法だと思う。

おわりに

こういう下回り的なところも気をつけていかないと。

github.com

RumaleでParallelを導入して高速化した

はじめに

Rumaleには、実行速度が遅いという問題があった。これに対して、version 0.11.0 では、one-vs-the-restやbaggingで並列化できる部分を、Parallelにより並列化してみた。

rumale | RubyGems.org | your community gem host

scikit-learnと同様に、n_jobsパラメータを持つ推定器では並列化が利用できる。ただし、試験的なモノなので、明確にParallelをrequireした場合にのみ、有効になるようにした。

使い方

Rumaleは、gemコマンドでインストールできる。データセットの読み込みでred-datasetsを使いたいので、あわせてインストールする。

$ gem install rumale red-datasets-numo-narray

そして、Parallelをインストールする。

$ gem install parallel

USPSデータセットを用いて、Random Forestによる分類器を例に、どれくらい速くなるか確認する。RandomForestClassifierでは、baggingで弱学習器を訓練する部分が並列化される。実行環境は、MacBook Early 2016で、 Intel Core m3 1.1Ghzである。

require 'datasets'
require 'datasets-numo-narray'
require 'rumale'

# Parallelがrequireされていない場合n_jobsパラメータは無視される.
require 'parallel' 

# データセットを読み込む.
datasets = Datasets::LIBSVM.new('usps').to_narray
labels = Numo::Int32.cast(datasets[true, 0])
samples = Numo::DFloat.cast(datasets[true, 1..-1])

# 訓練とテストデータセットに分割する.
ss = Rumale::ModelSelection::ShuffleSplit.new(n_splits: 1, test_size: 0.2, random_seed: 1)
train_ids, test_ids = ss.split(samples, labels).first
train_s = samples[train_ids, true]
train_l = labels[train_ids]
test_s = samples[test_ids, true]
test_l = labels[test_ids]

# Random Forestによる分類器を学習する.
start = Time.now
est = Rumale::Ensemble::RandomForestClassifier.new(
  n_estimators: 100,
  max_depth: 10,
  max_features: 8,
  # n_jobs: -1, # -1にするとcore数と同じだけプロセスを作る. 
  random_seed: 1
)
est.fit(train_s, train_l)

# テストデータで正確度を評価する.
puts("Accuracy: %.1f %%" % (100.0 * est.score(test_s, test_l)))
finish = Time.now
puts("Time: %.1f [sec]" % (finish - start))
~

まずは、並列化なし。

$ ruby hoge.rb
Accuracy: 95.1 %
Time: 31.0 [sec]

そして、n_jobsのオプションをコメントアウトを外して、有効にしたもの。

$ ruby hoge.rb
Accuracy: 95.3 %
Time: 15.9 [sec]

だいぶ速くなった。この他、one-vs-the-restの並列化では、線形SVMによる分類器SVCや勾配ブースティング木のGradientBoostingClassifierなどで実行速度が速くなるのを確認できた。

n_jobsを有効にする/しないで結果がわずかに変わった。これは、アルゴリズム中に乱数生成を含む場合、並列化することで乱数が生成されるタイミングが変わることに起因する。今後、乱数で初期化する部分など、乱数生成に関わる実装を再検討していく。

おわりに

Parallelを利用することで、Rumaleを一部高速化することができた。高速化の本質的な解決策としては、やはり、大部分をRuby Extensionsに書き換えるのが良いのだろうが、ある程度アルゴリズムがそろってからかな〜と考えている。

github.com

Rumaleにt-SNEを実装した

はじめに

Rumaleに高次元データの可視化として定番の一つになっている t-distributed Stochastic Neighbor Embedding（t-SNE）を実装した。「教師なし学習を増やさないとな〜少ないよな〜」と思っているところに、issueでリクエストを頂いたので実装してver. 0.10.0としてリリースした。

rumale | RubyGems.org | your community gem host

t-SNEが最初に提案された当時、Numpyで実装した経験があるが、勾配法により最適化するのもあって結構遅かった。t-SNEの高速化については、いくつか提案されているが、Rumaleでは不動点法による手法を採用した。

Z. Yang, I. King, Z. Xu, and E. Oja, “Heavy-Tailed Symmetric Stochastic Neighbor Embedding,” Proc. NIPS'09, pp. 2169–2177, 2009.

論文では、Symmetric SNE全般で使える手法として提案されている。t-SNEもSymmetric SNEの一種なので、これを応用できる。不動点法は、独立成分分析のFast ICAなどでも使用されており、一般に勾配法よりも速く収束する。また、不動点法は、勾配法と異なり、学習率やモーメンタム係数といったハイパーパラメータを必要としない。これらの利点から採用した。

使い方

Rumaleはgemコマンドでインストールできる。Numo::NArrayに依存している。

$ gem install rumale

データセットの読み込みでred-datasetsを使いたいので、これもインストールする。

$ gem install red-datasets-numo-narray

可視化したデータのplotにはNumo::Gnuplotを用いる。

$ brew install gnuplot
$ gem install numo-gnuplot

また、t-SNEはアルゴリズム上行列計算を繰り返すことになるので、Numo::Linalgのインストールすることをお薦めする。Intel MKLやOpenBLASと組み合わせると、行列計算が並列化されるので、計算が速くなる。詳しくはコチラを。

$ gem install numo-linalg

MNISTデータセットの可視化を試みる。

require 'rumale'
require 'numo/linalg/autoloader'
require 'datasets-numo-narray'
require 'numo/gnuplot'

# Numo::NArray形式でMNISTデータセットを読み込む.
dataset = Datasets::LIBSVM.new('mnist').to_narray
labels = Numo::Int32.cast(dataset[true, 0])
samples = Numo::DFloat.cast(dataset[true, 1..-1])

# データ数が多いと遅いので適当にサブサンプリングした.
rand_id = Array.new(samples.shape[0]) { |n| n }.sample(3000)
labels = labels[rand_id]
samples = samples[rand_id, true]

# データを[-1,1]の範囲に正規化する.
normalizer = Rumale::Preprocessing::MinMaxScaler.new(feature_range: [-1.0, 1.0])
samples = normalizer.fit_transform(samples)

# 元データの次元数が大きいと, 高次元空間の確率を計算する部分で計算時間を要する.
# MNISTの特徴量はほとんど背景部分の0なので, 主成分分析で1/10程度に削減する.
pca = Rumale::Decomposition::PCA.new(n_components: 80, random_seed: 1)
samples = pca.fit_transform(samples)

# t-SNEで2次元空間にマッピングすることで可視化する.
# verboseでtrueを指定すると, 最適化過程の数値を出力する.
tsne = Rumale::Manifold::TSNE.new(
  perplexity: 30.0, max_iter: 1000, verbose: true, random_seed: 1)
low_samples = tsne.fit_transform(samples)

# Numo::GnuplotでPNGファイルに書き出す.
x = low_samples[true, 0]
y = low_samples[true, 1]
plots = labels.to_a.uniq.sort.map { |l| [x[labels.eq(l)], y[labels.eq(l)], t: l.to_s] }

Numo.gnuplot do
  set(terminal: 'png')
  set(output: 'tsne.png')
  unset(:xtics)
  unset(:ytics)
  unset(:border)
  plot(*plots)
end

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

$ ruby rumale_tsne.py
[t-SNE] Computed conditional probabilities for sample 1000 / 3000
[t-SNE] Computed conditional probabilities for sample 2000 / 3000
[t-SNE] Computed conditional probabilities for sample 3000 / 3000
[t-SNE] Mean sigma: 4.30353336663795
[t-SNE] KL divergence after 100 iterations: 2.2765667341716527
[t-SNE] KL divergence after 200 iterations: 1.9882779434413718
[t-SNE] KL divergence after 300 iterations: 1.856052564481442
[t-SNE] KL divergence after 400 iterations: 1.7744255471415573
[t-SNE] KL divergence after 500 iterations: 1.7163771284160634
[t-SNE] KL divergence after 600 iterations: 1.6715272565792092
[t-SNE] KL divergence after 700 iterations: 1.6372943878281605
[t-SNE] KL divergence after 800 iterations: 1.6084045578469424
[t-SNE] KL divergence after 900 iterations: 1.5847659188134602
[t-SNE] KL divergence after 1000 iterations: 1.5633048653537873

そして、出力された結果は次のようになる。色のバリエーションがアレだったりして分かりにくいが、MNISTデータセットの各手書き数字が、種類別にある程度固まって可視化されている。

おわりに

可視化手法には、その他、多次元尺度構成法（Multi-dimensional scaling, MDS）やSammon nonlinear mappingなどがある。これらも順次、実装していきたい。t-SNEのように多様体学習ベースの次元削減手法という文脈では、Laplacian EigenmapsやLocally Linear Embeddingもある。これらは、k-近傍グラフの作成など、近傍探索が必要になる。データ数が大きくなると遅くなるので、近似最近傍探索の実装も考えないといけない。まだまだRumaleでできることあるな〜という感じ。

github.com