近似: 最小二乗・Chebyshev・ミニマックス近似

概要

「近似」とは、与えられた関数 $f$ あるいはデータ点列を、より扱いやすい簡単な関数 $p$ (多くは多項式) で置き換えることをいう。問題は「どの規範で $f$ と $p$ の近さを測るか」で性格が変わる。 sangi の approximation.hpp は次の 3 つの規範に対応する手法を提供する。

• 2 乗ノルム (最小二乗): 誤差の自乗の総和 $\sum (f-p)^2$ あるいは積分 $\int (f-p)^2\,dx$ を最小化する。計算は線形連立方程式 (正規方程式) を 1 回解くだけで、 外れ値やノイズの影響を平均化するためデータフィッティングに向く。
• $\infty$ ノルム (ミニマックス): 最大絶対誤差 $\max_x |f(x)-p(x)|$ を最小化する。 区間全体で誤差を最も均すため、関数を一定精度で評価したい数値ライブラリ実装に向くが、計算は反復的。
• Chebyshev (準ミニマックス): 第一種 Chebyshev 多項式で級数展開する。 係数を直交性から直接計算でき、打ち切り誤差がほぼ等振動になるため、 ミニマックスの良い代用 (準ミニマックス, near-minimax) になる。

関連 API: linearRegression, weightedLinearRegression, polynomialFit, ChebyshevApprox。

線形回帰と多項式フィッティング (最小二乗)

単回帰の正規方程式

データ点列 $(x_i, y_i)$ に直線 $y = a x + b$ を当てはめる。 残差平方和 $\sum_i (y_i - a x_i - b)^2$ を $a, b$ で偏微分して $0$ と置くと、$2 \times 2$ の正規方程式が出る。 sangi の linearRegression は平均 $\bar{x}, \bar{y}$ と中心化された積和 $S_{xx}, S_{xy}, S_{yy}$ から閉形式で求める:

$$a = \frac{S_{xy}}{S_{xx}}, \qquad b = \bar{y} - a\,\bar{x}, \qquad S_{xy} = \sum_i (x_i - \bar{x})(y_i - \bar{y})$$

多項式フィッティングの正規方程式

次数 $d$ の多項式 $p(x) = \sum_{j=0}^{d} a_j x^j$ を最小二乗で当てはめるときは、 計画行列 (design matrix) $X$ を $X_{ij} = x_i^{\,j}$ と置く。残差平方和を最小化する係数ベクトル $\beta = (a_0,\dots,a_d)^\top$ は 正規方程式を満たす:

$$X^\top X\,\beta = X^\top y$$

sangi の polynomialFit はこのグラム行列 $X^\top X$ とその右辺 $X^\top y$ を べき和 $\sum_k x_k^{\,p}$ から構成し、LU 分解 (algorithms::solve) で解く。

決定係数 $R^2$

当てはまりの良さは決定係数で測る。全平方和 $\mathrm{SS}_{\text{tot}}$ と残差平方和 $\mathrm{SS}_{\text{res}}$ から:

$$R^2 = 1 - \frac{\mathrm{SS}_{\text{res}}}{\mathrm{SS}_{\text{tot}}} = 1 - \frac{\sum_i (y_i - \hat{y}_i)^2}{\sum_i (y_i - \bar{y})^2}$$

$R^2 = 1$ で完全一致、$0$ で平均を予測するのと同程度。単回帰では $R^2 = S_{xy}^2 / (S_{xx} S_{yy})$ と等価で、LinearRegressionResult の r_squared に格納される。

重み付き回帰

標本ごとに信頼度が違う場合は、各残差に重み $w_i$ を掛けた $\sum_i w_i (y_i - a x_i - b)^2$ を最小化する。weightedLinearRegression は次の重み付き正規方程式を組み、$2\times 2$ 連立系を解く:

$$\begin{bmatrix} \sum w_i x_i^2 & \sum w_i x_i \\ \sum w_i x_i & \sum w_i \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum w_i x_i y_i \\ \sum w_i y_i \end{bmatrix}$$

