数値解析（01/XX） - 「大人の教養・知識・気付き」を伸ばすブログ

はじめに

　プログラミング云々ばかり言っていても、それがどのように計算されているかを知らないと仕様が無い。ということで数値解析を学んでいく。
　基本的には

数値解析入門 (サイエンスライブラリ現代数学への入門)

作者:山本哲朗
サイエンス社

Amazon

を参照しつつ、他書で補足する。

はじめに
0.　数値アルゴリズム
1.　数値計算における誤差
今日の復習
次回
参考文献

0.　数値アルゴリズム

　数値解析は実数上で定義された計算表現のアルゴリズムを研究する学問である。平方根がその例であり、今日ではこれらの評価を計算機上ないし計算プログラムにおいて二乗計算と同じくらいに簡単かのように計算する。これで可能になったのが数値解析であり、この数値解析がどのようになされるかを学んでいく。
　数値解析について説明するには2つの異なる段階がある：

アルゴリズムの発展
アルゴリズムの解析

原理的にはこれらは独立したものである。しかし現実にはアルゴリズムの発展はアルゴリズム解析もしくは同じものや同類の問題をより簡単に計算するアルゴリズムから導出されることもある。
　ある程度対立関係にある、実数を用いるアルゴリズムには3つの特徴がある：

アルゴリズムの精確さ(一貫性)
　アルゴリズムは適当な制約下において要求された精度で望まれた値に近似することを意味する。
アルゴリズムの安定性
　アルゴリズムの挙動がそのパラメータの観点で連続であることを意味する。
有限精度の計算による影響(すなわち丸め誤差)
　有限精度の計算のための良い数学モデルは存在せず、有限精度の計算の挙動に対しては荒い上限を為すことを甘受せざるを得ない。

　安定したアルゴリズムの精度ないし効率性を改善しようとすると、不安定、その結果として最低限度の値に関するアルゴリズムを検討することになることもままある。アルゴリズムの効率性はより複雑な概念であるが、一方のアルゴリズムを別のものから選別する基本になり得る。
　数値解析の別のテーマは適応性( $adaptivity$ )である。計算アルゴリズムを効率性および／ないし安定性を改善する方法として問題のデータに適用することを意味する。

1.　数値計算における誤差

　解の水準に基づいてアルゴリズムの精度は求められる。(物理)現象を数理モデリングする際、数値的に計算することが必要になり得るが、理論的なモデリングはその必要な精度を得るのには適当だとは限らない。必要に応じて変形したり別の手法を用いたりと、数値計算のために独自の手法やアイディアを考える必要がある。

1.1.　数値表現

　計算機内部では数値は2進数または16進数で表現されている。またひとつの数値を表現するのに一定の桁数（ビット数）が与えられているが、そのために実数をコンピュータを扱おうと思っても、コンピュータは有限桁の有限個の数値しか扱えず、表現可能な数値に制約がある。

$\begin{aligned} \pm 0.d_1 d_2\cdots d_N\times 10^{p} \end{aligned}$

　実数型の数値には32ビットを使用する単精度実数と64ビットを使用する倍精度実数の2種類がある。両者とも計算機内部では以下のような浮動小数点で記憶される。

$\begin{aligned} \pm 0.d_1 d_2\cdots d_t\cdot \beta^{m}=\pm \displaystyle{\sum_{i=1}^{t}d_i\beta^{m-i} },\\ \end{aligned}$

と表示される。ここで

$\beta\in\mathbb{Z}$	：	基数
$t+1$	：	桁数
$M$	：	最小指数
$N$	：	最大指数
$\pm0.d_1 d_2\cdots d_t$ , $d_i\in\mathbb{Z},d_1\neq0$ , $0\leq d_i\lt \beta$	：	仮数部
$m, -M\leq m\leq N$	：	指数部

と呼ぶ。この表記方法を $β$ 進 $t$ 桁浮動小数点数という。
　このため演算の結果は、たとえ無理数であっても、何らかの形（四捨五入、切り上げ、切り捨てなど）でこの形態に収められる。

