カハン加算アルゴリズム

数値解析において、補正和としても知られるカハン加算アルゴリズム^[1]は、単純な手法と比較して、有限精度浮動小数点数のシーケンスを加算することによって得られる合計の数値誤差を大幅に削減します。これは、別の実行補正（小さな誤差を蓄積する変数）を保持することで行われ、実際には補正変数の精度だけ合計の精度が拡張されます。

特に、単純に数字を順番に足し合わせると、最悪の場合の誤差は に比例して大きくなり、ランダム入力の場合は二乗平均平方根誤差が に比例して大きくなります（丸め誤差はランダムウォークを形成します）。^[2] 補正加算では、十分に高い精度の補正変数を使用して、最悪の場合の誤差境界は に実質的に依存しないため、結果の浮動小数点精度にのみ依存する誤差で多数の値を加算することができます。 ^[2]${\displaystyle n}$${\displaystyle n}$${\displaystyle {\sqrt {n}}}$${\displaystyle n}$

このアルゴリズムはウィリアム・カハンによるものとされている。^[3] イヴォ・バブシュカも独自に同様のアルゴリズムを考案したようだ（カハン・バブシュカ和と呼ばれる）。^[4]同様の先行技術としては、例えば、整数演算における累積誤差を追跡するブレゼンハムの直線アルゴリズム（同時期に初めて文書化されたが^[5]）やデルタシグマ変調などがある。^[6]

疑似コードでは、アルゴリズムは次のようになります。

関数KahanSum(入力)
    // アキュムレータを準備します。
    var合計 = 0.0
    // 失われた下位ビットの実行中の補償。
    変数c = 0.0
    // 配列input にはinput[1] から input[input.length] までのインデックスが付けられた要素があります。
     for i = 1 to input.length do 
        //最初は
        c は0 です。 var y = input[i] - c
        // 残念ながら、sumは大きく、y は小さいので、 yの下位桁は失われます。var         
         t = sum + y
        // (t - sum) はyの上位部分をキャンセルします。
        // y を減算すると負の値（ yの低い部分）が回復します
        c = (t - 合計) - y
        // 代数的には、cは常にゼロになるはずです。注意してください
        // 過度に積極的な最適化コンパイラ!
        合計 = t
    // 次回は、失われた低い部分が新たな試みで
    yに追加されます。 next i

    合計
を返す

このアルゴリズムはFast2Sumアルゴリズムを使って書き直すこともできる: ^[7]

関数KahanSum2(入力)
    // アキュムレータを準備します。
    var合計 = 0.0
    // 失われた下位ビットの実行中の補償。
    変数c = 0.0
    // 配列input には、 i = 1からinput.lengthまでのインデックスが付けられた要素があります 
    。 
        //最初は
        c は0 です。 var y = input[i] + c
        // sum + c は正確な合計の近似値です。
        (合計,c) = Fast2Sum(合計,y)
    // 次回は、失われた低い部分が新たな試みで
    yに追加されます。 next i
     return sum

このアルゴリズムでは、基数の特定の選択は必須ではなく、演算が「四捨五入または切り捨ての前に浮動小数点の合計を正規化」することだけを必須としています。^[3]コンピュータは通常、2進演算を使用しますが、例を読みやすくするために、10進数で示します。6桁の10進浮動小数点演算 を使用しているとしますsum。 が値 10000.0 に達し、 の次の2つの値がinput[i]3.14159 と 2.71828 であるとします。正確な結果は 10005.85987 で、10005.9 に切り捨てられます。単純な合計では、各入力値は に揃えられsum、多くの下位桁が失われます (切り捨てまたは四捨五入により)。四捨五入後の最初の結果は 10003.1 になります。 2 番目の結果は、四捨五入前は 10005.81828、四捨五入後は 10005.8 になります。これは正しくありません。

ただし、補正された合計を使用すると、正しく丸められた結果 10005.9 が得られます。

の初期値はゼロであると仮定しますc。6 桁の浮動小数点数で重要な場合は末尾のゼロが表示されます。

  y = 3.14159 - 0.00000              y = 入力[i] - c 
  t = 10000.0 + 3.14159              t = 合計 + y
    = 10003.14159 正規化が完了しました。次に 6 桁に切り捨てます。
    = 10003.1入力[i]の桁のうち、合計の桁と一致する桁はほとんどありません。多くの桁が失われています。
  c = (10003.1 - 10000.0) - 3.14159 c = (t - sum) - y   (注意: 括弧を最初に評価する必要があります!)= 3.10000 - 3.14159 y
    の同化された部分から元の完全なy を引いたもの。
    = -0.0415900 cはゼロに近いため、正規化により浮動小数点数の後の桁が多く保持されます。
合計 = 10003.1                       合計 = t

