数値計算の基礎

説明変数間の相関が高いとき、高次の多項式を当てはめるとき、平均が大きく分散が小さいデータを扱うときに、計算の途中で精度が失われることがあります。このページはその仕組みと MIDAS の対策を説明します。測定値そのものの精度ではなく、計算過程で失われる精度の話です。数値計算の精度 ページで使われている概念の背景でもあります。

浮動小数点数と有効桁数

コンピュータは実数を有限桁の近似値として扱います。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: 正解が分かっているデータで MIDAS の計算値と正解値が一致する有効桁数に相当する指標）がこの値に近いほど、計算が倍精度の限界まで正確であることを意味します。LRE の判断基準は同ページで扱います。

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

桁落ち

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

$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 で達成します。ただし平均が極端に大きく分散が極小なデータでは桁落ちの影響が残り、標準偏差の精度が低下する場合があります。その水準は 数値計算の精度 ページの NumAcc3・NumAcc4 で確認できます。

条件数

条件数は、入力の小さな摂動が出力にどれだけ増幅されるかを測る指標です。正則な行列 $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)$ です。

MIDAS は QR 分解によるモデルの当てはめで条件数の推定値を計算し、推定値が $10^{10}$ を超えた場合は結果に警告を表示します。この推定値は QR 分解の上三角因子 $R$ の Frobenius 条件数 $\|R\|_F \cdot \|R^{-1}\|_F$ で、真の条件数 $\kappa(X)$ に対して $\kappa(X) \le \|R\|_F \cdot \|R^{-1}\|_F \le p\,\kappa(X)$（$p$ は切片を含む計画行列の列数）を満たします。条件数を過小評価しないため、悪条件の計画行列を見逃しません。過大評価は最大で $p$ 倍にとどまり、有効桁数の見積もりでは $\log_{10} p$ 桁の差です。

警告が出ても係数は計算・表示されます。表示された推定値は失われた有効桁数を踏まえて解釈します。悪条件の主な原因は説明変数間の高い相関で（計画行列）、多項式回帰では直交多項式基底への切り替えで改善できます（多項式回帰への影響）。

計画行列

線形回帰モデル $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 になりやすくなります。回帰の前に説明変数間の相関は Statistics タブの相関行列で、当てはめ後の多重共線性は Linear Regression の係数テーブルの VIF で確認できます。

計画行列の定式化と 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)$ が約 $2 \times 10^{15}$、正規方程式に現れる $\kappa(X'X)$ はその二乗の約 $10^{30}$ です。見積もり式に $\kappa(X)$ を当てはめると期待できる有効桁数は 1 桁未満ですが、実測の LRE は 7.3 でこれを大きく上回ります。

見積もりが悲観的になるのは、QR 分解の丸め誤差が計画行列の各列のスケールに対して相対的に入るためです。各列を $\ell_2$ ノルムが 1 になるようスケールしてから条件数を計算すると、Filip では約 $5 \times 10^9$ となり、見積もりは約 6 桁で実測に近づきます。

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

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. ↩

このページの Markdown 版もあります。