1.1.1　近似方法

　実数をコンピュータで表示するのに近似する方法はいくつか存在する。

　1つが、切り捨てである。実数

$\begin{aligned} \tilde{x}&=\sigma\cdot\left(\displaystyle{\frac{d_{0}}{\beta^{0}}}+\cdots+\displaystyle{\frac{d_{n}}{\beta^{n}}}+\displaystyle{\frac{d_{n+1}}{\beta^{n+1}}}+\cdots\right)\cdot\beta^{m}\in\mathbb{R},\\ \sigma&=1\ \ \ \mathrm{or}\ \ \ -1 \end{aligned}$

において $-M\leq m\leq N$ が成立しているとする。このとき $\tilde{x}$ の近似値として

$\begin{aligned} x=\sigma\cdot\left(\displaystyle{\frac{d_{0}}{\beta^{0}}}+\cdots+\displaystyle{\frac{d_{n}}{\beta^{n}}}\right)\cdot\beta^{m} \end{aligned}$

を採用するものである。
　また

$\begin{aligned} x&=\sigma\cdot\left(\displaystyle{\frac{d_{0}}{\beta^{0}}}+\cdots+\displaystyle{\frac{d_{n}}{\hat{\beta}^{n}}}\right)\cdot\beta^{m},\\ \hat{d}_n&=\begin{cases}d_n,\ d_{n+1}\lt \displaystyle{\frac{\beta}{2}}\\d_n +1,\ d_{n+1}\geq \displaystyle{\frac{\beta}{2}}\end{cases} \end{aligned}$

を採用する方式を四捨五入（ $(\beta/2-1)$ 捨 $(\beta/2)$ 入）という。

1.2.　絶対誤差と相対誤差

　 $x$ を具体的に求めるべく計算しようとしても、何らかの原因から厳密解を求められるとは限らず、近似値を求めることも少なくない。そうなれば厳密解との差異を議論することでその精度を検討することが必要となる。
　そのために“誤差”という概念を導入する。 $a$ を $x$ の近似値とするとき、

$\begin{aligned} \varepsilon&=x-a,\\ \left|\varepsilon\right|&=\left|x-a\right|,\\ e&\geq \left|\varepsilon\right| \end{aligned}$

をそれぞれ $a$ の誤差、 $a$ の絶対誤差、そして誤差 $\varepsilon$ の限界という。
　単位（水準）が相違するときには相対誤差

$\begin{aligned} e_R=\displaystyle{\frac{|e|}{x}}=\displaystyle{\frac{|x-a|}{x}} \end{aligned}$

を用いる。
　また限界は

$\begin{aligned} \left|a\right|-\varepsilon\leq \left|a\right|-e\leq \left|a+e\right|\leq \left|a\right|+e\leq \left|a\right|+\varepsilon \end{aligned}$

に注意すると、 $\displaystyle{\frac{\varepsilon}{|a|}}$ が充分に小さい場合には

$\begin{aligned} e_R=\left|\displaystyle{\frac{|e|}{x}}\right|=\displaystyle{\frac{|e|}{|a+e|}}\leq\displaystyle{\frac{\varepsilon}{|a|-\varepsilon}}=\displaystyle{\frac{\varepsilon}{|a|}}+\left(\displaystyle{\frac{\varepsilon}{|a|}}\right)^2+⋯\approx \displaystyle{\frac{\varepsilon}{|a|}} \end{aligned}$

と近似できる。

1.3.　丸めの誤差

　現在の規格IEEE754-2008におけるbinary64（倍精度実数型）では $\beta=2,\$ $t=52,\$ $M=1022,\$ $N=1023$ としている。この規格で扱える最大および最小の数 $x_{\max},x_{\min}$ は

