2016/08/03(水)半精度浮動小数点数に関する思考実験

その他精度保証
半精度浮動小数点数というものがあります。よく使われている単精度(float, 32bit)、倍精度(double, 64bit)に対して、全長16bitと単精度の半分で浮動小数点数を表現するものです。IEEE754-2008でbinary16としてフォーマットが定められています。deep learningの隆盛とともに「精度が低くてもとにかく速く」計算するニーズが高まり、GPUでハードウェアサポートされるなど、最近注目を集めています(ような気がします)。

IEEE754-2008の半精度では、16bitを符号s(1bit)+指数部e(5bit)+仮数部m(10bit)に分割しています。指数部のオフセットは15で、従って正規化数は
x = (-1)s × 1.m × 2e-15
のように、非正規化数は
x = (-1)s × 0.m × 2-14
のように実数xと対応します。

ところで、URRという浮動小数点数の表現形式をご存知でしょうか。浜田穂積先生が80年代(IEEE754制定より前!)に提案された浮動小数点数の表現形式です。詳細は
を見ていただくとして、簡単に言えば、指数部と仮数部の区切りを可変にし、1に近い数(=指数部を表現するのに必要なbit数が少ない)ときには指数部を短くして仮数部を長くして精度を稼ぎ、非常に小さい数や非常に大きい数を表現するときには仮数部の長さを犠牲にして指数部に長いbitを割り当てる、というものです。このとき、何も考えずに指数部と仮数部を結合してしまうとその区切りが分からなくなってしまいますが、そこは指数部を表現するのに「bit列の末尾が分かるような自然数の表現方法」を用いることで解決します。例えば、Eliasのガンマ符号やデルタ符号といった符号化の方法がよく知られています。

さて、半精度浮動小数点数は、bit数が少ないこともあって表現できる数値の範囲が非常に狭く、簡単にアンダーフローやオーバーフローを起こしてしまいます。正の最大数は何と65504です。正の最小数は、精度を保っている正規化数で2-14≃6.1×10-5、非正規化数まで考えても2-24≃5.96×10-8にすぎません。

そこで、URR的な考え方を用いて16bit浮動小数点数を構成したらどうなるか考えてみました。URRは-infやNaNが無いなど、現代のIEEE754に慣れた我々には使いにくいところもあるので、指数部と仮数部の区切りを可変にするという思想はそのままで、適当にフォーマットを定めます。指数部は、Eliasのデルタ符号を用いることにします。デルタ符号は1,2,3,…の自然数しか表せないので、指数部とデルタ符号で表す数値を
デルタ符号12345…255256…508509510511
指数部0-11-22…127-128…-254±0±infNaN
デルタ符号の長さ13355…1415…15151515
のように対応させることにしました。指数部の最後の3つを特殊な数に割り当てています。仮数部は、IEEE754に倣って先頭の1を格納しない「ケチ表現」にします。このフォーマットとIEEE754-2008のbinary16で、指数部と仮数部の長さの関係を表にしてみます。
提案方式の仮数部長IEEE754-2008の仮数部長
2-2541-
⋮⋮⋮
2-1281-
2-1272-
⋮⋮⋮
2-2461
2-2362
2-2263
⋮⋮⋮
2-1669
2-15710
2-14711
⋮⋮⋮
2-8711
2-7811
⋮⋮⋮
2-4811
2-31111
2-21111
2-11311
201511
211311
221111
231111
24811
⋮⋮⋮
27811
28711
⋮⋮⋮
214711
215711
2166-
⋮⋮⋮
21272-
21281-
⋮⋮⋮
22531-
これを見ると、1付近では仮数部が長くなり、1から離れると徐々に仮数部が短くなっていき、(精度は低いものの)小さな数から大きな数まで表現できていることが分かります。全長16bitなどという極端に厳しい場面でこそ、このようなフォーマットが生きると思うのですがいかがでしょうか。

もちろんハードウェアのサポートが無くソフトウェアエミュレーションでは速度は絶望的ですが、将来このような優れたフォーマットが気軽に使えるようになればいいなと思っています。FPGAとかで作って遊んだりできないかなあ。
コメント（3件）

2016/06/20(月)方向付き丸めクイズ

精度保証
唐突ですがクイズです。a, b, c, dはIEEE 754に従う浮動小数点数とします。以下の計算をIEEE754の+∞方向への丸め(上向き丸め)で行い、計算値をx, 真値をx*とします。このとき、「必ずx ≥ x*が成立する」と言えるのはどれか。○と☓で答えて下さい。但し、計算中にゼロ除算が起きるケースは除外してよいものとします。
  1. x = a + b
  2. x = a - b
  3. x = a × b
  4. x = a / b
  5. x = (a + b) + c
  6. x = (a × b) × c
  7. x = (a - b) - c
  8. x = a - (b + c)
  9. x = -a + b
  10. x = -( (-a) + (-b) )
  11. x = a × b + c × d
  12. x = (a + b) × (c + d)
  13. x = a / (b + c)
  14. x = a - b × c
  15. x = a + (-b) × c
  16. x = sqrt(a)
  17. x = exp(a)
少し考えれば答えは分かると思いますが、これらの答えが全て○になると誤解してるのではないかという事例にたまに遭遇しますので、こういうクイズにも意味があるのではないかと思い、記事を書いてみました。
コメント（0件）

2016/05/09(月)kv-0.4.36

精度保証
というわけで、予告通りkvライブラリを0.4.36にアップデートしました。

