第14回 計算科学技術特論B (2024)

数値計算における丸め誤差 芝浦工業大学システム理工学部 数理科学科 尾崎 克久 1

自己紹介 • 芝浦工業大学 システム理工学部 • 数理科学科 教授 • 尾崎 克久（おざき かつひさ） • 専門 – 高速・高精度計算 – 丸め誤差解析 – 精度保証付き数値計算 – 数値線形代数 2

講義内容（14回目:2024/7/18） • 数値計算における誤差 – お手頃価格の電卓やスマホによる事例 – MATLABによる事例 • 浮動小数点数とその演算 • 誤差に関する話題 – 数値再現性 • MATLABとOctave，MATLAB内でCPUかGPUか • 組み込み関数の扱い – ソフトウェアによる結果の違い 3

講義内容（15回目:2024/7/25 ） • アルゴリズムと丸め誤差 –２次元凸包の問題 • 精度保証付き数値計算とは • 浮動小数点演算の丸めのモード –簡単な計算機援用証明 • 区間演算の紹介 • その他の誤差の問題の紹介 –三角形分割 • まとめと参考文献 予定です 4

数値計算の誤差（導入・実演） • 電卓・計算機で計算した結果はいつでも正しいと思っていますか? – そうではないことを実演します． – 電卓の事例は「計算機援用証明I」より紹介 • 計算結果があわない例 – 実数とは異なる性質 • ＭATLABを用いたデモですが，数値計算を扱う言語では同じ現象を確認 できるはずです． デモ 5

数値計算による誤差 • C言語やJavaにおいて – 整数を扱うならint型 – 小数点以下も扱うならfloatかdouble型 こっちの理解は微妙 を使うとよいと教わっていませんか? • C言語やJavaにおいて – 整数を扱うならint型 – 実数を扱うならfloatかdouble型 こっちの理解はダメ を使うとよいと教わっていませんか? 6

浮動小数点数 • IEEE 754によって定められる2進浮動小数点数 • 本日は数としての理解を進めるためにbinary32 （単精度浮動小数点数）とbinary64（倍 精度浮動小数点数）に限定し，かつ – 符号部・指数部・仮数部 – バイアス表記 – 正規化数・非正規化数 – 零・Inf・NaN などの説明は割愛します． 数の表現としてのビットパターンについては触れません． 7

浮動小数点数 • Binary64（倍精度浮動小数点数） • 2進数です • C言語やJavaではdouble型（MATLABではデフォルト） 21023 符号 20 21022 2-1 ・・・ 2-1073 2-1074 ・・・ ・ • 使えるのは連続する５３ビット 8

浮動小数点数 • １よりも大きい最小の倍精度浮動小数点数は? • 使えるのは連続する５３ビット 20 2-1 2-1 1 0 0 ・ ・・・ 2-51 2-52 0 1 2-53 これ以降は情報を保持できない 9

1 20 2-1 2-1 1 0 0 2-1 2-1 0 0 2-1 2-1 0 0 2-1 2-1 0 0 20 1+2-52 1 20 -51 1+2 1 20 -51 1+2 1 -52 +2 ・ ・ ・ ・ ・・・ ・・・ ・・・ ・・・ 2-51 2-52 0 0 2-51 2-52 0 1 2-51 2-52 1 0 2-51 2-52 1 1 1 1+2u 1+4u 1+6u u=2-53 10

浮動小数点数 １ １+2u １+4u １+6u １+8u １+10u １の後はしばらく2uの間隔で浮動小数点数は存在する 単精度ならu=2-24 倍精度ならu=2-53 11

例題 • 1+u+u+u+u+u+u+…+u-1-u • が数値計算ではマイナスになる?（デモの内容） • 1+uの浮動小数点演算の結果は丸められて1 • 1-1は0 • よって数値計算結果はマイナス • 計算順は結果に大きく影響 12

