高精度なインクリメンタル統計計算―平均値―

平均値の計算

ある時間関数の平均値を出したいとき，離散空間における数値計算的にはとりあえず算術平均してみることを思いつく．離散時間$t$に対してデータ数$N$の関数$f(t)$では

$$
average = \frac{1}{N}\sum_{t = 1}^N f(t)
$$

しかしそれでは計算時間に比例して丸め誤差が大きくなるため，誤差を減らすためのアルゴリズムを画策する．なお，以降は Fortran 準拠で説明する．また，分散については今回とは別にまとめる予定である．

コンパイラは

$ gfortran -v
...
gcc version 7.4.0 (Ubuntu 7.4.0-1ubuntu1~18.04.1)

単純な例題 (離散時間間隔DT = 1)

もっとも単純な問題の一つとして，ここでは，時間関数を

$$
f(t) = 0.01,\quad 1\le t\le N
$$

とする．すなわち，単精度実数に対して，(時間的には不変となっているが) 時間関数として常に 0.01 を取る関数を考える．この関数の平均値は解析的には任意の N に対して 0.01 である．

N = 10,000,000，離散時間間隔$\delta t = 1$ (DT = 1) で離散化した関数$f$に対して，算術平均でこの関数の平均を出すとすると，その Fortran コードおよび実行結果 (相対誤差) は，

stats_mean_arithmetic.f90

program stats_mean_arithmetic
    use,intrinsic :: iso_fortran_env
    implicit none

    integer(int32),parameter :: N = 10000000
    integer(int32),parameter :: DT = 1 !離散時間間隔

    real(real32) avg_arithmetic
    integer(int32) t

    avg_arithmetic = 0e0

    print*,'avg_arithmetic'
    do t=1,N,DT
        avg_arithmetic = avg_arithmetic+0.01e0
    enddo
    print*,0.01d0/(avg_arithmetic/(N/DT))-1d0 !倍精度相対誤差

end program stats_mean_arithmetic

$ gfortran stats_mean_arithmetic.f90
$ ./a.out
 avg_arithmetic
   4.5140203644990162E-002

なお，Fortran なら総和演算用組み込み関数としてsumが存在するが，ここではアルゴリズムによる丸め誤差の違いについて議論したいので用いない．また，表題にも書かれている $\delta t$ (離散時間間隔) をわざわざパラメータとして与えているのは，そのうち役に立つので実装にだけ気をつけてひとまず脇においておく．

