Fast-Mathに注意！浮動小数点演算の危険な最適化

Trapping mathがいい方法だよ。cosmopolitanのリポジトリに使い方を例として書いたから見てみて。URL: https://github.com/jart/cosmopolitan/blob/master/examples/tr…

なんかさ、本当の問題はIEEE規格自体にある気がするんだよね。多くのコードには関係ない細かいルールが多すぎるし、opt-outじゃなくてopt-inにすべきだったよ。ハードウェアの現実に合わせて規格も進化しないと。auto-vectorizeできないのは今のハードウェアじゃ致命的すぎ。でも決定論を犠牲にするのもキツいな。OpenCLとかCUDAはどうなんだろ？命令順序を保証して、しかもGPUでちゃんと並列化できる方法あるのかな？

IEEE 754ってどうやってauto-vectorizeを妨げるの？

浮動小数点演算は交換法則も結合法則も成り立たないから、やらない方がいいよ。

for x in array { sum += x }みたいなループを書いたらさ、それは要素を「この順番で」「一つずつ」足したいっていうプログラムの仕様なんだよ。Vectorizationするとその順番が変わっちゃうんだ。

なんで交換法則成り立たないの？

IEEE-754の細かい要件は個々のプログラムには重要じゃないかもしれないけど、ある程度ちゃんと定義されたセマンティクスを数字に持たせたいなら、どんな数値形式でもそういう欠点は出てくるよ。
GPUコードはね、基本的には各要素に個別作用するスカラーコードとして書かれるんだ。ハードウェアが全体としてループする感じ。マルチスレッドと同じで順序は保証されない。だから、操作の並べ替えが必要なVectorizeコードは物理的に書けないようになってる。「リストの各要素をアキュムレータに足す」みたいなコードは書けないんだよ（書けても非Vectorized CPUより遅くなるけどね）。

実際にはIEEE-754的には交換法則は成り立つんだ。ただ、NaNの結果になる場合は違うNaNの表現になるかもしれないけどね。

うん、浮動小数点数の不正確さがあるから、「(a + c) + (b + d)」が「a + b + c + d」と同じだと決めつけられないんだよね。これのせいで、いつかfixed pointとかbignumの実装の方が速くなったら皮肉だなあ。

Auto-vectorizeできないのはIEEE規格のせいじゃなくて、操作の順序が関係ないことを表現する方法がないプログラミング言語のせいだよ。人気の言語のほとんどは、必要ない場所でもシーケンシャルなセマンティクスを強制する欠陥があるね。
IEEEの浮動小数点標準は史上最も有用な規格の一つだよ。ハードウェア設計者や素人プログラマはベンチマークで良い結果を出すためにズルしたがるインセンティブがあるからさ。でも、何よりも正しい結果が必要なユーザーもいるし、それは生死に関わることさえある。正確さが重要じゃない用途、つまり主にグラフィックスやML/AI向けには、GPUやNPUみたいな専用アクセラレータを使うのがベター。汎用CPUでは、IEEE規格に完全準拠しないのは深刻な間違いだよ。
CUDA、OpenMPなんかについて言うと、定義上、並列化可能な操作なら実行順序は関係ない。順序が重要なら、どんなプラットフォームでも結果を保証するのは不可能だ。
Vectorizedコードが欲しいなら、CとかC++みたいな言語に頼るべきじゃないね。OpenMP、CUDA、OpenCLみたいな目的のために開発された拡張機能を使うべきだよ。そこではVectorizationは偶然に任されないからね。

CUDAの「カーネル」ってさ、1963年からある「parallel DO」とか「parallel FOR」と結局同じことだよね。NVIDIAはキーワードを変えたり、ループの本体と範囲を別々に書かせたりして分かりにくくしてるだけ。数十年前からある用語を新しい言葉で大量に作り直してて、GPUの仕組みを煙に巻こうとしてるみたいで超うざい。さらにAMDもNVIDIAを真似て、また違う独自の用語を作ってるし。

IEEEの仕様自体が問題って感じる人もいるけど、そうじゃないんだ。問題は浮動小数点コードが精度エラーに超敏感なこと。仕様に厳密に従っても精度エラーは直らないけど、エラーが出た時の挙動は決定的になる。これで9割方「チューニングの問題」として無視できちゃうんだけどね。でも精度エラーはバグなんだ。バグは直すべきで、決定論的なトリックで無視しちゃダメ。でもこれってデザインとか難しい数学（gimbal lockみたいに、直すにはquaternionsとか理解する必要があったりして大変！）が絡むから難しいんだよね。だからみんな対処してるだけ。個人的には–ffast-mathは良い面の方が大きいと思うし、プロジェクトでは絶対に有効にすべき。だってこれで見つかる「問題」は、どっちみち直すべきバグなんだから。

複数のNaNがあるのに、その振る舞いの仕様がちゃんと決まってないのって、なんか無駄なミスって感じだよね。

xargsもparallel forと同じことやってるよね。OFCourse、Forthの人たちなら、こんなの朝飯前だったかもね。

