数値計算の基礎

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

浮動小数点数と有効桁数

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

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

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

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

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

桁落ち

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

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

統計計算では、分散の計算がこの影響を受けやすい代表例です。定義式 $\sum(x_i - \bar{x})^2 / (n-1)$ でも、平均が大きく分散が小さいデータでは各偏差 $x_i - \bar{x}$ の計算で桁落ちが起きます。しかし二乗した偏差はすべて非負であり、それらの合計は引き算を含みません。そのため展開式のような1回の引き算で有効桁が壊滅する事態は起きません。ただし個々の偏差 $x_i - \bar{x}$ で失われた有効桁は回復しないため、展開式ほどではないものの精度は劣化します。展開式 $(\sum x_i^2 - n\bar{x}^2) / (n-1)$ では $\sum x_i^2$ と $n\bar{x}^2$ が近い値になる場合に、1回の引き算で大規模な桁落ちが発生します。MIDAS は Welford のオンラインアルゴリズム²を使用しています。平均と偏差平方和をデータ1件ごとに逐次更新するため、展開式のような大きな中間値同士の引き算が発生せず、定義式と同程度の精度を one-pass で達成します。

条件数

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

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

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

条件数が $\kappa$ の問題では、$\log_{10}(\kappa)$ のオーダーで有効桁数が失われます。倍精度の約 15.9 桁からこれを引くと、期待できる有効桁数のおおよその見積もりは次のようになります:

$$\text{有効桁数} \approx 15.9 - \log_{10}(\kappa)$$

これは厳密な下界ではなく、精度の桁数のオーダー見積もりです。実際の精度は問題の次元やアルゴリズムの詳細にも依存します。線形回帰で $\hat\beta$ を求める場合、QR 分解で直接 $X$ を分解するなら $\kappa(X)$ が関連し、正規方程式 $(X'X)^{-1}X'Y$ で解くなら $\kappa(X'X) = \kappa(X)^2$ が関連します。MIDAS は係数推定に QR 分解を使用しているため、精度の主要因は $\kappa(X)$ です。

計画行列

線形回帰モデル $Y = X\beta + \varepsilon$ における $X$ が計画行列です。$Y$ は応答変数のベクトル、$\beta$ は回帰係数のベクトル、$\varepsilon$ は誤差項です。計画行列は $n$ 行 $p$ 列の行列で、各行が 1 つの観測値に対応します。切片項がある場合はすべて 1 の定数列を含み、残りの列が各説明変数に対応します（$p$ は切片を含む総列数）。デザイン行列とも呼ばれます。

OLS（最小二乗法）推定量は残差平方和を最小化する解で、$X$ が列フルランクのとき $\hat\beta = (X'X)^{-1}X'Y$ と書けます（導出の詳細）。$X'$ は $X$ の転置行列です。計画行列の条件数が大きいと、丸め誤差が増幅されて $\hat\beta$ の計算精度が低下します。説明変数間の相関が高い場合に計画行列は ill-conditioned になりやすくなります。

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

多項式回帰への影響

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

数値計算の精度 ページの NIST StRD データセットでこの影響を確認できます。単回帰の Norris では係数の LRE が 12.3 ですが、10次多項式の Filip では 7.3 に低下します。Filip の $\kappa(X'X)$ は約 $10^{14}$（$\kappa(X)$ は約 $10^7$）です。MIDAS は係数推定に QR 分解を使用しているため $\kappa(X)$ が関連し、見積もりでは $15.9 - 7 \approx 8.9$ 桁の精度が期待されますが、実測の LRE は 7.3 です。この差は、見積もり式が問題の次元 $p$（Filip では $p = 11$）による後退代入での誤差蓄積などを反映していないことによるものです。

条件数を下げるには、単項式基底 $1, x, x^2, \ldots$ の代わりに直交多項式基底を使います。切片と直交多項式列のみで構成された計画行列では、列が互いに直交しノルムも等しいため条件数は理論的に 1 になります（浮動小数点の丸めにより数値的にはわずかに 1 からずれます）。MIDAS の Orthogonal Polynomials タブでこの変換を適用できます。

See also

• 数値計算の精度 - NIST データセットを使った MIDAS の精度検証結果
• OLS の基礎 - 線形回帰の数理的背景と計画行列の定式化
• 統計用語集 - 用語の定義

脚注

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

2. Welford, B. P. (1962). Note on a method for calculating corrected sums of squares and products. Technometrics, 4(3), 419–420. ↩