第15回 配信講義　計算科学技術特論B（2024）

精度保証付き数値計算入門 芝浦工業大学システム理工学部 数理科学科 尾崎 克久

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

講義内容（15回目:2024/7/25 ） • アルゴリズムと丸め誤差 –２次元凸包の問題 • 精度保証付き数値計算とは • 浮動小数点演算の丸めのモード –簡単な計算機援用証明 • 区間演算の紹介 • 事前誤差解析（時間の限り） • まとめと参考文献 3

凸包について • 凸包: すべての頂点を含む最小の凸多角形 4

凸包について • 凸包: すべての頂点を含む最小の凸多角形 5

凸包について • 凸包を求めるアルゴリズム – 逐次添加法（本日扱います） – 包装法 – Quick Hull – グラハムの方法 など 6

点と有向直線の位置関係 • 計算機は計算結果によって位置関係を判定する (bx,by) (cx,cy) (ax,ay) 7

点と有向直線の位置関係 • 計算機は計算結果によって位置関係を判定する (cx,cy) (ax,ay) ＋ ー (bx,by) ax ay 1 bx by 1 cx cy 1 8

丸め誤差の問題はここでも • 数値計算では結果の符号は正しいですか? • もしも結果（符号）が異なっていたら? if 符号が正 正の場合の処理 elseif 符号が負 負の場合の処理 else ゼロの場合の処理 処理内容を誤る可能性 9

逐次添加法 • Pk: k番目の点 • Ck: k番目までの点に対する凸包 • P1, P2, P3からC3（同一直線状にないと仮定して三角形）を作る • Ck からCk+1 を作る • PkがCkの内部にあるかどうかが大事 • Ck の点は反時計周りに並ぶと仮定 10

逐次添加法 ５点の例題 P3 P4 P5 P1 P2 11

逐次添加法 C3をP1, P2, P3で作る P3 P4 P5 P1 P2 P4が内部にあるかどかを判定する 12

逐次添加法 P4は内部にあるとわかるのでP4を除外 P3 ＋ P4 P1 ＋ ＋ P5 P2 C4はC3と同じ. 13

逐次添加法 P5は C4 の外側にある P3 ＋ P1 ー ＋ P5 P2 14

逐次添加法 P3 ＋ P5 P1 ＋ P2 15

逐次添加法 P5をくっつける P3 ＋ P5 P1 ＋ C5 P2 16

丸め誤差の問題? • もしも点と有向直線の位置関係の判定を間違えたら? （行列式の符号を間違えたら?） – 無視できるくらいの問題 – 問題を起こさない場合 – 大きな問題 – とても大きな問題 17

無視できるくらいの問題 18

無視できるくらいの問題 － ＋ ＋ 19

無視できるくらいの問題 20

問題を起こさない場合 P3 P3 P4 P4 P1 P2 P5 P1 P2 P5 21

問題を起こさない場合 P3 P3 P4 ＋ P4 P1 P2 P5 P1 ＋ － P5 P2 22

問題を起こさない場合 P3 ＋ P4 P1 ー P3 ＋ ー P5 P4 P2 P1 P2 P5 23

大きな問題 P3 P2 P1 P5 P4 24

大きな問題 P3 P2 P1 P5 P4 25

大きな問題 P3 ＋ ＋ P1 P2 ＋ P5 P4 26

大きな問題 P3 ＋ ＋ P1 P2 ー P5 27

大きな問題 P3 P2 P1 P5 28

大きな問題 P3 P2 P1 P5 P4 29

とても大きな問題 P5 P8 P9 P4 P7 P2 P3 P1 P6 30

とても大きな問題 P5 + P8 P9 P7 P2 P1 + P3 P6 - P4 31

とても大きな問題 P5 P8 P9 P4 P7 P2 P3 P1 P6 32

とても大きな問題 P5 P8 P9 P7 P2 P1 + P6 P3 + P4 33

とても大きな問題 P5 P8 P9 P4 P7 P2 P3 P1 P6 34

誤差の問題 • １つ１つの計算の多くは問題ない・・・が • １回の間違いでアルゴリズムの進行がグダグダになることもある – 符号を間違えればif文の選択を誤る • （他の計算幾何の問題で）最悪無限ループに陥る事例も報告されている • アルゴリズムの破綻も大きな問題である 35

メッシュ生成の例（その１） • MATLABの例 • 斜辺について節点を • x=linspace(0, 1, 11); • y=linspace(1, 0, 11); • で生成すると・・・・

メッシュ生成の例（その２） SciPy, delaunay関数の例 Sora Sawai, Kazuaki Tanaka, Katsuhisa Ozaki, Shin'ichi Oishi Polygonal Sequence-driven Triangulation Validator: An Incremental Approach to 2D Triangulation Verification, https://doi.org/10.48550/arXiv.2401.08242

アルゴリズムと数値誤差 正しい 正しい ＋ ＋ 厳密計算 アルゴリズム コーディング 正しい 正しい ＋ ＋ 数値計算 アルゴリズム コーディング 正しい結果

精度保証付き数値計算とは • 精度保証付き数値計算 – 数値計算を用いて，解の存在範囲を求めることができる 数値計算 精度保証付き 数値計算 数式処理 速 高速性 遅 低 信頼性? 高 多 対応する問題の多さ 少

精度保証付き数値計算とは • イメージ 真の解の存在範囲 • （例）真の解は近似解を中心に半径～以内に存在する – もちろん高速に実行したいし，小さい区間を得たい • 失敗すると – -InfからInfに存在 – そもそも，わからない

浮動小数点演算の丸めのモード • 数値計算では丸め誤差の影響で計算結果が正しくない可能性がある． • 例:２つのベクトルが直交しないことを確認したい – ベクトルの内積を計算して０ではない→２つのベクトルは直交しない 数値計算を使えば丸め誤差の影響は? • 例: デモ

