再帰級数のMemoization(メモ化)に関する覚書

Memoization (メモ化)

Fibonacci級数を例にMemoization(キャッシュ，メモ化)の覚書．
Mathematica(Wolfram言語)，python, juliaで実装比較してみる．
次回あたりにHermite多項式を例にMemoizationについて考えてみたい．

Memorizationかと思っていたがrは不要で，Memoizationと呼ぶらしい．

Fibonacci 級数

みんな大好きFibonacci級数

\begin{align}
&F_0 = 0,\\ 
&F_1 = 1,\\
&F_n = F_{n-1} + F_{n-2} \quad (n>1)
\end{align}

を素朴にFibonacci級数で$n=40$程度の項を評価してみる．Memoization(キャッシュ)をいくつか実装してみる．

python

pythonで素朴に書くと次の様だろうか．

fib.py

import time

def fib(n):
    if n < 2:
        return n
    else:
        return fib(n-1) + fib(n-2)

a = time.time()
n=40
res = fib(n)
print(time.time() -a, "[sec]", "f(%d)=%d" % (n, res))

手元のlaptopで

$ fmt="\nreal:%e[sec]\nuser:%U[sec]\nsys:%S[sec]\nMemory:%M[KB]"
$ /usr/bin/time -f ${fmt} python fib.py
53.77412390708923 [sec] f(40)=102334155

real:53.83[sec]
user:52.57[sec]
sys:0.04[sec]
Memory:6480[KB]

dictionaryを使って計算結果をキャッシュすればより早くなる

fib.py

import time

fib_d={0:0, 1:1}
def fib(n):
    if n in fib_d:
        return fib_d[n]
    else:
        res  = fib(n-1) + fib(n-2)
        fib_d[n] = res
        return res

a = time.time()
n=40
res = fib(n)
print(time.time() -a, "[sec]", "f(%d)=%d" % (n, res))

$ /usr/bin/time -f ${fmt} python fib.py
2.09808349609375e-05 [sec] f(40)=102334155

real:0.05[sec]
user:0.03[sec]
sys:0.03[sec]
Memory:6488[KB]
$ /usr/bin/time -f ${fmt} python fib.py

${F_n}$の計算結果を保存しているので少しメモリ使用量が上がる．同等のキャッシュは標準ライブラリーfunctools lru_cacheで低コストで実装できる．

fib.py

from functools import lru_cache
import time

@lru_cache(maxsize=1000)
def fib(n):
    if n < 2:
        return n
    else:
        return fib(n-1) + fib(n-2)

a = time.time()
n=40
res = fib(n)
print(time.time() -a, "[sec]", "f(%d)=%d" % (n, res))

$ /usr/bin/time -f ${fmt} python fib.py
2.0742416381835938e-05 [sec] f(40)=102334155

real:0.10[sec]
user:0.01[sec]
sys:0.03[sec]
Memory:6500[KB]

Mathematica(Wolfram言語)

素朴な実装は以下のようだろうか

fib.wl

F[n_]:=F[n-2] + F[n-1]
F[0] = 0
F[1] = 1

Print[Timing@F[40]]

$ /usr/bin/time -f ${fmt} wolframscript -script fib.wl
{336.265625, 102334155}

real:346.60[sec]
user:336.37[sec]
sys:1.28[sec]
Memory:114468[KB]

wolfram言語で$:=$は遅延評価(定義)で$=$は即時評価．組み合わせて
https://reference.wolfram.com/language/tutorial/FunctionsThatRememberValuesTheyHaveFound.html

fib.wl

F[n_]:=F[n]=F[n-2] + F[n-1]
F[0] = 0
F[1] = 1

Print[Timing@F[40]]

とすれば，計算結果をキャッシュできる

$ /usr/bin/time -f ${fmt} wolframscript -script fib.wl
{0., 102334155}

real:1.20[sec]
user:0.56[sec]
sys:0.71[sec]
Memory:114564[KB]

wolfram言語ではオーバーヘッドが結構あるのか，関数評価の時間と実測の時間で大きくずれる．また，メモリー使用量もpythonと比較して大きい傾向にあるようだ．

julia

素朴に書くと以下のようだろうか

fib.jl

function fib(n)
    if n<2
        return n
    else
        return fib(n-1) + fib(n-2)
    end
end

b = @time fib(40)
@show b

$ /usr/bin/time -f ${fmt} julia fib.jl
  0.783497 seconds (2.56 k allocations: 235.453 KiB)
b = 102334155

real:1.19[sec]
user:1.25[sec]
sys:0.32[sec]
Memory:125392[KB]

素朴な実装でもめっちゃ早い．pythonと同じように辞書(dictionary)を使ってキャッシュすると

fib.jl

known = Dict(0=>0, 1=>1)
function fibonacci(n)
    if n ∈ keys(known)
        return known[n]
    end
    res = fibonacci(n-1) + fibonacci(n-2)
    known[n] = res
    res
end

b = @time fibonacci(40)
@show b

$ /usr/bin/time -f ${fmt} julia fib.jl
  0.047068 seconds (7.81 k allocations: 508.287 KiB)
b = 102334155

real:0.81[sec]
user:0.43[sec]
sys:0.65[sec]
Memory:127560[KB]

さらに早い．標準ライブラリーではないが，Julia Memoizeをinstall とprecompileして使うと

fib.jl

using Memoize
@memoize function fib(n)
    if n<2
        return n
    else
        return fib(n-1) + fib(n-2)
    end
end

b = @time fibonacci(40)
@show b

$ /usr/bin/time -f ${fmt} julia fib.jl
  0.022669 seconds (33.39 k allocations: 1.714 MiB)
b = 102334155

real:2.81[sec]
user:2.70[sec]
sys:0.53[sec]
Memory:188900[KB]

と関数評価にかかる時間はさらに早くなる．一方パッケージをロードしてコンパイルする時間がかかるのか，プログラムの実行時間は少し長くなるようだ．