# 条件数と安定性：前進誤差と後退誤差で数値計算の精度を読む

> 問題の感度を表す相対条件数と、アルゴリズムの良さを表す後退誤差を区別し、前進誤差が条件数と後退誤差の積で抑えられることを平均値不等式から証明する。単根の摂動と Wilkinson の多項式、log1p、2 次方程式、分散の計算法を数値で確かめる。
> https://rikai.mugen-giken.com/mathematics/numerical-analysis/conditioning-and-stability

## 0. この記事の要点

計算機で計算した答えが期待とずれたとき、原因は 2 通り考えられます。計算のやり方が悪いのか、それとも問題そのものが入力のわずかな違いに敏感すぎるのか、です。この記事では、この 2 つを切り分けるための物差しを作ります。読み終えると、答えの誤差が「問題のせい」か「計算法のせい」かを見分け、後者なら式の書き換えで直せる理由を説明できるようになります。

- 数値計算の誤差を、**問題の感度**（条件数。入力の誤差が答えで何倍に拡大されるかを表す数）と**アルゴリズムの良さ**（後退誤差。入力をどれだけずらしたことにすれば計算値が正確な答えになるか）の 2 つに分けて考えます。前者は数学の問題そのものの性質で、どんなアルゴリズムでも変えられません。
- 微分可能な問題 $f$ の相対条件数は $\kappa(x)=\|J(x)\|\|x\|/\|f(x)\|$（$J$ は[ヤコビ行列](/mathematics/calculus/multivariable-differentiation)）、1 変数なら $|xf'(x)/f(x)|$ です。成分ごとの相対摂動に対する条件数は $\sum_i|x_i\,\partial_if(x)|/|f(x)|$ です。
- 計算値 $\hat y$ が「少しずらした入力 $x+\Delta x$ に対する正確な答え」と書けるとき、その $\Delta x$ の大きさが**後退誤差**です。**前進誤差 ≲ 条件数 × 後退誤差** が成り立ち、1 次の評価だけでなく、平均値不等式による非漸近の評価も証明できます。
- 浮動小数点数の四則・内積・Horner 法（多項式の値を、かけ算とたし算の繰り返しで求める方法）は**後退安定**（計算値が、少しずらした入力の正確な答えになっていること）です。多項式の単根は係数の摂動に対して $-\varepsilon q(x^\ast)/p'(x^\ast)$ だけ動き（陰関数定理。方程式の解が係数の変化に応じてなめらかに動くことを保証する定理）、Wilkinson の多項式では係数を $2^{-23}$ 動かすだけで 10 個の根が複素数になります。
- 2 次方程式の解の公式、$\log(1+x)$、分散の 1 パス公式は、条件のよい問題を不安定なアルゴリズムで解いている例です。式を書き換えると誤差は単位丸め程度に戻ります。

<Aside type="tip" title="この記事で使う事前知識">

- **単位丸め $u$**: 浮動小数点数で 1 回の演算をしたときの、丸めによる相対誤差の上限。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)
- **Horner 法**: 多項式を $(\cdots(a_nx+a_{n-1})x+\cdots)x+a_0$ の順に計算する方法。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)
- **桁落ち**: ほぼ等しい 2 数を引いたとき、有効数字が大きく減ること。→ [浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)
- **ヤコビ行列**: 多変数関数の各出力を各入力で偏微分して並べた行列で、局所的な一次近似を表す。→ [多変数関数の微分と偏微分](/mathematics/calculus/multivariable-differentiation)
- **全微分可能**: 点の近くで関数が一次式でうまく近似できる、という性質。→ [全微分可能](/mathematics/calculus/multivariable-differentiation#def-differentiable-n)

</Aside>

## 1. 動機

[浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)の章の最後で、Horner 法の計算値は「各係数を相対 $\gamma_{2n}$ だけ動かした多項式の正確な値」であることを示しました（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#thm-horner-backward" />）。それでも $(x-1)^{10}$ を展開した多項式を $x=1.01$ で計算すると、答えは 1 桁も信頼できません。Horner 法は悪くないのに、答えは悪い。この食い違いを説明するのがこの記事の主題です。

答えの誤差（**前進誤差**）は、2 つの量の積で決まります。

1. アルゴリズムが入力をどれだけ「ずらしたことにしたか」（**後退誤差**）。これはアルゴリズムの性質です。
2. 入力のずれが答えにどれだけ拡大されるか（**条件数**）。これは問題の性質です。

この分け方の利点は、それぞれを別々に、しかも定理として評価できることです。後退誤差は浮動小数点数の章のような丸め誤差の解析で、条件数は微分の計算で求まります。そして、条件数が大きい問題で答えが悪いのはアルゴリズムのせいではなく、条件数が小さいのに答えが悪いならアルゴリズムを替えるべきだ、と判断できます。

## 2. 条件数

### 2.1. 定義

以下、$\mathbb{R}^n$ のノルムはユークリッドノルム $\|\cdot\|$、行列のノルムはそれから誘導される作用素ノルム $\|J\|=\sup_{\boldsymbol{h}\ne\boldsymbol{0}}\|J\boldsymbol{h}\|/\|\boldsymbol{h}\|$ とします。成分ごとの相対摂動の大きさを測るため、$\boldsymbol{x}\in\mathbb{R}^n$ と $\Delta\boldsymbol{x}$ について
$$
\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})=\min\{\eta\ge0:|\Delta x_i|\le\eta|x_i|\ (i=1,\ldots,n)\}
$$
（そのような $\eta$ がなければ $\infty$。$x_i=0$ の成分は動かせない）とおきます。

$\eta_{\boldsymbol{x}}$ を使うのは、浮動小数点数では各成分が互いに独立に相対 $u$ 程度ずつ丸められるからです。入力全体を 1 つのノルムで測る条件数と、成分ごとに測る条件数の 2 通りを、次の定義で並べて定めます。どちらも「入力を小さくずらしたとき、出力の相対的な変化が入力の相対的な変化の何倍か」の最悪値（$\sup$）を、ずらしを 0 に近づけた極限で見ています。

<Definition id="def-condition-number" title="相対条件数">
$U\subset\mathbb{R}^n$ を開集合、$f:U\to\mathbb{R}^m$、$\boldsymbol{x}\in U$ で $\boldsymbol{x}\ne\boldsymbol{0}$, $f(\boldsymbol{x})\ne\boldsymbol{0}$ とする。

1. **（ノルムによる）相対条件数**を
   $$
   \kappa(\boldsymbol{x})=\lim_{\varepsilon\to+0}\ \sup\left\{\frac{\|f(\boldsymbol{x}+\Delta\boldsymbol{x})-f(\boldsymbol{x})\|/\|f(\boldsymbol{x})\|}{\|\Delta\boldsymbol{x}\|/\|\boldsymbol{x}\|}:0\lt\|\Delta\boldsymbol{x}\|\le\varepsilon\|\boldsymbol{x}\|\right\}
   $$
   と定める。
