カハンの加算アルゴリズム

提供: miniwiki
移動先:案内、検索

カハンの加算アルゴリズム（Kahan summation algorithm）とは、有限精度の浮動小数点数列の総和を計算する際の誤差を改善する計算手法・アルゴリズム。基本的にコンピュータ上で使用される。Compensated Summation（補正加算）とも呼ぶ[1]。

単純に n 個の数値の総和を計算すると、n に比例して誤差が増えていくという最悪のケースがありうる。また、無作為な入力では二乗平均平方根の誤差すなわち [math]\sqrt{n}[/math] に比例する誤差が生じる（丸め誤差はランダムウォークを形成する）[2]。補正加算では最悪の場合の誤り限界 (error bound) は n とは独立なので、多数の数値を合計しても、誤差は使用する浮動小数点数の精度に依存するだけとなる[2]。

このアルゴリズムは考案したウィリアム・カハンの名を取り、こう呼ばれる[3]。似たようなそれ以前の技法として、例えばブレゼンハムのアルゴリズムがあり、整数演算での誤差の蓄積を保持する（文書化されたのはカハンとほぼ同時期である[4]）。また、ΔΣ変調では誤差（雑音）を加算するのみでなく積分する[5]。

擬似コードと解説

このアルゴリズムの擬似コードは以下の通り: 

function kahanSum(input)
    var sum = 0.0
    var c = 0.0             // 処理中に失われる下位ビット群の補償用変数
    for i = 1 to input.length do
        y = input[i] - c          // 問題なければ、c はゼロ
        t = sum + y               // sum が大きく y は小さいとすると、y の下位ビット群が失われる
        c = (t - sum) - y         // (t - sum) は y の上位ビット群に相当するので y を引くと下位ビット群が得られる（符号は逆転）
        sum = t                   // 数学的には c は常にゼロのはず。積極的な最適化に注意
    next i                  // 次の繰り返しで y の失われた下位ビット群が考慮される
    return sum

6桁の十進浮動小数点演算を例として動作を見てみよう。コンピュータは普通二進演算だが、基本原理は同じである。sum の値が 10000.0 で、次の input(i) が 3.14159 と 2.71828とする。正確な加算結果は 10005.85987 であり、6桁に丸めると10005.9となる。単純に加算すると1回目の加算で 10003.1、2回目で 10005.8 となり、正しくない。

c の初期値はゼロとする。 

  y = 3.14159 - 0                    y = input[i] - c
  t = 10000.0 + 3.14159
    = 10003.1                        大部分の桁が失われた
  c = (10003.1 - 10000.0) - 3.14159  これは書かれた通りに評価される必要がある
    = 3.10000 - 3.14159              y の失われた部分を得るため、本来の y との差分を得る
    = -.0415900                      6桁なので、最後にゼロが付与される
sum = 10003.1                        このように input(i) の一部の桁しか加算されていない

合計が非常に大きいため、入力数値の一部の桁しか反映されない。しかし、ここで次の intput(i) の値が 2.71828 とすると、今回は c がゼロでなはないため次のようになる。 

  y = 2.71828 - -.0415900            前回加算できなかった部分をここで反映する
    = 2.75987                        大きな変化はなく、ほとんどの桁が有効に計算される
  t = 10003.1 + 2.75987              しかし、総和へ加算しようとすると一部の桁しか考慮されない
    = 10005.9                        丸めが発生している
  c = (10005.9 - 10003.1) - 2.75987  反映されなかった分を計算する
    = 2.80000 - 2.75987              今回は加算された値が大きすぎる
    = .040130                        しかし、次の繰り返しで反映するので問題ない
sum = 10005.9                       正確な値は 10005.85987 であり、正しく6桁に丸められている

ここで、総和は2つの部分に分けられて実行されると考えられる。すなわち、sum は総和を保持し、c は sum に反映されなかった部分を保持する。そして次の繰り返しの際に sum の下位桁への補正を試みる。これは、何もしないよりはよいが、精度（桁数）を倍にした方がずっと効果があるのも事実である。ただし、単純に桁数を増やすことが現実的とは限らない。input が倍精度だった場合、四倍精度をサポートしているシステムは少ないし、四倍精度を採用するなら input 内のデータも四倍精度にしなければならない。

精度

補正加算における誤差を注意深く分析することで、その精度の特性がわかる。単純な総和の計算よりも正確だが、悪条件の総和では相対誤差が大きくなる。

i=1,...,n の n 個の数値 xi の合計を計算するとする。その計算は次の式で表される。

[math]S_n = \sum_{i=1}^n x_i[/math] （無限の精度で計算する場合）

補正加算では、[math]S_n + E_n[/math] で総和が表され、誤差 [math]E_n[/math] について次が成り立つ[2]。

[math]|E_n| \leq \left[ 2\varepsilon + O(n\varepsilon^2) \right] \sum_{i=1}^n |x_i| [/math]

