RumaleにFastICAによる独立成分分析を追加した

はじめに

RumaleにFastICAによる独立成分分析（Independent Component Analysis, ICA）を追加した。ICAは、与えられたデータを、統計的に独立な元データが混合されてできたものと考え、元データを再構成したもの（独立な成分）を求める手法である。信号処理の信号分離から発展したもので、いまでは教師なし学習のベーシックな手法の一つとして定着している。これをRumaleにも実装した。

rumale | RubyGems.org | your community gem host

使い方

信号分離を例にFastICAの使い方を示す。データのplotにNumo::Gnuplotを使いたいので、Rumaleとあわせてインストールする。FastICAの実装では、Numo::Linalgを用いたのでこれもインストールする。

$ gem install rumale numo-linalg numo-gnuplot
$ brew install openblas gnuplot # 必要に応じて

サンプルコードは以下のようになる。2つの原信号（本来は未知で事前に知ることができない）がなんやかんやで混ざり、3つのセンサで観測されたとする（これが入力データとなる）。この3つの観測信号から、FastICAにより、2つの原信号を復元（推定）する。

require 'numo/linalg/autoloader'
require 'numo/gnuplot'
require 'rumale'

# 原信号を作成する.
# サイン波と矩形波にちょっとノイズを加えたもの（本来は未知の信号）.
n_samples = 1000
t = Numo::DFloat.linspace(0, 10, n_samples)
s1 = Numo::NMath.sin(2 * t)
s2 = Numo::NMath.sin(3 * t).sign
source = Numo::NArray.vstack([s1, s2]).transpose
source += 0.03 * Rumale::Utils.rand_normal([n_samples, 2], Random.new(1))
source -= source.mean(0)
source /= source.stddev(0)

# 観測信号（原信号が混ざってマイクとかで観測された感じの信号）
# 混合された2つの原信号は、3つのセンサで、3信号として観測されたとする.
mixing_mat = 0.2 * Rumale::Utils.rand_normal([2, 3], Random.new(1))
observed = source.dot(mixing_mat)

# 独立成分分析により観測信号から原信号を復元する.
fica = Rumale::Decomposition::FastICA.new(n_components: 2, max_iter: 500, tol: 1e-8)
reconstructed = fica.fit_transform(observed)

# 各信号をプロットする.
# 原信号
Numo.gnuplot do
  set(terminal: 'png')
  set(output: 'source.png')
  plot([source[true, 0], with: 'lines', title: '1'],
       [source[true, 1], with: 'lines', title: '2'])
end
# 観測信号
Numo.gnuplot do
  set(terminal: 'png')
  set(output: 'observed.png')
  plot([observed[true, 0], with: 'lines', title: '1'],
       [observed[true, 1], with: 'lines', title: '2'],
       [observed[true, 2], with: 'lines', title: '3'])
end
# 復元信号
Numo.gnuplot do
  set(terminal: 'png')
  set(output: 'reconst.png')
  plot([reconstructed[true, 0], with: 'lines', title: '1'],
       [reconstructed[true, 1], with: 'lines', title: '2'])
end

これを実行すると以下の三種類の信号の画像が得られる。

原信号のサイン波と矩形波が混合され、3つのセンサで観測された感じになっている。

原信号. ちょっとノイズが加わったサイン波と矩形波.

観測信号. 2つの原信号が混合され3つのセンサで観測されたイメージ.

この観測信号から、FastICAにより原信号に類似した信号を復元できている。

復元信号. もとのサイン波と矩形波がある程度復元できている.

ICAでは、独立性のみを頼りに推定を行う。注意点として、順序や大きさが原信号と異なるものが得られる場合があることがある（今回はわりと上手くいった）。

おわりに

ICAは、上記のサンプルで3信号から2信号を復元したように（3次元のデータから2次元の成分データを抽出したように）、次元削減として使うことも可能である。ICAには、主成分分析のような直交条件がないため、よりデータ分布にそった主軸をえることができる。ICA自体の研究は、まだまだ続いていて、様々な発展型が提案されている。

github.com

Rumaleのガウス混合モデルで共分散の種類を選べるようにした

はじめに

Rumaleのガウス混合モデル（Gaussian Mixture Model, GMM）によるクラスタリングでは、各クラスタの共分散行列には、対角要素のみの対角行列を使用していた。