2. $m=1$ のとき、**成分ごとの相対条件数**を
   $$
   \kappa_{\mathrm{c}}(\boldsymbol{x})=\lim_{\varepsilon\to+0}\ \sup\left\{\frac{|f(\boldsymbol{x}+\Delta\boldsymbol{x})-f(\boldsymbol{x})|}{|f(\boldsymbol{x})|\,\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})}:0\lt\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})\le\varepsilon\right\}
   $$
   と定める。

（$\varepsilon$ を小さくすると $\sup$ をとる集合が小さくなるので、$\sup$ は $\varepsilon$ について単調非減少で、極限は $[0,\infty]$ に存在する。$\Delta\boldsymbol{x}$ は $\boldsymbol{x}+\Delta\boldsymbol{x}\in U$ の範囲を動く。）
</Definition>

条件数は「入力の相対誤差を何倍に拡大して出力の相対誤差にするか」の、摂動が小さい極限での最悪値です。$\kappa\approx10^k$ なら、入力に相対 $u$ の誤差があるだけで、出力は $k$ 桁ほど失われうる、と読みます。浮動小数点数の入力は最初から相対 $u$ の誤差をもちうるので、成分ごとの条件数は計算機の上の問題の感度を直接表します。

### 2.2. 微分可能な問題の条件数

定義の $\sup$ と極限は、そのままでは計算できません。次の命題は、$f$ が微分できれば条件数が微分（ヤコビ行列）の大きさで書けることを示します。証明では、$f$ を一次近似 $J\boldsymbol{h}$ と残り $r(\boldsymbol{h})$ に分け、残りが十分小さいことを使って、条件数を上と下の両側から $\|J\|\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|$ ではさみます。

<Proposition id="prop-condition-differentiable" title="微分可能な問題の条件数">
<Ref to="def-condition-number" /> の状況で、$f$ が $\boldsymbol{x}$ で全微分可能（<Ref to="mathematics/calculus/multivariable-differentiation#def-differentiable-n" />）で、そのヤコビ行列を $J=J_f(\boldsymbol{x})$ とする。

1. $\kappa(\boldsymbol{x})=\dfrac{\|J\|\,\|\boldsymbol{x}\|}{\|f(\boldsymbol{x})\|}$。とくに $n=m=1$ なら $\kappa(x)=\left|\dfrac{xf'(x)}{f(x)}\right|$。
2. $m=1$ なら $\kappa_{\mathrm{c}}(\boldsymbol{x})=\dfrac{1}{|f(\boldsymbol{x})|}\displaystyle\sum_{i=1}^n|x_i|\left|\frac{\partial f}{\partial x_i}(\boldsymbol{x})\right|$。
</Proposition>

<Proof of="prop-condition-differentiable">
$r(\boldsymbol{h})=f(\boldsymbol{x}+\boldsymbol{h})-f(\boldsymbol{x})-J\boldsymbol{h}$ とおくと、全微分可能性から $\|r(\boldsymbol{h})\|/\|\boldsymbol{h}\|\to0$（$\boldsymbol{h}\to\boldsymbol{0}$）です。$\rho(\varepsilon)=\sup\{\|r(\boldsymbol{h})\|/\|\boldsymbol{h}\|:0\lt\|\boldsymbol{h}\|\le\varepsilon\|\boldsymbol{x}\|\}$ とおくと $\rho(\varepsilon)\to0$（$\varepsilon\to+0$）です。

1. **上からの評価.** $0\lt\|\boldsymbol{h}\|\le\varepsilon\|\boldsymbol{x}\|$ なら、三角不等式と作用素ノルムの定義から $\|f(\boldsymbol{x}+\boldsymbol{h})-f(\boldsymbol{x})\|\le(\|J\|+\rho(\varepsilon))\|\boldsymbol{h}\|$ なので、定義の $\sup$ は $(\|J\|+\rho(\varepsilon))\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|$ 以下です。$\varepsilon\to+0$ として $\kappa(\boldsymbol{x})\le\|J\|\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|$。

   **下からの評価.** $\tau\gt0$ を任意にとり、$\|\boldsymbol{v}\|=1$, $\|J\boldsymbol{v}\|\ge\|J\|-\tau$ となる $\boldsymbol{v}$ を選びます（作用素ノルムは上限なので選べます）。$\boldsymbol{h}=\varepsilon\|\boldsymbol{x}\|\boldsymbol{v}$ とすると $\|f(\boldsymbol{x}+\boldsymbol{h})-f(\boldsymbol{x})\|\ge\|J\boldsymbol{h}\|-\|r(\boldsymbol{h})\|\ge(\|J\|-\tau-\rho(\varepsilon))\|\boldsymbol{h}\|$ なので、$\sup$ は $(\|J\|-\tau-\rho(\varepsilon))\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|$ 以上です。$\varepsilon\to+0$, $\tau\to+0$ として逆向きの不等式を得ます。1 変数では $\|J\|=|f'(x)|$ です。

2. $\boldsymbol{g}=(\partial_1f(\boldsymbol{x}),\ldots,\partial_nf(\boldsymbol{x}))$ とすると $J\boldsymbol{h}=\boldsymbol{g}\cdot\boldsymbol{h}$ です。$S=\sum_i|x_i||g_i|$ とおきます。$0\lt\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})=\eta\le\varepsilon$ なら $\|\Delta\boldsymbol{x}\|\le\eta\|\boldsymbol{x}\|$ なので
   $$
   |f(\boldsymbol{x}+\Delta\boldsymbol{x})-f(\boldsymbol{x})|\le\sum_i|g_i||\Delta x_i|+|r(\Delta\boldsymbol{x})|\le\eta S+\eta\|\boldsymbol{x}\|\rho(\varepsilon)
   $$
   で、$\sup$ は $(S+\|\boldsymbol{x}\|\rho(\varepsilon))/|f(\boldsymbol{x})|$ 以下です。逆に、$g_i\ne0$ の成分は $\Delta x_i=\eta|x_i|\operatorname{sign}(g_i)$、$g_i=0$ の成分は $\Delta x_i=\eta x_i$ とおくと、$\boldsymbol{x}\ne\boldsymbol{0}$ から $\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})=\eta$ で、$\boldsymbol{g}\cdot\Delta\boldsymbol{x}=\eta S$ です。よって $|f(\boldsymbol{x}+\Delta\boldsymbol{x})-f(\boldsymbol{x})|\ge\eta S-\eta\|\boldsymbol{x}\|\rho(\varepsilon)$ で、$\eta=\varepsilon\to+0$ とすれば $\kappa_{\mathrm{c}}(\boldsymbol{x})\ge S/|f(\boldsymbol{x})|$ です。
</Proof>

<Example id="ex-elementary-conditions" title="基本的な関数の条件数">
<Ref to="prop-condition-differentiable" /> から、次の表が得られます（$x,x_1,x_2$ は条件数の分母が $0$ にならない範囲）。

