Aitken加速[programming][lisp][python][rust]

概要

数列のAitken加速というのを説明している http://cympfh.cc/aiura/_/aitken.html を見て、真似た。昔授業で習ったはずだが、すでに記憶の彼方だったので。あと最近Rustについて調べてみてるのでその練習も兼ねて幾つかの言語で書き比べてみようとした(python, rust, common lisp)。

数列の加速

http://cympfh.cc/aiura/_/aitken.html の通り、
sに収束する数列 $\{s_n\}$ が $s_n\sim s + \beta \lambda^n$ (|λ|<1)というふうに収束していくとき、隣り合う3項を使ってβとλを消去してやるとsに非常に近い値が得られ、結果として元の数列より早く極限値に収束する数列となるというもの。
解いた結果としては
$s_n'=s_n - \frac{{s_{n+1}-s_n}^2}{s_{n+2}-2s_{n+1}+s_{n}}$
となる。こうやって数列 $\{s_n\}$ から数列 $\{s'_n\}$ を得る手続きをAitken加速という。さらにこの手続きを $s_n'$ にもう一度行うことで、2段のAitken加速となる。k回繰り返したものをk段のAitken加速列 $\{s_n^{k}\}$ という。特に $\{s_n^0\}=\{s_n\}$ 。 $\{s_n\}$ が有限個しかないと、Aitken加速で数列の長さが2ずつ小さくなるので、N項の数列ではN/2段程度までしか行えない。

参考リンク2つ目の内容を読んだメモ:
数列 $\{s_n\}$ がsにp次収束するとは、
$\lim_{n\rightarrow \infty} \frac{|s_{n+1}-s_{n}|}{|s_n-s|^p}=C, (C>0)$
のときいう。またsに収束する数列 $\{s_n\}$ が極限値
$\lim_{n\rightarrow \infty} \frac{s_{n+1}-s}{s_{n}-s}=\lambda$
が存在するときは|λ|<1になるのだが、λ=1の場合を対数収束、λ=0を超線形収束、 $-1\le \lambda \lt 1$ を線形収束という。
超線形収束の場合は加速法は普通はなくても実用上問題ない。対数収束する数列の加速は一般に難しい。

Aitken加速は任意の線形収束する数列を加速する。数値積分とかに使われる。

実装

python(2.7), rust, common lisp (roswell script)で多段Aitken加速を実装してみた。
参考リンク1つ目と同じく、 $s_n = \sum_{i=1}^n\frac{1}{i^2}$ という数列の最初の20個をとってきて、順次加速している。
どれも実行すると

xs^0:
1.000000, 1.250000, 1.361111, 1.423611, 1.463611, 1.491389, 1.511797, 1.527422, 1.539768, 1.549768, 1.558032, 1.564977, 1.570894, 1.575996, 1.580440, 1.584347, 1.587807, 1.590893, 1.593663, 1.596163
xs^1:
1.450000, 1.503968, 1.534722, 1.554520, 1.568312, 1.578464, 1.586246, 1.592399, 1.597387, 1.601510, 1.604977, 1.607931, 1.610479, 1.612698, 1.614650, 1.616378, 1.617920, 1.619304
xs^2:
1.575465, 1.590296, 1.599981, 1.606776, 1.611798, 1.615658, 1.618716, 1.621197, 1.623250, 1.624977, 1.626449, 1.627720, 1.628827, 1.629801, 1.630664, 1.631434
xs^3:
1.618209, 1.622752, 1.626025, 1.628473, 1.630368, 1.631876, 1.633103, 1.634120, 1.634977, 1.635709, 1.636341, 1.636892, 1.637377, 1.637807
xs^4:
1.634458, 1.635750, 1.636852, 1.637743, 1.638467, 1.639062, 1.639559, 1.639980, 1.640341, 1.640653, 1.640925, 1.641165
xs^5:
1.643242, 1.641498, 1.641588, 1.641833, 1.642077, 1.642299, 1.642495, 1.642660, 1.642825, 1.642911
xs^6:
1.641584, 1.641444, 1.444996, 1.644407, 1.644065, 1.643488, 1.618429, 1.643004
xs^7:
1.641584, 1.543955, 1.644065, 1.644906, 1.644078, 1.630836
xs^8:
1.593382, 1.644913, 1.644489, 1.644961
xs^9:
1.644492, 1.644712

という出力になる。この数列はpi^2/6 = 1.644934...に収束する。

実行速度は Rust > Common lisp > Python という感じ。しかしRustはコンパイラの警告が厳しくて*1、あんまり気楽に書けないなあ。pythonは全体的に簡潔に書けてやっぱりうまくできた言語だと感じる。でも少しでも変則的な accumulation になるとlisp(というかloopマクロか)の方がずっと楽。ここのところPythonをずっと書いていたけどそこはずっと不便に感じていた。lispはlispで配列参照と数式を書くのが僕にはちょっと気が重いんだけど*2。

Rust

Common Lisp:

Python:

あとPython版については ./python aitken.py plot として起動すると加速した数列を重ね書きplotするようにしてみた。

加速1段ごとに最初から収束値(点線)に近いところに行っているのが見える。6段くらいで不安定になっている。要するに ${s_{n+2}-2s_{n+1}+s_{n}}$ が小さくなっているのだが、多段で加速しているうちに誤差が積もるかというとまだdoubleだったら余裕な気がするので、あんまり線形収束的な振る舞いをしていない部分ではそうなるということに思われるが。

参考

*1:その割にはindexが範囲外とかのエラーでは落ちるし、いやセキュリティを気にする用途なら落ちて正解なんだろうが

*2:http://qiita.com/y2q_actionman/items/0a7737710ba647697832 あたりを参考になんか作るといいのかもしれない