合計が非常に大きいため、入力数値の上位桁のみが累積されます。しかし、次のステップでは、c実行中の誤差の近似値である が問題に対処します。

  y = 2.71828 - (-0.0415900) c はyとほぼ同じ大きさなので、ほとんどの数字が一致します。
    = 2.75987 前回の反復の不足分（下位桁の損失）が正常に復元されました。t = 10003.1 + 2.75987 しかし、合計
  の桁を満たすものはまだわずかです。
    = 10005.85987 正規化が完了しました。次は 6 桁に丸めます。
    = 10005.9 ここでも多くの桁が失われていますが、c によって四捨五入が調整されています。c = (10005.9 - 10003.1) - 2.75987 調整されたy
  に基づいて累積誤差を推定します。
    = 2.80000 - 2.75987 予想どおり、下位部分はcで保持され、丸めの影響はまったくないか、またはわずかです。
    = 0.0401300 この反復ではtが少し高すぎたため、超過分は次の反復で減算されます。
合計 = 10005.9 正確な結果は 10005.85987 で、合計は正しく、6 桁に丸められています。

このアルゴリズムは、2 つのアキュムレータを使用して合計を実行します。 はsum合計を保持し、cに同化されない部分を累算して、次回のsumの下位部分を少しずつ動かします。したがって、 の「ガード桁」を使用して合計が実行されます。これは、ガード桁がないよりはましですが、入力の 2 倍の精度で計算を実行するほど良くはありません。ただし、計算の精度を単純に上げることは一般に実用的ではありません。 がすでに倍精度である場合、 4 倍精度 を提供するシステムはほとんどなく、提供されていたとしても は4 倍精度になる可能性があります。 sumcinputinput

補正加算の精度特性を理解するには、補正加算の誤差を注意深く分析する必要があります。補正加算は単純な加算よりも精度は高くなりますが、条件の悪い加算では相対誤差が大きくなる可能性があります。

の値を合計すると仮定します。正確な合計は ${\displaystyle n}$${\displaystyle x_{i}}$${\displaystyle i=1,\,\ldots ,\,n}$

${\displaystyle S_{n}=\sum _{i=1}^{n}x_{i}}$(無限の精度で計算されます)。

補正和を用いると、代わりに が得られる。ここで、誤差は^[2]で制限される。${\displaystyle S_{n}+E_{n}}$${\displaystyle E_{n}}$

${\displaystyle |E_{n}|\leq {\big [}2\varepsilon +O(n\varepsilon ^{2}){\big ]}\sum _{i=1}^{n}|x_{i}|,}$

ここで は使用される演算の機械精度です（例：IEEE標準倍精度浮動小数点）。通常、関心のある量は相対誤差であり、したがって、その上限は次のようになります。 ${\displaystyle \epsilon }$${\displaystyle \varepsilon \approx 10^{-16}}$ ${\displaystyle |E_{n}|/|S_{n}|}$

${\displaystyle {\frac {|E_{n}|}{|S_{n}|}}\leq {\big [}2\varepsilon +O(n\varepsilon ^{2}){\big ]}{\frac {\sum \limits _{i=1}^{n}|x_{i}|}{\left|\sum \limits _{i=1}^{n}x_{i}\right|}}.}$

相対誤差境界の式では、分数は総和問題の条件数です。本質的に、条件数は、計算方法に関係なく、誤差に対する総和問題の固有の感度を表します。 ^{[8]固定精度の固定アルゴリズム（}任意精度の演算を使用するものや、データに基づいてメモリと時間要件が変化するアルゴリズムではないもの）によるすべての（後方安定）総和方法 の相対誤差境界は、この条件数に比例します。 ^[2]悪条件の総和 問題とは、この比率が大きい問題であり、この場合は補正された総和でも相対誤差が大きくなる可能性があります。たとえば、加数が平均ゼロの相関のない乱数である場合、合計はランダムウォークであり、条件数は に比例して増加します。一方、平均がゼロでないランダム入力の場合、条件数は として有限定数に漸近します。入力がすべて非負の場合、条件数は 1 になります。 ${\displaystyle \Sigma |x_{i}|/|\Sigma x_{i}|}$${\displaystyle x_{i}}$${\displaystyle {\sqrt {n}}}$${\displaystyle n\to \infty }$

条件数が与えられると、補正された合計の相対誤差は実質的に に依存しません。原理的にはとともに線形に増加する がありますが、実際にはこの項は実質的にゼロです。最終結果は精度 に丸められるため、がおよそまたはそれより大きい場合を除き、項はゼロに丸められます。 ^[2] 倍精度では、これはおよそ の に相当し、ほとんどの合計よりもはるかに大きくなります。したがって、条件数が固定されている場合、補正された合計の誤差は実質的に であり、 に依存しません。 ${\displaystyle n}$${\displaystyle O(n\varepsilon ^{2})}$${\displaystyle n}$${\displaystyle \epsilon }$${\displaystyle n\varepsilon ^{2}}$${\displaystyle n}$${\displaystyle 1/\varepsilon }$${\displaystyle n}$${\displaystyle 10^{16}}$${\displaystyle O(\varepsilon )}$${\displaystyle n}$