数学的な用途なら、NaNのペイロードなんてどうでもいいし、挙動も同じ（quiet vs. signaling NaNsは除く）。等値比較でもNaNは常に等しくないから関係ないよ。

ほとんどのプログラミング言語は、必要ないのに逐次実行のセマンティクスを強制する欠点があるよね。Occamみたいにそうじゃない言語もあったけど、広まらなかった。メモリ順序付け関数でプログラマーがアトミック操作の優先順位を指定できるみたいに、数学とか一般的な式の集合に対しても同じことできないのかな？

それはそうだよ。486SXの固定小数点と486DXの浮動小数点を比べてみてよ。浮動小数点よりも、スケール変換とか合計とか必要な精度で表示するとかの方が固定小数点の方が速いんだぜ。

（gimbal lockについて）Euler anglesでも月に行ったんだし、他の方法でgimbal lockを避けることもできるんじゃない？本物の機械式ジンバルを使ってる時は、あまり選択肢がないだろうけどね。

順序を変えても違いは小さいって言うのは技術的には正しいけど、それは部分的な話だよ。matrix multiplicationみたいに算術が非可換な場合、もっと大きな違いが生じることがあるんだ。

それは今のアーキテクチャでも同じなの？例えばSIMDでも固定小数点ってあるの？

加算の順序を変えると計算結果が変わるかも。だから自動で並列処理する（オートベクタライゼーション）のが難しいんだよね。

順序変更で”小さな差”しか生まない、なんて思っちゃダメ。大きい値に小さい値を何度も足すと、とんでもない差になることもあるからね。例えば、huge_float_value に .01 を10億回足しても、結果が huge_float_value と同じになる、なんてことはないよ。具体的な状況を知らないと、その差が許容できるかは言えないね。

浮動小数点演算は結合法則は成り立たないけど、加算や乗算は交換法則（x + y == y + x）は成り立つよ。NaNの場合は例外的に引数の順序で戻り値のNaNが変わることがあるけど、仕様として保証されてないから誰も使わない。IEEE 754で決まったNaN値を用意して、それ以外のNaN（ペイロード付き）をちゃんと伝搬させるようにすれば、デバッグ情報とか入れられて便利なのにね、って思う。

そうそう、IEEE-754の足し算と掛け算は交換法則は成り立つんだよね。でも分配法則は成り立たないよ。

こういうのって、特別な型（例えばC++なら associative::float みたいに）とか、変数に属性を付けて「これは結合法則が成り立つと思っていいよ」ってコンパイラに教えられるようにすればいいんじゃない？

もっとデカい問題は、プログラミング言語が「この計算の順序は大事！」とか「順序はどうでもいいよ」って作者の意図を伝えられないことなんだよね。IEEE specでも言語でちゃんと定義してねって言ってるらしいし。ほとんどの人は順序なんて気にしないし、コンパイラが勝手に最適化して結果が変わっても驚く（し、むしろ最適化してほしいと思う）。でも、順序が超重要なケースは、すごく少ないけど本当にあるんだよ。

IEEE specは自動並列化（オートベクタライゼーション）を禁止してないんだよ。ただ、「再現性のある浮動小数点の結果」を出すって決めた言語は、自動並列化を避けるべきだって言ってるだけ（IEEE 754-2019の11章）。どう避けるかは実装次第。C++みたいにデバッグモードでシリアルに計算するなら、ユーザーが明示的に許可しない限り、自動並列化で結果が変わっちゃうのは避けるべきだね。

ユーザーから見ればNaNペイロードの扱いはそんなに問題ないんだけど、コンパイラには複雑なんだよ。LLVMだとNaNを比較するのすら未定義動作（UB）にしてて、ほとんどの浮動小数点演算が非決定論的になっちゃう。CPUのハードウェアにも余計な回路が必要になるしね。Positっていう新しい数値形式は、NaNが一つしかないとか、順序付けのルールがちゃんとしてるところは絶対正しいと思うな。

逐次実行が必要ない計算を、言語構文でどう表現するか。例えば、セミコロン区切りは順に実行、カンマ区切りは順不同か並列でもOK、みたいにね。多くの言語で関数の引数はカンマ区切りで評価順は任意だよね。C言語のカンマ演算子も最初はそうなりかけたけど、for文で使うためにセミコロンと同じにしちゃった。式の場所にブロックを書けるようにして、最後の式の値をそのブロックの値にする、みたいにした方が絶対良かったのにね。

精度にこだわるなら、そもそも浮動小数点数は使っちゃダメでしょ。用途によるけどね。標準規格が悪いか言語が悪いかの話はどっちでもいい。標準の目的は計算順序を定めて予測可能にすることだけど、それだと最適化できないのが問題なんだよね。精度が超大事なケースはまれだから、むしろ最適化をデフォルトにして、精度モードはオプトインにすべきだと思うよ。IEEE 754モードは必要だけどね。でも、今の大量並列処理時代には丸め誤差は受け入れるのが普通。それなのにIEEE 754にベッタリな言語が多いのはイケてないと思うな。浮動小数点数の型を増やすとか、コード内で指定するとか、エラー許容度を書くとか、改善策はいろいろありそう。