$\begin{aligned} x_{\max}&=\left\{\left(\displaystyle{\frac{1}{2}}\right)^0+\left(\displaystyle{\frac{1}{2}}\right)^1+\left(\displaystyle{\frac{1}{2}}\right)^2+\cdots+\left(\displaystyle{\frac{1}{2}}\right)^{52}\right\}\cdot2^{1023}\\ &=(2-2^{-52})\cdot2^{1023}\\ &\approx 1.7969\times 10^{308},\\ x_{\min}&=\displaystyle{\frac{1}{2}}\cdot2^{-1022}\approx2.2250\times 10^{-308} \end{aligned}$

である。
　「1よりも大きい最小の数値と1との差」で計算機イプシロンを定義し、 $\varepsilon_M$ と書くことにすると、

$\begin{aligned} 1=\left\{\left(\displaystyle{\frac{1}{2}}\right)^0+\left(\displaystyle{\frac{0}{2}}\right)^1+\left(\displaystyle{\frac{0}{2}}\right)^2 +\cdots+\left(\displaystyle{\frac{0}{2}}\right)^52\right\}\cdot2^0 \end{aligned}$

であり、これよりも大きく最小の数値は

$\begin{aligned} \left\{\left(\displaystyle{\frac{1}{2}}\right)^0+\left(\displaystyle{\frac{0}{2}}\right)+\left(\displaystyle{\frac{0}{2}}\right)^2+\cdots+\left(\displaystyle{\frac{1}{2}}\right)^{52}\right\}⋅2^{0} \end{aligned}$

であるから、

$\begin{aligned} \varepsilon_{M}=2^{-52}\approx 2.2204\times 10^{-16} \end{aligned}$

である。
　この規格では制約として

$\begin{aligned} 0\leq d_i\lt \beta,\ d_1\neq0,\ -M\leq m\leq N \end{aligned}$

があり、しかもその範囲内であってもすべての値を表すことができない。
　たとえば $\beta=2,\ t=2,\ M=2,\ N=1$ とすると

仮数部	$\displaystyle{\frac{1}{1}}+\displaystyle{\frac{0}{2}}+\displaystyle{\frac{0}{4}}=\displaystyle{\frac{4}{4}}$	指数部	$2^{-2}$
	$\displaystyle{\frac{1}{1}}+\displaystyle{\frac{0}{2}}+\displaystyle{\frac{1}{4}}=\displaystyle{\frac{5}{4}}$		$2^{-1}$
	$\displaystyle{\frac{1}{1}}+\displaystyle{\frac{1}{2}}+\displaystyle{\frac{0}{4}}=\displaystyle{\frac{6}{4}}$		$2^0$
	$\displaystyle{\frac{1}{1}}+\displaystyle{\frac{1}{2}}+\displaystyle{\frac{1}{4}}=\displaystyle{\frac{7}{4}}$		$2^1$

であるから、このときには

$\begin{aligned} \displaystyle{\frac{3}{2}}=\displaystyle{\frac{6}{4}}\cdot 1^{0},\displaystyle{\frac{7}{8}}=\displaystyle{\frac{7}{4}}\cdot 2^{-1} \end{aligned}$

とこの体系で表現できるものの、

$\begin{aligned} \displaystyle{\frac{3}{2}}+\displaystyle{\frac{7}{8}}=\displaystyle{\frac{19}{8}},\ \displaystyle{\frac{3}{2}}\times\displaystyle{\frac{7}{8}}=\displaystyle{\frac{21}{16}} \end{aligned}$

はいずれもこの体系では表現できない。したがって、計算のために入力した値や計算途中で算出される中間値、そして最終的な答えはすべてこの体系で表現できる値に近似される。このように「計算において入力値、中間値および最終的な解の値がコンピュータで利用できる規格に変換されることで生じる厳密値との乖離」を実数 $\tilde{x}$ を浮動小数点 $x$ に丸めるといい、誤差

$\begin{aligned} \tilde{x}-x \end{aligned}$

を丸め誤差という。1.1.1節で述べた切り捨てや四捨五入による誤差がまさに丸め誤差である。