分散の逆数 $w_i = 1/\sigma_i^2$ を重みにすると、ガウス・マルコフの意味で最良線形不偏推定 (一般化最小二乗) になる。

Vandermonde の悪条件に注意

計画行列 $X$ は本質的に Vandermonde 行列であり、次数が上がると グラム行列 $X^\top X$ の条件数が急激に悪化する。等間隔標本でべき基底を直接使うと、高次の係数が桁落ちして無意味な値になりうる。

sangi の polynomialFit は標本の平均 $\mu_x, \mu_y$ を引いて中心化した空間で解き、 最後に二項展開で元の係数へ戻すことで条件数をいくらか改善する。 それでも次数が高い (おおむね $d \gtrsim 6$) と限界があるため、高次では次節の Chebyshev 基底に切り替えるのが本筋である。

Chebyshev 近似

第一種 Chebyshev 多項式

第一種 Chebyshev 多項式 $T_n$ は $x = \cos\theta$ の置換で次のように定義される:

$$T_n(\cos\theta) = \cos(n\theta), \qquad n = 0, 1, 2, \dots$$

三角関数の加法定理から、三項漸化式が得られる:

$$T_0(x) = 1, \quad T_1(x) = x, \quad T_{n+1}(x) = 2x\,T_n(x) - T_{n-1}(x)$$

$T_n$ は区間 $[-1, 1]$ 上で値が $\pm 1$ の間を $n$ 回振動し、 重み $1/\sqrt{1-x^2}$ に関して直交する。この直交性が係数計算の鍵になる。

Chebyshev 係数 (離散直交による計算)

関数 $f$ を Chebyshev 級数 $f(x) \approx \tfrac{1}{2} c_0 + \sum_{j\ge 1} c_j T_j(x)$ で表すとき、係数 $c_j$ は直交性から積分で書ける。実際の計算では、$N$ 個の Chebyshev 節点 (第一種 Chebyshev–Gauss 点) で $f$ を標本化した離散直交和を使う:

$$x_k = \cos\!\left(\frac{\pi (k + \tfrac{1}{2})}{N}\right), \qquad c_j = \frac{2}{N} \sum_{k=0}^{N-1} f(x_k)\,\cos\!\left(\frac{\pi j (k + \tfrac{1}{2})}{N}\right)$$

これは離散コサイン変換 (DCT-II) そのものであり、節点が $f$ の零点に当たらないよう半整数 $k+\tfrac{1}{2}$ をずらしてある点に注意。 ChebyshevApprox はコンストラクタでこの $c_j$ を計算して保持する。

区間 $[a, b]$ への線形変換

Chebyshev 多項式の定義域は $[-1, 1]$ なので、一般の区間 $[a, b]$ には線形変換で写す:

$$y = \frac{2x - a - b}{b - a} \;\in [-1, 1], \qquad x = \frac{(b-a)\,y + (a+b)}{2}$$

ChebyshevApprox(f, a, b, n) は内部でこの変換を行うため、利用者は元の区間 $[a, b]$ の $x$ をそのまま渡せる。

評価 (Clenshaw 漸化) と打ち切り

近似値 $\sum_j c_j T_j(y)$ を計算するとき、$T_j$ を陽に展開すると不安定になる。 ChebyshevApprox は Clenshaw 漸化で安定に評価する。 $d_{m} = d_{m+1} = 0$ から始め、$j = m-1, \dots, 1$ について

$$d_j = 2y\,d_{j+1} - d_{j+2} + c_j, \qquad p(x) = y\,d_1 - d_2 + \tfrac{1}{2} c_0$$

係数 $c_j$ は次数とともに急速に小さくなる ($f$ が滑らかなら指数的) ので、適当な項数 $m$ で打ち切り (economization) しても精度がほとんど落ちない。 evaluate(x, m) は先頭 $m$ 項だけで評価でき、用途に応じて精度とコストを引き換えにできる。

