洋食の日記

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

Numo::LinalgのサブセットのNumo::TinyLinalgを作った

はじめに

Numo::NArrayで逆行列計算や固有値分解などの線形代数計算を担うNumo::LinalgのサブセットであるNumo::TinyLinalgを少しずつ作っていた。どのようなサブセットかというと、機械学習で使う(Rumaleで使う)メソッドだけを実装したものになる。

github.com

インストール

Numo::Linalgは、BLAS/LAPACKのためのライブラリを選択できる方式をとっているが、Numo::TinyLinalgではOpenBLASだけに対応している。必要に応じて、事前にOpenBLASをインストールする。

$ brew install openblas

その上で、OpenBLASのインストールディレクトリを指定してインストールする(/usr/libとか一般的なディレクトリにインストールしてあるならインストールディレクトリの指定は必要ない)。

$ gem install numo-tiny_linalg -- --with-opt-dir=/opt/homebrew/Cellar/openblas/0.3.23/

インストールでnative extensionsのビルド時に、OpenBLASが見つからない場合は、OpenBLASをダウンロードしてビルドするようにも作ってあるので、なにもせずgem installでもよい。

$ gem install numo-tiny_linalg

使い方

サブセットなのでNumo::Linalgと同様になる。Numo::NArrayは、Numo::Linalgがロードされている場合、実数配列同士のdotメソッドでNumo::Linalgのdotメソッドを呼ぶようにできている。これを再現するために、Numo::TinyLinalgでもdotメソッドを実装してある。なので、以下のように、Numo::Linalgがロードされていない場合に、Numo::TinyLinalgで置き換えることができる。

require 'numo/tiny_linalg'

Numo::Linalg = Numo::TinyLinalg unless defined?(Numo::Linalg)

あとはNumo::Linalgと同様となる。例えば、対称行列の固有値分解は以下のようになる。

x = Numo::DFloat.new(5, 3).rand - 0.5
c = x.dot(x.transpose)

vals, vecs = Numo::TinyLinalg.eigh(c, vals_range: [2, 4])

pp vals
# Numo::DFloat#shape=[3]
# [0.118795, 0.434252, 0.903245]

pp vecs
# Numo::DFloat#shape=[5,3]
# [[0.154178, 0.60661, -0.382961],
#  [-0.349761, -0.141726, -0.513178],
#  [0.739633, -0.468202, 0.105933],
#  [0.0519655, -0.471436, -0.701507],
#  [-0.551488, -0.412883, 0.294371]]

pp (c - vecs.dot(vals.diag).dot(vecs.transpose)).abs.max
# 3.3306690738754696e-16

おわりに

開発のきっかけはNumo::Linalgの開発が長らく停止しているように見えることにある。OSSなので停止していることは全く問題ない(私自身が腱鞘炎やら多忙やらで2年くらい大した活動ができなくなったことがあり、開発者に様々な事情があるのは理解できるし、OSSのために人生を捧げる必要はないと考える)。ただ、今後のRumaleの拡張を考えると、Numo::Linalgの発展が必要になったので、サブセットを作ることにした。forkすることも考えたが、Numo::Linalgは(Numo::NArrayも)erbによりC言語でGeneric programmingを実現していて、これ自体は大発明なのだが、常にコードに触れていないと概念というか開発のコツを忘れてしまう。(再びヒドい腱鞘炎になるとか)開発が断続的になってもスムーズに再開できるように、C++でサブセットを作ることにした。

今後はRumaleなど、機械学習ライブラリの開発に戻り、また必要に応じて拡張したいと思う。おそらく、ガウス過程回帰を実装するとかで、choloeskyを追加すると思う。