Rubyで連立一次方程式

Rubyで連立一次方程式を解く方法を２つ紹介する。

Ruby標準ライブラリ Matrix を使う

Rubyの標準ライブラリにあるMatrixクラスを使う方法。標準ライブラリにあるので追加でgemをインストールすることなく使える方法。
（参考: Ruby2.2 ではアレが死ぬほど使いやすくなるの！ )

require 'matrix'
m = Matrix[ [1.0, 2.0], [1.0, 3.0] ]
# => Matrix[[1.0, 2.0], [1.0, 3.0]]
b = Vector[1.0,0.0]
# => Vector[1.0, 0.0]
m.lup.solve(b)
# => Vector[3.0, -1.0]

Nmatrix gemを使う

gem install nmatrix でインストールする。

require 'nmatrix'
m = N[ [1.0, 2.0], [1.0, 3.0] ]    # 初期化
# =>
# [
#  [1.0, 2.0]
#  [1.0, 3.0]
# ]
b = N[[1.0],[0.0]]
pp m.solve(b)      # => [[ 3.0],[-1.0]]

NMatrixの行列要素へのアクセスの方法

m = N[ [1.0, 2.0], [1.0, 3.0] ]
m[0,1]   # => 2.0
m[0..-1,1]   # => [ [2.0], [3.0] ]  column vector
m[0,0..-1]   # => [ [1.0, 2.0] ]    row vector
m.shape      # => [2, 2]

速度比較

一般的にnmatrixの方がより数値計算に特化しているので速い。nmatrixの方は1000x1000くらいの行列なら一瞬で解ける。

require 'benchmark'
require 'matrix'
reuqire 'nmatrix'

m = Matrix.build(500) { rand }    # 100x100 random matrix
a = 500.times.map { rand }
b = Vector[*a]                # 100x1 random vector
m.lup.solve(b)

a2 = NMatrix.random([1000, 1000], dtype: :float64)
b2 = NMatrix.random([1000, 1], dtype: :float64)

Benchmark.benchmark(Benchmark::CAPTION) do |x|
  x.report("matrix:") { m.lup.solve(b) }
  x.report("nmatrix:") {  a2.solve(b2) }
end

#               user     system      total        real
#  matrix:  2.897694   0.002253   2.899947 (  2.901112)
# nmatrix:  0.129778   0.001833   0.131611 (  0.131675)