それは「チューニング」って感じの解決策だね。つまり「それはやっちゃダメ」って制限するか、ダメなら復旧する（リセットとかバックアップとか）方法。これで「まあうまくいく」けど、バグ自体は残ったままなんだよね。ソフト開発ならもっと頑丈なやり方が良いと思う。アポロ13号事故でジャイロ校正失って、太陽見て復旧した話がまさにそうだったのを思い出したよ。NASAは今思えばIEEE doublesとかクォータニオン選べばよかったんだろうね。

「浮動小数点数ほど私の心に恐怖をもたらすものはない。」ってGerald Jay Sussmanが言ってたらしいよ。
IEEEの標準化委員会は、UnumとかPositみたいな浮動小数点数の代替案についても検討してるのかな？[1],[2]
[1] Unum & Posit:https://posithub.org/about
[2] The End of Error:https://www.oreilly.com/library/view/the-end-of/978148223986…

UnumとかPositみたいな新しいやり方がCPUでハードウェアサポートされないのは、なんでだろうね？って思ってる。

Gustavsonの最後のプレゼンでは、Positをサポートする実際のハードウェアを手に持ってるところから始まるよ。知らんけど。
https://m.youtube.com/watch?v=vzVlQhaAZtQ

HPは1970年代からちゃんとした決定論的な10進数算術を持ってたんだぜ。

これ皮肉？そうじゃないなら、提案されてるposit標準、IEEE P3109だよ。

今のP3109ドラフトにはposit入ってないってさ。

PyTorchでもApple MPSでdeep learningとかでFast Mathがデフォルトで有効になってて問題起こるらしいよ。
autoregressive image generation modelの学習中にゴミみたいな出力になったとか。ここに議論があるよ: https://github.com/pytorch/pytorch/issues/84936

この話題、以前にもここで議論されてるよ: https://news.ycombinator.com/item?id=29201473 (記事の最後にもリンクあるけどね)

ForthやLispでは固定小数点とか有理数を使う哲学があるんだ。Forthの記事やSchemeでのexact-＞inexact変換の例を示してるよ。有理数でほとんどの計算は足りるし、円周率だって113/355で十分正確だったりするって話。Forthの記事リンク: https://www.forth.com/starting-forth/5-fixed-point-arithmeti…

Floatsって対数空間での固定小数点みたいなもんだよ。指数部があるから調整が必要になるけどね。

いやいや、全然違うよ。整数をスケールした方がずっと速いし、64bitあればどんな分野でも計算できるってば。

有理数なんて三角関数とか平方根が出てきたらもう使えなくなるじゃん…

浮動小数点数は有理数に戻せるんだぜ。(rational (sqrt 2d0))みたいにね。あるいは、自分で好きな精度で計算する処理を書くこともできるよ。

俺の記事でも不正確→正確変換は扱ったぜ。Schemeによってはできるけど、s9 Schemeは無理。Forthだと自分でやる必要があるんだ。三角関数ならCORDICってのを確認してみてくれよ。https://en.wikipedia.org/wiki/CORDIC

有理数で高精度な三角関数が必要なら、テイラー級数ライブラリを書けばいいじゃん。多項式近似とかルックアップテーブル、CORDICとかも使えるし。

CORDICをチェックしてくれよ。https://en.wikipedia.org/wiki/CORDIC
sqrt関数ならForthにもfsqrtがあるんだぜ。EForthとかPFE Forthの例を見てくれ。EForthの精度は低いけど、80年代のエンジニアが使ってたレベルからは遠くないんだ。

この記事のFTZの実装部分は面白かったな。“FTZがどう実装されてるかって問題は…有効にすると同じスレッドの他の操作全部に影響するんだ。” “GCCで-funsafe-math-optimizationsを使うとFTZが有効になる…共有ライブラリをロードするだけで結果が変わるってのは、楽しいデバッグ経験だね。”

この記事、科学技術ソフトでも問題の重要性を過大評価してると思うな。俺が書いた科学コードだと、ノイズの方が大きいし、fast-mathで性能が10%以上上がることもよくあるんだ。-fassociative-mathの話もそう。最適化で結果が変わるなら、シミュレートする数式自体を見直すべきでしょ。でも、ライブラリとか数学シミュレーションだと違うかもしれないけどね。

「計算の順番重要、この順番でやってくれ」みたいな構文があればいいのにね。そうすれば、指定した所以外は全部fast-mathにできるじゃん。

ほとんどの言語にあるんじゃない？普通に書けば言語仕様通りの順番だし、違う順番にしたいなら括弧を使えばいいでしょ。fast mathとどう関係するかはよく分かんないけど、俺はCを使わないからね。

それは別の種類の順序の話だよ。Pythonのsum()みたいな関数を想像してみて。本来ならリストの中身はどんな順番で足してもいいはずだよね？ でも、IEEE-754浮動小数点数だと、単純に1つずつ足した時と同じ結果は期待できないんだ。足す順番によって結果が変わる可能性がある。ここで話してるのはそういう順序のことだよ。

記事でもGCCとClangに拡張があるって言ってたね。でも、言語自体にあると便利だよね。Zigはそういうアプローチをとってるよ。