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

はじめに

前回Fibonacci級数を例題にMemoizationを考えた．今回はHermite多項式を例題にMemoizationの実装と言語による実行速度について考える．扱う言語はpython,wolfram言語(Mathematica)とjulia(準備中)．

先に断っておくと，Hermite多項式のような計算機科学において超基本的な直交多項式系の評価は，自身で実装しなくともパッケージとして提供されているのでそれを利用したほうが圧倒的によい．

pythonでは

• numpy.polynomial
• scipy.special
• mpmath (直交多項式系)
• sympy.Special

Wolfram言語では

• HermiteH

実はjuliaではよくわからんというのが本音(2020年)．実際に使用していないが調べた限りでは，

• PolyChaos
• ApproxFun

で直交多項式関係のパッケージが利用できそう．特殊関数に関しては

• SpecialFunctions

が利用できそう．juliaではnumpyやscipyほどのpythonにおける計算機科学のデファクトスタンダードとなっているものがあるのか理解していないが，少なくともGithubを見るとかなりアクティブなので将来が楽しみ．

動機

前回Fibonacci級数，Fibonacci級数では一度評価した項をキャッシュ(memoization)し，無駄な重複評価をなくすことで再帰級数の計算速度を改善した．同様にHermite多項式を例題に再帰定義関数のMemoizationを考える．

https://mathworld.wolfram.com/HermitePolynomial.htmlにあるように，Hermite多項式は様々な方法で定義することが可能だが，ここでは再帰関数として比較的簡単に実装可能な次の定義を評価する．

$$
H_0(x)= 1, \\
H_1(x)= 2x, \\
H_n(x) = 2xH_{n-1}(x) -2(n-1)H_{n-2}(x) \quad (n >1)\\
$$
wolfram言語でn=10までsymbolicに評価すると

H[0, x_] = 1;
H[1, x_] = 2x;
H[n_, x_] := H[n, x] = 2x*H[n-1, x] - 2(n-1)H[n-2, x];

res = Table["H_"<>ToString[n]<>"(x)="<>ToString@TeXForm@Expand@H[n, x], {n, 0, 10}]//MatrixForm;
Print@res

$$
\begin{array}{l}
H_0(x)=1 \\
H_1(x)=2 x \\
H_2(x)=4 x^2-2 \\
H_3(x)=8 x^3-12 x \\
H_4(x)=16 x^4-48 x^2 + 12\\
H_5(x)=32 x^5-160 x^3+120 x \\
H_6(x)=64 x^6-480 x^4+720 x^2-120\\
H_7(x)=128 x^7-1344 x^5+3360 x^3-1680 x \\
H_8(x)=256 x^8-3584 x^6+13440 x^4-13440 x^2+1680\\
H_9(x)=512 x^9-9216 x^7+48384 x^5-80640 x^3+30240 x\\
H_{10}(x)=1024 x^{10}-23040 x^8+161280 x^6-403200 x^4+302400 x^2-30240\\
\end{array}
$$

すぐにわかるが，n次のHermite多項式は$H_n(x)$は$x$に依存するため，前回Fibonacci級数のようなMemoizationでは意味がない．
また，偶数次の原点での値は$|H_n(0)|\sim 2(n-1)!!!!$ (4つ飛ばしの階乗，４度ビックリ!!!!)で増えていることがわかるので，$\exp(-x^2)$を掛けて直交関係を数値的に評価しようとすると，(n-1)!!!!はn!よりは遅いとはいえ(おそらく)指数関数より早く増大するので，$n$増大は指数関数より大きな値と指数関数的に小さな値の計算となり精度が失われる(と思う)．故に多倍長評価でHermite多項式の評価をしたいと言うのが動機．

pythonによる素朴実装

referenceとして素朴なHermite多項式をpython(numpy, mpmath)で評価すると次のようだろうか．

import numpy as np
import time

def H(n,x):
    if n==0:
        return np.array([1]*len(x))
    elif n==1:
        return 2*x
    else:
        return 2*x*H(n-1, x) - 2*(n-1)* H(n-2,x)