ここで ε は使用する算術体系の計算機イプシロンである（例えば、IEEEの倍精度浮動小数点数の場合は ε≈10−16）。通常、関心のある量は相対誤差（English版） [math]|E_n|/|S_n|[/math] であり、相対誤差は上の式から次のような条件となる。

[math]\frac{|E_n|}{|S_n|} \leq \left[ 2\varepsilon + O(n\varepsilon^2) \right] \frac{\sum_{i=1}^n |x_i|}{\left| \sum_{i=1}^n x_i \right|}. [/math]

この相対誤差の境界条件式において、分数 Σ|xi|/|Σxi| が総和問題の条件数である。計算方法がどうであれ、この条件数が総和問題の本質的な誤差への敏感さを表している[6]。固定精度を使った固定の（すなわち、任意精度演算のようにデータによって時間および領域の計算量が変化するアルゴリズムではない）アルゴリズムによる全ての（後方安定な）総和計算技法の相対誤差条件は、その条件数に比例する[2]。「悪条件」の総和問題では、その比率が大きくなり、補正加算であっても相対誤差が大きくなる。例えば被加数 xi が平均値がゼロの無相関の乱数列の場合、その総和はランダムウォークとなり、条件数は [math]\sqrt{n}[/math] に比例して成長する。一方、入力が無作為であっても平均がゼロでなければ、[math]n\to\infty[/math] に伴って条件数は有限の定数に漸近することになる。入力が全て負でない場合、条件数は1となる。

固定の条件数が与えられると、補正加算の相対誤差は事実上 n とは独立となる。原理的に n よって線型に成長する O(nε2) があるが、この項は実質的にゼロと見なせる。というのも最終結果が精度 ε に丸められるので、n がおよそ 1/ε かそれ以上でない限り nε2 という項はゼロに丸められる[2]。倍精度の場合、その項が無視できなくなる n の値は 1016 ぐらいであり、多くの場合それほどの数値の総和を求めるのは現実的でない。したがって、条件数が固定なら補正加算の誤差は事実上 O(ε) となって、n とは独立である。

それに比べ、加算のたびに丸めが発生する単純な総和計算では、相対誤差は [math]O(\varepsilon n)[/math] と条件数をかけた値として成長していく[2]。ただし、この最悪ケースは丸め方向が毎回同じ場合のみ発生するので、実際にはめったに観察されない。実際には丸め方向は毎回無作為に変化し、その平均はゼロに近づくことが多い。その場合の単純な総和の相対誤差は二乗平均平方根となり、[math]O(\varepsilon \sqrt{n})[/math] と条件数をかけた値として成長していく[7]。その場合でも補正加算より誤差が大きくなる。ただし、精度を2倍にすれば ε が ε2 となるので、単純総和の誤差は O(nε2) となって、元の精度での補正加算に匹敵するようになる。

代替手法

カハンのアルゴリズムでの誤差成長は入力数 n に対して [math]O(1)[/math] だが、pairwise summation では若干悪い [math]O(\log n)[/math] となる。これは入力を再帰的に半分に分けていって、再帰的に加算を行う方式である[2]。この技法は単純な総和計算に比べて加算回数が増えないという利点があり（カハンのアルゴリズムでは演算回数が4倍に増える）、並列計算も可能という利点がある。通常、再帰の行き着いた基本ケースでは1回または0回の加算になるが、再帰のオーバーヘッドを低減させるため n がある程度小さくなったらそれ以上再帰させないという方式も考えられる。pairwise summation と同様の技法は高速フーリエ変換 (FFT) アルゴリズムでよく使われており、そのためFFTでは誤差が対数的に成長することが多い[8]。実際には丸め誤差の符号は無作為に変化するので、pairwise summation の二乗平均平方根誤差の成長は [math]O(\sqrt{\log n})[/math] となる[7]。

もう1つの代替手法として、何よりも丸め誤差を生じないことが優先されるなら、任意精度演算を使うことが考えられる。Shewchukは、特に高精度ではないが正確に丸められた総和を求める技法として、任意精度演算を使ってそれなりのコストで計算する方法を示した[9][10]。KirchnerとKulischはビット幅の大きいアキュムレータを使い、整数演算だけで総和を求める技法を示した[11]。ハードウェアの実装を示した例として Müller、Rüb、Rülling の論文がある[12]。

10の20乗といったような絶対値が極端に大きい数を扱うことが無い、会計などには、浮動小数点方式ではなく、必要十分な桁数以上を確保した固定小数点方式（また、場合によっては十進演算）を、可能であれば使う方がよい。

最適化コンパイラにおける注意

コンパイラが生成するコードは書かれたプログラムと正確に対応しているのが普通である。コンパイラはプログラムの内容を数式処理的に解析して最適化を施すことがあるが、コンパイラの設計によってはカハンの加算アルゴリズムについて間違った最適化を行なうことがある。例えば、以下のコードについて 

  t:=sum + y;
  c:=(t - sum) - y;