| 問題 | 条件数 | 読み方 |
|---|---|---|
| $x_1+x_2$ | $\kappa_{\mathrm{c}}=\dfrac{\vert x_1\vert +\vert x_2\vert }{\vert x_1+x_2\vert }$ | 同符号なら $1$、打ち消し合うと大きい（桁落ち） |
| $x_1x_2$, $x_1/x_2$ | $\kappa_{\mathrm{c}}=2$ | つねに条件がよい |
| $\sqrt{x}$ | $\kappa=1/2$ | 誤差を半分にする |
| $e^x$ | $\kappa=\vert x\vert $ | $\vert x\vert $ が大きいと悪い |
| $\log x$ | $\kappa=1/\vert \log x\vert $ | $x\approx1$ で悪い |
| $\log(1+x)$（$x$ を入力とみる） | $\kappa=\dfrac{\vert x\vert }{(1+x)\vert \log(1+x)\vert }$ | $x\to0$ で $1$ に近づく |

最後の 2 行は、同じ値 $\log(1+x)$ でも「何を入力とみるか」で条件数が違うことを示しています。$y=1+x$ を入力とすると $\kappa=1/|\log y|\approx1/|x|$ ですが、$x$ そのものを入力とすると $\kappa\approx1$ です。$x$ が与えられているのに $1+x$ を先に計算すると、条件の悪い問題を自分で作ってしまいます（6 節）。
</Example>

## 3. 多項式の根の摂動

### 3.1. 単根の 1 次の変化

<Proposition id="prop-simple-root-perturbation" title="単根の摂動">
$p$, $q$ を実係数の多項式とし、$x^\ast$ を $p$ の単根（$p(x^\ast)=0$, $p'(x^\ast)\ne0$）とする。このとき $\delta\gt0$ と $x^\ast$ の開近傍 $V$、$C^1$ 級の関数 $x:(-\delta,\delta)\to V$ があって、$|\varepsilon|\lt\delta$ のとき $p+\varepsilon q$ の $V$ の中の根は $x(\varepsilon)$ だけであり、
$$
x(0)=x^\ast,\qquad x'(0)=-\frac{q(x^\ast)}{p'(x^\ast)}
$$
が成り立つ。
</Proposition>

<Proof of="prop-simple-root-perturbation">
$F(\varepsilon,x)=p(x)+\varepsilon q(x)$ は $\mathbb{R}^2$ 上の多項式なので $C^1$ 級で、$F(0,x^\ast)=0$, $\partial F/\partial x(0,x^\ast)=p'(x^\ast)\ne0$ です。陰関数定理（<Ref to="mathematics/calculus/inverse-implicit-functions#thm-implicit-function" />）を $n=m=1$ で使うと、$0$ の開近傍 $X$ と $x^\ast$ の開近傍 $V$、$C^1$ 級の $g:X\to V$ があって、$X\times V$ の中の $F=0$ の点はちょうど $(\varepsilon,g(\varepsilon))$ で、$g(0)=x^\ast$ です。$X\supset(-\delta,\delta)$ となる $\delta$ をとり $x=g$ とします。同じ定理の導関数の公式から
$$
x'(0)=-\left(\frac{\partial F}{\partial x}(0,x^\ast)\right)^{-1}\frac{\partial F}{\partial\varepsilon}(0,x^\ast)=-\frac{q(x^\ast)}{p'(x^\ast)}
$$
です。
</Proof>

係数 $a_j$ だけを相対 $\varepsilon$ だけ動かすのは $q(x)=a_jx^j$ の場合で、根の相対変化は 1 次で $\varepsilon\,|a_j(x^\ast)^{j-1}/p'(x^\ast)|$ です。すべての係数を相対 $\varepsilon$ 以内で動かすと（$q=\sum_j\theta_ja_jx^j$, $|\theta_j|\le1$）、根の変化は 1 次で $\varepsilon\sum_j|a_j||x^\ast|^j/|p'(x^\ast)|$ 以下で、この上限は $\theta_j$ の符号を選べば達成されます。

<Example id="ex-wilkinson-polynomial" title="Wilkinson の多項式">
$w(x)=\prod_{k=1}^{20}(x-k)=x^{20}-210x^{19}+\cdots+20!$ の根は $1,\ldots,20$ で、互いに 1 ずつ離れた「よく分かれた」根です。$x^{19}$ の係数を $-210-2^{-23}$ に変える（$q(x)=-x^{19}$, $\varepsilon=2^{-23}$）と、<Ref to="prop-simple-root-perturbation" /> から根 $k$ は 1 次で $\varepsilon k^{19}/w'(k)$ 動きます。$w'(k)=\prod_{j\ne k}(k-j)=(-1)^{20-k}(k-1)!\,(20-k)!$ です。

次のコードは、(i) 1 次の予測 $\varepsilon k^{19}/w'(k)$、(ii) 摂動を $10^{-12}$ 倍に小さくしたときの実際の変化と 1 次の予測の比（mpmath の 50 桁で根を求める）、(iii) 係数を binary64 に丸めて `np.roots` で解いた根の誤差、(iv) すべての係数に相対 $u$ の誤差があるときの 1 次の評価 $u\sum_j|a_j|k^j/|w'(k)|$ を並べます。最後の行は $\varepsilon=2^{-23}$ そのもので摂動したときの結果です。

```python
# Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0（x86-64）
from math import factorial as fac
import numpy as np
import mpmath as mp

mp.mp.dps = 50
c = [1]                                # prod (x - k) の係数（整数, 降べき）
for k in range(1, 21):
    c = [a - k * b for a, b in zip(c + [0], [0] + c)]
def roots(e):                          # x^19 の係数から e を引いた根
    cp = [mp.mpf(a) for a in c]
    cp[1] -= e
    return mp.polyroots(cp, maxsteps=300, extraprec=300)
e0 = mp.mpf(2) ** -23
small = roots(e0 * mp.mpf(10) ** -12)
big = roots(e0)
rn = np.roots(np.array(c, dtype=float))   # 係数を binary64 に丸めて解く
u = 2.0 ** -53
for k in (1, 4, 8, 12, 16, 20):
    dp = (-1) ** (20 - k) * fac(k - 1) * fac(20 - k)  # p'(k)
    pred = e0 * mp.mpf(k) ** 19 / dp     # 1 次の予測 e k^19/p'(k)
    z = min(small, key=lambda z: abs(z - k)) - k
    ratio = mp.re(z) / (pred * mp.mpf(10) ** -12)
    cnd = sum(abs(a) * k ** (20 - i) for i, a in enumerate(c)) / abs(dp)
    err = np.min(np.abs(rn - k))
    print(f"{k:>2} {float(pred):9.2e} {float(ratio):.4f}"
          f" {err:.2e} {u * cnd:.2e}")
nonreal = [z for z in big if abs(mp.im(z)) > 1e-30]
print("2^-23 で実でない根", len(nonreal),
      "虚部の最大", mp.nstr(max(abs(mp.im(z)) for z in big), 3))
```

```text
 1 -9.80e-25 1.0000 1.92e-13 4.66e-14
 4  2.61e-10 1.0000 1.62e-08 2.29e-08
 8  7.12e-03 1.0000 7.21e-04 3.48e-04
12  2.37e+01 1.0000 5.17e-02 3.79e-02
16  2.87e+02 1.0000 5.76e-02 6.29e-02
20  5.14e+00 1.0000 1.96e-04 3.06e-04
2^-23 で実でない根 10 虚部の最大 2.81
```