time0 = time.time()
x = np.linspace(-1, 1, 100)
n=30
y = H(n,x)
y0 = y[len(x)//2]
print("numpy:", time.time()-time0, "sec", "H_%d(0)=%.18e" % (n, y0))

n=20
time0 = time.time()
x = np.linspace(-1, 1, 100)
y = H(n,x)
y0 =y[len(x)//2]
print("numpy:", time.time()-time0, "sec", "H_%d(0)=%.18e" % (n, y0))

import mpmath
mpmath.mp.dps = 100
time0 = time.time()
x = np.array(mpmath.linspace(-1, 1, 100), dtype="object")
y = H(n, x)
y0 = float(y[len(x)//2])
print("mpmath:", time.time()-time0, "sec", "H_%d(0)=%.18e" % (n,y0))

$ python hermite-poly.py 
numpy: 10.97638750076294 sec H_30(0)=-2.022226152780266537e+20
numpy: 0.09095430374145508 sec H_20(0)=6.690748809755917969e+11
mpmath: 6.780889987945557 sec H_20(0)=6.690748809755917969e+11

Memoizationの方針

$H_n(x)$は次数$n$の$x$多項式なので，$x$を未評価で，係数だけを評価すれば良いと気づく．lambda式を使えば$x$を未評価のままキャッシュすることが出来きるので，計算速度が改善される・・・と考えたがこれは上述素朴実装と効率は変わらない．なぜなら，lambda式で，計算演算の順序がキャッシュされるだけで，計算結果は未評価のままなだからである(具体例はmpmath.polynomial後述する)．

さてlambda式のmemoizationでは計算速度の改善がもたらされない．故にlamda式を使わずに$H_n(x)$の係数だけを評価する．
つまり，$m$次の多項式

$$
p_m(x) = a_0 + a_1x + a_2x^2 + \cdots + a_m x^m
$$
が与えられたとする．長さ$m$の係数リストを

$$
L(p_m(x)) = [a_0, a_1, \ldots, a_m]\
L(q_m(x)) = [b_0, b_1, \ldots, b_m]
$$

とすると$p_m(x) \pm q_m(x)$は

$$
L(p_m(x)) \pm L(q_m(x)) = [a_0\pm b_0,\ldots, a_m\pm b_m]
$$

乗法$x\times p_m$は

$$
L(x) \times L(p_m(x)) = [0, a_0, a_1, a_2, \cdots, a_m]
$$

として長さ$m+1$のリストを作る相当する．この演算の操作によって，$p_m(x)$ の加減乗(除)，積分，微分を実現することができるがHermite多項式の評価にはこれで十分．除法がかっこ付きなのは特別な$p_m$を除き，一般に演算を厳密に評価することができないためである．係数リストL(p_m(x))であれば具体的に評価した値をキャッシュすることができるのでMemoizationによる計算改善が見込めるだろうという見込み．

numpyでmemoization

python numpyを使っていくつかの素朴な実装を試す．

import numpy as np
import time
from functools import lru_cache

def HermiteCoeff(n):
    if n == 0:
        return np.array([1])
    elif n == 1:
        return np.array([0, 2])
    else: 
        coeff = np.zeros(n+1)        
        coeff[:n-1] += -2*(n-1)*HermiteCoeff(n-2)
        coeff[1:] += 2*HermiteCoeff(n-1)
        return coeff
    
known0={0:np.array([1], dtype=np.int64),
        1:np.array([0, 2], dtype=np.int64)}
def HermiteCoeffMemo0(n):
    if n in known0:
        return known0[n]
    else: 
        coeff = np.zeros(n+1, dtype=np.int64)        
        coeff[:n-1] += -2*(n-1)*HermiteCoeffMemo0(n-2)
        coeff[1:] += 2*HermiteCoeffMemo0(n-1)
        known0[n] = coeff
        return coeff

known1={0:np.array([1], dtype=o),
        1:np.array([0, 2], dtype=object)}
def HermiteCoeffMemo1(n):
    if n in known1:
        return known1[n]
    else: 
        coeff = np.zeros(n+1, dtype=object)        
        coeff[:n-1] += -2*(n-1)*HermiteCoeffMemo1(n-2)
        coeff[1:] += 2*HermiteCoeffMemo1(n-1)
        known1[n] = coeff
        return coeff

@lru_cache()
def HermiteCoeffMemo2(n):
    if n == 0:
        return np.array([1], dtype=object)
    elif n == 1:
        return np.array([0, 2], dtype=object)
    else: 
        coeff = np.zeros(n+1, dtype=object)        
        coeff[:n-1] += -2*(n-1)*HermiteCoeffMemo2(n-2)
        coeff[1:] += 2*HermiteCoeffMemo2(n-1)
        return coeff
    
def Hermite(coeff):
    return lambda x:sum(c*x**i for i, c in enumerate(coeff))  

import time

funcs = [HermiteCoeff, HermiteCoeffMemo0, HermiteCoeffMemo1, HermiteCoeffMemo2]
x = np.linspace(-1, 1, 100)
n = 30
for f in funcs:
    time0 = time.time()    
    coeff = f(n)
    H = Hermite(coeff)
    y = H(x)
    y0 = y[len(x)//2]
    print("eval. time:%10f[sec]" % (time.time()-time0), "H_%d(0)=%5e" % (n, y0), "by", f.__name__, coeff.dtype)    
 

実行結果は

eval. time: 11.646824[sec] H_30(0)=-2.022226e+20 by HermiteCoeff float64
eval. time:  0.000543[sec] H_30(0)=7.076253e+16 by HermiteCoeffMemo0 int64
eval. time:  0.000932[sec] H_30(0)=-2.022226e+20 by HermiteCoeffMemo1 object
eval. time:  0.000933[sec] H_30(0)=-2.022226e+20 by HermiteCoeffMemo2 object

Memoizationによって速度が改善されることは期待通りであるがいくつか注意しないといけない．
numpy.zerosはゼロで初期化された長さ$n$のarrayを返すが，dtypeを指定していない場合，floatのzeroで埋められる．エルミート多項式の評価では係数リストの要素はintなので無駄である．そこで(Memoization)+dtype=np.intにした結果が2番目．速度が早くなっているのは良いが，$H_{30}(0)$の結果が異なる．これは以前の記事がoverflowしている可能性がある．

そこで3番目はmemoizationとarrayで多倍長整数を扱えるようにするため，dtype=objectにした結果intの場合より若干遅くなるが最初の定義と同じ値を返すことが確認できる．もう少し，次数$n=100$にした結果が次の通り

eval. time:  0.001842[sec] H_100(0)=-2.410424e+18 by HermiteCoeffMemo0 int64
eval. time:  0.003694[sec] H_100(0)=3.037263e+93 by HermiteCoeffMemo1 object
eval. time:  0.003563[sec] H_100(0)=3.037263e+93 by HermiteCoeffMemo2 object

dtypeをnp.intにした場合は全くでたらめな値となっているので，やはりoverflowしているであろうと想像できる．またlru_cacheでもお手製のキャッシュ機構も殆ど速度が変わらないので，既存のlru_cacheを使った方が良いだろう．

次で比較するために$n=250$の結果も載せておく(仮数部100桁なので新の値から0.1程度ずれているが・・・)．

$ python hermite-poly-coeff.py
eval. time:  0.005363[sec] H_250(0)=0.000000e+00 by HermiteCoeffMemo0 int64
eval. time:  0.012851[sec] H_250(0)=-1.673543e+283 by HermiteCoeffMemo1 object
eval. time:  0.013394[sec] H_250(0)=-1.673543e+283 by HermiteCoeffMemo2 object

polynomial manipulation packages

一つ前の章の素朴実装では手で実装したが，汎用性は皆無であり多項式の代数演算の需要は極めて大きいので，自ら実装するのは車輪の再開発であろう．実際google先生に相談したところ

pythonでは

• numpy.polynomial
• mpmath.polynomials
  juliaでは
• Polynomials
  のパッケージを使えば自ら実装することなく，多項式代数演算を実現することが可能．

との啓示をうける．素朴実装を改良することができのではないかと期待できるわけだ(が上手く行かない)．
また，Symbolicな計算が可能なsympyやwolfram言語ではよりsimpleに実装可能(後に考察する)．

numpy.polynomial

Document numpy.polynomial.polynomial.Polynomialを読んでもよくわからないのだが，

>>> import numpy as np 
>>> from numpy.polynomial.polynomial import Polynomial
>>> f = Polynomial([1,2,3])                                                                                         
Polynomial([1., 2., 3.], domain=[-1,  1], window=[-1,  1])
>>> f.coef                                                                                                      
array([1., 2., 3.])

とすれば$P(x) = 1 + 2x + 3x^2$の多項式クラスが返される．係数リストの要素は整数として[1,2,3]を代入したが，各要素は浮動小数点となっている．これは潜在的にoverflowの可能性があるので多倍長整数で評価したい．ドキュメントには書いていないが

>>> f = Polynomial(np.array([1,2,3], dtype=object))                                                             
>>> f.coef                                                                                                      
array([1, 2, 3], dtype=object)

とすれば良さそう．memoizationはお手製のキャッシュ機構でも同速度なのでlru_cacheをつかう．

import numpy as np
from numpy.polynomial.polynomial import Polynomial
from functools import lru_cache
import time
import mpmath
mpmath.mp.dps = 100

@lru_cache(maxsize=1000)
def HermitePoly(n, dtype=None):
    if n == 0:
        return Polynomial(np.array([1], dtype=dtype))
    elif n == 1:
        return Polynomial(np.array([0,2], dtype=dtype))
    else:
        xc = Polynomial(np.array([0,1], dtype=dtype) )
        return 2*xc*HermitePoly(n-1) -2*(n-1)*HermitePoly(n-2)


mpmath.mp.dps = 100
x = np.array(mpmath.linspace(1, 1, 100)) 
n = 250
for dtype in [np.int, np.float, object]: 
    time0 = time.time()
    H = HermitePoly(n, dtype=dtype)
    y = H(x)
    c = H.coef
    print(time.time()-time0, "[sec]", H(0), c.max(), c.dtype)
    HermitePoly.cache_clear()

$ python numpy-poly.py 
0.14997124671936035 [sec] 5.922810347657652e+296 2.46775722331116e+305 float64
0.15068888664245605 [sec] 5.922810347657652e+296 2.46775722331116e+305 float64
0.14925169944763184 [sec] 5.922810347657652e+296 2.46775722331116e+305 object

dtypeでintを指定した場合はfloatにオーバーライドされているが，dtype=objectの係数はobjectoのままだがよくよく考えていみると，係数にはintしか現れないはずなので，c.max()が浮動小数点となっているのはdtypeはobjectだが実際の要素はfloatに丸められているということだ．実際n=260まで上げると

$ python numpy-poly.py 
/home/hanada/Applications/anaconda3/lib/python3.7/site-packages/numpy/polynomial/polynomial.py:788: RuntimeWarning: invalid value encountered in double_scalars
  c0 = c[-i] + c0*x
0.14511513710021973 [sec] nan nan float64
0.14671945571899414 [sec] nan nan float64
0.14644980430603027 [sec] nan inf object

となり倍精度浮動小数点のoverflowが生じる．一応この問題はmpmathを使えば解決できて

mpmath.mp.dps = 10000   
@lru_cache(maxsize=1000)
def HermitePoly1(n, dtype=None):
    if n == 0:
        return Polynomial(np.array([mpmath.mpf("1")]))
    elif n == 1:
        return Polynomial(np.array([mpmath.mpf("0"), mpmath.mpf("2")]))
    else:
        xc = Polynomial(np.array([mpmath.mpf("0"), mpmath.mpf("1")]))
        return 2*xc*HermitePoly1(n-1) -2*(n-1)*HermitePoly1(n-2)


x = np.array(mpmath.linspace(1, 1, 100)) 
n = 250

time0 = time.time()
H = HermitePoly1(n, dtype=dtype)
y = H(x)
c = H.coef
print(time.time()-time0, "[sec]", float(H(0)), c.max(), c.dtype)

とすればちょっと長いが

0.30547070503234863 [sec] -1.7171591075700597e+283 4742832273990616849979871677859751938138722866700194098725540674413432950130755339448602832297348736154189959265033164596368647282052560604201151158911400057598325345216981770764974144951365648905162496309065906209718571075651658676878296726900331786598665420800000000000000000000000000000000.0 object

となることがわかる．ただこれも係数が浮動小数点表記になっているので，おそらくnumpy.polynomial.polynomialの親クラスであるABCPolyBaseでの内部実装の影響だろうか．そこまでは追っていない．

また計算速度は自分で実装したほうが明らかに早いので，多倍長精度計算でnumpy.polynomial.polynomial.Polynomial(最後は大文字ではじまる!!)を使うのは避けた方が良いだろう．

mpmath.polynomial

そもそもnumpyはc言語baseなので無理はできない．餅は餅屋ということでmpmathでも実装してみる(がこれも問題がある)．
mpmathではmpmath (直交多項式系)でHermite多項式が提供されているのでそれを使えば良いのだが，一般的な場合での拡張を考えてmpmath.polynomialを使ってmemoizationできるか考える．

numpyの場合と順番が逆で，係数list, $[c_0, c_1, c_2]$を与えると$P(x) = c_0x^2 + c_1x^1 + c_2$ので注意．

仮数部1000桁で評価してみる．

import mpmath as mp
mp.mp.dps = 1000
import time

素朴に実装すると

def HermitePoly1(n,x):
    if n == 0:
        return mp.polyval([1], x)
    elif n==1:
        return mp.polyval([2,0], x)
    else:
        xc = mp.polyval([1,0], x)        
    return 2*xc* HermitePoly1(n-1, x) - 2*(n-1)*HermitePoly1(n-2,x)

であろうか．polyvalのxに配列(リスト)を代入できないようなので素朴にmemoizationすると

known={0:lambda x: mp.polyval([1], x),
       1:lambda x: mp.polyval([2,0], x)}
def HermitePoly2(n):
    if n in known:
        return known[n]
    else:
        H = lambda x: 2*mp.polyval([1,0], x) * HermitePoly2(n-1)(x) - 2*(n-1)*HermitePoly2(n-2)(x)
        known[n] = H
        return H

であろうか

x = mp.linspace(-5,5,100,endpoint=False)

n = 20
time0 = time.time()
y = [HermitePoly1(n, xx) for xx in x]
print(time.time()-time0, "[sec]", HermitePoly1(n, 0))

time0 = time.time()
H = HermitePoly2(n)
y = [H(xx) for xx in x]
print(time.time()-time0, "[sec]", H(0))

として実行する．結果は

$ python mpmath-poly.py 
17.37535572052002 [sec] 670442572800.0
17.98806405067444 [sec] 670442572800.0

高速化できない．なぜならlamda式のキャッシュはあくまで，評価する順番のキャッシュであり，
実際の評価結果はキャッシュされていないからである．あまりmpmath.Polynomialを使ったことが無いのでもっと良い方法があるのかもしれないが，mpmath.Polynomialだででmemoizationする方法が思いつかなかった．mp.hermiteの内部実装が気になるところ．

Hermite多項式の場合は係数だけ多倍長評価すれば良いのでその場しのぎであれば，

from functools import lru_cache
import numpy as np

@lru_cache()
def HermiteCoeff(n):
    if n == 0:
        return np.array([1], dtype=object )
    elif n == 1:
        return np.array([0, 2], dtype=object)
    else: 
        coeff = np.zeros(n+1, dtype=object)        
        #2xH(n-1) - 2*(n-1)*H(n-2)
        coeff[:n-1] += -2*(n-1)*HermiteCoeff(n-2)
        coeff[1:] += 2*HermiteCoeff(n-1)
        return coeff

def Herm(coeff):
    return lambda x:sum(c*x**i for i, c in enumerate(coeff))  
    

time0 = time.time()
coeff = HermiteCoeff(250)
H = Herm(coeff)
coeff = coeff.tolist()[::-1]
y = [mp.polyval(coeff, xx) for xx in x]
print(time.time()-time0, "[sec]", float(mp.polyval(coeff, 0)))

のようにすれば

$ python mpmath-poly.py
0.19260430335998535 [sec] -1.7171591075700597e+283

memoizationできないわけではない．

sympy

次にsympyを使ってみることにする．

import sympy, mpmath
import time
x = sympy.symbols('x')

def HermitePoly1(n):
    if n==0:
        return 1
    elif n==1:
        return x
    else:
        return 2*x*HermitePoly1(n-1) - 2*(n-1)*HermitePoly1(n-2)


known={0:1, 1:x}
def HermitePoly2(n):
    if n in known:
        return known[n]
    else:
        Hn = 2*x*HermitePoly2(n-1) - 2*(n-1)*HermitePoly2(n-2)
        known[n] = Hn
        return Hn

known1={0:1, 1:x}
def HermitePoly3(n):
    if n in known1:
        return known1[n]
    else:
        Hn = 2*x*HermitePoly3(n-1) - 2*(n-1)*HermitePoly3(n-2)
        known1[n] = sympy.expand(Hn)        
        return known1[n]

mpmath.mp.dps = 100
x1 = mpmath.linspace(0, 1, 100)

n = 20
time0 = time.time()
H = HermitePoly1(n)
time1 = time.time()
y = [H.subs(x, i) for i in x1]
print(time.time() - time1, time1 - time0, H.subs(x, 0))


time0 = time.time()
H = HermitePoly2(n)
time1 = time.time()
y = [H.subs(x, i) for i in x1]
print(time.time() - time1, time1 - time0, H.subs(x, 0))


time0 = time.time()
H = HermitePoly3(n)
time1 = time.time()
y = [H.subs(x, i) for i in x1]
print(time.time() - time1, time1 - time0, H.subs(x, 0))

として実行すれば

$ python sympy-poly.py 
1.612032413482666 0.1680293083190918 670442572800
1.5776398181915283 0.012787342071533203 670442572800
0.2140791416168213 0.07543253898620605 670442572800

memoizationできるがsympyもあまりつかったことがないのでこれが最適なのかよくわからない．

Mathematica(wolfram言語)

Mathematica(wolfram言語)であれば，Fibonacci 級数と同じようにHermite多項式の漸化式を定義できる．
素朴な実装では約14秒

He[n_, x_] := 2x*He[n - 1, x] - 2 (n - 1)*He[n - 2, x] 
He[0, x_] = 1;
He[1, x_] = 2x; 

{First@Timing@Expand@He[30, x], "sec"}
Out[4]= {14.214, sec}

計算結果をキャッシュした場合

H[n_, x_] := H[n, x] = 2 x*H[n - 1, x] - 2 (n - 1)*H[n - 2, x] 
H[0, x_] = 1;
H[1, x_] = 2x; 

{First@Timing@Expand@H[30, x], "sec"}
Out[8]= {7.89031, sec}

どうもExpandに時間がかかっていそうなので

h[n_, x_] := h[n, x] = 2x*h[n - 1, x] - 2 (n - 1)*h[n - 2, x] // Expand 
h[0, x_] = 1;
h[1, x_] = 2x; 

{First@Timing@Expand@h[30, x], "sec"}
Out[12]= {0.007158, sec}

Mathematicaで定義されたHermiteHを使うともっと早い．

In[14]:= {First@Timing@HermiteH[30, x], "sec"}
Out[14]= {0.000216, sec}

一体内部で何をやっているんだろう．

Julia

pythonでいうところのpolynomial 演算に関する操作はPolynomials.jlを使えば実現できる．

(疲れてきたのでプログラムは準備中)