数値計算の基礎

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

浮動小数点数と有効桁数

MIDAS は実数を IEEE 754 倍精度浮動小数点数で近似します。この形式は、実数を 2 進の仮数と指数の組で表し、仮数部として 52 ビットを格納します。0 でない値は、仮数の先頭桁が 1 になるように指数を調整した形で格納されます。先頭桁は必ず 1 になるため格納の必要がなく、暗黙に補われます。格納する 52 ビットと暗黙の 1 ビットを合わせた有効ビット数は 53 で、10 進で約 15.9 桁に相当します:

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

これが単一の倍精度浮動小数点数で表せる有効桁数の上限です。

数値計算の精度 ページの精度検証では、正解が分かっているデータで MIDAS の計算値が正解値と何桁一致するかを、LRE（Log Relative Error）という指標で報告しています。LRE がこの上限に近いほど、計算の過程で失われた有効桁が少ないことを意味します。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^2 - n\bar{x}^2) / (n-1)$ は、平均が大きく分散が小さいデータで $\sum x_i^2$ と $n\bar{x}^2$ が近い値になり、1 回の引き算で有効桁の大半を失います。定義式 $\sum(x_i - \bar{x})^2 / (n-1)$ でも各偏差 $x_i - \bar{x}$ の計算で桁落ちは起きますが、二乗した偏差の合計は引き算を含まないため、失われる有効桁は個々の偏差の計算で失われた分に留まり、展開式のように大半を一度に失うことはありません。

MIDAS は Welford のオンラインアルゴリズム（Welford, 1962）で分散を計算します。平均と偏差平方和をデータ 1 件ごとに更新するため、展開式のような大きな中間値同士の引き算が発生しません。さらに先頭の値でシフトしてから逐次更新し、大きなオフセットによる逐次平均の精度劣化も抑えます。これにより、データを 1 回走査するだけで定義式と同程度の精度を達成します。

Welford のアルゴリズムが防ぐのは、展開式の 1 回の引き算と逐次平均の劣化までです。個々の偏差での桁落ちは定義式と同じように残るため、平均が極端に大きく分散が極小なデータでは標準偏差の精度が低下します。どの程度低下するかは 数値計算の精度 ページの NumAcc3・NumAcc4 で確認できます。

計画行列

線形回帰モデル $Y = X\beta + \varepsilon$ における $X$ が計画行列（design matrix）です。$Y$ は応答変数のベクトル、$\beta$ は回帰係数のベクトル、$\varepsilon$ は誤差項です。計画行列は $n$ 行 $p$ 列の行列で、各行が 1 つの観測値に対応します。切片項がある場合はすべて 1 の定数列を含み、残りの列が各説明変数に対応します（$p$ は切片を含む総列数）。定式化と OLS（最小二乗法）推定量の導出は OLS の基礎 で扱っています。

回帰係数 $\hat\beta$ は計画行列から計算されるため、その計算精度は計画行列の数値的性質に依存します。説明変数間の相関が高いほど、次節で説明する条件数が大きくなり、丸め誤差が増幅されます。回帰の前に説明変数間の相関は Statistics タブの Relationships セクション で、当てはめ後の多重共線性は Linear Regression と GLM の係数テーブルの VIF で確認できます。

条件数

条件数は、入力の小さな摂動が出力にどれだけ増幅されるかを測る指標です。正則な行列 $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$ を求める場合、$X$ を直接 QR 分解するなら $\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$ 桁の差です。

Linear Regression・二元配置 ANOVA・DoE では、解く問題が重みなしの最小二乗なので、この推定値は $\kappa(X)$ に対応します。GLM は IRLS（反復再重み付け最小二乗）で当てはめるため、各反復の重み付き計画行列 $\sqrt{W}\,X$ の条件数を推定し、反復を通じた最大値で警告を判定します（$W$ は反復ごとの重み行列）。Gaussian・identity リンクの組み合わせでは重みが一定になり推定値は $\kappa(X)$ に一致しますが、ロジスティック回帰やポアソン回帰では重みが観測ごとに大きく異なるため、$X$ の列相関だけでは警告の原因を説明できないことがあります。

警告が出ても係数は計算・表示されます。表示された係数は、失われた有効桁数を踏まえて解釈します。

悪条件は、説明変数間の高い相関と、説明変数の値の大きさという、性質の異なる 2 つの要因から生じます。相関は 計画行列 の節で述べた Relationships セクションと VIF で確認できます。

値の大きさが条件数を押し上げる形は 2 つあります。説明変数の平均が 0 から大きく離れていると、その列は切片の列と疑似的に相関します。また、説明変数どうしで値の分散が大きく異なると、その差だけで条件数が上がります。

どちらの要因が原因かは、各説明変数を平均 0・分散 1 に標準化して当てはめ直すと切り分けられます。標準化は値の大きさの影響だけを取り除くため、標準化後も警告が出るなら原因は相関で、警告が出なくなるなら原因は値の大きさです。

原因が値の大きさなら、条件数からの見積もりほどの精度低下は実際には起きておらず、元のモデルをそのまま使えることが多いです¹。理由と具体例は 多項式回帰への影響 で説明します。

多項式回帰への影響

多項式回帰では計画行列の列が $1, x, x^2, \ldots, x^d$ になります。次数 $d$ が高くなると、累乗列の間が強く相関し、列ごとの値の大きさも大きく開くため、前節の 2 つの要因がともに働いて計画行列は悪条件になります。

数値計算の精度 ページの NIST StRD データセットでこの影響を確認できます。単回帰の Norris では係数の LRE が 12.3 ですが、10次多項式の Filip では 7.3 に低下します。

条件数からの有効桁数の見積もりは最悪ケースを想定したオーダーの目安で、実際の精度はこれを上回ることも多いです。Filip の計画行列は $\kappa(X)$ が約 $2 \times 10^{15}$ です。見積もり式に当てはめると期待できる有効桁数は 1 桁未満ですが、実測の LRE は 7.3 でこれを大きく上回ります。

見積もりが悲観的になるのは、MIDAS が係数推定に使う QR 分解の丸め誤差が、計画行列の各列のスケールに対して相対的に生じるためです。列の間のスケール差は条件数を押し上げますが、その増分は実際の精度低下には効きません。スケールを揃える前の条件数は、スケール差の分だけ計算の難しさを過大評価します。各列を $\ell_2$ ノルムが 1 になるようスケールしてから条件数を計算すると、Filip では約 $5 \times 10^9$ となり、見積もりは約 6 桁で実測に近づきます。相関だけによる悪条件では列のスケールがもともとそろっているため、この過大評価は起こりません。

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

参考文献

• Welford, B. P. (1962). Note on a method for calculating corrected sums of squares and products. Technometrics, 4(3), 419-420. https://www.jstor.org/stable/1266577

See also

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

脚注

1. 例外は、平均が標準偏差に比べて極端に大きい説明変数が原因の場合です。桁落ち の節で述べた偏差の情報の損失が既に起きており、標準化しても失われた情報は戻りません。 ↩

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