Approximation — 近似

概要

sangi の近似モジュールは、データや関数を 多項式・有理関数で近似する アルゴリズム群をまとめている。 用途別に次の 5 つに分かれる。

• 回帰 — 観測点 $(x_i, y_i)$ に直線 $y = ax + b$ を最小二乗で当てはめる (linearRegression 系)
• 多項式フィッティング — 観測点に次数 $d$ の多項式を最小二乗で当てはめる (polynomialFit)
• Chebyshev 近似 — 既知の関数 $f$ を区間 $[a, b]$ 上で Chebyshev 多項式級数により準ミニマックス近似する (ChebyshevApprox)
• Padé 近似 — Taylor 係数列から有理関数 $p(x)/q(x)$ を構成し、 単純な打ち切り Taylor より広い範囲で精度を改善する (padeApprox 系)
• 連分数・有理補間 — Taylor 級数を連分数に変換する qd アルゴリズム (taylorToSFraction / taylorToJFraction)、 相異なる標本点での Thiele 型有理補間 (rationalInterpolation)

全テンプレートは namespace sangi に属し、 要素型 T は concept concepts::OrderedField を満たす必要がある (double、 多倍長 Float 等)。

ビルド

#include <math/approx/approximation.hpp>

ヘッダオンリー。 リンクするライブラリは不要。 内部で線形系を解くために線形代数モジュール (algorithms::solve・algorithms::expm) に依存するが、 これらも同梱ヘッダで完結する。

線形回帰

観測点に直線 $y = ax + b$ を最小二乗で当てはめる。 結果は傾き・切片・決定係数 $R^2$ を保持する LinearRegressionResult<T> で返る。

LinearRegressionResult<T>

template<concepts::OrderedField T>
struct LinearRegressionResult {
    T slope;       // 傾き a
    T intercept;   // 切片 b   (y = a*x + b)
    T r_squared;   // 決定係数 R²
};

関数詳細

linearRegression

// (1) x, y の両方を与える版
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(std::span<const T> x, std::span<const T> y);
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(const std::vector<T>& x, const std::vector<T>& y);

// (2) y のみを与える版 (x = 0, 1, ..., n-1 を自動生成)
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(std::span<const T> y);
template<concepts::OrderedField T>
LinearRegressionResult<T> linearRegression(const std::vector<T>& y);

動作: 標準的な最小二乗で $y = a x + b$ の $a, b$ を求める。 平均 $\bar x, \bar y$ を引いてから $S_{xx} = \sum (x_i - \bar x)^2$、 $S_{xy} = \sum (x_i - \bar x)(y_i - \bar y)$ を計算し、 $a = S_{xy}/S_{xx}$、 $b = \bar y - a\,\bar x$ とする。 $y$ のみ版は独立変数を添字 $i = 0, 1, \dots, n-1$ として等間隔データの傾き推定に使う。

パラメータ:

注意: $x$ と $y$ のサイズ不一致、 要素数が 2 未満、 全 $x$ が同一 ($S_{xx} = 0$) の場合は std::invalid_argument を投げる。

// 例: 完全な直線 y = 2x + 1 にフィット
std::vector<double> x = {0, 1, 2, 3, 4};
std::vector<double> y = {1, 3, 5, 7, 9};
auto r = sangi::linearRegression(x, y);
std::cout << r.slope << ", " << r.intercept << ", " << r.r_squared;
// 実行結果: 2, 1, 1   (傾き=2, 切片=1, R²=1)

weightedLinearRegression

template<concepts::OrderedField T>
LinearRegressionResult<T> weightedLinearRegression(
    std::span<const T> x, std::span<const T> y, std::span<const T> weights);
template<concepts::OrderedField T>
LinearRegressionResult<T> weightedLinearRegression(
    const std::vector<T>& x, const std::vector<T>& y, const std::vector<T>& weights);

動作: 各点に重み $w_i$ を付けた最小二乗 $\min \sum w_i (y_i - a x_i - b)^2$。 正規方程式 $\begin{bmatrix} \sum w x^2 & \sum w x \\ \sum w x & \sum w \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} \sum w x y \\ \sum w y \end{bmatrix}$ を $2 \times 2$ の直接解法 (LU) で安定に解く。 $R^2$ は重み付き残差平方和 $1 - \mathrm{SS}_{\mathrm{res}} / \mathrm{SS}_{\mathrm{tot}}$ で算出する。