**検算.** (a) 既知値: 摂動前の根は整数 $k$ で、$\varepsilon=2^{-23}$ のとき 10 個の根が 5 組の共役複素数になり、虚部が最大の組は $16.7307\pm2.8126i$ です（Wilkinson が報告した値と一致します）。(b) 別精度・別手法: 50 桁の mpmath と、係数を binary64 に丸めた `np.roots` の誤差は、係数の丸めだけから予測される (iv) と同じ桁で、比は 0.6〜4 倍の範囲です（`np.roots` 自身の丸め誤差が加わるため完全には一致しません）。(c) 次数: 摂動を $10^{-12}$ 倍にすると、実際の変化と 1 次の予測の比はどの根でも $1.0000$ で、変化が摂動に比例する（1 次の項が支配する）ことがわかります。一方 $\varepsilon=2^{-23}$ では予測される変化が $k=16$ で $287$ と根の間隔 $1$ を大きく超え、1 次の近似は成り立たず、根は複素数に移ります。
</Example>

Wilkinson の多項式は、根がよく分かれていても、**係数で表した**多項式の根は極端に条件が悪くなりうることを示しています。根 $16$ の (iv) の評価は $6.3\times10^{-2}$ で、係数を倍精度で持つだけで小数第 2 位が怪しくなります。固有値を[特性多項式](/mathematics/linear-algebra/eigenvalues#def-charpoly)の根として計算しないのはこのためです（数値線形代数のテーマで、行列のまま固有値を計算する方法を扱います）。

## 4. 前進誤差と後退誤差

### 4.1. 定義

問題 $f$ を解くアルゴリズムが、入力 $\boldsymbol{x}$ に対して計算値 $\hat{\boldsymbol{y}}$ を返したとします。

後退誤差の考え方は、「計算値 $\hat{\boldsymbol{y}}$ が、入力を少しだけずらした問題の正確な答えになっていれば、そのずらし幅がアルゴリズムの誤差だ」というものです。次の定義の $\inf$ は、そのようなずらし方のうち最小のものを選ぶ、という意味です。

<Definition id="def-forward-backward-error" title="前進誤差と後退誤差">
<Ref to="def-condition-number" /> の状況で、計算値 $\hat{\boldsymbol{y}}\in\mathbb{R}^m$ について、

1. **（相対）前進誤差**を $\|\hat{\boldsymbol{y}}-f(\boldsymbol{x})\|/\|f(\boldsymbol{x})\|$ と定める。
2. **（ノルムによる相対）後退誤差**を
   $$
   \eta(\hat{\boldsymbol{y}})=\inf\left\{\frac{\|\Delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|}:\boldsymbol{x}+\Delta\boldsymbol{x}\in U,\ f(\boldsymbol{x}+\Delta\boldsymbol{x})=\hat{\boldsymbol{y}}\right\}
   $$
   と定める（集合が空なら $\infty$）。$\|\Delta\boldsymbol{x}\|/\|\boldsymbol{x}\|$ の代わりに $\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})$ を使ったものを**成分ごとの後退誤差** $\eta_{\mathrm{c}}(\hat{\boldsymbol{y}})$ という。
</Definition>

前進誤差は「答えがどれだけ違うか」、後退誤差は「どれだけ違う問題の正確な答えか」です。後者は、真の答え $f(\boldsymbol{x})$ を知らなくても、$\hat{\boldsymbol{y}}$ を代入して確かめられることがあります（線形方程式の残差など）。

<Definition id="def-backward-stable" title="後退安定と混合安定">
問題 $f$ の定義域の部分集合 $D$ の入力に対し、単位丸め $u$ の浮動小数点演算で計算値 $\hat f(\boldsymbol{x})$ を返すアルゴリズムを考える。$\boldsymbol{x}$ と $u$ によらない定数 $C$（入力の次元 $n$ などには依ってよい）があって、

1. すべての $\boldsymbol{x}\in D$ で後退誤差（ノルムによる、または成分ごと）が $Cu$ 以下のとき、アルゴリズムは（その意味で）**後退安定**であるという。
2. すべての $\boldsymbol{x}\in D$ で $\hat f(\boldsymbol{x})+\Delta\boldsymbol{y}=f(\boldsymbol{x}+\Delta\boldsymbol{x})$, $\|\Delta\boldsymbol{y}\|\le Cu\|\hat f(\boldsymbol{x})\|$, $\|\Delta\boldsymbol{x}\|\le Cu\|\boldsymbol{x}\|$ となる $\Delta\boldsymbol{x}$, $\Delta\boldsymbol{y}$ があるとき、**混合安定**（前進後退安定）であるという。
</Definition>

後退安定なら混合安定です（$\Delta\boldsymbol{y}=\boldsymbol{0}$）。混合安定は「少しずらした入力の正確な答えに近い」ことで、多くのアルゴリズムで期待できる最良の性質です。

<Figure caption="前進誤差と後退誤差。計算値 ŷ を、ずらした入力 x + Δx の正確な答えとみなす。答えのずれ（前進誤差）は、入力のずれ（後退誤差）を条件数倍したもので抑えられる。">
<svg viewBox="0 0 640 220" width="100%" role="img" aria-label="入力空間の x と x+Δx、出力空間の f(x) と ŷ を矢印で結んだ図">
  <ellipse cx="150" cy="110" rx="120" ry="85" fill="none" stroke="currentColor" stroke-width="1.5" />
  <ellipse cx="490" cy="110" rx="120" ry="85" fill="none" stroke="currentColor" stroke-width="1.5" />
  <text x="150" y="214" fill="currentColor" font-size="13" text-anchor="middle">入力</text>
  <text x="490" y="214" fill="currentColor" font-size="13" text-anchor="middle">出力</text>
  <circle cx="120" cy="80" r="4" fill="currentColor" />
  <text x="104" y="72" fill="currentColor" font-size="13" text-anchor="end">x</text>
  <circle cx="150" cy="140" r="4" fill="var(--sl-color-accent)" />
  <text x="140" y="160" fill="var(--sl-color-accent)" font-size="13" text-anchor="end">x + Δx</text>
  <circle cx="460" cy="80" r="4" fill="currentColor" />
  <text x="476" y="72" fill="currentColor" font-size="13">f(x)</text>
  <circle cx="500" cy="140" r="4" fill="var(--sl-color-accent)" />
  <text x="514" y="156" fill="var(--sl-color-accent)" font-size="13">ŷ = f(x + Δx)</text>
  <line x1="126" y1="80" x2="452" y2="80" stroke="currentColor" stroke-width="1.5" />
  <polygon points="452,76 460,80 452,84" fill="currentColor" />
  <text x="290" y="72" fill="currentColor" font-size="12" text-anchor="middle">f（正確な計算）</text>
  <line x1="124" y1="86" x2="492" y2="136" stroke="var(--sl-color-accent)" stroke-width="1.5" stroke-dasharray="6 4" />
  <polygon points="484,138 494,138 486,130" fill="var(--sl-color-accent)" />
  <text x="300" y="132" fill="var(--sl-color-accent)" font-size="12" text-anchor="middle">アルゴリズム</text>
  <line x1="122" y1="86" x2="148" y2="134" stroke="currentColor" stroke-width="1" />
  <text x="100" y="118" fill="currentColor" font-size="12" text-anchor="end">後退誤差</text>
  <line x1="463" y1="86" x2="497" y2="134" stroke="currentColor" stroke-width="1" />
  <text x="560" y="112" fill="currentColor" font-size="12">前進誤差</text>