それに比べて、単純な加算（単純に数字を順番に追加し、各ステップで丸めを行う）の相対誤差境界は、条件数で乗じた分だけ大きくなります。^[2] ただし、この最悪のケースの誤差は実際にはめったに見られません。丸め誤差がすべて同じ方向にある場合にのみ発生するためです。実際には、丸め誤差がランダムな符号を持ち、平均がゼロであるためランダムウォークを形成する可能性の方がはるかに高く、この場合、単純な加算の二乗平均平方根相対誤差は、条件数で乗じた分だけ大きくなります。^[9] ただし、これは補正加算よりもはるかに悪いです。ただし、合計を 2 倍の精度で実行できる場合は、 がに置き換えられ、単純な加算の最悪のケースの誤差は、元の精度での補正加算の項に匹敵します。 ${\displaystyle O(\varepsilon n)}$${\displaystyle O\left(\varepsilon {\sqrt {n}}\right)}$${\displaystyle \epsilon }$${\displaystyle \varepsilon^{2}}$${\displaystyle O(n\varepsilon ^{2})}$

同様に、上記に現れる は、すべての丸め誤差が同じ符号を持つ（かつ最大限可能な大きさを持つ）場合にのみ発生する最悪の境界です。^[2] 実際には、誤差がランダムな符号を持つ可能性の方が高く、その場合 の項はランダムウォークに置き換えられ、その場合、平均がゼロのランダム入力であっても、誤差は（項を無視すると）合計の増加率と同じ割合でのみ増加し、相対誤差が計算されるときに要因がキャンセルされます。そのため、漸近的に悪条件の合計であっても、補正された合計の相対誤差は、最悪のケースの分析が示唆するよりもはるかに小さくなることがよくあります。 ${\displaystyle \Sigma |x_{i}|}$${\displaystyle E_{n}}$${\displaystyle \Sigma |x_{i}|}$${\displaystyle E_{n}}$${\displaystyle O\left(\varepsilon {\sqrt {n}}\right)}$${\displaystyle n\varepsilon ^{2}}$${\displaystyle S_{n}}$${\displaystyle {\sqrt {n}}}$

Neumaier ^[10]は、 Kahanアルゴリズムの改良版を導入しました。彼はこれを「改良Kahan-Babuškaアルゴリズム」と呼んでいます。これは、次に加算される項の絶対値が累積和よりも大きい場合もカバーし、大きいものと小さいものの役割を効果的に入れ替えます。疑似コードでは、アルゴリズムは次のようになります。

function KahanBabushkaNeumaierSum(input)
     var sum = 0.0
     var c = 0.0 // 失われた下位ビットの実行中の補正。

    for i = 1 to input.length do 
        var t = sum + input[i]
         if |sum| >= |input[i]| then 
            c += (sum - t) + input[i] // sumの方が大きい場合、 input[i]の下位桁は失われます。
         else c += (input[i] - t) + sum // それ以外の場合、 sum
            の下位桁は失われます。
         endif
        合計 = t
    次に私

    return sum + c // 修正は最後に一度だけ適用されます。

この機能強化は、Fast2Sum で書き直された Kahan のアルゴリズムで Fast2Sum を2Sumに置き換えることに似ています。

多くの数列では、両方のアルゴリズムは一致しますが、Peters ^[11]による簡単な例では、それらがどのように異なるかを示しています。倍精度で合計する場合、Kahanのアルゴリズムは0.0を生成しますが、Neumaierのアルゴリズムは正しい値2.0を生成します。 ${\displaystyle [1.0,+10^{100},1.0,-10^{100}]}$

より精度の高い高次の修正も可能です。例えば、Klein ^[12]が提案した変形は、 2次の「反復Kahan-Babuškaアルゴリズム」と呼ばれています。疑似コードでは、アルゴリズムは次のようになります。

関数KahanBabushkaKleinSum(入力)
     var sum = 0.0
     var cs = 0.0
     var ccs = 0.0
     var c = 0.0
     var cc = 0.0

    i = 1 の場合、 input.lengthに対してvar t = sum + input[i]
        を実行し
        、 |sum| >= |input[i]|の場合、
            c = (合計 - t) + 入力[i]
        それ以外
            c = (入力[i] - t) + 合計
        終了
        合計 = t
        t = cs + c
        |cs| >= |c |ならば
            cc = (cs - t) + c
        それ以外
            cc = (c - t) + cs
        終了
        cs = t
        ccs = ccs + cc
    ループ終了 

    合計 + cs + ccs を