浮動小数点演算における 丸めのモード • 丸めのモード（一部を紹介） – （多くの環境でデフォルト）最近点への丸め（偶数丸め） – 上向き丸め – 下向き丸め a+b 浮動小数点数 浮動小数点数

浮動小数点演算における 丸めのモード • 真値が包含された区間を求める→区間演算 • :括弧内の数式を下向き丸めで計算した結果 • :括弧内の数式を最近点丸めで計算した結果 • ∆ :括弧内の数式を上向き丸めで計算した結果 が成立する．

丸めモードを使ってみる • Ｃ言語でよくある例 fenv.hをインクルードして int fesetround(int mode) を使用する 使えるモードは FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO • CUDAでは演算毎に丸めの方向を指定することが可能 • MATLABはあとでデモします．

浮動小数点演算における 丸めのモード • に関する計算 ∆ はダメ ∆ はよい 単純に ・下向き丸めを用いれば下限が求まる ・上向き丸めを用いれば上限が求まる などと思わないこと ∆

ベクトルが直交しない とすると以下が成立する 下限が正または上限が負であればベクトルが直交しない のときはわからないと判断する デモ

行列積では? が成立する? • をどう評価するかによる． • 例:Strassenの方法 ＝

区間という考え方 • 実数のすべてを浮動小数点数で表現できない – • 実数を包含する区間を浮動小数点数で表現することは 可能 – 下端・上端型 – 中心・半径型

区間という考え方 浮動小数点数 浮動小数点数

区間演算という考え方 • • と に関して，この和を考える = が成立するが，数値計算では? ⊆ ノート

区間演算という考え方 • と に関して，この差を考える = が成立するが，数値計算では? ⊆

区間演算という考え方 • と に関して，この積を考える = が成立するが，数値計算では? ⊆ 入力に制約があれば，もっと簡単に計算ができることがある

区間演算という考え方 区間演算すると，真の値を包含する区間を求めることが可能 四則演算および平方根の計算のみを用いる場合， 真の値を包含する区間を求めることが可能． 指数関数・対数関数・三角関数などは注意 デモ

精度保証付き数値計算 ある値を計算する if 符号が正 正の場合の処理 elseif 符号が負 負の場合の処理 else ゼロの場合の処理 精度保証化 ある値を区間演算で計算 if 区間の下端が正 正の場合の処理 elseif 区間の上端が負 負の場合の処理 わからないとき else ?? • 単に警告として捉える • 高精度計算・有利数演算に移行

精度保証法の例 • 点と有向直線の位置関係の判定問題 (cx,cy) (ax,ay) ＋ ー (bx,by) ax ay 1 bx by 1 cx cy 1

精度保証法の例 • 計算値（数値計算で評価） • • ある値（数値計算で評価） • 倍精度:2-53 • if文の条件として 正の正規化数の最小数 倍精度:2-1022 が真ならば の符号は正しい

精度保証法の例 • 精度保証は計算コストと区間幅にトレードオフの関係がある（図の幅はコスト） 計算時間 手法A 手法B 手法C（完全精度計算） • 適応的（adaptive）な手法を構築して，なるべく計算時間がかからないようにする • 手法Aはフィルタリングとも呼ぶ（static, semi-static, dynamic）

事前誤差解析（足し算の例） • fl（）:括弧内の計算を浮動小数点演算で評価した結果 ≦ ≦ ≦ ≦ はunit roundoff 単精度浮動小数点数であれば =2-24 倍精度浮動小数点数であれば =2-53

事前誤差解析（足し算の例） ≦ ≦ 真のa+b fl(a+b) 誤差 解析デモ

本日のまとめ • アルゴリズムによる数値誤差の影響 • 精度保証付き数値計算の紹介 • 区間および区間演算 • 丸め誤差解析について

本日の参考文献 • 点と有向直線の位置関係に関する浮動小数点フィルタ • K. Ozaki, T. Ogita, F. Bunger, S. Oishi, S.M. Rump: • Simple floating-point filters for the two-dimensional orientation problem, • BIT Numerical Mathematics,Vol. 56(2): 729-749, 2016. • 連立一次方程式の精度保証 • S. Oishi and S.M. Rump, • Fast verification of solutions of matrix equations, • Numer. Math., 90(4):755–773, 2002. • 大規模計算と精度保証付き数値計算 • K. Ozaki, T. Terao, T. Ogita, T. Katagiri, • Verified numerical computations for large-scale linear systems, • Applications of Mathematics, 66(2): 269-285, 2021.

精度保証付き数値計算の文献 • 大石進一編: • 精度保証付き数値計算の基礎， • コロナ社，2018. • S.M. Rump: • Verification methods: Rigorous results using floating-point arithmetic, • Acta Numerica, 19:287–449, 2010. • 凸包の問題点 Classroom examples of robustness problems in geometric computations: https://www.sciencedirect.com/science/article/pii/S0925772107000697

• https://www.sciencedirect.com/science/article/pii/S0925772107000697

精度保証付き数値計算の ライブラリ INTLAB • ハンブルク工科大学のSiegfried M. Rump先生が開発されているライブラリ • MATLAB / Octave 上で動く • https://www.tuhh.de/ti3/rump/intlab/ kv - C++による精度保証付き数値計算ライブラリ • 早稲田大学柏木雅英先生が開発されているライブラリ • 精度保証付き数値計算を行うためにC++で作成したライブラリ群 • http://verifiedby.me/kv/ • 加えて柏木先生のページは情報の宝庫

• https://www.tuhh.de/ti3/rump/intlab/
• http://verifiedby.me/kv/