</svg>
</Figure>

### 4.2. 前進誤差 ≲ 条件数 × 後退誤差

<Theorem id="thm-forward-backward" title="前進誤差と後退誤差の関係">
$U\subset\mathbb{R}^n$ を開集合、$f:U\to\mathbb{R}^m$、$\boldsymbol{x}\in U$, $\boldsymbol{x}\ne\boldsymbol{0}$, $f(\boldsymbol{x})\ne\boldsymbol{0}$ とし、$\hat{\boldsymbol{y}}=f(\boldsymbol{x}+\Delta\boldsymbol{x})$（$\boldsymbol{x}+\Delta\boldsymbol{x}\in U$）とする。

1. **（1 次の評価）** $f$ が $\boldsymbol{x}$ で全微分可能なら
   $$
   \frac{\|\hat{\boldsymbol{y}}-f(\boldsymbol{x})\|}{\|f(\boldsymbol{x})\|}\le\kappa(\boldsymbol{x})\frac{\|\Delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|}+R(\Delta\boldsymbol{x}),\qquad\lim_{\Delta\boldsymbol{x}\to\boldsymbol{0}}\frac{R(\Delta\boldsymbol{x})}{\|\Delta\boldsymbol{x}\|}=0
   $$
   を満たす $R\ge0$（$\boldsymbol{x}$ と $f$ だけで決まり、$\Delta\boldsymbol{x}\ne\boldsymbol{0}$ で定義される）がある。これを $o(\|\Delta\boldsymbol{x}\|)$ の項と呼ぶ。
2. **（非漸近の評価）** 線分 $L=\{\boldsymbol{x}+t\Delta\boldsymbol{x}:0\le t\le1\}$ が $U$ に含まれ、$f$ が $L$ の各点で全微分可能なら
   $$
   \frac{\|\hat{\boldsymbol{y}}-f(\boldsymbol{x})\|}{\|f(\boldsymbol{x})\|}\le\left(\sup_{\boldsymbol{z}\in L}\frac{\|J_f(\boldsymbol{z})\|\,\|\boldsymbol{x}\|}{\|f(\boldsymbol{x})\|}\right)\frac{\|\Delta\boldsymbol{x}\|}{\|\boldsymbol{x}\|}.
   $$
3. **（成分ごと）** $m=1$ で、2 と同じ仮定のもとで
   $$
   \frac{|\hat y-f(\boldsymbol{x})|}{|f(\boldsymbol{x})|}\le\left(\sup_{\boldsymbol{z}\in L}\frac{1}{|f(\boldsymbol{x})|}\sum_{i=1}^n|x_i|\left|\frac{\partial f}{\partial x_i}(\boldsymbol{z})\right|\right)\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x}).
   $$
</Theorem>

<Proof of="thm-forward-backward">
1. $r(\boldsymbol{h})=f(\boldsymbol{x}+\boldsymbol{h})-f(\boldsymbol{x})-J_f(\boldsymbol{x})\boldsymbol{h}$ とすると、全微分可能性（<Ref to="mathematics/calculus/multivariable-differentiation#def-differentiable-n" />）から $\|r(\Delta\boldsymbol{x})\|/\|\Delta\boldsymbol{x}\|\to0$（$\Delta\boldsymbol{x}\to\boldsymbol{0}$）です。三角不等式から $\|\hat{\boldsymbol{y}}-f(\boldsymbol{x})\|\le\|J_f(\boldsymbol{x})\|\|\Delta\boldsymbol{x}\|+\|r(\Delta\boldsymbol{x})\|$ で、両辺を $\|f(\boldsymbol{x})\|$ で割り、<Ref to="prop-condition-differentiable" /> の 1 の $\kappa(\boldsymbol{x})=\|J_f(\boldsymbol{x})\|\|\boldsymbol{x}\|/\|f(\boldsymbol{x})\|$ を使えば、$R(\Delta\boldsymbol{x})=\|r(\Delta\boldsymbol{x})\|/\|f(\boldsymbol{x})\|$ として主張です（$\|f(\boldsymbol{x})\|$ は定数なので $R(\Delta\boldsymbol{x})/\|\Delta\boldsymbol{x}\|\to0$）。
2. 平均値不等式（<Ref to="mathematics/calculus/multivariable-differentiation#thm-mean-value-inequality" />）から、ある $\theta\in(0,1)$ で $\|f(\boldsymbol{x}+\Delta\boldsymbol{x})-f(\boldsymbol{x})\|\le\|J_f(\boldsymbol{x}+\theta\Delta\boldsymbol{x})\|\|\Delta\boldsymbol{x}\|$ です。右辺の $\|J_f\|$ を $L$ 上の上限で置き換え、$\|f(\boldsymbol{x})\|$ で割って $\|\boldsymbol{x}\|$ を掛け・割りすれば主張です。
3. $\eta=\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})$ とし、$\varphi(t)=f(\boldsymbol{x}+t\Delta\boldsymbol{x})$ とおきます。$U$ は開集合なので $\{t:\boldsymbol{x}+t\Delta\boldsymbol{x}\in U\}$ は $[0,1]$ を含む開集合で、連鎖律（<Ref to="mathematics/calculus/multivariable-differentiation#thm-chain-rule-n" />）から $\varphi$ は $[0,1]$ の各点で微分可能、$\varphi'(t)=\sum_i\partial_if(\boldsymbol{x}+t\Delta\boldsymbol{x})\Delta x_i$ です。平均値の定理（<Ref to="mathematics/calculus/mean-value-and-taylor#thm-mvt" />）から、ある $\xi\in(0,1)$ で
   $$
   |\hat y-f(\boldsymbol{x})|=|\varphi(1)-\varphi(0)|=|\varphi'(\xi)|\le\sum_i|\partial_if(\boldsymbol{z})||\Delta x_i|\le\eta\sum_i|x_i||\partial_if(\boldsymbol{z})|
   $$
   （$\boldsymbol{z}=\boldsymbol{x}+\xi\Delta\boldsymbol{x}\in L$）です。上限で置き換えて $|f(\boldsymbol{x})|$ で割れば主張です。
</Proof>

1 は「前進誤差 ≈ 条件数 × 後退誤差」の正確な意味です。2 と 3 は $o(\cdot)$ の項を含まない評価で、誤差の上界を保証したいとき（計算機援用証明など）に使えます。とくに $f$ が入力について 1 次（線形）なら、2・3 の上限は $\boldsymbol{x}$ での値と等しく、評価は $\kappa\times$ 後退誤差そのものになります。

