Aitken 加速

があるときに、その収束値 \(\alpha = a_\infty\) を知りたい. 加速法 は、元の数列よりもその収束が速い数列 \[a_1', a_2', \ldots, a_N', \ldots\]

\(a_n\) の値の計算が定数時間でしかないのなら、 \(a_1, a_2, a_3, \ldots\) と計算する代わりに (例えば) \(a_1, a_4, a_9, \ldots\) と計算すれば、加速するじゃんってことなのだけど、例えば、 \(a_n\) は漸化式で計算できない等の理由によって、\(a_n\) の値の計算に \(O(n)\) の時間が掛かる等、何かそういう理由で、無限列の中の有限の断片 \[a_1, a_2, \ldots, a_N\] しか手に入らないと仮定してるのだと思う. たぶん.

Richardson 加速

ここで \(\alpha\) が実際の収束値で、\(\gamma\) は大きさが \(1\) 未満の値であって、 \(n \to \infty\) に対してこれが \(0\) に収束する.

今、都合よく右辺の \(\gamma\) だけが既知だとすると (これは結構強い仮定だが)、

を得る. この右辺は \(a_n\) よりもずっと \(\alpha\) そのものに近い. というわけで

また、\(a_n'\) の計算コストは \(a_n\) を計算する時間の二倍でしかないし、というか \(a_1, \ldots, a_N\) を初めに全て計算してメモしておけば、計算時間は、差を取って割り算する分しか増えない.

Aitken 加速

Aitken 加速は、Richardson 加速で既知だと仮定していた \(\gamma\) も推定する. 2つの未知変数 (\(\alpha, \beta\)) を消すために3つ用意して弄る. すなわち、

\[\begin{align} a_n' & = \frac{a_{n+1} - \gamma a_n}{1 - \gamma} \\ & = a_n - \frac{(a_{n+1} - a_n)^2}{a_{n+2} - 2a_{n+1} + a_n} \end{align}\]

を得る. この方法によって、数列 \(\{a_n\}\) から \(\{a_n'\}\) を得る方法を Aitken 加速 という. 特にこの方法は繰り返し適用できる. すなわち、\(\{a_n'\}\) から \(\{a_n''\}\) を構成し、更なる加速が期待できる. \(k\) 回 Aitken 加速を繰り返し適用したものを、\(k\) 段のAitken加速列 \(\{a_n^k\}\) という. 特に \(\{a_n\} = \{a_n^0\}\) である.

Aitken 加速の実装

初めに述べたように、有限の断片 \(a_1, a_2, \ldots, a_N\) しか手に入らない場合を考える.

また、一回の Aitken 加速は、元の数列の連続した3つの成分から一つの成分を作るので、数列の長さは 2 ずつ減る. 従って、\(N/2\) 段階程度しか繰り返しの加速は適用できない.

; リストとして与えた有限長の数列を Aitken 加速する
(define (accell ls)

  (define return (list ls))
  (define n (length ls))

  ; k 段の第 n 成分 を (memo-get (cons k n)) に置くメモ
  ; k = 0 には元の数列の値を入れて初期化する
  (define memo (make-hash-table 'equal?))
  (define (memo-get key) (hash-table-get memo key #f))
  (define (memo-put! key val) (hash-table-put! memo key val))

  ; 数列 a の k 段 第 n 成分 の計算
  (define (accell-get k n)

    (let ((i   (memo-get (cons (- k 1) n)))
          (i+1 (memo-get (cons (- k 1) (+ n 1))))
          (i+2 (memo-get (cons (- k 1) (+ n 2)))))

      (- i
         (/
           (expt (- i+1 i) 2)
           (+ i+2 (- (* 2 i+1)) i)))))

  ; 初期化
  (for-each (lambda (i val) (memo-put! (cons 0 i) val))
            (iota n) ls)

  ; k 段加速列を求める
  ; len は k-1 段加速列の長さ (3以上必要)
  (let loop-k ((k 1) (len n))
    (if (>= len 3)
      (begin
        (let ((k-list (map (lambda (i) (accell-get k i)) (iota (- len 2)))))
            ; memo
            (for-each (lambda (i val) (memo-put! (cons k i) val))
                      (iota (- len 2)) k-list)
            ; append to return
            (set! return (cons k-list return))
            (loop-k (+ k 1) (- len 2))))
      (reverse return))))

Aitken 加速のテスト

という数列の初めの 20 個だけがある時にこの数列の収束を加速する.

; 元の数列 a_n
(define (a n) (exact->inexact (apply + (map (lambda (i) (/ 1 i i)) (iota n 1)))))

(define (main)
  (let* ((ls (map a (iota 20)))
         (accelled (accell ls)))

    (for-each print accelled)))

(main)

出力は次のようになった. 第 \(k\) 行の \(i\) 個目が \(k\) 段階の \(i\) 成分の値である (\(i=0,1,\ldots\))

0 行目が元の数列であるが、これを見ただけではまさか収束値が 1.6 を超えるようには見えない. 最後の行を見ると、1.64 くらいを定めそうである.

概要・加速法とは何か

Richardson 加速

Aitken 加速

Aitken 加速の実装

Aitken 加速のテスト