1.4.　桁落ちと情報落ち

1.4.1　浮動小数点数の演算

　浮動小数点数での四則演算は以下のようになされる。
　 $x=\xi\cdot\beta^m,y=\eta\cdot\beta^n$ が与えられているとする。このとき

加減算
	(1)	指数を一致させるように仮数部を調整する。
	(2)	仮数部の加減算を実行する。
	(3)	$t+1$ 桁の浮動小数点数になるように丸める。
乗除算
	(1)	仮数の乗除算を行う。
	(2)	乗算のときは $m+n$ を、除算のときは $m-n$ を新たな指数とする。
	(3)	$t+1$ 桁の浮動小数点数になるように丸める。

1.4.2　桁落ち

　ある実数 $x,y$ の近似値として $a=2.62347,$ $b=2.62315$ が得られており、それらが小数第5位まで正しいとする。このとき

$\begin{aligned} a-b=2.62347-2.62315=0.00032=0.32\times 10^{-3} \end{aligned}$

であるため、

$\begin{aligned} \left|e_R\right|&\approx\left|\displaystyle{\frac{x-y-(a-b)}{a-b}}\right|\leq\displaystyle{\frac{|x-a|+|y-b|}{a-b}}\\ &\leq\displaystyle{\frac{10^{-5}}{0.32×10^{-3}}}\\ &\approx 0.31\times10^{-1} \end{aligned}$

が得られる。この誤差限界は $a,b$ の相対誤差限界である $0.5×10^{-5}$ に比べてかなり大きい。
　なぜこのように大きな誤差が生じたのか。これは途中で丸めが発生し有効数字が減ったためである。このように「途中で丸めが発生し有効数字が減った中で値がほぼ等しい数値差を求めた時に発生する誤差」を桁落ちといい、その発生に注意しなければ相当の誤差が生じ得る。

1.4.3　情報落ち

　たとえば2次方程式 $ax^2+2bx+c=0,\ a\neq0$ の解を、解の公式

$\begin{aligned} x_1=\displaystyle{\frac{-b+\sqrt{b^2-ac}}{a}},\ x_2=\displaystyle{\frac{-b-\sqrt{b^2-ac}}{a}} \end{aligned}$

で求めることを考える。
　 $b\gt0,b^2\gg ac$ ならば、 $x_1$ の分子の計算において桁落ちが生じる。そのため

$\begin{aligned} x_1=\displaystyle{\frac{(-b+\sqrt{b^2-ac})(-b-\sqrt{b^2-ac})}{a(-b-\sqrt{b^2-ac})}}=\displaystyle{\frac{-c}{b+\sqrt{b^2-ac}}} \end{aligned}$

と変形すべきである。他方で $b\lt0$ かつ $b^2\gg ac$ ならば

$\begin{aligned} x_2=\displaystyle{\frac{c}{b+\sqrt{b^2-ac}}} \end{aligned}$

と変形すべきである。
　実際に $a=0.0001,b=300,c=0.00005$ とおいたとき、

(1)	解の公式	$-8.35598257253878\times10^{-8}$
(2)	変形式	$-8.33333333333345\times10^{-8}$
(3)	絶対誤差	$2.26492392053279\times10^{-10}$
(4)	相対誤差	$2.71790870463931\times10^{-3}$

と相対誤差ベースでみると誤差の大きさを無視し得ないことが推察される。
　また指数部に大きな差のある2つの数による和や差でも同様に大きな誤差が生じる。これを情報落ちという。

1.5.　打ち切り誤差

　関数 $f(x)=e^x$ の無限級数展開

$\begin{aligned} f(x)=1+\displaystyle{\frac{x}{1!}}+\displaystyle{\frac{x^2}{2!}}+\displaystyle{\frac{x^3}{3!}}+\cdots,x\gt0 \end{aligned}$

を有限項で打ち切って