この定理を後退安定なアルゴリズムに当てはめると、**後退安定なアルゴリズムの前進誤差は、およそ $\kappa\cdot Cu$ 以下**です。したがって、

- $\kappa$ が小さい問題を後退安定なアルゴリズムで解けば、答えは正確です。
- $\kappa u\gtrsim1$ の問題は、どんなアルゴリズムでも（入力の丸めだけで）答えが保証されません。精度を上げるか、問題の定式化を変えます。
- 前進誤差が $\kappa u$ より桁違いに大きいなら、そのアルゴリズムは後退安定ではありません（6 節の分散の 1 パス公式）。

## 5. 基本演算の後退安定性

<Proposition id="prop-basic-backward-stable" title="四則・内積・Horner 法の後退安定性">
オーバーフローとアンダーフローは起きないとする。

1. $x,y\in\mathbb{F}$ について、$\mathrm{fl}(x\pm y)=x(1+\delta)\pm y(1+\delta)$, $\mathrm{fl}(xy)=x(1+\delta)\,y$, $\mathrm{fl}(x/y)=x(1+\delta)/y$（$|\delta|\le u$）。したがって四則は成分ごとの後退誤差 $u$ 以下で後退安定である。
2. $\boldsymbol{x},\boldsymbol{y}\in\mathbb{F}^n$、$nu\lt1$ とし、内積を前から順に計算した値を $\hat s$ とする。$\hat s=(\boldsymbol{x}+\Delta\boldsymbol{x})^{\mathsf{T}}\boldsymbol{y}$, $|\Delta\boldsymbol{x}|\le\gamma_n|\boldsymbol{x}|$（成分ごと）となる $\Delta\boldsymbol{x}$ がある。すなわち内積は $\boldsymbol{x}$ について成分ごとの後退誤差 $\gamma_n$ 以下である。
3. Horner 法で計算した $\sum_{i=0}^na_ix^i$ の値は、係数 $(a_0,\ldots,a_n)$ について成分ごとの後退誤差 $\gamma_{2n}$ 以下である（$2nu\lt1$）。
</Proposition>

<Proof of="prop-basic-backward-stable">
1. 標準モデル（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#cor-standard-model-ops" />）から $\mathrm{fl}(x\circ y)=(x\circ y)(1+\delta)$, $|\delta|\le u$ で、右辺を書き直したものです。
2. 内積の丸め誤差の定理（<Ref to="mathematics/computer-assisted-proofs/floating-point-interval-arithmetic#thm-dot-product" />）から $\hat s=\sum_ix_iy_i(1+\theta^{(i)})$, $|\theta^{(i)}|\le\gamma_n$ なので、$\Delta x_i=x_i\theta^{(i)}$ とおけばよいです。
3. Horner 法の後退誤差（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#thm-horner-backward" />）そのものです。
</Proof>

<Proposition id="prop-horner-forward" title="多項式の値の条件数と Horner 法の前進誤差">
$p(x)=\sum_{i=0}^na_ix^i$ の値を係数 $\boldsymbol{a}=(a_0,\ldots,a_n)$ の関数とみたとき、$p(x)\ne0$ ならその成分ごとの条件数は
$$
\kappa_p(x)=\frac{\sum_{i=0}^n|a_i||x|^i}{|p(x)|}
$$
であり、$x,a_i\in\mathbb{F}$, $2nu\lt1$ で Horner 法の途中でオーバーフローもアンダーフローも起きなければ、計算値 $\hat p$ は
$$
\frac{|\hat p-p(x)|}{|p(x)|}\le\gamma_{2n}\,\kappa_p(x)
$$
を満たす。
</Proposition>

<Proof of="prop-horner-forward">
$\boldsymbol{a}\mapsto\sum_ia_ix^i$ は 1 次式で $\partial/\partial a_i=x^i$ なので、<Ref to="prop-condition-differentiable" /> の 2 から $\kappa_p(x)$ は上の形です。<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#thm-horner-backward" /> から $\hat p=\sum_ia_i(1+\theta_i)x^i$, $|\theta_i|\le\gamma_{2n}$ なので、$|\hat p-p(x)|=|\sum_ia_i\theta_ix^i|\le\gamma_{2n}\sum_i|a_i||x|^i$ で、$|p(x)|$ で割れば主張です（<Ref to="thm-forward-backward" /> の 3 を 1 次式に使ったものと同じです）。
</Proof>

[浮動小数点数と丸め誤差](/mathematics/numerical-analysis/floating-point-arithmetic)の演習で扱った $(x-1)^{10}$ の展開は $\kappa_p(1.01)\approx10^{23}$ で、$\gamma_{20}\kappa_p\gg1$ でした。Horner 法は後退安定でも、問題（展開した係数から値を求めること）の条件が悪すぎたわけです。

<Example id="ex-outer-product" title="外積の計算は後退安定でない">
$\boldsymbol{x},\boldsymbol{y}\in\mathbb{F}^n$ の外積 $A=\boldsymbol{x}\boldsymbol{y}^{\mathsf{T}}$ を $\hat a_{ij}=\mathrm{fl}(x_iy_j)=x_iy_j(1+\delta_{ij})$ で計算すると、各成分の相対誤差は $u$ 以下で、前進誤差は小さく混合安定です。しかし後退安定ではありません。もし $\hat A=(\boldsymbol{x}+\Delta\boldsymbol{x})(\boldsymbol{y}+\Delta\boldsymbol{y})^{\mathsf{T}}$ と書けるなら $\hat A$ の階数は 1 以下ですが、$\delta_{ij}$ はそろわないので、$\hat A$ は一般に階数 2 以上になるからです。たとえば binary64 で $\boldsymbol{x}=(1,\mathrm{fl}(1/3))$, $\boldsymbol{y}=(\mathrm{fl}(0.1),\mathrm{fl}(0.7))$ とすると、第 1 行は厳密で、第 2 行の 2 成分の相対誤差は $-1.4\times10^{-17}$ と $1.6\times10^{-17}$ と異なり、有理数で計算した $\det\hat A$ は約 $6.9\times10^{-19}\ne0$ です。後退安定性は「出力が入力の集合の像に入るか」という性質なので、出力の次元が入力より大きい問題では成り立たないことがよくあります。
</Example>

## 6. 条件のよい問題を安定に解く書き換え

ここでは、問題の条件数は小さいのに、素朴な式が不安定なアルゴリズムになっている例を 3 つ見ます。どれも、桁落ち（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#prop-cancellation" />）が丸め誤差を拡大しています。

<Example id="ex-quadratic-formula" title="2 次方程式の解の公式">
$x^2+bx+c=0$（$b^2\gt4c$, $c\ne0$）の解を $x_\pm=(-b\pm\sqrt{b^2-4c})/2$ で計算します。$b^2\gg4|c|$ だと $\sqrt{b^2-4c}\approx|b|$ なので、$-b$ と同じ符号の平方根を足す方の解は問題ありませんが、他方は桁落ちを起こします。binary64 で $b=10^8$, $c=1$ とすると、小さい方の解の真値は約 $-1.0\times10^{-8}$ ですが、$(-b+\sqrt{b^2-4c})/2$ は $-7.45\times10^{-9}$ になり、相対誤差は $0.25$ です。