yoshoku.hatenablog.com

これは、逆行列や行列式の計算に特別なアルゴリズムを使用しないためであった。Rumaleのver. 0.13系では、Numo::Linalgへの対応を進めており、逆行列や行列式の計算が可能になったため、共分散行列の全要素を用いる方法を実装した。

rumale | RubyGems.org | your community gem host

使い方

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

$ gem install rumale

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

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

Numo::Linalgもインストールする。Numo::Linalgがロードされていなければ、共分散行列の全要素を使ったバージョンは使用できない。

$ brew install openblas
$ gem install numo-linalg

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

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

# 楕円の大きさを求める
def ellipse_size_full(covar)
  v, = Numo::Linalg.eigh(covar)
  nv = 2 * Math.sqrt(2) * Numo::NMath.sqrt(v)
  nv.to_a
end

# 楕円の回転角を求める
def ellipse_angle_full(covar)
  _v, w = Numo::Linalg.eigh(covar)
  u = w[0, true] / Numo::Linalg.norm(w[0, true])
  angle = Numo::NMath.atan2(u[1], u[0])
  (180.0 * angle[0] / Math::PI).to_i
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でクラスタリングする.
# covariance_typeで共分散行列の要素をどう使うかを指定できる.
gmm = Rumale::Clustering::GaussianMixture.new(n_clusters: 3, covariance_type: 'full', 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_full(gmm.covariances[2, 0..1, 0..1]),
      angle: ellipse_angle_full(gmm.covariances[2, 0..1, 0..1 ]),
      front: true, fs: true, empty: true, bo: 1)
  set(:object, 2, 'ellipse',
      center: [gmm.means[0, 0], gmm.means[0, 1]],
      size: ellipse_size_full(gmm.covariances[0, 0..1, 0..1]),
      angle: ellipse_angle_full(gmm.covariances[0, 0..1, 0..1]),
      front: true, fs: true, empty: true, bo: 2)
  set(:object, 3, 'ellipse',
      center: [gmm.means[1, 0], gmm.means[1, 1]],
      size: ellipse_size_full(gmm.covariances[1, 0..1, 0..1]),
      angle: ellipse_angle_full(gmm.covariances[1, 0..1, 0..1]),
      front: true, fs: true, empty: true, bo: 3)
  plot(*plots)
end

これを実行すると、以下のような画像を得られる。以前の対角要素のみを使う方法（covariance_type: 'diag'）と比較して、楕円の角度が考慮されており、よりクラスタの形状を捉えることができている。

共分散行列の全要素を使用する場合

共分散行列の対角要素のみを使用する場合

共分散行列の全要素を使用する場合、各クラスタの共分散行列をすべて保持する必要があるためメモリをより多く使用する。また、逆行列や行列式の計算のため実行時間も遅くなる。

おわりに

その他、因子分析（Factor Analysis）を追加したり、スペクトルクラスタリング（Spectral Clustering）を追加したり、HDBSCANを追加したりしている。これで、ひととおり、Rumaleに実装したアルゴリズムで、Numo::Linalgを使用してパワーアップできるものは対応できた気がする。

ちなみに、会社ではマネ〜ジャ〜兼エンジニアで、タスクフォースな仕事が忙しくなりそうで、これから平日のOSS開発は難しくなる。これまで、月に4回ほどのペースでRumaleをバージョンアップしてきたが、月に1〜2回とかに絞ろうと思う。Rumaleをはじめ、Rubyでデータ分析・機械学習なライブラリの開発が、昼間の仕事としてできれば良いのだけど。幸いベーシックなアルゴリズムは実装できているので、ガウス過程か、計量学習か、あらたにModuleを切るようなものを、じっくり作ろうかな〜といったところ。

github.com

Rumaleにカーネルリッジ回帰を追加した

はじめに

Rumaleのversion 0.13系では、固有値分解や逆行列計算といった、Numo::Linalgにある線形代数のテクニックを利用する機械学習アルゴリズムの実装を進めている。version 0.13.3では、カーネルリッジ回帰とカーネル主成分分析を追加した。

rumale | RubyGems.org | your community gem host

使い方

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

$ gem install rumale red-datasets-numo-narray