パラメータ:

注意: $x, y, \mathrm{weights}$ のサイズ不一致・2 点未満・重み和ゼロで std::invalid_argument。 信頼度の異なる観測 (分散重み等) を扱う場合に使う。

多項式フィッティング

観測点 $(x_i, y_i)$ に次数 $d$ の多項式 $p(x) = a_0 + a_1 x + \dots + a_d x^d$ を最小二乗で当てはめる。

PolynomialFitResult<T>

template<concepts::OrderedField T>
struct PolynomialFitResult {
    std::vector<T> coefficients;   // 多項式係数 [a0, a1, ..., aN] (昇べき順)
    T             residual;       // 残差平方和 Σ (y_i - p(x_i))²
};

関数詳細

polynomialFit

template<concepts::OrderedField T>
PolynomialFitResult<T> polynomialFit(
    std::span<const T> x, std::span<const T> y, std::size_t degree);
template<concepts::OrderedField T>
PolynomialFitResult<T> polynomialFit(
    const std::vector<T>& x, const std::vector<T>& y, std::size_t degree);

動作: 正規方程式 $A c = b$ ($A_{ij} = \sum_k x_k^{i+j}$、 $b_i = \sum_k x_k^i y_k$) を組んで LU 直接解法で解く。 数値安定化のため $x, y$ から平均を引いてシフト空間でフィットし、 二項展開で元の空間の係数へ戻す。 残差平方和も同時に返す。

パラメータ:

注意: サイズ不一致・データ点が $\mathrm{degree} + 1$ 未満で std::invalid_argument。 高次にすると正規方程式の条件数が悪化する (Vandermonde 由来) ため、 高次が必要なときは ChebyshevApprox による近似を検討する。

// 例: 放物線 y = x^2 にフィット (degree=2)
std::vector<double> x = {-2, -1, 0, 1, 2};
std::vector<double> y = { 4,  1, 0, 1, 4};
auto r = sangi::polynomialFit(x, y, 2);
// r.coefficients ≈ {0, 0, 1}  (a0=0, a1=0, a2=1 → p(x)=x²)
// r.residual ≈ 0

Chebyshev 近似

既知の関数 $f$ を区間 $[a, b]$ 上で次数 $n$ の Chebyshev 多項式級数 $f(x) \approx \tfrac{c_0}{2} + \sum_{k=1}^{n} c_k\, T_k(y)$ ($y$ は $[a,b] \to [-1,1]$ の線形写像) で展開する。 Chebyshev 節点 (極値点) で関数を標本化し DCT で係数を求めるため、 多項式近似の中でも誤差が区間全体に均等に分布する 準ミニマックス (near-minimax) 近似になる。 単純な Taylor 打ち切りより最大誤差がはるかに小さい。

class ChebyshevApprox<T>

template<concepts::OrderedField T>
class ChebyshevApprox {
public:
    ChebyshevApprox(const std::function<T(T)>& f, T a, T b, std::size_t n);

    T operator()(T x) const;                       // 全項で評価
    T evaluate(T x, std::size_t m) const;          // 先頭 m 項のみで評価 (打ち切り)

    const std::vector<T>& coefficients() const;     // Chebyshev 係数 c_k
    std::size_t          size() const;             // 項数 (= n + 1)
    T                    lower() const;            // 区間下端 a
    T                    upper() const;            // 区間上端 b

    std::vector<T> toPolynomialCoefficients() const; // 通常の多項式係数 [a0,...,aN] へ変換
};

コンストラクタ ChebyshevApprox(f, a, b, n)

動作: 構築時に Chebyshev 節点 $y_k = \cos\!\big(\pi (k + \tfrac12)/N\big)$ ($N = n + 1$) で $f$ を評価し、 DCT-II により Chebyshev 係数 $c_0, \dots, c_n$ を一括計算する。 以後の評価は係数の再利用で高速。

パラメータ:

注意: $a \ge b$ の場合は std::invalid_argument。

operator() / evaluate

動作: operator()(x) は全 $n + 1$ 項を Clenshaw 漸化式で評価する。 evaluate(x, m) は先頭 $m$ 項のみ ($m \le n + 1$) を使う打ち切り評価で、 精度と速度を引き換えにできる ($m = 0$ なら $0$ を返す)。 評価点 $x$ は区間 $[a, b]$ 内が前提 (外挿は誤差が急増)。

パラメータ:

toPolynomialCoefficients

動作: Chebyshev 係数 $c_k$ から通常の単項基底の多項式係数 $[a_0, a_1, \dots, a_n]$ ($p(x) = \sum_j a_j x^j$、 昇べき順) を構成する。 $[a, b] \to [-1, 1]$ の変数変換も込みで展開するため、 結果は元の $x$ についての多項式。 他モジュール (Polynomial・求根等) に渡したいときに使う。

注意: 高次では Chebyshev → 単項基底の変換で条件数が悪化するため、 評価だけが目的なら operator() (Clenshaw) の方が数値的に安定。

// 例: exp(x) を [-1, 1] で 8 次 Chebyshev 近似
sangi::ChebyshevApprox<double> approx(
    [](double x){ return std::exp(x); }, -1.0, 1.0, 8);
double v = approx(0.5);             // ≈ exp(0.5) = 1.6487212707...
std::cout << v;
// 実行結果: 1.64872127  (真値 exp(0.5) と高精度で一致)

Padé 近似

Taylor 係数列から有理関数 $[M/N]$ Padé 近似 $p(x)/q(x)$ ($\deg p = M$, $\deg q = N$, $q(0) = 1$ に正規化) を構成する。 同じ項数の Taylor 打ち切りより収束半径外まで精度が伸び、 極を持つ関数の近似に強い。

PadeResult<T>

template<concepts::OrderedField T>
struct PadeResult {
    std::vector<T> numerator;     // 分子係数 [p0, p1, ..., pM] (昇べき順)
    std::vector<T> denominator;   // 分母係数 [1, q1, ..., qN]  (q0 = 1 に正規化)
    bool          valid;         // 解が得られたか
};

関数詳細

padeApprox

template<concepts::OrderedField T>
PadeResult<T> padeApprox(std::span<const T> taylor, std::size_t M, std::size_t N);
template<concepts::OrderedField T>
PadeResult<T> padeApprox(const std::vector<T>& taylor, std::size_t M, std::size_t N);

動作: Taylor 係数 $a_0, a_1, \dots, a_{M+N}$ から $[M/N]$ Padé を構成する。 まず分母の係数 $q_1, \dots, q_N$ を $N \times N$ の線形系 (LU 直接解法) で求め、 続いて分子 $p_m = a_m + \sum_{l=1}^{\min(m,N)} q_l\, a_{m-l}$ を計算する。 $N = 0$ のときは分母 $1$ で Taylor 係数をそのまま分子とする。

パラメータ:

注意: 係数が $M + N + 1$ 個未満なら std::invalid_argument。 線形系が解けない (病条件) 場合は valid = false を返す。

// 例: exp(x) の Taylor から [2/2] Padé
std::vector<double> taylor = {1.0, 1.0, 0.5, 1.0/6, 1.0/24};  // 1, x, x²/2, ...
auto pa = sangi::padeApprox(taylor, 2, 2);
double v = sangi::evaluatePade(pa, 0.5);   // ≈ exp(0.5)
std::cout << v;
// 実行結果: 1.64872...  (真値 exp(0.5)=1.6487212707 に近い)

evaluatePade

template<concepts::OrderedField T>
T evaluatePade(const PadeResult<T>& result, T x);

動作: padeApprox の結果について、 分子・分母をそれぞれ Horner 法で評価し $p(x)/q(x)$ を返す。

パラメータ:

padeTable

template<concepts::OrderedField T>
std::vector<std::vector<PadeResult<T>>>
padeTable(std::span<const T> taylor, std::size_t Mmax, std::size_t Nmax);
template<concepts::OrderedField T>
std::vector<std::vector<PadeResult<T>>>
padeTable(const std::vector<T>& taylor, std::size_t Mmax, std::size_t Nmax);

動作: $0 \le M \le \mathrm{Mmax}$, $0 \le N \le \mathrm{Nmax}$ の全 $(M, N)$ について $[M/N]$ Padé を計算し、 行 $M$・列 $N$ の二次元表 ($\mathrm{table}[M][N]$) として返す。 対角 $[M/M]$ や反対角 $[M + N = \text{const}]$ を見比べて最良近似を選ぶのに使う。 退化セルは valid = false でそのまま格納される。

パラメータ:

注意: 係数不足で std::invalid_argument。 各セルが $N \times N$ 系を解くため計算量は $O(\mathrm{Mmax} \cdot \mathrm{Nmax} \cdot N^2)$。 対角だけ欲しいなら padeApprox を直接呼ぶ方が安い。