安定な計算法は、桁落ちしない方の解 $x_1=(-b-\operatorname{sign}(b)\sqrt{b^2-4c})/2$ を先に求め、もう一方を解と係数の関係 $x_1x_2=c$ から $x_2=c/x_1$ とするものです。$x_2$ の計算は掛け算と割り算だけで、相対誤差は数 $u$ に収まります（上の例で $7.9\times10^{-17}$）。一方、根そのものの条件数は、$c$ に関して $|\partial x_2/\partial c|\cdot|c/x_2|=|c/(x_2(x_2-x_1))|=|x_1/(x_2-x_1)|\approx1$、$b$ に関して $|\partial x_2/\partial b|\cdot|b/x_2|=|b/(x_2-x_1)|\approx1$（$|x_2|\ll|x_1|$ のとき。$\partial x_2/\partial c=-1/(x_2-x_1)$, $\partial x_2/\partial b=-x_2/(x_2-x_1)$）で、問題は条件がよいことがわかります。
</Example>

<Example id="ex-log1p" title="log1p と expm1">
<Ref to="ex-elementary-conditions" /> のとおり、$x$ を入力とする $\log(1+x)$ の条件数は $x\to0$ で $1$ に近づきます。ところが `np.log(1 + x)` は、まず $\mathrm{fl}(1+x)=(1+x)(1+\delta)$ を作り、条件数 $\approx1/|x|$ の問題 $\log y$ を解いてしまうので、相対誤差は $u/|x|$ 程度になります。`np.log1p(x)` は $x$ から直接 $\log(1+x)$ を計算する関数です。$e^x-1$ と `np.expm1` も同じ関係です。

```python
# Python 3.12.10 / NumPy 2.5.2 / mpmath 1.3.0（x86-64）
import numpy as np
import mpmath as mp

mp.mp.dps = 50
for k in (1, 3, 5, 7, 9, 11, 13, 15, 17):
    x = 10.0 ** -k
    ref = mp.log1p(mp.mpf(x))            # x は binary64 の値
    e_log = abs(np.log(1 + x) - ref) / ref
    e_l1p = abs(np.log1p(x) - ref) / ref
    ref2 = mp.expm1(mp.mpf(x))
    e_exp = abs((np.exp(x) - 1) - ref2) / ref2
    e_em1 = abs(np.expm1(x) - ref2) / ref2
    print(f"1e-{k:<2} {float(e_log):.2e} {float(e_l1p):.2e}"
          f" {float(e_exp):.2e} {float(e_em1):.2e}")
```

```text
1e-1  7.32e-16 4.00e-18 7.75e-16 1.68e-17
1e-3  1.10e-13 3.90e-17 4.29e-14 2.60e-17
1e-5  6.55e-12 1.46e-16 9.70e-12 3.11e-17
1e-7  5.84e-10 7.36e-17 5.66e-10 2.95e-17
1e-9  8.27e-08 7.49e-17 8.22e-08 7.44e-17
1e-11 8.27e-08 8.06e-17 8.27e-08 8.06e-17
1e-13 7.99e-04 1.08e-16 7.99e-04 1.78e-17
1e-15 1.10e-01 1.06e-16 1.10e-01 1.06e-16
1e-17 1.00e+00 5.00e-18 1.00e+00 5.00e-18
```

列は $x$、`log(1+x)`、`log1p`、`exp(x)-1`、`expm1` の相対誤差です。

**検算.** (a) 既知値: $\log(1+x)=x-x^2/2+\cdots$ なので参照値は $x\to0$ で $x$ に近づきます。参照値は同じ binary64 の $x$ に対する mpmath の 50 桁の値です。(b) 別手法: `log1p` と `expm1` の誤差は全行で $1.5u\approx1.7\times10^{-16}$ 以下（$u\approx1.1\times10^{-16}$ を超えるのは $x=10^{-5}$ の `log1p` だけ）で、条件数 $\approx1$ と後退安定性から予想される大きさです。(c) 次数: 素朴な式の誤差は $u/x$ の程度で、$x$ を $10^{-2}$ 倍すると誤差はおおむね $10^2$ 倍になります（$10^{-3}$ から $10^{-7}$ まで）。$x=10^{-9},10^{-11}$ で同じ値になり、$x=10^{-13}$ で跳ねるのは、$\mathrm{fl}(1+x)-1$ が $2u$ の整数倍にしかならず、誤差が $u/x$ の上限の中で $x$ の 2 進展開に応じて揺れるためです。$x=10^{-17}\lt u$ では $\mathrm{fl}(1+x)=1$ で、答えは $0$ になります。
</Example>

<Example id="ex-variance" title="分散の 3 つの計算法">
データ $x_1,\ldots,x_n$ の標本分散 $s^2=\frac1{n-1}\sum_i(x_i-\bar x)^2$ を、次の 3 通りで計算します。

- **1 パス公式**: $s^2=\frac1{n-1}\bigl(\sum x_i^2-\frac1n(\sum x_i)^2\bigr)$。データを 1 回読むだけで済みます。
- **2 パス公式**: 平均 $\bar x$ を先に求め、$\frac1{n-1}\sum(x_i-\bar x)^2$ を計算します。
- **Welford 法**: $M_k=M_{k-1}+(x_k-M_{k-1})/k$、$S_k=S_{k-1}+(x_k-M_{k-1})(x_k-M_k)$ を更新し、$s^2=S_n/(n-1)$ とします。1 パスで済みます。

$\partial s^2/\partial x_i=2(x_i-\bar x)/(n-1)$（$\sum_j(x_j-\bar x)=0$ を使う）なので、<Ref to="prop-condition-differentiable" /> の 2 から成分ごとの条件数は $\kappa=2\sum_i|x_i||x_i-\bar x|/((n-1)s^2)$ で、平均 $m$ が標準偏差 $\sigma$ より大きいと $\kappa\approx2|m|\cdot\mathrm{mean}|x_i-\bar x|/s^2\sim|m|/\sigma$ 程度になります。1 パス公式は $\sum x_i^2\approx nm^2$ と $(\sum x_i)^2/n$ の差をとるので、丸め誤差が $m^2/\sigma^2$ 倍に拡大されます。$\kappa u$ より桁違いに大きい誤差は、<Ref to="thm-forward-backward" /> により後退誤差が $u$ 程度では説明できず、1 パス公式は後退安定ではありません。

乱数種を固定した標準正規乱数 $10^4$ 個に $m=10^3,10^5,10^7$ を足したデータ（binary64）について、3 つの相対誤差と $\kappa u$ を並べます。真の分散はデータの浮動小数点数を有理数として厳密に計算します。

