オイラー変換

はじめに

オイラー変換なるものはいくつかあるようですが、ここでは交代級数の収束を早くする手法であるオイラー変換について話します。これによりコンピュータでの交代級数の計算を高速にできます。
オイラー変換という難しそうな名前がついていますが高校数学程度の知識があれば理解できると思います。

交代級数

級数とは数を無限に足し合わせたもののことです。交代級数はその一種で、正の数、負の数を交互に無限に足し合わせたもののことです。$${a_n\geq0}$$として、

$$
\sum_{n=0}^{\infty} (-1)^na_n=a_0-a_1+a_2-a_3+\cdots
$$

具体例:

$$
\begin{align*}
\ln{2}=\log_e{2} = 1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}+\frac{1}{5}-\cdots\\
\frac{\pi}{4} = 1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\frac{1}{9}-\cdots
\end{align*}
$$

オイラー変換のイメージ

具体例として先ほど挙げた$${\ln{2}}$$を考えます。つまり$${a_n=\frac{1}{n+1}}$$です。どのように$${\ln{2}}$$に収束するのか少し計算してみましょう。第$${n}$$項までの部分和を$${S_n}$$とします。

$$
\begin{align*}
S_3&=1-\frac{1}{2}+\frac{1}{3}=0.8333\ldots\\
S_4&=1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}=0.5888\ldots\\
S_5&=1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}+\frac{1}{5}=0.78333\ldots\\
S_6&=1-\frac{1}{2}+\frac{1}{3}-\frac{1}{4}+\frac{1}{5}-\frac{1}{6}=0.61666\ldots\\
\ln{2}&=0.693147\ldots
\end{align*}
$$

なかなか収束に時間がかかりそうです。収束のイメージをするために、$${a_n(0\leq n\leq8)}$$の値をグラフにしましょう。

線と$${x}$$軸との間に色付けをしました。$${x}$$軸より下にできる長方形は負の面積だと考えます。この長方形の横の長さは1なので、面積を足し合わせることで部分和となります。例えば左側三つの長方形の面積は$${1,-0.5,0.333\ldots}$$であり、足し合わせると$${S_3=0.8333\ldots}$$です。つまり、この長方形達の面積を無限に足し合わせたものが$${\ln{2}}$$というわけです。逆にいえば、真の値($${\ln{2}}$$)に対する部分和の打ち切り誤差は足し合わせなかった長方形の面積の和です。
さて、次の図のように長方形を二つに分割して、隣同士の長方形の面積を足し合わせることを考えます。

隣同士は符号が異なるので、打ち消し合います。結果として下図の橙色の線のようになります。

橙色の線は先頭の部分を除いて、青色の線と同様に符号が反転しているため交代級数となっています。なので、その部分に対して先ほどと同じように長方形の分割、隣同士の長方形の面積を足し合わせることをしましょう。

さらにもう4回繰り返すと

交代級数の符号の異なる値同士が相殺されて、先頭の+符号のみになりその分早く収束します。注意点としては、隣同士の面積を足しているのでその項よりも先の項も計算利用していることです。例えば、先ほどの図の緑色の線の区間(1,2)での値$${0.08333\ldots}$$は$${(a_0-2a_1+a_2)/2^2}$$で、$${a_2}$$の値が必要です。しかし、それを考慮しても長方形の分割&隣同士の加算を行った方が真の値に近い値を得られます。

オイラー変換による加速

それでは隣同士の項を足し合わせるのを式で考えてみましょう。$${\binom{n}{r}}$$は二項係数で、$${\binom{n}{r}={}_n\mathrm{C}_r}$$です。