$\begin{aligned} f_n(x)=1+\displaystyle{\frac{x}{1!}}+\displaystyle{\frac{x^2}{2!}}+\displaystyle{\frac{x^3}{3!}}+\cdots+\displaystyle{\frac{x^n}{n!}} \end{aligned}$

とおくとき誤差

$\begin{aligned} e_T (x)=f(x)-f_{n}(x)=\displaystyle{\frac{x^{n+1} e^{\xi}}{n+1}},\ 0\lt\xi\lt x \end{aligned}$

が生じる。このように無限項の和の操作を有限項で打ち切ったことで生じた誤差を打ち切り誤差という。
　コンピュータは有限界の四則演算しかできないため、ある関数 $f(x)$ の値を求める際、有限項からなるそれの近似値 $f_n(x)$ を用いて近似する。このとき $f_n (x)$ 自体でも代入される $x$ でも丸めの誤差が発生するため、コンピュータには $\tilde{x},\ \tilde{f_n}(\tilde{x})$ が表示される。その誤差は

$\begin{aligned} f(x)-\tilde{f}_n(\tilde{x})=(f(x)-f(\tilde{x}))+(f(\tilde{x})-f_n (\tilde{x}))+(f_n(\tilde{x})-\tilde{f_n}(\tilde{x})) \end{aligned}$

となる。このとき右辺第1項を代入誤差（伝搬誤差）、右辺第2項を打ち切り誤差、右辺第3項を生成された誤差という。

1.6.　オーバーフロー・アンダーフローと非正規化数

　計算の途中で数値が表現可能な上限値を値が超えてしまい、その値を表現できなくなることをオーバーフローという。IEEE754-2008では

$\begin{aligned} x_{\max}\approx 1.7969\times10^{308} \end{aligned}$

であるから、これよりも大きくなればオーバーフローを起こす。
　逆に絶対値が小さすぎるために $0$ でない実数が $0$ に丸められてしまうことをアンダーフローという。
　IEEEでは、逆に $x_{\min}\approx 2.2250\times10^{-308}$ よりも小さい正数を用いることもできる。実際、

$\begin{aligned} y=\pm\left\{\left(\displaystyle{\frac{d_0}{2}}\right)^0 +\cdots+\left(\displaystyle{\frac{d_{51}}{2}}\right)^{51}\right\}\cdot 2^{-1023},\ d_i=0,1,\ 0\leq i\leq51 \end{aligned}$

が利用でき、これを非正規化数という。非正規化数のうち絶対値が最も小さい値 $y_{\min}$ は

$\begin{aligned} y_{\min}=2^{-51-1023}\approx 4.9406\times10^{-324} \end{aligned}$

であり、絶対値が最大の正数は

$\begin{aligned} y_{\min}=(2-2^{-51})\cdot2^{-1023}\lt x_{\min} \end{aligned}$

である。

1.7　誤差評価

　アルゴリズムを評価するには、2つの観点からの評価が必要となる。

求めたい解の近くにある合理的な初期値を選んだと仮定したときのアルゴリズムの局所的な挙動を知りたい
アルゴリズムの大域的な挙動(任意の初期値を与えたときの挙動)を知りたい

1.8　初期値をどこにすべきか

　呼び出しアルゴリズムでは呼び出しの初期値をどこかに設定しなければならず、その選択が興味のある問題であることがある。これがアルゴリズム全体のパフォーマンスに影響し得る。

1.9.　計算量と誤差

　あらゆる計算は有限の資源を用いて有限の時間内に終了しなければならない。そこで計算さの複雑さ（計算量）を考慮する必要がある。計算の開始から終了までに要する数値演算の総量で計算量を見積もりたいものの、計算機の計算作業は複雑であるため、負荷の大部分を占める特定の演算の実行回数を数えることで計算量を大雑把に把握しておく。
　なお計算量には別の制約としてメモリの量的な問題がある。

1.9.1　Landauの記号

　 $x\rightarrow a$ のとき

