逆数の総和を取る方法 strikes back

# これは何？

以前

https://qiita.com/Nabetani/items/6d328f626776c2531c52

という記事を書いたんだけど、その後知見が増えたので再挑戦。

多倍長整数、四倍八倍精度浮動小数点数、有理数を使うと精度は上がるけど、浮動小数点数の誤差の話をしたいので逆に単精度縛りで。

# この案件の性質

正の数だけを加算しているが、合計は 1〜100 のどこかに着地する。という点が肝だと思う。

合計がある程度見えていることを利用できる。

単精度縛りなのに $2^{-30}$ みたいな小さな値を扱うことになるので、先の方に行くと普通に足したら負け確定。

# 単精度縛りの中で、加算の精度を上げる方法

精度を上げる方法はいくつかある。

## Kahan-Babuska を使う

詳しくは

https://qiita.com/Nabetani/items/59a8ec6160cd409c97f0

を参照のこと。

## 大きな値に対して足さない

逆数の総和というこの文脈では 10 のような大きさの数に対して 百万分の1 のような相対的に小さな数を加算するからたくさん情報落ちする。

浮動小数点数として誤差がなく表現できる値を避難しておいて、端数に対して加算をすれば情報落ちが起きにくくなる。

避難する数値としては、2の整数乗の整数倍がよい。2の-8乗の整数倍とか。

## 両方やる

前述の「Kahan-Babuska」と「大きな値に対して足さない」を両方やればもっと良くなるかもしれない。

Kahan-Babuska は、最終的に 普通の合計と誤差の累計を足したものを返す。
なので、普通の合計と誤差の累計の両者に対して「大きな値を足さない」を実施するのが良さそうだ。

# 測定方法

まず、正解の値の代用として、逆数を float にした値を double の変数に加算していった。単に double だと不安なので、さらに「大きな値に対して足さない」（後述）を仕込んで精度を高める努力をしてみた。

今回は除算で生じる誤差には今回は興味がないので、正解の値を計算するにあたっても逆数を一旦単精度にキャストしてから倍精度の合計に積んでいった。

つまり
$$\sum_{1}^{N}{ float(\frac{1}{n}) }$$
を計算するんだけど、N の値を 2の6乗から2の31乗まで二倍回しで変化させて、正解の代用とする値との差を見た。

ちなみに、 $N=2^{6}$ の場合、総和は 4.7 ぐらい。$N=2^{31}$ の場合、総和は 22 ぐらいになる。

差を見るに当たり、例えば 64件の場合

* 1回目から64回目まで、いちいち正解の代用との差が正解の値に対してどれだけ大きいかを見て、差の2乗の平均の平方根を指標とする
* 最終的な値と正解の代用との差を、正解の代用で割ったもの だけを見る

の２つを調べた。前者を平均誤差。後者を最終誤差と呼ぶ。

# 参加アルゴリズム

## 単純加算

単に足す。
単精度なので、分母が $2^{-23}$ とかになると絶望的。

名前を simple_sum とした。

## 大きな値に対して足さない

合計から基準値の整数倍を減じつつ積算することで誤差を抑える。
基準値として $2^{0}$ を選んだものを keep_small 0、
$2^{-8}$ を選んだものを keep_small 8、
$2^{-16}$ を選んだものを keep_small 16 と呼ぶ。

## Kahan-Babuska

ruby の sum で使われている Kahan-Babuska を使う。

で。
単に Kahan-Babuska を使う、ruby の sum と同等のものを kb と呼ぶ。

Kahan-Babuska の誤差項に Kahan-Babuska を使うものを、kb2 と呼ぶ。

Kahan-Babuska の誤差項に Kahan-Babuska を使い、その誤差項にもKahan-Babuska を使うものを kb3 と呼ぶ。

## 両方使う

こちらの Kahan-Babuska は ruby の sum と同様のシンプルなやつ。

で。

Kahan-Babuska でおこなう２つの累積和の両者に対して「大きな値に対して足さない」の処理を入れた。

「大きな値に対して足さない」で使う基準値として $2^{0}$ を選んだものを mix 0、
$2^{-8}$ を選んだものを mix 8、
$2^{-16}$ を選んだものを mix 16 と呼ぶ。

# 結果

対数目盛にしないと何も見えないんだけど、誤差ゼロだとグラフが書けない。
仕方ないので、誤差ゼロは 2の-40乗 のところにプロットした。

みんなこういうときどうしてるんだろう。

## 大きな値に対して足さない

* 平均誤差
![average error (keep_small).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/fe4c9e33-b1ae-759d-c42b-8ead88c3a0cf.png)

* 最終誤差
![filal error (keep_small).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/1e659ccd-d9a9-e915-de55-098791ea65fc.png)

keep_small 0 でも多少効果はあるけど、やはり基準値を小さくすると強い。

2の-16乗の整数倍を退避すれば、最後 ($ N = 2^{31} $) までほぼ誤差なく加算できる。
退避の基準値は、退避した値で情報落ちしない範囲でできるだけ小さくするのが良い。
今回は合計が高々 2の5乗 ぐらいなので、単精度なら $ 2^{-18} $ ぐらいが限界かな。

## Kahan-Babuska

* 平均誤差

![average error (kb).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/a9292a01-2a39-d2c8-d600-5264a381693b.png)

* 最終誤差

![filal error (kb).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/24017ad0-83b5-054f-ac08-03058576973e.png)

kb3 は、三段目をやっても意味がないよ、ということを言いたくて入れたんだけど、どうも役に立つ場面があるもよう。

単純加算に対する強さは圧倒的だけど「大きな値に対して足さない」に比べると弱い。

kb3 でも 件数が $2^{25}$ を超えてくると情報落ちを避けられない模様。

## 混合戦略

* 平均誤差

![average error (mix).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/8d2791a0-542c-2e2b-930d-6e2c2ecb2292.png)

* 最終誤差

![filal error (mix).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/3245355b-e372-b66d-cc02-29f281f0e539.png)

強い。

## ３つの戦略の頂上決戦

* 平均誤差

![average error (pickup).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/60082604-0d17-8c31-aa34-d4ef78cd8027.png)

* 最終誤差

![filal error (pickup).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/6854/77be55ad-2915-63a2-1723-cbbc5b204413.png)

三段 Kahan-Babuska は $2^{24}$ までは他と互角に戦っているけど、そのあたりから負けてくる。

平均を見ると、混合戦略は他の戦略より優れているようだけど「大きな値に対して足さない」との差はあまり大きくない。

# まとめ

今回は単精度縛りで、倍精度で求めた値との差を見た。

[もともとの課題](https://qiita.com/ttf1998seiya/items/a60ff99c9fdffa09c28f) はそんな縛りはないので、真の値と、コンピュータで現実的な時間で求められる値との差が気になるところではある。

その場合、除算の結果を浮動小数点数に入れた場合に発生する誤差をいかに刈り取るかが勝負になると思うんだけど、今はアイディアがない。

今回導入した
「大きな値に対して足さない」というアルゴリズムは合計（厳密には、合計を求めるにあたって途中で到達する値のピーク）が見えているときには常に有効なんじゃないかと思う。

Kahan-Babuska より計算量が多いけど、わりとおすすめ。

仮平均との差を使って誤差を減らすという作戦と似ているけど、それよりは気楽な感じがする。計算量は多いけど。

混合戦略は精度面ではそうとう強いけど、これが必要になる場面は稀なんじゃないかとも思う。