返す

カハンのアルゴリズムはn個の数を合計すると誤差が増大しますが、ペアワイズ加算では誤差増大はそれよりわずかに悪くなります。ペアワイズ加算では、数のセットを再帰的に2つに分割し、各半分を合計してから、2つの合計を加算します。 ^[2] この方法には、単純な加算と同じ数の算術演算しか必要とせず（4倍の算術演算を必要とし、単純な加算の4倍のレイテンシを持つカハンのアルゴリズムとは異なります）、並列に計算できるという利点があります。再帰の基本ケースは、原理的には1つ（または0）の数の合計にすることができますが、再帰のオーバーヘッドを償却するには、通常、より大きな基本ケースを使用します。ペアワイズ加算に相当するものは、多くの高速フーリエ変換（FFT）アルゴリズムで使用されており、それらのFFTでの丸め誤差の対数的増大の原因となっています。^[13]実際には、ランダムな符号の丸め誤差により、ペアワイズ和の二乗平均平方根誤差は実際には として増大します。^[9]${\displaystyle O(1)}$${\displaystyle O(\log n)}$${\displaystyle O\left({\sqrt {\log n}}\right)}$

もう一つの選択肢は、任意精度演算を使用することです。これは原理的に丸めをまったく必要としませんが、計算コストがはるかに高くなります。任意精度を使用して正確に丸められた合計を実行する方法は、複数の浮動小数点コンポーネントを使用して適応的に拡張することです。これにより、高精度が必要とされない一般的なケースで計算コストが最小限に抑えられます。^{[14] [11]} 整数演算のみを使用し、大きなアキュムレータを使用する別の方法は、KirchnerとKulischによって説明されました。^[15]ハードウェア実装は、Müller、Rüb、およびRüllingによって説明されました。^[16]

原理的には、十分に積極的な最適化コンパイラは、カハン加算の有効性を破壊する可能性がある。例えば、コンパイラが実数演算の結合規則に従って式を簡略化すると、シーケンスの2番目のステップが「簡略化」される可能性がある。

t = sum + y;

c = (t - sum) - y;

に

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

そして

c = 0;

したがって、誤差補正が不要になります。^[17] 実際には、多くのコンパイラは、コンパイラオプションで「安全でない」最適化を有効にするように明示的に指示されない限り、結合規則（浮動小数点演算では近似値に過ぎない）を簡略化に使用しません。^{[18] [19] [20] [21]}ただし、Intel C++コンパイラは、結合規則に基づく変換をデフォルトで許可する一例です。^[22] Cプログラミング言語の 元のK&R Cバージョンでは、コンパイラが実数演算の結合規則に従って浮動小数点式の順序を並べ替えることができました。しかし、その後のANSI C標準では、Cを数値アプリケーションにより適したものにするために（並べ替えを禁止しているFortranに似せるために）並べ替えが禁止されました。 ^[23]ただし、実際には、前述のように、コンパイラオプションで並べ替えを再度有効にすることができます。

このような最適化をローカルで抑制する移植可能な方法は、元の定式化の行の 1 つを 2 つのステートメントに分割し、中間生成物の 2 つをvolatile にすることです。

関数KahanSum(入力)
     var sum = 0.0
     var c = 0.0

    i = 1からinput.lengthまで、 var y = input[i] - c
         volatile var t = sum + y
         volatile var z = t - sum を実行します。
        
        c = z - y
        合計 = t
    次に私

    合計
を返す

一般に、コンピュータ言語に組み込まれている「合計」関数は、特定の合計アルゴリズムが使用されることを保証するものではなく、カハン合計はなおさらです。^{[引用が必要]}線形代数サブルーチンの BLAS標準では、パフォーマンス上の理由から、特定の計算順序を明示的に強制しないようにしており、^[24] BLAS実装では通常、カハン合計は使用されません。

Pythonコンピュータ言語の標準ライブラリでは、正確な合計を求めるためにfsum関数が指定されています。Python 3.12以降では、組み込みの「sum()」関数はノイマイヤー合計を使用します。^[25]

Julia言語では、関数のデフォルトの実装は、良好なパフォーマンスで高精度のペアワイズ合計sumを行いますが^[26] 、外部ライブラリは、より高い精度が必要な場合のために名付けられたNeumaierの変種の実装を提供します。 ^[27]sum_kbn

C#言語では、HPCsharp NuGetパッケージは、スカラー、 SIMDプロセッサ命令を使用したデータ並列、並列マルチコアの両方として、Neumaierバリアントとペアワイズ合計を実装します。 ^[28]

• 分散を計算するアルゴリズム（安定和を含む）

• 浮動小数点加算、ドブ博士のジャーナル、1996 年 9 月