漸近アルゴリズム (Welford's Online algorithm)

この相対誤差4.5140203644990162E-002を改善していく (0に近づける)．といっても，ほとんどこちらの記事に紹介されている方法を実装していくだけ (なんなら「Fortran!?」と思う人達にとっては元記事のほうが優良物件)．まず，某研究室では「漸化式」と呼んでいる漸近平均アルゴリズム (元記事ではさらに誤差を抑えるアルゴリズム，WikipediaではWelford's online algorithmの中で紹介) を実装する．証明は元記事に譲るとして，その数式は関数$x_i$に対して現在の平均値を$m_1$，一つサンプルが増えたときの平均値を$m_2$とすると，

$$
m_2 = m_1+\left(\frac{1}{N+1}\right)(x_{N+1}-m_1)
= \left(1-\frac{1}{N+1}\right)m_1+\frac{1}{N+1}x_{N+1}
$$

なぜ最初にこの漸近アルゴリズムを？というと，今あるデータから常にそのときまでの平均値を得られると，時間方向に値を全て保持する必要がなくなるのでメモリ運用に優しいから．また，これなら後でデータを追加しても大丈夫なので，ある問題の平均値を得るのに十分な時間が分かっていないときに強い．これをインクリメンタルなアルゴリズムという．

stats_mean_recurrence.f90

program stats_mean_recurrence
    use,intrinsic :: iso_fortran_env
    implicit none

    integer(int32),parameter :: N = 10000000
    integer(int32),parameter :: DT = 1 !離散時間間隔

    real(real32) avg_recurrence
    integer(int32) t

    avg_recurrence = 0e0

    print*,'avg_recurrence'
    do t=1,N
        if(mod(t,DT)==0) avg_recurrence = avg_recurrence*(1e0-real(DT)/t)+0.01e0*real(DT)/t
    enddo
    print*,0.01d0/avg_recurrence-1d0

end program stats_mean_recurrence

$ gfortran stats_mean_recurrence.f90
$ ./a.out
 avg_recurrence
  -1.7902774972853197E-002

まぁ，小さくなった．オーダーすら変わらなかったのはちょっと寂しい気もするが...．

漸近+マージアルゴリズム

次に，上記漸近アルゴリズムコードにマージ平均のアルゴリズム (元記事ではマージする) を追加する．マージ平均の利用場面としては，別プロセスで出したそれぞれの平均値を最後にまとめて全体の平均値を得たい場合が考えられる．しかしそれ以外にも，上記漸近アルゴリズムを複数に分割して一定間隔ごとにマージ平均を計算するようにすれば，漸化式に現れる分母tの最大値を小さくすることができ，漸化式の和算におけるオーダーを離れづらくできる．例えば時間間隔$A$ごとにマージ平均することを考えると，

$t$	average
$1\ -\ A$	$avg_\mathrm{merge1} = avg_A$
$A+1\ -\ 2A$	$avg_\mathrm{merge2} = \dfrac{A}{2A}avg_A+\dfrac{A}{2A}avg_{2A} = \dfrac{1}{2}(avg_\mathrm{merge1}+avg_{2A})$
$2A+1\ -\ 3A$	$avg_\mathrm{merge3} = \dfrac{A}{3A}avg_A+\dfrac{A}{3A}avg_{2A}+\dfrac{A}{3A}avg_{3A} = \dfrac{1}{3}(2avg_\mathrm{merge2}+avg_{3A})$
...	...
$(\mathrm{ct}-1)A+1\ -\ \mathrm{ct}A$	$avg_\mathrm{merge,ct} = \dfrac{1}{\mathrm{ct}}((\mathrm{ct}-1)avg_\mathrm{merge,ct-1}+avg_{\mathrm{ct}A})$

したがって，マージ間隔を固定すれば，それぞれの区間で計算される平均 (上記表の$avg_A,avg_{2A},...$) では分母が最大でも$A$に限定される (何もしなければ最大は$\mathrm{ct}A$)．以上の考え方をアルゴリズムとして実装すると，

stats_mean_recurrence_merge.f90

program stats_mean_reccurence_merge
    use,intrinsic :: iso_fortran_env
    implicit none

    integer(int32),parameter :: N = 10000000
    integer(int32),parameter :: DT = 1 !離散時間間隔
    integer(int32),parameter :: M_DT = 50000 !マージ時間間隔

    real(real32) avg_part,avg_merge
    integer(int32) t,t_part,ct

    avg_part = 0e0
    avg_merge = 0e0
    ct = 0

    print*,'avg_merge'
    do t=1,N
        if(mod(t,DT)==0)then
            t_part = t-ct*M_DT !ここでマージごとに次式の分母をリセットできるので和算のオーダーが離れづらくなる
            avg_part = avg_part*(1e0-real(DT)/t_part)+0.01e0*real(DT)/t_part !漸近平均アルゴリズム
        endif

        if(mod(t,M_DT)==0)then
            ct = t/M_DT
            avg_merge = 1e0/ct*((ct-1)*avg_merge+avg_part)
            avg_part = 0e0 !マージごとに個別の平均値はリセット
        endif
    enddo
    print*,0.01d0/avg_merge-1d0

end program stats_mean_reccurence_merge

M_DT = 50000とした場合，

$ gfortran stats_mean_recurrence_merge.f90
$ ./a.out
 avg_merge
   9.8944689321811552E-006

かなり0に近づいた．

漸近+マージ+カハンアルゴリズム

さらにここにカハンのアルゴリズム (元記事では補正加算，Wikipedia ではカハンの加算アルゴリズム) を加えてみる．その全てはWikipediaに譲るとして，えいやと実装すると，

stats_mean_recurrence_merge_Kahan.f90

program stats_mean_recurrence_merge_Kahan
    use,intrinsic :: iso_fortran_env
    implicit none

    integer(int32),parameter :: N = 10000000
    integer(int32),parameter :: DT = 1 !離散時間間隔
    integer(int32),parameter :: M_DT = 50000 !マージ時間間隔

    real(real32) avg_merge
    real(real32) avg_kahan,yk,tk,ck,t_inv
    integer(int32) t,t_part,ct

    avg_kahan = 0e0
    avg_merge = 0e0
    ck = 0e0
    ct = 0

    print*,'avg_kahan'
    do t=1,N
        if(mod(t,DT)==0)then
            t_part = t-ct*M_DT !ここでマージごとに次式の分母をリセットできるので和算のオーダーが離れづらくなる
            t_inv = real(DT)/t_part

            !avg = (avg*(1-t_inv)-ck)+f(t)*t_inv
            !    = avg*(1-t_inv)+(f(t)*t_inv-ck)
            !    = avg*(1-t_inv)+yk
            !    = tk
            yk = 0.01e0*t_inv-ck
            tk = avg_kahan*(1e0-t_inv)+yk
            ck = (tk-avg_kahan*(1e0-t_inv))-yk
            avg_kahan = tk
        endif

        if(mod(t,M_DT)==0)then
            ct = t/M_DT !ct回目のマージ
            avg_merge = 1e0/ct*((ct-1)*avg_merge+avg_kahan)

            avg_kahan = 0e0 !マージごとに個別の平均値はリセット
            ck = 0e0
        endif
    enddo
    print*,0.01d0/avg_merge-1d0

end program stats_mean_recurrence_merge_Kahan

$ gfortran stats_mean_recurrence_merge_Kahan.f90
$ ./a.out
 avg_kahan
   9.8944689321811552E-006

なんと，マージ平均と完全に一致．

一旦まとめ

以上をまとめると，

 avg_arithmetic
   4.5140203644990162E-002
 avg_recurrence
  -1.7902774972853197E-002
 avg_merge
   9.8944689321811552E-006
 avg_kahan
   9.8944689321811552E-006

マージ平均が非常に強力だと分かる．強力すぎてカハンのアルゴリズムでの改善結果が見えなくなっている？

おまけ：漸近+カハンアルゴリズム

試しに漸近+カハン平均を計算してみると，

stats_mean_recurrence_Kahan.f90

program stats_mean_recurrence_Kahan
    use,intrinsic :: iso_fortran_env
    implicit none

    integer(int32),parameter :: N = 10000000
    integer(int32),parameter :: DT = 1 !離散時間間隔

    real(real32) avg_kahan,yk,tk,ck,t_inv
    integer(int32) t,t_part,ct

    avg_kahan = 0e0
    ck = 0e0

    print*,'extra: avg_recurrence+kahan'
    do t=1,N
        if(mod(t,DT)==0)then
            t_inv = real(DT)/t

            yk = 0.01e0*t_inv-ck
            tk = avg_kahan*(1e0-t_inv)+yk
            ck = (tk-avg_kahan*(1e0-t_inv))-yk

            if(t==N) print*,'  ck(t = N) =',ck

            avg_kahan = tk
        endif
    enddo
    print*,0.01d0/avg_kahan-1d0

end program stats_mean_recurrence_Kahan

$ gfortran stats_mean_recurrence_Kahan.f90
$ ./a.out
 extra: avg_recurrence+kahan
   ck(t = N) =  -1.75169546E-10
  -4.6186828883503783E-002

！？算術平均よりもわずかにだが相対誤差が大きくなっている．ここで，t = Nの下位ビット群ckを出力してみると，ck(i = N) = -1.75169546E-10となっていることから，関数のオーダー (1E-2) と比較しても非常に小さくなっており，ckが (少なくとも改善の方向には) 作用しないと考えられる．そうはいっても改悪した理由はイマイチ分からない...漸近とカハンは相性が悪いのだろうか．問題依存の可能性もある．

単純な例題 (離散時間間隔DT = 1000)

さてここで，これまでわざわざパラメータとして与えていた DT を動かしてみる．これは離散時間間隔に相当し，つまりはサンプリング間隔である ($t = 1-N$ですでに元の関数は離散化されているが，それをさらに離散化する)．この時間間隔が例えばある信号に対して広すぎるとサンプリング定理を満たさなくなり，エイリアシング誤差が生じる．今は常に一定値をとる関数なので特に間隔に制限はないが，ひとまずDT = 1000にしてみる．適当にシェルスクリプトを作って実行すると，

$ ./shell_stats.sh
 avg_arithmetic
  -2.9499703615321060E-005
 avg_recurrence
   1.1399401307343737E-006
 avg_merge
   3.9488092751227555E-007
 avg_kahan
   3.9488092751227555E-007
 extra: avg_recurrence+kahan
   ck(t = N) =   3.24916982E-10
   1.1399401307343737E-006

得も言えぬ．算術平均avg_arithmeticの相対誤差が減少したのは，総和のサンプル数が減ったために丸め誤差が生じにくくなったことと，現問題の関数は一定値をとるため離散間隔が広がった影響は平均に現れないことによる．そのせいもあって誤差の改善が微妙な結果になっているが，少なくとも算術平均 → 漸近平均でさらに 1 桁，漸近+マージ平均で 1 桁，相対誤差が小さくなることは分かった．漸近+マージ+カハン平均ではやはり変わらない．また，漸近+カハンは先程と異なり，漸近とまったく同じ結果になっているので，少なくとも算術平均よりかは改善されている．ただ，やはり今回の問題に対しては漸近とカハンは相性が悪いのか，計算コストが増えるだけでカハンの利点が生まれない．

単純な例題を倍精度で

試しに倍精度実数で全ての平均を再計算してみると，$\delta t = 1$ (DT = 1) のとき，

 avg_arithmetic
   1.3690315547876253E-010
 avg_recurrence
   3.4061642395499803E-013
 avg_merge
  -8.8817841970012523E-016
 avg_kahan
   1.5765166949677223E-014
 extra: avg_recurrence+kahan
   ck(t = N) =   1.7732622343059589E-019
  -1.6320278461989801E-014

算術平均に対して相対誤差は，それぞれ 3 桁 (漸近)，6 桁 (漸近+マージ)，4 桁 (漸近+マージ+カハン) だけ改善された．むしろカハンを追加することで相対誤差は改悪された．しかし漸近平均に比べると 1 桁だけ改善されており，漸近とカハンの組み合わせはやはり問題の設定にかなり依存するようである．計算時間とのトレードオフを考えると，カハンを組み合わせないのが無難かもしれない．

$\delta t = 1000$ (DT = 1000) のとき，

 avg_arithmetic
  -1.4255263636187010E-013
 avg_recurrence
   1.1102230246251565E-015
 avg_merge
  -8.8817841970012523E-016
 avg_kahan
  -8.8817841970012523E-016
 extra: avg_recurrence+kahan
   ck(t = N) =  -5.5480658045156672E-020
   5.9952043329758453E-015

なお，これらの改善により多少計算コストは増加するが，DT = 1000の場合は 10,000 回程度の do ループなので (この例題では) さしたる問題ではない．1 ms を争うような計算ならともかく，各アルゴリズムはループを追加したりせず，数個の追加計算や if 文を一つ追加する (if 内の計算頻度は do ループに比べて極小) 程度のコスト増加なので，気にする必要はひとまずない．

CPU時間計測

参考程度に，ある基準時点からのCPU時間を取得する組み込みサブルーチンcpu_timeを使ってループの計算時間を計ってみる．コードは算術平均の場合，

stats_mean_arithmetic.f90

    ...
    real(real32) clock_start,clock_finish
    ...
    print*,'avg_arithmetic'

    call cpu_time(clock_start)
    do t=1,N,DT
        avg_arithmetic = avg_arithmetic+0.01e0
    enddo
    call cpu_time(clock_finish)
    print*,clock_finish-clock_start,'seconds'

    print*,0.01d0/(avg_arithmetic/(N/DT))-1d0 !倍精度相対誤差
    ...

計算条件は，倍精度 (selected_real_kind(15))，N = 100000000，DT = 1，M_DT = 50000として，3回計算した平均をとると (単位は全て秒)，

算術平均	漸近アルゴリズム	漸近+マージアルゴリズム	漸近+マージ+カハンアルゴリズム	extra: avg_recurrence+kahan
0.239583333	1.177083333	1.177083333	1.213541667	1.260416667

まとめ

上記サンプルでは一次モーメント (関数の平均値) の統計計算アルゴリズムの改善を行った．結果としては，単純な算術平均よりも数桁だけ相対誤差を小さくすることができた．もっとも安定的に有効なアルゴリズムは漸近+マージ平均であり，漸近+カハン平均は問題に依存する可能性が得られた．

今後は二次モーメント (関数の分散値) に対する統計計算アルゴリズムの改善を行う．

追記: 三角関数の検証結果について

上記では単純な問題として$0.01$を定数関数として与えていたが，三角関数$1+\sin(x)$を与えて検証した結果を高精度なインクリメンタル統計計算―分散値―に載せている．ただし今回のアルゴリズムを網羅しているわけではない．