Irugo's Learning Record

# 数値計算の常識

# Info

  • 書名: 数値計算の常識
  • 著者: 伊理正夫・藤野和建
  • 読んだ日付: 2021.12

# Motivation

授業や趣味を通して数値計算手法(数値積分、ODEの数値解法、勾配法、ニュートン法、補間など)の基本的な部分には触れてきた。よく知られている本で、常識として知って置かなければならないことがたくさん書かれているらしいので、ある程度手法を知った今、読むべきだと思った。

# Contents

# 1. 数の表現と誤差

  • 加減に対しては絶対誤差が和に、乗除に対しては相対誤差が和になって伝播する。
  • 浮動小数点数の表現
    • 仮数部と指数部
  • 具体的な精度の話はあまりに昔のものだった

# 2. 桁落ちに気をつけよう(その1)

  • 絶対値のごく近い2つの数を加減すると、絶対値が小さくなった分だけ相対誤差が大きくなり、有効数字が減る(桁落ち)
  • 有利化・倍角/半角の変換などでどうにかして桁落ちしづらい形に変形してから計算する。
  • 有名な例として、ほとんど1次方程式な系に小さい2次項を付けたような2次方程式に対して解の公式を使うと桁落ちでハマるというのがある(これは聞いたこと有り)
    • m重根の有効数字の桁数は計算桁数の1/m程度になるらしい

# 3. 桁落ちに気をつけよう(その2)

  • 収束判定で差の絶対値をεと比較するやつ
    • f(x)の計算方法とxの値に強く依存している
  • 多数の数を加える
    • 現在の和に比べて被加数が非常に小さい時、被加数の情報は抜け落ちてしまう。丸め誤差。
    • 数値積分やODEの数値計算などで、刻みをいたずらに小さくすればいいというものでもない！
    • ↑の理由でかえって結果が悪くなることもある。

# 4. たった1回だけの計算なんて...

  • だった1回の計算だけでは結果の正しさについてほとんど何もわからない
  • 系統的に計算方法を変えて複数回計算すると、非常に有用な情報を得られる
    • 刻み幅、計算桁数を変えて実験
    • 刻み幅の変更で打ち切り誤差の検討・計算桁数の変更で丸め誤差の検討
  • 結果に置いて正しい数字の桁数は計算桁数より少ない(桁落ち)が、落ちた桁数はほぼ一定。

# 5. 逆行列よさようなら

  • Gauss-Jordanの掃出し方
    • 乗算n^3回
  • LU分解
    • 乗算n^3/3回 (対称行列ならさらに1/2、疎行列ならより効率的！)
    • AからGauss Jordanで直接Ax = bにいく、あるいは逆行列をn^3回乗算で計算してからAx=bにいくよりも、AからLU分解へn^3/3回で言って、n^2回でAx=bの解に行ったほうが良い。
    • 正定値対称行列ならコレスキー分解
    • Aが疎行列のとき、A^{-1}は往々にして密行列。一方、LUは疎行列性を保つことが多いので、うれしい！
  • 数値計算的には、逆行列ではなくLU分解を使うべき。

# 6.単位と次元

  • 単位のあるパラメータは変形して無次元量に規格化してからSolveする
    • 特にODEなど。

# 7. 数値積分法-台形則を使いこなすには

  • 台形則は刻み幅を1/2にすると誤差が1/4になる。この傾向からはずれるということは、f(x)が扱いづらかったり、丸め誤差の影響が出ているということ。これを逆手に取れば、計算環境の丸め誤差やf(x)の性質について観察できる。
  • Richardson補外により、系統的にいくらでも精度のよい公式を作れる。
  • 不連続点があればそこで区間を分割する。逆にRichardson補外の式で刻み幅を倍にしたときの精度の挙動を調べれば、不連続点をまたいでいるかもしれないということになる。
  • 端点の特異性は変数変換で除く。
  • 周期関数を周期に渡って数値積分するときは、台形則は非常に精度が良い。区間両端での導関数の値が一致するため。
  • 無限積分に対しては、刻み幅hを指定して、f(nh) (n=0, 1, 2, ...)と計算していき、例えば連続するいくつかの和がこれまでの累積和のε倍になれば打ち切る、などとする。
  • 万能な変数変換として、二重指数型変換が知られている。[-1, 1]に正規化してx=tanh(pi/2 sinht)とすれば端点特異性を飛ばせる。
  • 多重積分は難しい。次元の増加が課題で、結局はモンテカルロ積分になりがち。