浮動小数点数演算 • double a, b, c; • //a,bに何かしらの値をセット • c = a + b; • 浮動小数点数と浮動小数点数の演算結果は浮動小数点数か? – そうでない例がある→丸め誤差が発生 – a = 1, b=2-53, a+bは浮動小数点数ではない • 丸め誤差が発生する場合，最も近い浮動小数点数を返す – デフォルトの丸めモード – 偶数丸め方式 13

非正規化数の計算 • 倍精度浮動小数点数においては2-1022よりも小さな値を計算する ときには計算時間に気をつけることが必要な場合あり • MATLABの例: • n=1000; • A=2^-530*ones(n); • B=2^-530*ones(n); • tic; C=A*B; toc %ここでは計算時間を測る デモ 14

例外 • NaN (Not a number)が発生する – inf – inf – -inf + inf など • NaNとの比較は，すべて「偽」になる 15

浮動小数点数の基礎はここまで • デモによる?な例はなぜ起こるのか，わかったことと思います． • 浮動小数点数とその演算は – 扱えない実数があること – 実数のときにいえる性質は使えないものがあること • 例:結合則や分配則 を押さえておいてください． 16

余談:整数演算 • 浮動小数点数の基礎のお話をいたしました． • 丸め誤差には注意する必要性があります． • 整数演算はどうでしょうか? デモ • ラップアラウンドに注意 • （フォーマットは異なるが）ドラゴンクエスト４ カジノのお話が有名です 17

数値再現性 • 先ほどのデモで，数値計算（浮動小数点演算）は計算順によって結果が変わることは 理解できたと思います． ２スレッドで総和を計算 + + + + + + + + + + + + + + + スレッド２が計算 スレッド１が計算 各スレッド内では 逐次計算を仮定 ４スレッドで総和を計算 + + + スレッド１が計算 + + + + スレッド２が計算 + + + + スレッド３が計算 + + + + スレッド４が計算 18

数値再現性 • AVX （Advance Vector Extension）を用いる（倍精度の例） a b c d 256ビット e f g h 256ビット • a+e, b+f, c+g, d+hが同時に計算できる． • レジスタによって256ビットか512ビットの違いがある． 19

数値再現性 • 数値再現性に関連する要素 – スレッド数 – レジスタのサイズ（256ビットや512ビット） – Fused Multiply-Add • a*b+c – コンパイラの違い – コンパイルオプション（最適化） – 動的スケジューリング • ユーザの環境の違いが大きく数値計算結果に影響する 20

数値再現性 • ある手法の有効性が論文で示されていた． • 自分もその手法を実装して実験してみた!→有効性が確認できなかった? • その論文で書かれていたことが怪しいのか? • （著者にクレームをする前に）数値再現性の問題を思い出してください． • 数値再現性（Reproducible Numerical Algorithms）自体が研究のキーワードの １つ 21

連立一次方程式 • 連立一次方程式の近似解の精度を検証したい • 標準的な方法（密行列）:LU分解＋前進代入＋後退代入 – MATLAB:A¥b – LAPACK: dgesv • 疎行列に対しては反復解法を用いることも可 22

連立一次方程式 • 連立一次方程式のソルバの性能（精度）を検証したい． Ax=b, Aは係数行列，bは右辺ベクトル • 行列Aを何かしらの方法で作る • 解xを自分で事前設定する • 右辺ベクトルbをb:=A*xで生成する • xが解であるので，数値解（近似解）の精度を検証できる!? 丸め誤差の問題を忘れていませんか? デモ 23

連立一次方程式 A(x+x’)=b + b’ • 右辺ベクトルがb’だけずれたら，解はx’だけずれる Ax’=b’ x’=A-1b’ Aの正則性を仮定 || x’ || = ||A-1b’|| xの⾧さは0ではないと仮定 || x’ || ≦ ||A-1|| ||b’|| || x’ || / || x || ≦ ||A-1|| ||b’|| / ||x|| || x’ || / || x || ≦ ||A-1|| ||A|| ||b’|| / ||b|| 条件数 ||Ax|| = ||b|| ||A|| ||x||≧ ||b|| ||x||≧ ||b|| / ||A|| 1 / ||x||≦ ||A|| / ||b|| 24

