kenju's blog

About Programming and Mathematics

Introduction to Matrix API in Ruby for Matrix Operation

広告配信サーバーを書いていると、最適化問題によくぶち当たる。

Optimization problem - Wikipedia

最適化問題は、純粋な数理モデルに落とし込んで考えることができる。広告は、業務上のタスクを、純粋な数学上の問題に置き換えてコードを書く珍しいドメインの一つの気がする。

数理モデルを解くとき、大抵の場合行列演算に落とし込んで考えると綺麗にコードが書けることが多い。 最適化問題に限らず、単純にループを書くと O(n^2) にも O(n^3) にもなるオーダーが、ギュッとまとまって書けることが多い(計算量が減るかどうかは行列演算 API の内部実装によるけれど)

Ruby では標準 API に Matrix class があり、基本的なメソッドは用意されている。

ruby-doc.org

"400x faster Matrix multiplication for Ruby" で検証されているように、SciRuby で開発されている NMatrix API のほうがパフォーマンスが良い場合が多い。要件上パフォーマンスを達成できない場合は、標準 API 以外も見て見る必要はある。

github.com

個人の感想としては、API は直感的で、軽く試すのには適している。column/row ごとの enumerator が無いのが少し驚いたが、標準メソッドを組み合わせば以下のように書ける。

x = Matrix[
  [100, 100, 100, 100, 100],
  [200, 200, 200, 200, 200],
  [300, 300, 300, 300, 300],
  [400, 400, 400, 400, 400],
  [500, 500, 500, 500, 500],
]

puts "total value of x"
puts x.reduce(0) {|acc, e| acc + e } #=> 7500

puts "total value per columns of x"
puts x.column_count.times.map {|n| x.column(n).sum } #=> [1500, 1500, 1500, 1500, 1500]

puts "total value per rows of x"
puts x.row_count.times.map {|n| x.row(n).sum } #=> [500, 1000, 1500, 2000, 2500]

独自のデータ構造を定義して計算をラッピングするときに、もし行列演算に置き換えて考えられるなら、Matrix API を wrap するときれいに設計できることもある。

例えば、行列に対してインデックス(ラベル)を持つ独自のデータ構造を以下のように定義してみる。エクセルの Worksheet のようなイメージ。ラベル別に行や列も検索できるとする。

require 'matrix'
require 'forwardable'

class MatrixWithIndex
  include Enumerable

  extend Forwardable
  def_delegators :@matrix, :each

  def initialize(matrix, column_index, row_index)
    @matrix = matrix
    @column_index = column_index
    @row_index = row_index
  end

  def column_by(val)
    val_index = @column_index.find_index(val)
    @matrix.column(val_index)
  end

  def row_by(val)
    val_index = @row_index.find_index(val)
    @matrix.row(val_index)
  end

  def pretty_print
    lines = []
    lines << @row_index.unshift('*').join(' | ')
    @matrix.column_count.times.each_with_index.to_a.map {|c, i|
      lines << "#{@column_index[i]} | #{@matrix.column(c).to_a.join(' | ')}"
    }
    puts lines.join("\n")
  end
end

x = Matrix[
  [100, 100, 100, 100, 100, 100, 100],
  [200, 200, 200, 200, 200, 200, 200],
  [300, 300, 300, 300, 300, 300, 300],
  [400, 400, 400, 400, 400, 400, 400],
  [500, 500, 500, 500, 500, 500, 500],
]
column_index = %w(Monday Tuesday Wednesday Thursday Friyday Saturday Sunday)
row_index = %w(USD JPY HKD EUR CNY)

mwi = MatrixWithIndex.new(x, column_index, row_index)
mwi.pretty_print
#=> * | USD | JPY | HKD | EUR | CNY
#=> Monday | 100 | 200 | 300 | 400 | 500
#=> Tuesday | 100 | 200 | 300 | 400 | 500
#=> Wednesday | 100 | 200 | 300 | 400 | 500
#=> Thursday | 100 | 200 | 300 | 400 | 500
#=> Friyday | 100 | 200 | 300 | 400 | 500
#=> Saturday | 100 | 200 | 300 | 400 | 500
#=> Sunday | 100 | 200 | 300 | 400 | 500

mwi.sum #=> 10500

mwi.column_by("Monday") #=> Vector[100, 200, 300, 400, 500]
mwi.column_by("Monday").sum #=> 1500

mwi.row_by("JPY") #=> Vector[200, 200, 200, 200, 200, 200, 200]
mwi.row_by("JPY").sum #=> 1400

行列演算ベースで考えると設計の幅も実装の選択肢も広がるので、覚えておいて損はない。