$$
\begin{align*}
S&=a_0-a_1+a_2-a_3+a_4-\cdots\\
&=\frac{a_0}{2}+\frac{1}{2}\left\{\left(a_0-a_1\right)-\left(a_1-a_2\right)+\left(a_2-a_3\right)-\left(a_3-a_4\right)+\cdots\right\}\\
&=\frac{a_0}{2}+\frac{a_0-a_1}{2^2}+\frac{1}{2^2}\left\{\left(a_0-2a_1+a_2\right)-\left(a_1-2a_2+a_3\right)+\left(a_2-2a_3+a_4\right)-\left(a_3-2a_4+a_5\right)+\cdots\right\}\\
&=\frac{a_0}{2}+\frac{a_0-a_1}{2^2}+\frac{a_0-2a_1+a_2}{2^3}+\frac{1}{2^3}\left\{\left(a_0-3a_1+3a_2-a_3\right)-\left(a_1-3a_2+3a_3-a_4\right)+\left(a_2-3a_3+3a_4-a_5\right)-\left(a_3-3a_4+3a_5-a_6\right)+\cdots\right\}\\
&\vdots\\
&=\frac{a_0}{2}+\frac{a_0-a_1}{2^2}+\frac{a_0-2a_1+a_2}{2^3}+\cdots+\frac{1}{2^m}\sum_{i=0}^{m-1}(-1)^i\binom{m-1}{i}a_i+\frac{1}{2^m}\left(\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i+\sum_{i=0}^m(-1)^i\binom{m}{i}a_{i+1}+\cdots\right)
\end{align*}
$$

$${a_{m-1}}$$まで使って級数を計算する時、第一項から第$${m-1}$$項までが部分和$${S'_m}$$、その後が剰余項$${R'_m}$$です。

$$
\begin{align*}
S'_m&=\frac{a_0}{2}+\frac{a_0-a_1}{2^2}+\frac{a_0-2a_1+a_2}{2^3}+\cdots+\frac{1}{2^m}\sum_{i=0}^{m-1}(-1)^i\binom{m-1}{i}a_i\\
R'_m&=\frac{1}{2^m}\left(\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i+\sum_{i=0}^m(-1)^i\binom{m}{i}a_{i+1}+\cdots\right)
\end{align*}
$$

剰余項$${R'_m}$$は部分和$${S'_m}$$の誤差を表します。$${R'_m}$$が速く0に収束するなら$${S'_m}$$も速く収束するわけです。
$${a_0>a_1>a_2>\cdots}$$であれば、数列$${{a_n}}$$は単調減少数列です。さらに$${m}$$回まで隣同士の項を足し合わせてできる数列も単調減少数列とします。つまり

$$
a_0>a_1>a_2>a_3>\cdots\\
a_0-a_1>a_1-a_2>a_2-a_3>\cdots\\
a_0-2a_1+a_2>a_1-2a_2+a_3>a_2-2a_3+a_4>\cdots\\
a_0-3a_1+3a_2-a_3>a_1-3a_2+3a_3-a_4>a_2-3a_3+3a_4-a_5>\cdots\\
\vdots\\
\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i>\sum_{i=0}^m(-1)^i\binom{m}{i}a_{i+1}>\sum_{i=0}^m(-1)^i\binom{m}{i}a_{i+2}>\cdots
$$

任意の$${m\in\mathbb{N}}$$についてこれが成り立つ場合、数列$${{a_n}}$$は完全に単調な減少数列といいます。また、完全に単調な減少数列であれば明らかに

$$
a_0>a_0-a_1>a_0-2a_1+a_2>\cdots>\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i>\cdots
$$

よって

$$
\begin{align*}
R'_m&=\frac{1}{2^m}\left(\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i+\sum_{i=0}^m(-1)^i\binom{m}{i}a_{i+1}+\cdots\right)\\
&=\frac{1}{2^m}\left(\frac{1}{2^1}\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i+\frac{1}{2^2}\sum_{i=0}^{m+1}(-1)^i\binom{m+1}{i}a_i+\frac{1}{2^3}\sum_{i=0}^{m+2}(-1)^i\binom{m+2}{i}a_i+\cdots\right)\\
&<\frac{1}{2^m}\left(\frac{1}{2^1}+\frac{1}{2^2}+\frac{1}{2^3}+\cdots\right)\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i\\
&=\frac{1}{2^m}\sum_{i=0}^{m}(-1)^i\binom{m}{i}a_i\\
&<\frac{a_0}{2^m}
\end{align*}
$$