$\begin{aligned} \left|\displaystyle{\frac{f(x)}{g(x}}\right|\rightarrow k\in\mathbb{R} \end{aligned}$

と有界ならば、 $f(x)$ は高々 $g(x)$ の位数（order）にあるといい、

$\begin{aligned} f(x)=O(g(x)) \end{aligned}$

と書く。

　例えば、 $f(x)=2x+3x^2+4x^3,\ g(x)=x$ に対して

$\begin{aligned} \left|\displaystyle{\frac{f(x)}{g(x)}}\right|=\left|\displaystyle{\frac{2x+3x^2+4x^3}{x}}\right|=|2+3x+4x^2|\rightarrow2(x\rightarrow0) \end{aligned}$

であるから、

$\begin{aligned} x\rightarrow 0\Rightarrow f(x)=O(g(x)) \end{aligned}$

これは $x\rightarrow 0$ で $f(x)$ が $g(x)$ にほぼ比例することを意味する。
　特に $x\rightarrow a$ のとき

$\begin{aligned} \left|\displaystyle{\frac{f(x)}{g(x)}}\right|→0 \end{aligned}$

ならば(すなわち $k=0$ の場合)、

$\begin{aligned} f(x)=o(g(x)) \end{aligned}$

と書く。

1.10　実装（計算機イプシロン）

1.10.1　Python

########################
### 計算機イプシロン ###
########################

import sys
info_eps = sys.float_info.epsilon

print(info_eps)

# 定義通りの計算方法
EPS = 0.0
tmp = 1.0

while 1.0 + tmp > 1.0:
    EPS = tmp
    tmp /= 2
print(EPS)

# もう1つの計算方法
a = 4 / 3
b = a - 1
c = b + b + b
EPS = 1 - c
print(EPS)

1.10.2　Julia

########################
### 計算機イプシロン ###
########################

# 計算機イプシロン
ϵ_m = eps()

println(ϵ_m)

# 定義による計算
EPS = 0.0
tmp = 1.0

while 1.0 + tmp > 1.0
    EPS = tmp
    tmp /= 2
end

println(EPS)

今日の復習

実数 $x$ をあるコンピュータに入力すると浮動小数点で記憶されるとき、その実数 $x$ を $β$ 進 $t$ 桁浮動小数点数で表示せよ。
$a$ を実数 $x$ の近似値とするとき、誤差、絶対誤差、限界および相対誤差を定義せよ。
丸め誤差を説明せよ。
桁落ちを説明せよ。
情報落ちを説明せよ。
打切り誤差を説明せよ。
オーバーフローを説明せよ。

次回

power-of-awareness.com

参考文献

伊理正夫・藤野和建（1985）「数値計算の常識」（共立出版）
菊地文雄・齊藤宣一（2016）「数値解析の原理　現象の解明をめざして」（岩波書店）
齊藤宣一（2012）「数値解析入門」（東京大学出版）
齊藤宣一（2017）「数値解析」（共立出版）
高橋大輔（1996）「理工系の基礎数学8　数値計算」（岩波書店）
山本哲朗（2003）「数値解析入門[増訂版]」（サイエンス社）
日本応用数理学会監修・櫻井鉄也編・松尾宇泰編・片桐孝洋編（2018）「数値線形代数の数理とHPC」（共立出版）
Ridgway Scott, Larkin (2011), "Numerical Analysis", (Princeton University Press)

はじめに

0. 数値アルゴリズム

1. 数値計算における誤差

1.1. 数値表現

1.1.1 近似方法

1.2. 絶対誤差と相対誤差

1.3. 丸めの誤差

1.4. 桁落ちと情報落ち

1.4.1 浮動小数点数の演算

1.4.2 桁落ち

1.4.3 情報落ち

1.5. 打ち切り誤差

1.6. オーバーフロー・アンダーフローと非正規化数

1.7 誤差評価

1.8 初期値をどこにすべきか

1.9. 計算量と誤差

1.9.1 Landauの記号

1.10 実装（計算機イプシロン）

1.10.1 Python

1.10.2 Julia