kashiの日記

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 × 2^e-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,…の自然数しか表せないので、指数部とデルタ符号で表す数値を

デルタ符号	1	2	3	4	5	…	255	256	…	508	509	510	511
指数部	0	-1	1	-2	2	…	127	-128	…	-254	±0	±inf	NaN
デルタ符号の長さ	1	3	3	5	5	…	14	15	…	15	15	15	15

のように対応させることにしました。指数部の最後の3つを特殊な数に割り当てています。仮数部は、IEEE754に倣って先頭の1を格納しない「ケチ表現」にします。このフォーマットとIEEE754-2008のbinary16で、指数部と仮数部の長さの関係を表にしてみます。

	提案方式の仮数部長	IEEE754-2008の仮数部長
2^-254	1	-
⋮	⋮	⋮
2^-128	1	-
2^-127	2	-
⋮	⋮	⋮
2^-24	6	1
2^-23	6	2
2^-22	6	3
⋮	⋮	⋮
2^-16	6	9
2^-15	7	10
2^-14	7	11
⋮	⋮	⋮
2^-8	7	11
2^-7	8	11
⋮	⋮	⋮
2^-4	8	11
2^-3	11	11
2^-2	11	11
2^-1	13	11
2⁰	15	11
2¹	13	11
2²	11	11
2³	11	11
2⁴	8	11
⋮	⋮	⋮
2⁷	8	11
2⁸	7	11
⋮	⋮	⋮
2¹⁴	7	11
2¹⁵	7	11
2¹⁶	6	-
⋮	⋮	⋮
2¹²⁷	2	-
2¹²⁸	1	-
⋮	⋮	⋮
2²⁵³	1	-

これを見ると、1付近では仮数部が長くなり、1から離れると徐々に仮数部が短くなっていき、(精度は低いものの)小さな数から大きな数まで表現できていることが分かります。全長16bitなどという極端に厳しい場面でこそ、このようなフォーマットが生きると思うのですがいかがでしょうか。

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

コメント（3件）

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

精度保証

唐突ですがクイズです。a, b, c, dはIEEE 754に従う浮動小数点数とします。以下の計算をIEEE754の+∞方向への丸め(上向き丸め)で行い、計算値をx, 真値をx*とします。このとき、「必ずx ≥ x*が成立する」と言えるのはどれか。○と☓で答えて下さい。但し、計算中にゼロ除算が起きるケースは除外してよいものとします。

x = a + b
x = a - b
x = a × b
x = a / b
x = (a + b) + c
x = (a × b) × c
x = (a - b) - c
x = a - (b + c)
x = -a + b
x = -( (-a) + (-b) )
x = a × b + c × d
x = (a + b) × (c + d)
x = a / (b + c)
x = a - b × c
x = a + (-b) × c
x = sqrt(a)
x = exp(a)

少し考えれば答えは分かると思いますが、これらの答えが全て○になると誤解してるのではないかという事例にたまに遭遇しますので、こういうクイズにも意味があるのではないかと思い、記事を書いてみました。

コメント（0件）

2016/05/21(土)カシオの関数電卓の奇妙な挙動

その他

カシオの関数電卓をいじっていて、とても奇妙な挙動に気づいたので記録を残しておきます。使用したのはfx-5800Pという機種で、内部演算精度は10進数15桁です。結論から書くと、次のような現象を発見しました。

123456789123456 - 123456789121111 = 2345
123456789123456 - 123456789123111 = 345
123456789123456 - 123456789123411 = 0 (45でない!)

非常に近い2つの数同士の引き算を行なうと、「桁落ち」と呼ばれる有効数字の減少が発生します。これは有限精度の計算を行っている以上仕方のないことですが、どうやらこの関数電卓では、

「減算で桁落ちが発生した結果有効数字が2桁以下になったとき、その数を強制的に0に書き換える」

という仕様になっているようです。桁落ちが起きるような厳しい計算こそ高精度な電卓の助けが欲しいところなのに、そこでわざわざかろうじて残った情報を捨て去ってしまうというのは極めて不可解な仕様だと感じます。

この仕様で困ってしまう場面はいくらでもありそうですが、一つ例を作ってみました。2次方程式

0.9999996 x² - 2x + 1.0000004 = 0

の解をいわゆる2次方程式の解の公式で計算してみます。判別式の値は、