# 8. 数値微分法-打切り誤差と丸め誤差の闘い

  • hを小さくしていくとはじめの内は理論通り誤差が小さくなるが、ある時点からは桁落ちの影響で誤差が大きくなる。
  • 普通の片側差分の場合O(h^2)の誤差。これと丸め誤差の影響を考えると、演算桁数の半分くらいがよい。
  • 中心差分の場合はO(h^4)の打ち切り誤差だが、丸め誤差の限界は大差ない。
  • Richardson補外を繰り返してRomberg1段公式などの精度のよい公式をいくらでも作れる。
    • 片側差分なら、hを半分にしたら差が1/2になる傾向が正常
    • 両側差分なら、hを半分にしたら差が1/4になる傾向が正常
    • これらの傾向は、打ち切り誤差が丸め誤差より卓越的で、打ち切り誤差の漸近理論が適用できることの実際的な判定法となる。

# 9. Newton法

  • 非線形方程式の数値解法
  • 2次収束
  • 解x=aの近くで関数値を計算する際に予想される丸め誤差などの影響をだいたい見積もり、関数値がそれより小さくなったら終える、とするしかない。
  • m重解の精度は計算精度の1/mの桁の精度くらいしかない。

# 10. 複素数の計算

  • 複素数の形式として、直角座標表現と極形式とがある。後者への変換のコストは結構高い。
  • 絶対値の大きさの精密な値が必要？大きさがだいたいわかればOK？後者のほうが当然雑でよいが、前者なら気をつける必要あり。
  • 偏角は{n, +-, 90度以内のズレ}の情報だけ持っておいて、明示的に角度を計算せずに、極形式から直交座標に戻すときに式変形でなんとかするべき。

# 11. 代数方程式

  • 計算誤差の上界を勘弁に見積もれる。打ち切りεはそれでOK。
  • 問題は初期値の選び方。
    • 解の重心b=-a_1/nのまわりで多項式p(z)を展開し、その係数に対して1次の項だけ抜かして係数を負にした実係数方程式を考える。これにはただ一つの正の解がありこれを半径として解の重心を中心とする円盤のなかにある。ここを初期値としてNewton法を走らせると良い。
  • 平野法
    • ニュートン法の課題である大域的収束性を保証
  • n次方程式のn個の解すべてを同時に求めるという連立法
    • 2次法
    • 3次法
  • Smithの誤差評価法
    • 解の存在円を逐次的に求めていって、円がそれぞれ分離されていれば十分な計算桁数で分離できている。そうでなければその円が重なっている連結成分に重解あるいは今の計算精度では分離できないくらい近い解があることになり、まだ続けるか？という話になる。

# 12 常微分方程式の初期値問題

  • Euler法で満足しない
    • Euler法はO(h)
    • Euler法で求めた仮のt+1ステップの点との平均増加率を用いるという修正Euler法はO(h^2)
    • より高精度な公式として常識的なのは
      • Runge-Kutta法 O(h^4)
      • 予測子・修正子法 O(h^4)

# 13. 数値計算の手間の理論と実際

  • 多項式の計算
    • Horner法が最良
  • 単項式
    • 繰り返し二乗法(Doublingの応用)
    • O(log N)
  • 線形計算
    • 内積、行列ベクトル積は愚直が最良
    • 行列行列積はO(n^3)よりオーダーの良い方法がいろいろ考えられているが、定数倍がクソデカだったり、密行列にしか有用でなかったりする。
    • たいていは疎行列なので、疎行列の性質を活かす方法を考えるのが良い。
  • フーリエ変換
    • FFT O(N^2)をO(N log N)に削減できる。

# 14. 数値計算の中の非数値計算

  • 疎行列において非零要素だけを格納するような効率的なデータ構造の方法など。
  • CSの観点から工夫したり。
  • 数値計算は数学だけではない。

# 15. 補間(内挿)

  • ラグランジュ補間が代表的
    • 少し性質が悪い関数だと、次数を上げても振動するだけでうまく近似できない。
  • 分点の選び方を工夫

# 16. 雲形定規とスプライン-続補間

  • スプライン補間
    • 区分的に低次の多項式で近似しつつ、その区間の境界では2階の導関数までが連続につながるようにする。
    • 曲率の2乗積分が最小になるような曲線という意味合いを持つ。