コンパイラは、結合法則を適用して以下のように推論することがある。 

  c = 0

すると補正がなされなくなる[13]。ただし、一般にコンパイラは浮動小数点演算では近似的にしか結合法則が成り立たないとして、明示的に「安全でない」最適化を指示されないかぎり、このような最適化を行わない[14][15][16]。なお、Intel C++ Compiler のようにデフォルトで結合法則を適用した最適化を行うコンパイラもある[17]。K&R版の最初のC言語では、結合法則による浮動小数点演算の項の並べ替えを許していたが、ANSI C 規格ではそれが禁止され、C言語を数値解析分野で使いやすくした（FORTRANでも並べ替えが禁止されている）[18]。それでもオプションを指定すれば並べ替えできるようにしてあるコンパイラが多い。

一方、一部の言語では総和機能を提供している（APLの +/input、FORTRANのSUMなど）。これらは最良の手法で総和を計算するような実装であると期待される。しかし、FORTRANのマニュアルを見ても単に入力と同じ精度で計算するとあるだけで詳細は不明である。線型代数学の標準的サブルーチン集であるBLASでは、高速化のための演算順序の並べ替えを明確に避けているが[19]、BLASの実装ではカハンのアルゴリズムを採用していないことが多い。Pythonの標準ライブラリには fsum という総和関数があり、Shewchukのアルゴリズムを使い丸め誤差の蓄積を防いでいる。

脚注

  1. ↑ 「補正加算」と呼ばれる技法は他にもある。詳しくは Higham, Nicholas (2002). Accuracy and Stability of Numerical Algorithms (2 ed). SIAM, 110–123. 
  2. ↑ 2.0 2.1 2.2 2.3 2.4 2.5 2.6 Higham, Nicholas J. (1993), “The accuracy of floating point summation”, SIAM Journal on Scientific Computing 14 (4): 783–799, doi:10.1137/0914050 
  3. ↑ Kahan, William (January 1965), “Further remarks on reducing truncation errors”, Communications of the ACM 8 (1): 40, doi:10.1145/363707.363723 
  4. ↑ Jack E. Bresenham, "Algorithm for computer control of a digital plotter", IBM Systems Journal, Vol. 4, No.1, January 1965, pp. 25–30
  5. ↑ H. Inose, Y. Yasuda, J. Murakami, "A Telemetering System by Code Manipulation – ΔΣ Modulation," IRE Trans on Space Electronics and Telemetry, Sep. 1962, pp. 204–209.
  6. ↑ L. N. Trefethen and D. Bau, Numerical Linear Algebra (SIAM: Philadelphia, 1997).
  7. ↑ 7.0 7.1 Manfred Tasche and Hansmartin Zeuner Handbook of Analytic-Computational Methods in Applied Mathematics Boca Raton, FL: CRC Press, 2000).
  8. ↑ S. G. Johnson and M. Frigo, "Implementing FFTs in practice, in Fast Fourier Transforms, edited by C. Sidney Burrus (2008).
  9. ↑ Jonathan R. Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, Discrete and Computational Geometry, vol. 18, pp. 305–363 (October 1997).
  10. ↑ Raymond Hettinger, Recipe 393090: Binary floating point summation accurate to full precision, Python implementation of algorithm from Shewchuk (1997) paper (28 March 2005).
  11. ↑ R. Kirchner, U. W. Kulisch, Accurate arithmetic for vector processors, Journal of Parallel and Distributed Computing 5 (1988) 250-270
  12. ↑ M. Muller, C. Rub, W. Rulling Exact accumulation of floating-point numbers, Proceedings 10th IEEE Symposium on Computer Arithmetic (Jun 1991), doi 10.1109/ARITH.1991.145535
  13. ↑ Goldberg, David (March 1991), “What every computer scientist should know about floating-point arithmetic” (PDF), ACM Computing Surveys 23 (1): 5–48, doi:10.1145/103162.103163, http://www.validlab.com/goldberg/paper.pdf 
  14. ↑ GNU Compiler Collection manual, version 4.4.3: 3.10 Options That Control Optimization, -fassociative-math (Jan. 21, 2010).
  15. ↑ Compaq Fortran User Manual for Tru64 UNIX and Linux Alpha Systems, section 5.9.7 Arithmetic Reordering Optimizations (retrieved March 2010).
  16. ↑ Eric Fleegal, "Microsoft Visual C++ Floating-Point Optimization", Microsoft Visual Studio Technical Articles (June 2004).
  17. ↑ Martyn J. Corden, "Consistency of floating-point results using the Intel compiler," Intel technical report (Sep. 18, 2009).
  18. ↑ Tom Macdonald, "C for Numerical Computing", Journal of Supercomputing vol. 5, pp. 31–48 (1991).
  19. ↑ BLAS Technical Forum, section 2.7 (August 21, 2001).

外部リンク



「http:///http:///mymemo.xyz/mymemo.xyz/wiki/index.php?title=カハンの加算アルゴリズム&oldid=112323」から取得