数値計算の基礎

数値計算の精度 ページで使われている概念の背景です。

浮動小数点数と有効桁数

コンピュータは実数を有限桁の近似値として扱います。MIDAS を含むすべてのブラウザアプリケーションは IEEE 754 倍精度浮動小数点数を使用します。

倍精度浮動小数点数は実数を (1)s×(1+f)×2e(-1)^s \times (1 + f) \times 2^e の形式で表現します。ss は符号 1 ビット、ee は指数部 11 ビットで値のスケールを決め、ff は仮数部 52 ビットで 0f<10 \le f < 1 の小数部分を表します。正規化数1では (1+f)(1 + f) の整数部分が必ず 1 になるため、この 1 は格納せず暗黙に補います。格納する 52 ビットと暗黙の 1 ビットを合わせた 53 ビットが有効精度です。10進で約 15.9 桁に相当します:

log10(253)15.95\log_{10}(2^{53}) \approx 15.95

この制約から、単一の倍精度浮動小数点数で表せる有効桁数は最大で約 15.9 桁です。数値計算の精度 ページの LRE(Log Relative Error: 一致する有効桁数に相当する指標)がこの値に近いほど、計算が倍精度の限界まで正確であることを意味します。

実数を浮動小数点数に変換するとき、有限ビットで表現できない端数は丸められます。この丸め誤差は個々の演算では微小ですが、計算を重ねると蓄積する場合があります。蓄積の程度はアルゴリズムとデータの性質に依存し、後述する桁落ちや条件数がその代表的な要因です。

桁落ち

桁落ちは、近い値の浮動小数点数同士を引き算したときに有効桁数が大幅に減少する現象です。

a=1.23456789012345a = 1.234567890123\mathbf{45}b=1.23456789012300b = 1.234567890123\mathbf{00} はどちらも 15 桁の有効数字を持ちますが、差 ab=0.00000000000045a - b = 0.00000000000045 の有効数字は 2 桁です。先頭の 0 は桁の位置を示しているだけで有効数字に含まれないため、意味のある情報は末尾の 45 だけです。aabb で一致していた 13 桁分の情報は相殺で消え、もともと不確かだった下位桁だけが残ります。

統計計算では、分散の計算がこの影響を受けやすい代表例です。定義式 (xixˉ)2/(n1)\sum(x_i - \bar{x})^2 / (n-1) は各観測値と平均の偏差を計算するため、近い値同士の引き算にはなりません。一方、展開式 (xi2nxˉ2)/(n1)(\sum x_i^2 - n\bar{x}^2) / (n-1) では xi2\sum x_i^2nxˉ2n\bar{x}^2 が近い値になる場合に桁落ちが発生します。平均が大きく分散が小さいデータでこの問題が顕著になります。

条件数

条件数は、入力の小さな摂動が出力にどれだけ増幅されるかを測る指標です。正則な行列 AA の条件数は次のように定義されます:

κ(A)=AA1\kappa(A) = \|A\| \cdot \|A^{-1}\|

A\|A\| は行列の 2-ノルム(最大特異値)です。AA1AA1=I=1\|A\| \cdot \|A^{-1}\| \ge \|AA^{-1}\| = \|I\| = 1 であるため、条件数の下限は 1 です。条件数が 1 に近い行列を well-conditioned、大きい行列を ill-conditioned と呼びます。明確な閾値はなく、求める精度との兼ね合いで判断します。

線形方程式 Ax=bAx = b を浮動小数点演算で解くとき、bb の相対摂動が δ\delta なら、解 xx の相対誤差は最大で κ(A)δ\kappa(A) \cdot \delta に増幅されます。条件数が κ\kappa なら最大で log10(κ)\log_{10}(\kappa) 桁の精度が失われるため、倍精度の約 15.9 桁から引いて、残る有効桁数の下限はおおよそ次のようになります:

有効桁数15.9log10(κ)\text{有効桁数} \gtrsim 15.9 - \log_{10}(\kappa)

この式は数値的に安定なアルゴリズムを使った場合の下限です。不安定なアルゴリズムではこれより多くの桁を失います。

計画行列

線形回帰モデル Y=Xβ+εY = X\beta + \varepsilon における XX が計画行列です。YY は応答変数のベクトル、β\beta は回帰係数のベクトル、ε\varepsilon は誤差項です。計画行列は nnpp 列の行列で、各行が 1 つの観測値、各列が 1 つの説明変数に対応します。切片項がある場合、すべて 1 の定数列が加わります。デザイン行列とも呼ばれます。

OLS(最小二乗法)推定量は β^=(XX)1XY\hat\beta = (X'X)^{-1}X'Y と定義されます。XX'XX の転置行列です。計画行列の条件数が大きいと、丸め誤差が増幅されて β^\hat\beta の計算精度が低下します。これは推定量の統計的な分散が増大する多重共線性の問題とは別の、数値計算上の問題です。説明変数間の相関が高い場合に計画行列は ill-conditioned になりやすく、OLS の基礎 で扱っている VIF は関連する現象を個々の変数について測る指標です。

計画行列の定式化と OLS 推定量の性質は OLS の基礎 で詳しく扱っています。

多項式回帰への影響

多項式回帰では計画行列の列が 1,x,x2,,xd1, x, x^2, \ldots, x^d になります。次数 dd が高くなると累乗列の間が強く相関し、計画行列は ill-conditioned になります。

数値計算の精度 ページの NIST StRD データセットでこの影響を確認できます。単回帰の Norris では係数の LRE が 12.3 ですが、10次多項式の Filip では 7.3 に低下します。Filip の計画行列の条件数は約 101410^{14} で、多項式の次数に伴う条件数の増大が精度を制限しています。

条件数を下げるには、単項式基底 1,x,x2,1, x, x^2, \ldots の代わりに直交多項式基底を使います。直交多項式基底では計画行列の列が互いに直交し、条件数が約 1 に低下します。MIDAS の Orthogonal Polynomials タブでこの変換を適用できます。

See also

脚注

  1. 正規化とは、整数部分が 1 になるように指数を調整して表現を一意にする操作です。たとえば 0.50.51.0×211.0 \times 2^{-1} と表現されます。ただし指数部には下限があるため、0 に十分近い値では指数が下限に達し、整数部分を 1 にできません。このような数を非正規化数と呼び、暗黙の 1 が適用されないため有効精度が下がります。統計計算で扱う値の範囲では問題になりません。