最後の行の式は公比$${\frac{1}{2}}$$の等比級数の打ち切り誤差と同じです。したがって、数列$${{a_n}}$$が完全に単調な減少数列であれば、$${S'_m}$$は公比$${\frac{1}{2}}$$の等比級数と同じかそれより速く収束します。

オイラー変換を用いて実際に計算

$${S'_m}$$を変形し、$${a_j}$$にかかる係数を求めましょう。

$$
\begin{align*}
S'_m&=\frac{a_0}{2}+\frac{a_0-a_1}{2^2}+\frac{a_0-2a_1+a_2}{2^3}+\cdots+\frac{1}{2^m}\sum_{i=0}^{m-1}(-1)^i\binom{m-1}{i}a_i\\
&=a_0\sum_{i=0}^{m-1}\frac{1}{2^{i+1}}\binom{i}{0}-a_1\sum_{i=1}^{m-1}\frac{1}{2^{i+1}}\binom{i}{1}+a_2\sum_{i=2}^{m-1}\frac{1}{2^{i+1}}\binom{i}{2}-\cdots+(-1)^{m-1}a_{m-1}\sum_{i=m-1}^{m-1}\frac{1}{2^{i+1}}\binom{i}{m-1}\\
&=\frac{a_0}{2^{m}}\sum_{i=0}^{m-1}\binom{m}{i+1}-\frac{a_1}{2^{m}}\sum_{i=1}^{m-1}\binom{m}{i+1}+\frac{a_2}{2^{m}}\sum_{i=2}^{m-1}\binom{m}{i+1}-\cdots+(-1)^{m-1}\frac{a_{m-1}}{2^{m}}\sum_{i=m-1}^{m-1}\binom{m}{i+1}\\
&=\frac{1}{2^m}\sum_{j=0}^{m-1}\left\{(-1)^ja_j\sum_{i=j}^{m-1}\binom{m}{i+1}\right\}
\end{align*}
$$

なので、$${S'_m}$$を計算する際の$${a_j}$$にかかる係数を$${c_j^{(m)}}$$とすると

$$
\begin{align*}
c_j^{(m)}&=\frac{(-1)^j}{2^m}\sum_{i=j}^{m-1}\binom{m}{i+1}\\
S'_m&=\sum_{j=0}^{m-1}c_j^{(m)}a_j
\end{align*}
$$

$${c_j^{(m)}}$$は二項係数の和を$${2^m}$$で割った形になります。 $${2^m=(1+1)^m=\sum_{i=0}^m\binom{m}{i}}$$ なので $${\left|c_j^{(m)}\right|<1}$$です。
オイラー変換を使った和の計算をPythonで実装してみましょう。

def euler_transform_sum(a):
    c_sum = 0
    s = 0
    l = len(a)
    c = 1
    for i, e in enumerate(reversed(a)):
        c_sum += c
        s += e*c_sum
        c = c*(l-i)//(i+1)
    return s/2**l

$${\ln{2}}$$を計算してみます。

print((v:=math.log(2)))
print('10項', '通常計算の和', (s:=sum(reversed([(-1 if i%2==1 else 1)/(i+1) for i in range(0,10)]))), '%.1e'%((s-v)/v))
print('10項', 'オイラー変換', (s:=euler_transform_sum([(-1 if i%2==1 else 1)/(i+1) for i in range(0,10)])), '%.1e'%((s-v)/v))
print('30項', '通常計算の和', (s:=sum(reversed([(-1 if i%2==1 else 1)/(i+1) for i in range(0,30)]))), '%.1e'%((s-v)/v))
print('30項', 'オイラー変換', (s:=euler_transform_sum([(-1 if i%2==1 else 1)/(i+1) for i in range(0,30)])), '%.1e'%((s-v)/v))
print('50項', '通常計算の和', (s:=sum(reversed([(-1 if i%2==1 else 1)/(i+1) for i in range(0,50)]))), '%.1e'%((s-v)/v))
print('50項', 'オイラー変換', (s:=euler_transform_sum([(-1 if i%2==1 else 1)/(i+1) for i in range(0,50)])), '%.1e'%((s-v)/v))
# <結果>
# 0.6931471805599453
# 10項 通常計算の和 0.6456349206349206 -6.9e-02
# 10項 オイラー変換 0.6930648561507936 -1.2e-04
# 30項 通常計算の和 0.6767581376913978 -2.4e-02
# 30項 オイラー変換 0.6931471805307892 -4.2e-11
# 50項 通常計算の和 0.6832471605759182 -1.4e-02
# 50項 オイラー変換 0.6931471805599453 0.0e+00

$${1-\frac{1}{2}+\frac{1}{3}-\cdots}$$をそのまま計算すると50項計算しても一桁しか合っていないのに対して、オイラー変換を行った場合は64ビット浮動小数点数の限界まで精度よく計算できてます。元の数列に対して$${c_j^{(m)}}$$をかけて足し合わせていくので、$${c_j^{(m)}}$$の分だけ計算量が増えますが、オーダーとしてはどちらも$${O(m)}$$なのでそこまで問題にはならないと思います。
ただし、64ビット浮動小数点数を使う場合、指数部は1023までしか表せないので、上記のeuler_transform_sumに1024項以上のfloat型のリストを渡すとオーバーフローします。Fraction型やbigfloat等を使えば計算できますが、計算量は増えると思います。
例えば、ライプニッツの公式をオイラー変換して、1500項使って円周率を求めると

s:Fraction = 4*euler_transform_sum([Fraction(-1 if i%2==1 else 1, 2*i+1) for i in range(0,1500)])
decimal.getcontext().prec = int(math.log10(max(s.numerator, s.denominator)))+1
print(decimal.Decimal(s.numerator)/decimal.Decimal(s.denominator))
# 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196442881097566593344612847564823378678316527120190914564856692346034861045432664821339360726024914127372458700660631558817488152092096282925409171536436789259036001133053054882046652138414695194151160943305727036575959195309218611738193261179310511854807185434891688891802598771695551206361267651616711300629551820661115453817768559884996151338228562171990223575929243228051633422815051225276149431043169718873808789239926163956390073401729253816347735931420301067391909220473923512280331985912800475874521369993563822651206638054351554265827242626165881986103016240714586668144351578344584920968293719928321898998308974205773694449647244121769479054862068241812299671557156295655997179007872065149660914957919074658122706373000636853057372594605640431914194194408663719492625892771139319940520480950333865346443582128922458306633037602222700408913475259494974687882113783641230108016563161344083585925604641962325584909202900292897916213292613132232548041311845470410451927704879572624750210920882845022668760361052660709340660521938472990790981597900888723941490293069888413309419762166553859470529

453桁合っています。(そもそもライプニッツの公式で円周率を求めるのは実用的ではないですが、、)

補足

他の文献では演算子$${\Delta}$$を用いて、オイラー変換を

$$
a_0-a_1+a_2-a_3+\cdots=\frac{a_0}{2}-\frac{\Delta a_0}{2^2}+\frac{\Delta^2 a_0}{2^3}-\frac{\Delta^3 a_0}{2^4}+\cdots
$$

のように書いていたりします。演算子$${\Delta}$$の定義は

$$
\begin{align*}
\Delta a_n&=a_{n+1}-a_n\\
\Delta^m a_n&=\Delta^{m-1}a_{n+1}-\Delta^{m-1}a_n \quad(m=2,3,4,\ldots)
\end{align*}
$$

この記法の方が見た目はすっきりしますね。

おわりに

オイラー変換を使うことで交代級数の収束を加速することができます。何かしらの数列の収束を加速する方法のことを加速法といいます。オイラー変換のほかにもリチャードソン補外やエイトケン$${\Delta^2}$$法などがあります。こちらを級数に適用する場合は部分和の数列$${{S_n}}$$の収束を加速させることになると思います(多分)。
数値計算というと連続的な関数に対して適用する補間法や数値積分を使うことが現実的には多いと思います。しかも値を求める効率的な方法が無い交代級数というのも少なそうなので、オイラー変換を使うような場面はほとんど無いかと思いますが、こういうものもあるんだよということで。
私は授業で高速逆ラプラス変換(FILT)を計算するためにオイラー変換を習ったのですが、ネット上にオイラー変換についての記述が少なかったので授業のメモがてらnoteに残しておくことにしました。専門外の内容なので間違っている部分があるかもしれませんが許してください。

参考文献

「お話：数値解析」(理系への数学連載)サポートページ
森口繁一, 計算数学夜話 -数値で学ぶ高等数学-. 日本評論社, 1978