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

RumaleにGradient Tree Boostingによる分類と回帰を実装した

はじめに

Rumaleでは決定木系のアルゴリズムの高速化と追加を進めている。ついに人気のGradient Tree Boosting（Gradient Boosting MachineやGradient Boosted Regression Treeなどとも呼ばれる）を実装して、ver. 0.9.2としてリリースした。

rumale | RubyGems.org | your community gem host

Gradient Tree Boosting (GTB)は、PythonではXGBoost、LightGBM、CatBoostなどで有名なアルゴリズムである。 GTBを実装しようと思ったのは、Scikit-Learnでも実装されているのと、ver. 0.2.1でLightGBMを参考にした実装であるHistGradientBoostingClassifier/Regressorが追加されるのを知って、Rubyで実装してみてみるか〜という気持ちになったのもあって。

使い方

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

$ gem install rumale

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

$ gem install red-datasets-numo-narray

LIBSVM Dataのcpusamllデータセット（データ数 8,192、次元数12）を用いた回帰を試してみる。

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

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

# ランダムに訓練とテストに分ける.
ss = Rumale::ModelSelection::ShuffleSplit.new(n_splits: 1, test_size: 0.1, random_seed: 1)
train_ids, test_ids = ss.split(samples, values).first
train_s = samples[train_ids, true]
train_v = values[train_ids]
test_s = samples[test_ids, true]
test_v = values[test_ids]

# Gradient Boosting Treeによる回帰を訓練する.
# ※ハイパーパラメータは勘で決めている.
est = Rumale::Ensemble::GradientBoostingRegressor.new(
  n_estimators: 100,   # 生成する決定木の数、Boostingにおける繰り返し数でもある
  learning_rate: 0.1,  # 学習率
  reg_lambda: 0.001,   # L2正則化の係数
  subsample: 0.8,      # ランダムサンプリングする際のデータの割合
  max_depth: 4,        # 決定木の深さ
  max_features: 8,     # 使用する特徴数、column samplingとも呼ばれる
  random_seed: 1       # 乱数のシード
)
est.fit(train_s, train_v)

# テストセットの決定係数（1に近づくほどよい）確認する.
puts("GTB R2-Score: %.4f" % est.score(test_s, test_v))

# 比較のためにRandom Forestでも同様のことを行う.
est = Rumale::Ensemble::RandomForestRegressor.new(
  n_estimators: 100,
  max_depth: 4,
  max_features: 8,
  random_seed: 1
)
est.fit(train_s, train_v)
puts("RF R2-Score: %.4f" % est.score(test_s, test_v))

これを実行すると以下のようになる。GTBの方が良い値を得られている。

GTB R2-Score: 0.9703
RF R2-Score: 0.9433

分類器のRumale::Ensemble::GradientBoostingClassifierも、同様のパラメータと手順で利用できる。

特徴量の離散化

LightGBMなどでは、特徴量を離散値に変換すること（離散化）で高速な計算を実現している。決定木では特徴量の値（特徴ベクトルをuniqして残る値）が、木を分割する際の閾値の候補となる。離散化することで候補値が少なくなると、分割の評価計算の回数が少なくなるので、そのぶん速くなる。この離散化込みでアルゴリズムを考え、全体的に高速化しているものもあるが、Rumaleでは、離散化したくない場合もあると思い、アルゴリズムとは別で、特徴量を離散化する BinDiscretizer クラスを用意した。実行例は以下のようになる。[-1, 1]な実数値を4段階に離散化している。

irb(main):001:0> require 'rumale'
=> true
irb(main):002:0> t = Rumale::Preprocessing::BinDiscretizer.new(n_bins: 4)
=> #<Rumale::Preprocessing::BinDiscretizer... 省略
irb(main):003:0> x=Numo::DFloat.new(5, 3).rand - 0.5
=> Numo::DFloat#shape=[5,3]
[[-0.438246, -0.126933, 0.294815],
 [-0.298958, -0.383959, -0.155968],
 [0.039948, 0.237815, -0.334911],
 [-0.449117, -0.391935, -0.431292],
 [0.404121, -0.0213559, -0.157031]]
irb(main):004:0> t.fit_transform(x)
=> Numo::DFloat#shape=[5,3]
[[0, 1, 3],
 [0, 0, 1],
 [2, 3, 0],
 [0, 0, 0],
 [3, 2, 1]]

離散化することでGTB高速になるかを以下のコードで確認した。

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

datasets = Datasets::LIBSVM.new('cpusmall').to_narray
values = Numo::DFloat.cast(datasets[true, 0])
samples = Numo::DFloat.cast(datasets[true, 1..-1])

ss = Rumale::ModelSelection::ShuffleSplit.new(n_splits: 1, test_size: 0.1, random_seed: 1)
train_ids, test_ids = ss.split(samples, values).first
train_s = samples[train_ids, true]
train_v = values[train_ids]
test_s = samples[test_ids, true]
test_v = values[test_ids]

est = Rumale::Ensemble::GradientBoostingRegressor.new(
  n_estimators: 100,
  learning_rate: 0.1,
  reg_lambda: 0.001,
  subsample: 0.8,
  max_depth: 4,
  max_features: 8,
  random_seed: 1
)

Benchmark.bm 10 do |r|
  r.report 'non-transform' do
    est.fit(train_s, train_v)
    puts(" (R2-Score: %.4f)" % est.score(test_s, test_v))
  end

  r.report 'discretized' do
    # 4段階の離散値に変換する（ちょっと極端な例で実用では32段階以上が良いと思われる）
    t = Rumale::Preprocessing::BinDiscretizer.new(n_bins: 4)
    dis_train_s = t.fit_transform(train_s)
    dis_test_s = t.transform(test_s)
    est.fit(dis_train_s, train_v)
    puts(" (R2-Score: %.4f)" % est.score(dis_test_s, test_v))
  end
end

実行結果は以下のとおり。高速化されるが、予測精度は落ちるようだ。このあたりは、データセットの大きさや特徴量の次元数、パラメータの兼ね合いで変わってくると思われる。

                 user     system      total        real
non-transform (R2-Score: 0.9703)
  8.030000   0.810000   8.840000 (  8.866683)
discretized (R2-Score: 0.9203)
  6.330000   0.920000   7.250000 (  7.275490)

おわりに

Rumaleではベーシックなアルゴリズムの追加を進めてきたが、これでちょっとRumaleもモダンになったな？今後は、もう少し決定木関連の実装を進めつつ、教師なし学習のアルゴリズムを増やしたいと考えている。

github.com