ruby 2.6.0p0 (2018-12-25 revision 66547) [x86_64-linux]で下記コードを実行し、速度を検証しました。
# frozen_string_literal: true
require 'benchmark'
require 'parallel'
require 'numo/narray'
require 'nmatrix'
list_size = 1_000_000
range = (-100.0..100.0)
out = Array.new(list_size) do
rand(range) + rand(range) * 1i
end
resiSiPower = Array.new(list_size) do
rand(range)**2
end
# Numo::NArray
out_na = Numo::DComplex.new(list_size)
out_na.store(out)
resiSiPower_na = Numo::DFloat.new(list_size)
resiSiPower_na.store(resiSiPower)
# NMatrix
out_nm = N[out, dtype: :complex128]
resiSiPower_nm = N[resiSiPower, dtype: :float64]
Benchmark.bm(23) do |x|
x.report('orginal') do
(0...out.length).each do |j|
resiSiPower[j] += out[j].abs**2
end
end
x.report('r^2 + i^2') do
(0...out.length).each do |j|
resiSiPower[j] += out[j].real**2 + out[j].imag**2
end
end
x.report('abs2') do
(0...out.length).each do |j|
resiSiPower[j] += out[j].abs2
end
end
x.report('times abs2') do
out.length.times do |j|
resiSiPower[j] += out[j].abs2
end
end
x.report('upto abs2') do
0.upto(out.length - 1) do |j|
resiSiPower[j] += out[j].abs2
end
end
x.report('zip map abs2 map sum') do
resiSiPower = resiSiPower.zip(out.map(&:abs2)).map(&:sum)
end
x.report('map abs2 zip map sum') do
resiSiPower = out.map(&:abs2).zip(resiSiPower).map(&:sum)
end
x.report('each_with_index map') do
resiSiPower = resiSiPower.each_with_index.map do |r, i|
r + out[i].abs2
end
end
x.report('map! with i') do
i = 0
resiSiPower.map! do |r|
r + out[i].abs2
i += 1
end
end
x.report('map with each') do
out_e = out.each
resiSiPower.map! do |r|
r + out_e.next.abs2
end
end
x.report('parallel') do
resiSiPower = Parallel.map(0...out.length, in_threads: 4) do |j|
resiSiPower[j] + out[j].abs2
end
end
x.report('Numo::NArray') do
resiSiPower_na += out_na.abs**2
end
x.report('NMatrix') do
resiSiPower_nm += out_nm.abs**2
end
end
実行前にgem install parallel numo-narray nmatrix
として必要なライブラリをインストールしておいてください。
out
を浮動小数点数複素数の配列、resiSiPower
を浮動小数点数の配列という前提にしています。型が異なる場合は、違いが出るかも知れませんので、それに合わせて修正してください。また、OSやRubyのバージョンによる違いも出る場合があります。
絶対値の2乗の求め方
絶対値の2乗の場合はComplex#abs2が利用できます。と言っても、x.real**2 + x.imag**2
と速度はほとんど変わりませんでした。どちらか好みを使うときでしょう。
ループの回し方
(0...n).each
はn.times
や0.upto(n-1)
とも書けます。これらも速度はほとんど変わりませんでしたので、好みを使うときでしょう。
代入ではなくArray#map!を使う方法もありますが、インデックスを別途増加させる必要があります。こちらも(0...n).each
とほとんど変わりませんでした。インデックスではなくeach
でEnumeratorを作る方法は逆に遅くなりました。
その他、関数型プログラミング風にzip
やmap
を組み合わせて使う方法は配列の生成があるため遅くなります。
並列化
CRubyにはGVLがあるため単純なマルチスレッドにしても遅くなるだけでした。ただ、JRubyではマルチスレッドが並列化するため、速くなる可能性はあります(未検証)。
マルチプロセスはそれぞれのプロセスをまとめる処理があるせいか、遅すぎて計測できませんでした。ただ、配列を適当な数に分割してマルチプロセスで処理した場合は早くなる可能性はあります(未検証)。
行列ライブラリの使用
Pythonでいうnumpyのようなライブラリとして、Numo::NArrayとNMatrixがあります。これらはマトリックス(ベクトル含む)をマトリックスのまま計算することが出来ます。内部処理はCレベルで実装されており、それぞれの要素に対して同じ処理を高速に行うことが出来ます。
Numo::NArrayは大幅に速度が向上しました。高速を目指すのであれば、Numo::NArrayを使うのが一番いいかもしれません。
NMatrixは逆に遅くなりました。何故かはわかりません。(私のコードの書き方が悪いのでしょうか?)
以上の事から、今のからあまり変えずに少しでも速くするならabs2
を使う。さらに速度の向上を目指すならNumo::NArrayを使うというのがいいのではないでしょうか。
out[j]
は複素数(complex number)でしょうか?(out[j].abs)**2
は絶対値の2乗なので、sqrt(out[j].real ** 2 + out[j].imag ** 2) ** 2
ということになります。つまり、out[j].real ** 2 + out[j].imag ** 2
に書き替えることができます(高速化に寄与するかどうかは分かりませんが…)。