```python
# Python 3.12.10 / NumPy 2.5.2（x86-64）／乱数種 20261006
from fractions import Fraction
import numpy as np

rng = np.random.default_rng(20261006)
z = rng.standard_normal(10**4)
u, n = 2.0 ** -53, len(z)
for m in (1e3, 1e5, 1e7):
    x = m + z                                   # 平均 m, 分散約 1
    xs = x.tolist()
    q = [Fraction(t) for t in xs]               # 厳密な分散（有理数）
    mq = sum(q) / n
    v = sum((t - mq) ** 2 for t in q) / (n - 1)
    one = (np.sum(x * x) - np.sum(x) ** 2 / n) / (n - 1)
    two = np.sum((x - np.mean(x)) ** 2) / (n - 1)
    mean, m2 = 0.0, 0.0                          # Welford 法
    for k, t in enumerate(xs, 1):
        d = t - mean
        mean += d / k
        m2 += d * (t - mean)
    wel = m2 / (n - 1)
    kap = 2 * sum(abs(t) * abs(t - mq) for t in q) / ((n - 1) * v)
    e = [abs(Fraction(t) - v) / v for t in (one, two, wel)]
    print(f"{m:.0e} " + " ".join(f"{float(t):.2e}" for t in e)
          + f" {float(kap) * u:.2e}")
```

```text
1e+03 4.30e-10 1.03e-16 2.69e-14 1.78e-13
1e+05 1.12e-06 6.05e-17 1.93e-12 1.78e-11
1e+07 7.57e-03 6.72e-17 2.39e-11 1.78e-09
```

列は $m$、1 パス公式、2 パス公式、Welford 法の相対誤差、$\kappa u$ です。

**検算.** (a) 既知値: 参照値はデータの有理数による厳密な分散です。(b) 別手法: 2 パス公式と Welford 法の誤差はどの $m$ でも $\kappa u$ 以下で、後退安定なアルゴリズムに期待される範囲です。2 パス公式がとくに正確なのは、$x_i$ と $\bar x$ の比が $1/2$ と $2$ の間にあるので $x_i-\bar x$ が Sterbenz の補題（<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#lem-sterbenz" />）により厳密に計算されるからです（$\bar x$ の誤差は残ります）。(c) 次数: $m$ を $100$ 倍するごとに、$\kappa u$ は $100$ 倍（$m$ に比例）、1 パス公式の誤差は $2.6\times10^3$ 倍・$6.8\times10^3$ 倍で、$m^2$ の比 $10^4$ に近い増え方をし、$m=10^7$ では $\kappa u$ の $4\times10^6$ 倍です。Welford 法の誤差は $\kappa u$ の $1/7$ から $1/70$ の範囲で増えます。
</Example>

## 7. 演習

<Exercise id="exr-tan-condition" difficulty="易">
$f(x)=\tan x$ の相対条件数を求め、$x$ が $\pi/2$ に近いときと $0$ に近いときの値を見積もれ。
<Solution>
$f'(x)=1/\cos^2x$ なので $\kappa(x)=|x/(\sin x\cos x)|=|2x/\sin2x|$ です。$x\to0$ では $1$ に近づき、条件がよい問題です。$x=\pi/2-h$（$h$ は小さい正の数）では $\sin2x=\sin(\pi-2h)=\sin2h\approx2h$ なので $\kappa\approx\pi/(2h)$ です。たとえば $h=10^{-6}$ なら $\kappa\approx1.6\times10^6$ で、入力の丸め（相対 $u$）だけで約 6 桁が失われます。
</Solution>
</Exercise>

<Exercise id="exr-sum-backward" difficulty="標準">
3 つの浮動小数点数の和を $\hat s=\mathrm{fl}(\mathrm{fl}(a+b)+c)$ で計算する。成分ごとの後退誤差が $\gamma_2$ 以下であることを示し、前進誤差の評価 $|\hat s-s|/|s|\le\gamma_2(|a|+|b|+|c|)/|a+b+c|$ を <Ref to="thm-forward-backward" /> の 3 から導け（オーバーフロー・アンダーフローはなく、$a+b+c\ne0$ とする）。
<Solution>
標準モデルから $\hat s=((a+b)(1+\delta_1)+c)(1+\delta_2)=a(1+\delta_1)(1+\delta_2)+b(1+\delta_1)(1+\delta_2)+c(1+\delta_2)$ で、誤差因子の積の補題から各係数は $1+\theta$, $|\theta|\le\gamma_2$ の形です。よって $\hat s=f(\boldsymbol{x}+\Delta\boldsymbol{x})$（$f(a,b,c)=a+b+c$）で $\eta_{\boldsymbol{x}}(\Delta\boldsymbol{x})\le\gamma_2$ です。$f$ は 1 次式で $\partial_if=1$ なので、<Ref to="thm-forward-backward" /> の 3 の上限は $(|a|+|b|+|c|)/|a+b+c|$ で、主張の評価が出ます。これは $\kappa_{\mathrm{c}}\times$ 後退誤差の形です。
</Solution>
</Exercise>

<Exercise id="exr-sqrt-difference" difficulty="標準">
大きな $x\gt0$ に対して $f(x)=\sqrt{x+1}-\sqrt{x}$ を計算する。(1) $f$ の相対条件数を求め、問題の条件がよいことを示せ。(2) 素朴な式の相対誤差がおよそ $u\cdot4x$ になる理由を説明し、安定な式を与えよ。
<Solution>
(1) $f'(x)=\frac12(1/\sqrt{x+1}-1/\sqrt{x})=-f(x)/(2\sqrt{x}\sqrt{x+1})$ なので $\kappa(x)=x/(2\sqrt{x(x+1)})\lt1/2$ です。(2) $\sqrt{x+1}$ と $\sqrt{x}$ はそれぞれ相対 $u$ 程度の誤差をもち、<Ref to="mathematics/numerical-analysis/floating-point-arithmetic#prop-cancellation" /> の倍率 $(\sqrt{x+1}+\sqrt{x})/f(x)=(\sqrt{x+1}+\sqrt{x})^2\approx4x$ で拡大されます。有理化した $f(x)=1/(\sqrt{x+1}+\sqrt{x})$ は引き算がなく、どの段も相対誤差が $u$ 程度なので、全体の相対誤差も数 $u$ に収まります。
</Solution>
</Exercise>

## 参考文献

- N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — 第 1 章（前進誤差・後退誤差・条件数、分散の計算）、第 5 章（多項式）。
- L. N. Trefethen, D. Bau III, *Numerical Linear Algebra*, SIAM, 1997 — Lecture 12（条件数）、Lecture 14・15（安定性と後退安定性）。
- J. H. Wilkinson, *Rounding Errors in Algebraic Processes*, Prentice-Hall, 1963 — 多項式の根の感度と Wilkinson の多項式。
- B. P. Welford, "Note on a method for calculating corrected sums of squares and products", *Technometrics* 4 (1962), 419–420. [doi:10.1080/00401706.1962.10490022](https://doi.org/10.1080/00401706.1962.10490022)
- T. F. Chan, G. H. Golub, R. J. LeVeque, "Algorithms for computing the sample variance: analysis and recommendations", *The American Statistician* 37 (1983), 242–247. [doi:10.1080/00031305.1983.10483115](https://doi.org/10.1080/00031305.1983.10483115)