係数を通常の多項式 $\sum_j a_j x^j$ に戻したい場合は、toPolynomialCoefficients() が 三項漸化式で各 $T_k$ をべき基底に展開し、区間変換 $y = \alpha x + \beta$ を二項展開で合成して係数列を返す。

なぜ「準ミニマックス」なのか

滑らかな $f$ を $n$ 次で Chebyshev 級数打ち切りすると、主要誤差項はおおむね $c_{n+1} T_{n+1}(x)$ に比例する。 $T_{n+1}$ は区間全体で $\pm 1$ の間を等しく振動するので、近似誤差 $f - p$ もほぼ等振幅で振動する。 後述するように等振動こそが最良一様近似 (ミニマックス) の特徴なので、Chebyshev 近似は真のミニマックスに非常に近い (典型的に最大誤差で数 % 以内) 誤差を達成する。 これを真のミニマックスと区別して「準ミニマックス (near-minimax)」と呼ぶ。

ミニマックス近似

一様ノルムの最小化

連続関数 $f$ を次数 $n$ の多項式空間 $\mathcal{P}_n$ で近似するとき、最大絶対誤差 ($\infty$ ノルム, 一様ノルム) を最小にする多項式 $p^*$ を最良一様近似 (ミニマックス近似) という:

$$p^* = \arg\min_{p \in \mathcal{P}_n} \;\max_{x \in [a,b]} |f(x) - p(x)|$$

$f$ が連続なら最良近似 $p^*$ は存在し、しかも一意に定まる (Chebyshev の定理)。

Chebyshev の等振動定理

最良近似を特徴づけるのが等振動定理 (equioscillation theorem) である。 $p \in \mathcal{P}_n$ が $f$ の最良一様近似である必要十分条件は、誤差関数 $e(x) = f(x) - p(x)$ が 区間内の少なくとも $n+2$ 個の点 $x_0 < x_1 < \dots < x_{n+1}$ で、 最大絶対誤差 $E = \max_x |e(x)|$ を符号を交互に変えながら達成することである:

$$e(x_k) = (-1)^k\,\sigma\,E, \qquad \sigma = \pm 1, \qquad k = 0, 1, \dots, n+1$$

つまり誤差が振幅 $E$ で交互に「振り切れる」点が $n+2$ 個並ぶ。 前節で見たように Chebyshev 近似の誤差はこの等振動に近い形をしているため、ミニマックスの良い近似になる。

Remez 交換法 (概念)

真のミニマックス多項式を求める古典的な反復法が Remez 交換法 (exchange algorithm) である。概念的には次を繰り返す:

1. $n+2$ 個の基準点 (参照集合) を仮に選ぶ。Chebyshev 節点が良い初期点になる。
2. その点で誤差が等振動する条件 $e(x_k) = (-1)^k E$ を線形連立系として解き、係数と等振幅 $E$ を求める。
3. 得られた $p$ の誤差 $e(x)$ の極値点を探し、参照集合をその極値点へ「交換」する。
4. 参照集合が動かなくなる (等振動が達成される) まで繰り返す。

Remez 法は最良近似に二次的に収束するが、各反復で極値探索を要し、実装が重い。

比較表

大まかな指針: ノイズを含む観測データを当てはめるなら最小二乗 (linearRegression / polynomialFit)。 既知の滑らかな関数を区間全体で高精度に近似し、係数も扱いたいならChebyshev (ChebyshevApprox)。 最大誤差を厳密に最小化したいときだけミニマックス (Remez) を使い、その初期値には Chebyshev 近似を充てる、というのが実務的な流れである。

参考文献

• Trefethen, L. N. (2013). Approximation Theory and Approximation Practice. SIAM.
• Cheney, E. W. (1966). Introduction to Approximation Theory. McGraw-Hill.
• Powell, M. J. D. (1981). Approximation Theory and Methods. Cambridge University Press.
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical Recipes: The Art of Scientific Computing (3rd ed.). Cambridge University Press. (Chebyshev 近似と Clenshaw 漸化)
• Björck, Å. (1996). Numerical Methods for Least Squares Problems. SIAM.