内容は前回の記事の予告通りで、maffine3を正式アルゴリズムに昇格させました。
  • maffine → 一応maffine0として残したが消滅候補
  • maffine2 → そのまま
  • maffine3 → maffineに改名
使い勝手は0.4.34以前と変わらず、単に性能がよくなっているはずです。
コメント（0件）

2016/05/08(日)新しいODE solverの性能評価

精度保証
さて、前回の記事で書いた新しいアルゴリズムの性能を簡単な例で示します。例として選んだのはvan der Pol方程式です。
x''- μ(1-x2)x'+x = 0
この方程式はμ≥0の値によって計算のしやすさが変化し、μが1くらいまでなら計算しやすいnon-stiffな方程式ですが、μが大きくなると極めて計算しづらいstiffな方程式になります。ここでは、μ=1とμ=100でこの方程式の軌道を計算し、ステップ幅の変化を調べてみました。ステップ幅はある方法で局所誤差がmachine epsilon程度になるように計算していますが、stiffなODEの場合そこを変えて大きなステップ幅にしても精度保証の根幹となる不動点定理が成立しなくなってしまい小さなステップ幅への修正を強いられてしまいます。初期値は
x(0) = x'(0) = 1
とし、t=200まで計算しました。この軌道を相図(phase diagram)で示すと、次のような感じです。まずμ=1の場合:
dense1.png

周期解(limit cycle)に巻き付いている様子がよく分かります。次にμ=100の場合:
dense100.png

こちらは急激な変化を含んだ周期の長い軌道になります。

さて、これを、
  • ode-maffine (従来の標準的なアルゴリズム)
  • ode-maffine2 (従来の高速アルゴリズム。初期値に関する微分が出来ない欠点がある)
  • ode-maffine3 (新しいアルゴリズム)
  • awa (Lohnerによる有名な精度保証プログラム)
で計算して比較してみます。まずμ=1で:
step1.png

横軸は時刻tで、縦軸は上のグラフがx、下のグラフがODE Solverが選んだステップ幅です。non-stiffな方程式なので、どれでもあまり大きな違いがないことが分かります。計算時間は、
アルゴリズム計算時間(sec)
maffine2.276
maffine21.240
maffine31.798
awa3.815
でした。次に、μ=100の場合を示します。
step100.png

どの方法もμ=1の場合に比べると刻み幅が小さくなっています。変化の激しさに応じて刻み幅を頑張って適応させていることがよく分かります。刻み幅はawaが一番小さく、次にmaffine、maffine2とmaffine3はほぼ同じで比較的大きい刻み幅で計算できていることが分かります。結果として、計算時間にも次のように大きな差が出ました。
アルゴリズム計算時間(sec)
maffine53.820
maffine29.287
maffine313.478
awa176.78
maffine3はmaffineと完全に同じことが出来て性能が良さそうなので、maffineを残す意味はなさそうです。次期versionでは、
  • maffine → 消滅 (別名で残すかも)
  • maffine2 → 従来通り
  • maffine3 → maffineに改名
としたいと思います。これなら0.4.34以前と使い勝手は変わらずに単に性能が良くなったということになります。関数名の変更はユーザに混乱をもたらすので、現version (0.4.35) はなるべく短命にするつもりです。
コメント（0件）

2016/05/06(金)kv-0.4.35

精度保証
kvライブラリを0.4.35にアップデートしました。

まず、前回(0.4.34)に混ぜてしまったバグの修正です。1.00e-05みたいな数値が1.00e0-5のように表示されてしまうというバグを直しました。

また、常微分方程式のsolverに新しいアルゴリズム(maffine3)を追加しました。kv-0.4.30とstiffなODEで日記に書いたように、maffineはstiffなODEに弱く、maffine2はstiffに強いが初期値に関する精度保証された微分が得られないという状況でしたが、今回、「stiffに強くて初期値に関する微分もできる」maffine3を追加しました。今年のGWはほぼこの実装に費やしてしまいました。

実装したアルゴリズムのアイデアは極めて単純です。常微分方程式の初期値問題を精度よく解くためには、いわゆるWrapping Effectを防ぐために解の初期値に関する微分を利用することがほぼ必須です。それを計算するのに必要なのが「初期値に関する変分方程式」というものです。すなわち、
x'(t)= f(x(t),t)
x(t0) = x0
に対して、x*(t)をこの式の真の解としたy(t)に関する行列微分方程式
y'(t) = fx(x*(t), t)y(t)
y(t0) = I (単位行列)
を初期値に関する変分方程式といいます。ある時刻t1におけるyの値y(t1)は、x(t0)にx(t1)を対応させる写像のヤコビ行列と一致します。

このyを計算するのに、従来は上の2つの式を連立させて、
x'(t)= f(x(t),t)
y'(t) = fx(x(t), t)y(t)
x(t0) = x0
y(t0) = I (単位行列)
のようなxとyの組を未知関数とするn+n×n次元の新しい初期値問題を作成し、これを解くという方法を使っていました。その方が実装が楽だったという事情もあります。

それを変更して、まずxだけの式を解いて真の解x*を区間関数として求め、それを代入してyの式を解く、素直な2段階法が、新しいアルゴリズムです。stiffでない普通の方程式だとどちらを使っても大差ありませんが、stiffな場合だと後者の方が刻み幅が大きくなり、結果として高速になるようです。

アルゴリズムの詳細や性能評価はまた改めて書く予定です。問題が無ければ、将来のバージョンではこの新しいmaffine3をmaffineとし、旧maffineは抹消してしまうことも考えています。
コメント（0件）
OK キャンセル 確認 その他