# 17. 刻みは細ければ良いというものではない-偏微分方程式I(拡散方程式)

  • 陽的解法公式
    • r = (Δt / Δx^2)がほぼ1/2以下でないと陽的解法では求められない(安定性条件)
  • 陰的解法公式
    • 無条件に安定。例えばt->infとした定常解にだけ興味ある時とか

# 18. 偏微分方程式I(楕円形方程式)

  • ラプラス方程式、ポアソン方程式など。
  • SOR法、ガウスザイデル法
  • 大規模なPDEの数値解法は「しかるべき専門家に相談すべきである」

# 19. 数列の収束の速度

  • 1次収束
    • |a_n+1 - a| <= c|a_n - a|
    • nが大きいところでの差分の比λの大きさに収束速度が強く依存する
  • 2次収束
    • Newton法など
    • 加速が不要なくらい十分
  • p次収束
    • 一般の場合
  • O(n^(-a))の収束
    • |a_n - a| <= cn^(-a)
    • 早いとは言えない。
    • また、近似列の部分列ごとに規則的な収束性があり、全体としてはじわじわ収束していくみたいなパターンもあり。
    • 収束の加速法が必要

# 20. 数列の収束の加速

  • Romberg-Richardson法

    • a_n = a + b_1 * λ_1^n + b_2 * λ_2^n + ...の形で表され、係数λはすでにわかっている場合
    • a_n(1) = (a_n - λ_1 a_{n-1}) / (1 - λ_1)　とするとλ_1のべき乗のファクタを除くことができて、より早く収束させられる。

  • Aitken法

    • lambdaがわかっていない場合に有用。
    • https://ja.wikipedia.org/wiki/%E3%82%A8%E3%82%A4%E3%83%88%E3%82%B1%E3%83%B3%E3%81%AE%CE%942%E4%B9%97%E5%8A%A0%E9%80%9F%E6%B3%95#:~:text=%E3%82%A8%E3%82%A4%E3%83%88%E3%82%B1%E3%83%B3%E3%81%AE%CE%942%E4%B9%97%E5%8A%A0%E9%80%9F%E6%B3%95%EF%BC%88%E3%82%A8%E3%82%A4%E3%83%88%E3%82%B1%E3%83%B3%E3%81%AE%E3%83%87%E3%83%AB%E3%82%BF%E3%81%98,%E3%81%AE%E4%B8%80%E3%81%A4%E3%81%A7%E3%81%82%E3%82%8B%E3%80%82

  • ε-algorithm

# 21. 中間変数の活用

- 実装的な話。事前計算しておく、定数を使う、Cache活用、などなど。今ならCompilerの最適化でやってくれる部分も多いが、明示的にできるならすべき。

# 22. 固有値の計算

  • 入門者程度ではまともなプログラムを組むのは難しい。どのように既存手法をチョイスすればいいか？

  • 部品

    • 行列Aを扱いやすい形に相似変換
      • Householder変換が主流
    • 反復法で固有値固有ベクトルを求める
      • QR法が支流
    • 減次
    • 原点シフト
      • AとA-μIは同じ固有ベクトルをもち、固有値はμだけシフトすることを活用して収束を高速化
    • 逆反復
      • 固有値の近似値μと固有ベクトルの近似ベクトルvが得られたあとに(A-μI)x = vをとく。

  • 固有値問題の難易度

    • 実対称、エルミート行列
    • 固有値の絶対値が互いに離れている
    • 固有値が皆単純
    • 対角化可能行列

    なら一般的に解きやすい

# 23. モンテ・カルロ法

  • きちんとした方法で解けるならそちらのほうが良い。特に多重積分などはモンテカルロに頼りがち。
  • 数値積分
    • 一様乱数でサンプリングして、領域の面積・体積を求める。
    • サンプリング点数にのみ依存し、次元にはよらないため、次元が高いときに有用。
    • 分散の推定値もしっかり出しましょう。

# 感想

数値計算はCSの中でも歴史が深く、職人的な工夫も多いため、よほど専門として関わるのでない限りはすべてを網羅しようとするのは得策でないというのは非常に感じた。 「いざというときにしかるべき専門家に聞くことを常識とすべき」というのがよくわかった。 ユーザとして最低限知っておくべきことは理解できたと思う。基本的な方法はEigenなどを使って実装したい。

Last Updated: 12/30/2021, 8:49:18 PM

← Displayed Title →