連立一次方程式 • 誤差のチェックは残差で十分か?（近似解をx’とする） x’-x = x’ – A-1b = A-1 A x’ – A-1b = A-1 (A x’ – b) 両辺にノルムをとると ||x’-x|| = ||A-1 (A x’ – b)|| ||x’-x|| ≦ ||A-1 || || A x’ – b|| 誤差のノルム 残差のノルム 条件数が関連する式もあります 25

連立一次方程式 • 残差が良い（小さい）ほど近似解の精度がよいか? 26

組み込み関数の精度 • sin, cos, tan, log, expなどなど，組み込み関数を使 う方も多いのではないかと思います． • これらの関数，四則演算と平方根とは異なり， 精度に関する標準的な規定がない! • 計算環境によって異なる． 27

組み込み関数の精度 • Unit in the Last Place（ulp）という単位で精度は良く議論されます． だいたいの イメージ 1 ulp 誤差が0.5 ulp以内であれば，最近点への丸め（最高の近似）といえる 誤差が1.0 ulp以内であれば，faithful rounding といえる なお，ulpは定義によって異なる（Goldberg, Muller, Harrisonなどを参照） 28

組み込み関数の精度 • 私的な調査結果（単精度浮動小数点数，Harrison型） 環境 概ね 最悪の場合 注意 MATLAB 3.0以下 1000 ulp以上 asec, acoth, acsc Octave 3.0以下 すごい悪い sin, cos, csc Visual Studio (cl) 0.5より少し上 MATLABにおける例:a = 1 + 2^-13; double(asec(single(a))) asec(a) asec(mp(a)) の結果 0.015625158324838 0.01562420533397061 0.01562420533397061625680522414072565 ー optionによる? デモ 29

MATLABとOCTAVE • sum関数の違い • linspaceの違い • 初等関数の違い • 初等関数の精度について • MATLABにおけるGPU計算 デモ 30

同じコードを実行する • 先ほど紹介した事例は，関数がどう実装されているか に依存するので，特に驚きを感じなかった方もおられ ると思います． • では，同じコードをコンパイルした場合についてデモ を見ます． デモ 31

足し算の事例 • C言語で以下の計算をすると? res=0; • for (int i = 0; i < n; i++) { res += vector[i]; • • } • 配列vectorには順に1, u, u, u, u, u, u, u, u, …が入っていることにします • 1+u+u+u+u+u+u+u+u+… • を前から順に計算すればresには1が入りますね? デモ 32

Fused Multiply-Addの確認例 1+u • a=1+2u; • b=1-u; • c=u; １ 1+2u • とする． • a*b = (1+2u)(1-u)=1+u-2u2→1, 1+u →1 • a*b+c= 1+2u-2u2 →1+2u デモ 33

まとめ • 数値計算の誤差に関する言葉（ W. Kahan ） – 数値計算における誤差はいつでも心配する問題ではないが，無視 できるほど少なくはない問題である． • 自分が計算原理は知っていて損はない • 計算結果は – アルゴリズムの問題?計算誤差? – コンパイラ?コンパイルオプション? に依存する． • 数値計算結果の妥当性を知るには? → 次回のお話し 34

補足 • ここまでの講義で数値計算の誤差の問題ばかりを取り上げたの で，数値計算にネガティブな印象を持った方もおられるかもし れませんが，あくまで注意点の紹介です． • 紹介したツールに関しても，あくまで例として紹介しています． • IEEE 754の設計はしっかりとしたものだと思っています． • 数値計算は計算速度の観点で非常に重要なものです． （１秒間に～ギガ回の計算ができる） • 計算の原理がIEEE 754に規定されていない話も増えています． 35