4 - 4 × 0.9999996 × 1.0000004
= 4 - 3.99999999999936
= 0.00000000000064
= 6.4 × 10^-13

ですが、この値が

4 - 3.99999999999936 = 0

となってしまいます。3.99999999999936を正しく格納できる仮数部を持っているにも関わらず、です。実にもったいない仕様だとは思いませんか。素直な10進数15桁の演算ならば、

sqrt(6.4 × 10^-13) = 8 × 10^-7
x₁ = (2 - 8×10^-7)/2/0.9999996 = 1
x₂ = (2 + 8×10^-7)/2/0.9999996 ≈ 1.000000800000032

のようにほぼ正確に計算できたはずなのに、

x₁ ≈ 1.00000040000016
x₂ ≈ 1.00000040000016

のような重解になってしまいます。また、手作業で解の公式を使って計算した場合のみならず、関数電卓に組み込まれた2次方程式ソルバーを使っても全く同じ計算結果が得られました。すなわち、ユーザに直接見える部分のみならず、内部の演算ルーチンそのものがこの仕様に蝕まれていることは確実と思われます。

この仕様はカシオfx-5800Pに特有のものなのか、他の機種にも見られるものなのか、ヨドバシカメラの関数電卓売り場の展示機を触って調べてみました。すると、

機種名	仮数部の桁数	桁落ち時に捨てられる桁数
カシオ fx-5800P/JP900/FD10Pro/375ES/72F	15	2
キヤノン F-789SG	18	2
キヤノン F-715SA	16	3
シャープ EL-509M	14	1

のようになりました。試した限りの全ての関数電卓で、桁数の違いこそあれこのような「過桁落ち」処理がなされているようです。全く理解できないのですが、どなたか事情をご存知なら教えて欲しいものです。

コメント（6件）

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-x²)x'+x = 0

この方程式はμ≥0の値によって計算のしやすさが変化し、μが1くらいまでなら計算しやすいnon-stiffな方程式ですが、μが大きくなると極めて計算しづらいstiffな方程式になります。ここでは、μ=1とμ=100でこの方程式の軌道を計算し、ステップ幅の変化を調べてみました。ステップ幅はある方法で局所誤差がmachine epsilon程度になるように計算していますが、stiffなODEの場合そこを変えて大きなステップ幅にしても精度保証の根幹となる不動点定理が成立しなくなってしまい小さなステップ幅への修正を強いられてしまいます。初期値は

x(0) = x'(0) = 1

とし、t=200まで計算しました。この軌道を相図(phase diagram)で示すと、次のような感じです。まずμ=1の場合:

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

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

さて、これを、

ode-maffine (従来の標準的なアルゴリズム)
ode-maffine2 (従来の高速アルゴリズム。初期値に関する微分が出来ない欠点がある)
ode-maffine3 (新しいアルゴリズム)
awa (Lohnerによる有名な精度保証プログラム)

で計算して比較してみます。まずμ=1で:

横軸は時刻tで、縦軸は上のグラフがx、下のグラフがODE Solverが選んだステップ幅です。non-stiffな方程式なので、どれでもあまり大きな違いがないことが分かります。計算時間は、

アルゴリズム	計算時間(sec)
maffine	2.276
maffine2	1.240
maffine3	1.798
awa	3.815

でした。次に、μ=100の場合を示します。

どの方法もμ=1の場合に比べると刻み幅が小さくなっています。変化の激しさに応じて刻み幅を頑張って適応させていることがよく分かります。刻み幅はawaが一番小さく、次にmaffine、maffine2とmaffine3はほぼ同じで比較的大きい刻み幅で計算できていることが分かります。結果として、計算時間にも次のように大きな差が出ました。

アルゴリズム	計算時間(sec)
maffine	53.820
maffine2	9.287
maffine3	13.478
awa	176.78

maffine3はmaffineと完全に同じことが出来て性能が良さそうなので、maffineを残す意味はなさそうです。次期versionでは、

maffine → 消滅 (別名で残すかも)
maffine2 → 従来通り
maffine3 → maffineに改名

としたいと思います。これなら0.4.34以前と使い勝手は変わらずに単に性能が良くなったということになります。関数名の変更はユーザに混乱をもたらすので、現version (0.4.35) はなるべく短命にするつもりです。

コメント（0件）