matrixExpPade

template<typename T>
Matrix<T> matrixExpPade(const Matrix<T>& A, int m = 6);

動作: 行列指数 $e^A$ を Scaling-and-Squaring 法 + 対角 Padé 核で計算する (Higham 2005)。 ライブラリの行列指数ルーチン algorithms::expm への薄いラッパであり、 Padé 次数は内部で自動選択される。

パラメータ:

注意: 非正方行列で std::invalid_argument。 引数 m は精度に影響しない。

連分数 (qd アルゴリズム)

Taylor 級数 $f(z) = c_0 + c_1 z + c_2 z^2 + \dots$ を 連分数に変換する。 Rutishauser の qd (quotient-difference) アルゴリズムで係数を求める。 連分数は同じ項数の Padé 対角と等価で、 段ごとに評価を打ち切れるため、 関数値の安定した数値評価に向く。 Stieltjes 型 (S 分数) と Jacobi 型 (J 分数) の 2 形式を提供する。

StieltjesFraction<T> / JacobiFraction<T>

// Stieltjes 連分数:  f(z) = b0 / (1 + a1 z / (1 + a2 z / (1 + a3 z / ...)))
template<concepts::OrderedField T>
struct StieltjesFraction {
    T             b0;          // c_0 (先頭定数、 通常 taylor[0])
    std::vector<T> a;          // a_1, a_2, a_3, ...
    bool          terminated;  // 途中で qd 表が 0 に当たり有限長で完結したか
    int           order;       // 入力 Taylor 係数の個数 - 1
};

// Jacobi 連分数:  f(z) = c0 / (1 - β0 z - α1² z² / (1 - β1 z - α2² z² / ...))
template<concepts::OrderedField T>
struct JacobiFraction {
    T             c0;          // 先頭定数
    std::vector<T> beta;       // β_0, β_1, β_2, ...
    std::vector<T> alpha2;     // α_1², α_2², ...
    bool          terminated;
    int           order;
};

関数詳細

taylorToSFraction / taylorToJFraction

template<concepts::OrderedField T>
StieltjesFraction<T> taylorToSFraction(std::span<const T> taylor, int kmax = -1);
template<concepts::OrderedField T>
StieltjesFraction<T> taylorToSFraction(const std::vector<T>& taylor, int kmax = -1);

template<concepts::OrderedField T>
JacobiFraction<T> taylorToJFraction(std::span<const T> taylor, int kmax = -1);
template<concepts::OrderedField T>
JacobiFraction<T> taylorToJFraction(const std::vector<T>& taylor, int kmax = -1);

動作: Taylor 係数列を qd アルゴリズムにかけ、 S 分数係数 $a_k$ または J 分数係数 $\beta_k, \alpha_k^2$ を生成する。 qd 表が $0$ に当たった (Padé 退化) 時点で打ち切り、 terminated = true を立てる。 J 分数は S 分数 2 段を 1 段にまとめた形で、 直交多項式の三項漸化式に対応する。

パラメータ:

注意: 入力が空なら std::invalid_argument。 $c_0 = 0$ (先頭項ゼロ) では連分数が定義できず terminated = true で空を返す (呼び出し側で平行移動等を要する)。

evalSFraction / evalJFraction

template<concepts::OrderedField T>
T evalSFraction(const StieltjesFraction<T>& sf, T x);
template<concepts::OrderedField T>
T evalJFraction(const JacobiFraction<T>& jf, T x);

動作: 連分数を最下段から積み上げて (bottom-up) 評価する。 evalSFraction は $b_0 / (1 + a_1 x / (1 + a_2 x / \dots))$、 evalJFraction は $c_0 / (1 - \beta_0 x - \alpha_1^2 x^2 / (1 - \beta_1 x - \dots))$ を計算する。 係数列が空なら先頭定数 ($b_0$ / $c_0$) をそのまま返す。

パラメータ:

// 例: log(1+x) の Taylor を S 分数に変換して評価
std::vector<double> taylor = {0.0, 1.0, -0.5, 1.0/3, -0.25, 0.2};  // log(1+x)
auto sf = sangi::taylorToSFraction(taylor);
double v = sangi::evalSFraction(sf, 1.0);   // ≈ log(2) = 0.6931...
std::cout << v;
// 実行結果: 0.6931...  (打ち切り Taylor より広い x で安定)