そして、Numo::Linalgをインストールする。Numo::LinalgはBLAS/LAPACK系のライブラリに依存する。別途OpenBLASなどをインストールする必要がある。

$ gem install numo-linalg

Boston Housingデータセットを利用して、カーネルリッジ回帰を試してみる。Boston Housingデータセットは、地域の住宅価格を、犯罪発生率や平均部屋数などから推定するタスクである。リッジ回帰・カーネルリッジ回帰で住宅価格を推定し、決定係数により評価を行う。ハイパーパラメータはそれらしい値をいれた（※カーネルリッジ回帰は正則化パラメータを調整したほうがよい）。

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

# データセットを読み込む.
datasets = Datasets::LIBSVM.new('housing').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.2, 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]

# 特徴量を正規化する.
t = Rumale::Preprocessing::MinMaxScaler.new
train_s = t.fit_transform(train_s)
test_s = t.transform(test_s)

# リッジ回帰を学習し決定係数で評価する.
ridge = Rumale::LinearModel::Ridge.new(reg_param: 0.1, solver: 'svd')
ridge.fit(train_s, train_v)
score = ridge.score(test_s, test_v)
puts "linear ridge: #{score.round(4)}"

# カーネルリッジ回帰を学習し決定係数で評価する.
kridge = Rumale::KernelMachine::KernelRidge.new(reg_param: 0.1)
## 訓練データ間のRBFカーネル関数の値をもとめ訓練を行う.
gamma = 2.0
train_kernel_mat = Rumale::PairwiseMetric::rbf_kernel(train_s, nil, gamma)
kridge.fit(train_kernel_mat, train_v)
## テストデータと訓練データ間のRBFカーネル関数の値をもとめ決定係数を計算する.
test_kernel_mat = Rumale::PairwiseMetric.rbf_kernel(test_s, train_s, gamma)
score = kridge.score(test_kernel_mat, test_v)
puts "kernel ridge: #{score.round(4)}"

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

linear ridge: 0.6443
kernel ridge: 0.8655

決定係数の値が大きいほどよい。カーネルリッジ回帰のほうが、Boston Housingデータセットに関しては、回帰の推定性能がリッジ回帰よりも優れていると思われる（上記スクリプトは、シンプルなホールドアウト検証になってるので、本当は交差検証をちゃんとしたほうが良い）。

おわりに

Rumaleにカーネルリッジ回帰を実装した。カーネル法は、データ数が十分にない場合などでは、まだまだ有効であると思われる。また、ニューラルネットワークがReLUやDropoutなどの学習方法を発見して復活したように、とんでもないカーネル関数やすごいマルチカーネル手法などが発見されると、復活の可能性はあるかもしれない??

Rumaleでは、その他、カーネル主成分分析やissueでリクエストのあったShared Nearest Neighborクラスタリングなどを追加している。しばらくは、Numo::Linalgありきのアルゴリズムを追加していく。そろそろRumaleのユーザーガイドやチュートリアル的な文書が必要だな〜とかも考えている。

github.com

Rumaleの主成分分析でNumo::Linalgの固有値分解を利用できるようにした

はじめに

機械学習アルゴリズムでは、固有値分解が逆行列計算といった、線形代数のテクニックを利用する場合がある。これらをRubyで使用するには、Numo::Linalgが適している。一方で、Numo::Linalgは、Numo::NArrayと異なり、BLAS/LAPACK系の外部ライブラリに依存するため、Rumaleでは依存することを避けてきた。ただ、今後のRumaleの拡張を考えると必要となる。そこで、以前にParallelを導入したときのように、別でNumo::Linalgをrequireしている場合にのみ、機能が有効になるよう実装することにした。

まず、試験的に、主成分分析で固有値分解を使うものを実装してみた。主成分分析は、主成分ベクトルを、共分散行列の固有値分解により求める方法がメジャーである。他にも主成分ベクトルを求めるアルゴリズムはあり、Rumaleでは不動点法によるアルゴリズムを実装していた。今回、solverオプションを追加して、固有値分解によるアルゴリズムを選択できるようにした。これを version 0.13.0 としてリリースした。

rumale | RubyGems.org | your community gem host

使い方

Rumaleは、gemコマンドでインストールできる。データセットの読み込みでred-datasetsを、主成分分析したデータの可視化にNumo::Gnuplotを使いたいので、あわせてインストールする。

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

