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

106 Views

July 08, 24

スライド概要

第14回 7月18日 数値計算における丸め誤差
概要: 数値計算には有限桁の精度でよって表現される数が用いられている。 よ
って、実数演算とは異なることをまず理解しておく必要がある。 本講義では,数
値計算の問題点である誤差の問題点について解説し、 数値計算に用いる数の性質を
正しく理解することを目的とする。

シェア

またはPlayer版

埋め込む »CMSなどでJSが使えない場合

関連スライド

各ページのテキスト
1.

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

2.

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

3.

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

4.

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

5.

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

6.

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

7.

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

8.

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

9.

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

10.

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

11.

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

12.

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

13.

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

14.

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

15.

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

16.

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

17.

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

18.

数値再現性 • 先ほどのデモで,数値計算(浮動小数点演算)は計算順によって結果が変わることは 理解できたと思います. 2スレッドで総和を計算 + + + + + + + + + + + + + + + スレッド2が計算 スレッド1が計算 各スレッド内では 逐次計算を仮定 4スレッドで総和を計算 + + + スレッド1が計算 + + + + スレッド2が計算 + + + + スレッド3が計算 + + + + スレッド4が計算 18

19.

数値再現性 • 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

20.

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

21.

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

22.

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

23.

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

24.

連立一次方程式 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

25.

連立一次方程式 • 誤差のチェックは残差で十分か?(近似解を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.

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

27.

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

28.

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

29.

組み込み関数の精度 • 私的な調査結果(単精度浮動小数点数,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

30.

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

31.

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

32.

足し算の事例 • 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

33.

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

34.

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

35.

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