有理補間

相異なる標本点 $x_i$ での値 $y_i$ に有理関数 $p(x)/q(x)$ を当てはめる多点 Padé (Cauchy / Newton-Padé 型)。 padeApprox が 単一点での Taylor 係数から構成するのに対し、 rationalInterpolation は 複数の標本点 $(x_i, y_i)$ から構成する点が異なる (Taylor 係数版が必要なら padeApprox を使う)。

RationalInterpResult<T>

template<concepts::OrderedField T>
struct RationalInterpResult {
    std::vector<T> numerator;     // [p0, p1, ..., pM] (昇べき順)
    std::vector<T> denominator;   // [1, q1, ..., qN]  (q0 = 1 に正規化)
    bool          valid;
};

関数詳細

rationalInterpolation

template<concepts::OrderedField T>
RationalInterpResult<T> rationalInterpolation(
    const std::vector<T>& x, const std::vector<T>& y,
    std::size_t M, std::size_t N);

動作: 各標本点で $p(x_i) - y_i\, q(x_i) = 0$ となる線形系 ($M + N + 1$ 変数) を組み、 LU 直接解法で分子・分母の係数を一括で求める。 $[M/N]$ 有理関数が $f(x_i) = y_i$ を満たす ($M + N + 1$ がデータ点数)。 浮動小数では補間条件は厳密でないため、 必要なら呼び出し側で残差を確認する。

パラメータ:

注意: $x$ と $y$ のサイズ不一致、 または点数が $M + N + 1$ でないとき std::invalid_argument。 線形系が解けない場合は valid = false。 出典: Stoer & Bulirsch "Introduction to Numerical Analysis" §2.2.4。

evalRational

template<concepts::OrderedField T>
T evalRational(const RationalInterpResult<T>& r, T x);

動作: 分子・分母をそれぞれ Horner 法で評価し $p(x)/q(x)$ を返す。

パラメータ:

// 例: 5 点から [2/2] 有理補間
std::vector<double> x = {-2, -1, 0, 1, 2};
std::vector<double> y = {/* f(x_i) */};
auto r = sangi::rationalInterpolation(x, y, 2, 2);
if (r.valid) {
    double v = sangi::evalRational(r, 0.5);
    std::cout << v;
}

使用例

#include <math/approx/approximation.hpp>
#include <cmath>
#include <iostream>
#include <vector>
using namespace sangi;

int main() {
    // (1) Chebyshev 近似: cos(x) を [0, π] で 10 次近似して評価
    ChebyshevApprox<double> cheb(
        [](double x){ return std::cos(x); }, 0.0, std::acos(-1.0), 10);
    std::cout << "cos(1.0) ≈ " << cheb(1.0) << '\n';
    // 出力: cos(1.0) ≈ 0.540302...  (真値 cos(1) = 0.5403023059)

    // (2) Padé 近似: exp(x) の Taylor から [3/3] 有理近似
    std::vector<double> taylor = {
        1.0, 1.0, 0.5, 1.0/6, 1.0/24, 1.0/120, 1.0/720
    };  // exp(x) の Maclaurin 係数 (7 個 = 3+3+1)
    auto pade = padeApprox(taylor, 3, 3);
    if (pade.valid) {
        std::cout << "exp(0.5) ≈ " << evaluatePade(pade, 0.5) << '\n';
        // 出力: exp(0.5) ≈ 1.648721...  (真値 1.6487212707)
    }

    // (3) 線形回帰: 直線へのフィット
    std::vector<double> xs = {0, 1, 2, 3, 4};
    std::vector<double> ys = {1.1, 2.9, 5.2, 6.8, 9.1};
    auto fit = linearRegression(xs, ys);
    std::cout << "slope=" << fit.slope
              << " intercept=" << fit.intercept
              << " R²=" << fit.r_squared << '\n';
    // 出力例: slope≈2.0 intercept≈1.0 R²≈0.99
}

関連する数学的背景

以下の記事では、近似モジュールの基盤となる数学的概念を解説している。

• 近似理論 — 関数近似の枠組みと最良近似
• 最小二乗法 — 線形回帰・多項式フィットの基礎
• Chebyshev 近似 — Chebyshev 級数による準ミニマックス近似
• ミニマックス近似 — 一様ノルム最小化と等振動定理
• Padé 近似 — 有理関数による近似と連分数