そして、Numo::Linalgをインストールする。

$ gem install numo-linalg

USPSという手書き数字画像のデータセットを、主成分分析で二次元に射影して可視化する。

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

# Numo::Linalgを明示的にrequireしない限りは有効にならない.
require 'numo/linalg/autoloader' 

# Numo::NArray形式でUSPSデータセットを読み込む.
dataset = Datasets::LIBSVM.new('usps').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(1000)
labels = labels[rand_id]
samples = samples[rand_id, true]

# 共分散行列を固有値分解することによる主成分分析で、二次元空間にデータを射影する.
# solverオプションに 'evd' を指定する（指定しなければコレまでどおり不動点法を使う）.
pca = Rumale::Decomposition::PCA.new(n_components: 2, solver: 'evd', random_seed: 1)
low_samples = pca.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: 'pca.png')
  plot(*plots)
end

以上のスクリプトを実行した結果は次のようになる。USPSを主成分分析で可視化するとだいたいこんな感じになる。上手く動いていることがわかる。

おわりに

Numo::Linalgを導入したので、カーネル主成分分析やカーネルRidge回帰など、固有値分解や逆行列計算を必要とする機械学習アルゴリズムを実装していこうと思っている。

github.com

Numo::Pocketfftにフーリエ変換による畳み込み演算するメソッドを追加した

はじめに

畳み込み演算は、そのまま実装すると、データが大きくなると重くなる。一方、フーリエ変換により、畳み込み演算は単純な掛け算に変換される。これを応用して、畳み込み演算したい二つの配列をフーリエ変換し、乗算を行った後に、フーリエ逆変換する高速化手法が一般に利用される。これをNumo::Pocketfftに実装しversion 0.2.0とした。

numo-pocketfft | RubyGems.org | your community gem host

使い方

Numo::Pocketfftはgemコマンドでインストールできる。pocketfftのコードを同梱しているので、外部ライブラリを別でインストールする必要はない。

$ gem install numo-pocketfft

SciPyにあやかって、メソッド名はfftconvolveとした。二つのNumo::NArrayな配列を渡せば、その（離散）畳み込み演算の結果を返す。

irb(main):001:0> require 'numo/pocketfft'
=> true
irb(main):002:0> Numo::Pocketfft.fftconvolve(Numo::DFloat[1, 2, 3], Numo::DFloat[4, 5])
=> Numo::DFloat#shape=[4]
[4, 13, 22, 15]

画像処理への応用

画像のフィルタリングは、画像とフィルタの畳み込みで実現される。Magroも使って、アイコンに使ってるハンバーグの画像にSobelフィルタをかけて、エッジを抽出してみる。

require 'magro'
require 'numo/pocketfft'

# Lena画像を読み込む.
img = Magro::IO.imread('lena.png')

# グレイスケール化する.
gray = img.median(axis: 2)

# Sobelフィルタを定義する.
fx = Numo::DFloat[[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]
fy = Numo::DFloat[[-1, -2, -1], [0, 0, 0], [1, 2, 1]]

# 畳み込みによりSobelフィルタをかける.
gx = Numo::Pocketfft.fftconvolve(gray, fx)
gy = Numo::Pocketfft.fftconvolve(gray, fy)

# エッジ画像を得る.
g = Numo::NMath.sqrt(gx**2 + gy**2)
g[g>255] = 255
g[g<0] = 0

# 畳み込みによりできた画像の端を削る.
# ※
# 元の画像サイズは512x512であり, フィルタのサイズが3x3であることから, 
# 畳み込みにより (512 + 3 - 1) x (512 + 3 - 1) の画像ができる.
# 上下左右の端を1ピクセルずつ削り512x512の画像にする.
h, w = gray.shape
fh, fw = g.shape
sy = (fh - h) / 2
sx = (fw - w) / 2
g = Numo::UInt8.cast(g[sy..h, sx..w])

# 画像を保存する.
Magro::IO.imsave('lena_fil.png', g)

これを実行すると、以下のような画像が得られる。Sobelフィルタによりエッジが抽出されていることがわかる。

おわりに

Numo::NArrayな配列で畳み込み演算ができるようになった。使い方の例のように、画像のフィルタリングも簡単になったので、Magroの開発も進むかな？

github.com