Root Finding — 求根アルゴリズム

概要

sangi の求根モジュールは 3 つのカテゴリに分かれる。

• 1D 求根 — $f(x) = 0$ を満たす実数 $x$ を求める。二分法から Brent 法まで 10+ 種類
• 多項式求根 — $p(x) = 0$ の全根 (複素根含む) を求める。低次は公式解、高次は反復法
• N-D 求根 — $\mathbf{F}(\mathbf{x}) = \mathbf{0}$ を満たすベクトル $\mathbf{x}$ を求める

ビルド

// 全モジュールまとめてインクルード
#include <math/roots/root_finding.hpp>

// 個別インクルードも可能
#include <math/roots/root_finding_1d.hpp>   // 1D 求根
#include <math/roots/root_finding_nd.hpp>   // N-D 求根
#include <math/roots/polynomial_roots.hpp>  // 多項式求根

ヘッダオンリー。リンクするライブラリは不要。

結果型と収束判定

RootFindingResult<T>

1D 求根関数および多項式根反復関数の戻り値型。 math/roots/root_finding_base.hpp 定義。

template<typename T>
struct RootFindingResult {
    std::optional<T> root;     // 成功時: 根の近似値、失敗時: std::nullopt
    bool             converged; // 収束したかどうか
    size_t           iterations;// 実際に行った反復回数
    T                error_estimate; // 推定誤差 (絶対値、収束時は最終ステップ幅 等)
};

使い方:

auto r = brent_method(f, a, b);
if (r.converged && r.root) {
    double x = *r.root;            // または r.root.value()
    // ... x を使う
} else {
    // 失敗 (区間に符号変化なし / max iter 到達 / 数値特異 等)
    std::cerr << "no root: iter=" << r.iterations
              << " err=" << r.error_estimate << '\n';
}

N-D 求根関数は同じ構造の SystemRootFindingResult<V, T> を返す (root がベクトル型 std::optional<V> である点だけ異なる)。

ConvergenceCriteria<T>

収束判定の閾値を保持する struct。 全求根関数の最後の引数として渡せる (省略時はデフォルト)。

template<concepts::OrderedField T>
struct ConvergenceCriteria {
    T      abs_ftol = std::numeric_limits<T>::epsilon() * 100;   // |f(x)| 絶対許容誤差
    T      abs_xtol = std::numeric_limits<T>::epsilon() * 100;   // |Δx|   絶対許容誤差
    T      rel_ftol = std::numeric_limits<T>::epsilon() * 1000;  // |f|相対許容誤差
    T      rel_xtol = std::numeric_limits<T>::epsilon() * 1000;  // Δx 相対許容誤差
    size_t max_iterations = 100;                                  // 反復上限
};

使用例:

// 緩い許容誤差で素早く打ち切る
ConvergenceCriteria<double> loose;
loose.abs_ftol = 1e-6;
loose.max_iterations = 50;
auto r = newton_raphson(f, df, x0, loose);

// 多倍長 (Float) で高精度
ConvergenceCriteria<Float> tight;
tight.abs_ftol = Float("1e-50");
tight.max_iterations = 200;
auto r2 = brent_method<Float>(f, a, b, tight);

1D: 囲い込み法 (Bracketing)

区間 $[a, b]$ 内で $f(a)$ と $f(b)$ が異符号であることを前提とし、確実に収束する。

関数詳細

bisection

template<concepts::OrderedField T>
RootFindingResult<T> bisection(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 区間 $[a, b]$ を毎ステップ半分に絞り、$f$ の符号が変わる側を残す。 線形収束 (毎反復で誤差が 1/2)。

パラメータ: f 1 変数実関数 / a, b 初期区間端点 ($f(a) \cdot f(b) < 0$ が必須) / criteria 収束判定。

返り値: RootFindingResult<T>。 区間に符号変化がない場合は converged = false、 root = std::nullopt。

用途: ロバスト性が最優先で速度は二の次のとき。 区間幅が確実に半減するため最悪反復回数を予測可能 (log2((b-a)/tol))。

auto r = bisection<double>([](double x){ return x*x - 2.0; }, 1.0, 2.0);
if (r.converged) std::cout << *r.root;
// 実行結果 (build & run 検証済): iter=42  root=1.41421356237311
//   ※ 真値 sqrt(2) = 1.41421356237310 と 12 桁一致

false_position

template<concepts::OrderedField T>
RootFindingResult<T> false_position(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 区間端点を結ぶ直線と $x$ 軸の交点を次の試行値とする。 凹な関数では一方の端点が動かず線形になる弱点がある (Illinois / Anderson-Björk 等の改良法はここでは未実装)。

パラメータ: bisection と同じ。 $f(a) \cdot f(b) < 0$ が必須。

ridders_method

template<concepts::OrderedField T>
RootFindingResult<T> ridders_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 区間の中点 $m$ と $e^{Q(x-m)}$ で重み付けした false position を組合せる。 1 反復で $f$ を 2 回評価する代わりに $\sqrt{2} \approx 1.414$ 次の超線形収束を得る (実効収束次数は約 1.84)。

用途: bisection より速く、 Brent より実装が単純。 関数評価が安価な場合の中庸選択。

brent_method

template<concepts::OrderedField T>
RootFindingResult<T> brent_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 逆 2 次補間 / 割線法 / 二分法を内部状態に応じて切替える Brent (1973) のハイブリッド。 区間幅が縮まないステップは bisection で強制縮退するため、 速度と保証収束を両立する。

用途: 1D 求根のデフォルト推奨。 SciPy brentq 等の標準アルゴリズム。 区間端点で符号が異なれば常に確実に根を返す。

auto r = brent_method<double>([](double x){ return std::cos(x) - x; }, 0.0, 1.0);
if (r.converged) std::cout << *r.root;
// 実行結果 (build & run 検証済): iter=26  root=0.739085133215161
//   ※ Dottie 数 = 0.7390851332151607 と 14 桁一致

toms748

template<concepts::OrderedField T>
RootFindingResult<T> toms748(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Alefeld-Potra-Shi (1995) の TOMS Algorithm 748。 逆 3 次補間と double-length step を組合せて Brent より高い実効収束次数 ($\approx 1.84$ for smooth $f$) を持つ。

用途: Brent と同じく区間保証付きで、 滑らかな $f$ では Brent より少ない反復で収束。 SciPy scipy.optimize.brentq の代替に推奨。

1D: 非区間法 (Non-bracketing)

初期値 $x_0$ (と任意で導関数) から反復する。収束は速いが初期値依存。

関数詳細

newton_raphson

template<concepts::OrderedField T>
RootFindingResult<T> newton_raphson(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 各反復 $x_{n+1} = x_n - f(x_n) / f'(x_n)$。 単根近傍で 2 次収束 ($|x_{n+1} - x^*| \le C |x_n - x^*|^2$)。 重根では線形に劣化。

パラメータ: f 関数 / df 導関数 / x0 初期値 / criteria 収束判定。

注意: $|f'(x)| < \epsilon \cdot 10$ の near-stationary 点で発散判定 → converged=false。 初期値が悪いと別の根に収束、 または発散することがある。 区間保証がほしい場合は brent_method を推奨。

// 例: f(x) = x^2 - 2 の根を Newton 法で求める
auto r = newton_raphson<double>(
    [](double x) { return x*x - 2.0; },
    [](double x) { return 2.0 * x; },
    1.0);
if (r.converged) std::cout << *r.root;
// 実行結果 (build & run 検証済): iter=5  root=1.41421356237310
//   ※ 真値 sqrt(2) = 1.41421356237310 と 15 桁一致 (5 反復で収束、 2 次収束の威力)

secant_method

template<concepts::OrderedField T>
RootFindingResult<T> secant_method(
    const std::function<T(T)>& f,
    T x0, T x1,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 2 点 $(x_{n-1}, x_n)$ を結ぶ割線 (secant) と $x$ 軸の交点を次の試行値: $x_{n+1} = x_n - f(x_n) \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})}$。 導関数を 1 次差分で代用するため Newton よりわずかに遅い ($\approx 1.618$ 次収束、 黄金比)。

パラメータ: f 関数 / x0, x1 初期 2 点 / criteria 収束判定。 $x_0, x_1$ は近すぎず、 根を挟む必要はない。

注意: $|f(x_n) - f(x_{n-1})| < \epsilon$ (傾きほぼゼロ) で停滞検出 → converged=false。

// 例: f(x) = x - cos(x) の根 (= Dottie 数)
auto r = secant_method<double>(
    [](double x) { return x - std::cos(x); }, 0.0, 1.0);
if (r.converged) std::cout << *r.root;
// 実行結果 (build & run 検証済): iter=6  root=0.739085133215161
//   ※ Dottie 数 0.7390851332151607 と 14 桁一致

1D: 高次収束法

関数詳細

halley_method

template<concepts::OrderedField T>
RootFindingResult<T> halley_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 更新式 $x_{n+1} = x_n - \frac{2 f f'}{2 f'^2 - f f''}$。 Newton の 2 次収束を 3 次に引き上げる Householder 1 次。 2 階導関数が必要、 計算コストは Newton の約 2 倍だが反復回数は減る。

パラメータ: f, df, d2f 関数・1 階・2 階導関数 / x0 初期値 / criteria 収束判定。

注意: 分母 $2 f'^2 - f f''$ がほぼゼロで発散判定。

// 例: f(x) = x^3 - 2 の根 (= 2^(1/3))
auto r = halley_method<double>(
    [](double x) { return x*x*x - 2.0; },
    [](double x) { return 3.0 * x * x; },
    [](double x) { return 6.0 * x; },
    1.0);
if (r.converged) std::cout << *r.root;
// 実行結果 (build & run 検証済): iter=3  root=1.25992104989487
//   ※ 真値 2^(1/3) = 1.2599210498948732 と 14 桁一致 (3 反復、 3 次収束の威力)

schroder_method

template<concepts::OrderedField T>
RootFindingResult<T> schroder_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    const std::function<T(T)>& d2f,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 更新式 $x_{n+1} = x_n - \frac{f f'}{f'^2 - f f''}$。 重根に対して通常の Newton が 1 次に劣化するのに対し、 Schröder (1870) は 重根でも 2 次収束を維持する設計。

用途: 多項式の重根を扱う場合、 Halley より安定。 単根では Halley (3 次) の方が速い。

king_method

template<concepts::OrderedField T>
RootFindingResult<T> king_method(
    const std::function<T(T)>& f, T a, T b,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: King-Pegasus 法 (R.F.King, BIT 1973)。 Regula Falsi で補間点を求めた後、 同じ端点が連続停滞した場合に $f_1 \leftarrow f_1 \cdot f_2 / (f_2 + f)$ で修正して再補間。 Illinois 法より速い超線形収束。

パラメータ: f 関数 / a, b 端点 ($f(a) \cdot f(b) < 0$ 必須) / criteria 収束判定。

popovski_method

template<concepts::OrderedField T>
RootFindingResult<T> popovski_method(
    const std::function<T(T)>& f,
    const std::function<T(T)>& df,
    T x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 2 ステップ更新 (Newton ステップ + 補正項) で 1 階導関数のみから 4 次収束を実現。 King-Werner / Ostrowski 系の改良。 2 階導関数を計算しない代わり 1 反復あたり $f$ を 2 回評価する。

用途: 2 階導関数が高コスト (e.g. AD で計算) なケースで、 反復回数より関数評価回数の方が問題な時。

多項式求根

全関数は std::vector<Complex<T>> を返す (実根は虚部 = 0)。 根はベクトル空間の元ではなく単なる複素数の集合なので、 sangi::Vector (線形代数用) ではなく STL の std::vector を使う。 入力 Polynomial<T> は係数を 昇べき順で保持 (p[0] = 定数項)。

低次 (公式解)

高次 (反復法)

関数詳細

solveLinear / solveQuadratic / solveCubic / solveQuartic

template<typename T> std::vector<Complex<T>> solveLinear(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveQuadratic(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveCubic(const Polynomial<T>& p);
template<typename T> std::vector<Complex<T>> solveQuartic(const Polynomial<T>& p);

動作: それぞれ 1/2/3/4 次方程式に対する閉形式解。 次数不一致は assert で検出。

桁落ち対策: solveQuadratic は判別式が小さい場合に大きい絶対値の根を直接計算し、 もう一方を Vieta の公式で導出 (Numerical Recipes 推奨)。 出典: 奥村「アルゴリズム事典」p205。

// 例: x^2 - 5x + 6 = (x-2)(x-3) = 0
Polynomial<double> p({6.0, -5.0, 1.0});  // 昇べき: c0=6, c1=-5, c2=1
auto roots = solveQuadratic(p);
for (auto& r : roots) std::cout << r.re << ' ';
// 実行結果 (build & run 検証済): "3 2" (真値 {2, 3} と完全一致)

jenkinsTraub

template<typename T>
std::vector<Complex<T>> jenkinsTraub(const Polynomial<T>& poly);

動作: Jenkins-Traub (1970) の RPOLY アルゴリズム (3 段階 shift-shift-shift)。 各根を 1 つずつ求め、 デフレーション (因数除去) で次数を下げる。 同分野で数十年にわたり標準的な「最も信頼できる」多項式求根法。

用途: 中次数 (5-19 次) のデフォルト推奨。 solvePolynomial も 5-19 次でこれを内部呼出。

// 例: x^6 - 1 = 0 (1 の 6 乗根)
Polynomial<double> p({-1.0, 0, 0, 0, 0, 0, 1.0});
auto roots = jenkinsTraub(p);
// 実行結果 (build & run 検証済): 6 根 = {±1, ±exp(±iπ/3)}
//   |p(root)| 残差最大 4.97e-16 (16 桁精度)

laguerre

template<typename T>
std::vector<Complex<T>> laguerre(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: Laguerre 法 + デフレーション。 各反復で $L(z) = -n / \left( G \pm \sqrt{(n-1)(n H - G^2)} \right)$ 更新 ($G = p'/p$, $H = G^2 - p''/p$)。 任意の初期値から全根に大域収束する珍しい性質、 3 次収束。

パラメータ: eps 収束許容誤差 / maxIter 反復上限。

durandKernerAberth

template<typename T>
std::vector<Complex<T>> durandKernerAberth(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: Aberth-Ehrlich の全根同時反復 ($z_i \leftarrow z_i - p(z_i)/p'(z_i) / (1 - p(z_i)/p'(z_i) \cdot \sum_{j \ne i} 1/(z_i - z_j))$)。 収束後に元の多項式で各根を Newton polish 精密化する。 全根同時のため デフレーション誤差が蓄積しない。

用途: 高次 (20 次以上) のデフォルト推奨。 solvePolynomial も 20 次以上で内部呼出。

bairstow

template<typename T>
std::vector<Complex<T>> bairstow(
    const Polynomial<T>& poly,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

動作: 試行 2 次因子 $x^2 + rx + s$ で多項式を割り、 Newton 法で $r, s$ を修正して 2 次因子を抽出。 実数演算のみで複素根 (共役ペア) を得る。 出典: Press et al., Numerical Recipes §9.5。

用途: 複素数演算を避けたい組込み環境、 実係数のみで複素根の組をペア単位で扱いたい場合。

統合ディスパッチャ

template<typename T>
std::vector<Complex<T>> solvePolynomial(
    const Polynomial<T>& p,
    T eps = std::numeric_limits<T>::epsilon() * T(100),
    size_t maxIter = 1000);

solvePolynomial は次数に応じて自動的に最適なアルゴリズムを選択する:

• 1 次: solveLinear (閉形式)
• 2 次: solveQuadratic (閉形式、 桁落ち対策付)
• 3 次: solveCubic (Cardano + 三角関数)
• 4 次: solveQuartic (Ferrari)
• 5-19 次: jenkinsTraub (1 根ずつ + デフレーション、 高精度)
• 20 次以上: durandKernerAberth (全根同時 + Newton polish、 デフレーション誤差回避)

N-D: 多変数求根

$\mathbf{F}: \mathbb{R}^n \to \mathbb{R}^n$ の零点を求める。 戻り値は SystemRootFindingResult<V, T> (1D の RootFindingResult のベクトル版)。

関数詳細

newton_raphson_system

template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
SystemRootFindingResult<V, T> newton_raphson_system(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 各反復で $\mathbf{x}_{n+1} = \mathbf{x}_n - J(\mathbf{x}_n)^{-1} F(\mathbf{x}_n)$。 ヤコビアン特異時は対角要素を abs_xtol * 10 で正則化して solve をリトライ。

パラメータ: F 系を表すベクトル値関数 / J ヤコビアンを返す関数 / x0 初期ベクトル / criteria 収束判定。

返り値: SystemRootFindingResult<V, T>。 root は std::optional<V>。

収束判定: $\|F(\mathbf{x})\|_2 < \mathrm{abs\_ftol}$ または $\|\Delta\mathbf{x}\|_2$ が許容値以下。

例: 単位円 $x^2 + y^2 = 1$ と直線 $y = x$ の交点を求める。 sangi::Matrix<T>::solve() を経由して concept 要件を満たすため <math/linalg/solvers.hpp> も include する。

#include <math/roots/root_finding.hpp>
#include <math/linalg/solvers.hpp>   // Matrix::solve() の本体定義はここ

using namespace sangi;

auto F = [](const Vector<double>& v) {
    Vector<double> r(2);
    r[0] = v[0]*v[0] + v[1]*v[1] - 1.0;   // x^2 + y^2 - 1
    r[1] = v[0] - v[1];                   // x - y
    return r;
};
auto J = [](const Vector<double>& v) {
    Matrix<double> m(2, 2);
    m(0,0) = 2*v[0]; m(0,1) = 2*v[1];
    m(1,0) = 1.0;    m(1,1) = -1.0;
    return m;
};
Vector<double> x0(2); x0[0] = 0.5; x0[1] = 0.5;
auto r = newton_raphson_system<double, Vector<double>, Matrix<double>>(F, J, x0);
if (r.converged && r.root) {
    const auto& v = *r.root;
    std::cout << v[0] << ", " << v[1];
}
// 実行結果 (build & run 検証済): iter=5  root=(0.707106781186548, 0.707106781186548)
//   ※ 真値 1/√2 = 0.7071067811865476 と 15 桁一致

broyden_method

template<concepts::OrderedField T, typename V, typename M>
SystemRootFindingResult<V, T> broyden_method(
    const std::function<V(const V&)>& F,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: 「Good Broyden」式の準 Newton 法。 ヤコビアンを毎回計算せず、 ランク 1 更新で近似 ($B_{n+1} = B_n + \frac{(\Delta F - B_n \Delta \mathbf{x}) \Delta \mathbf{x}^\top}{\|\Delta \mathbf{x}\|^2}$)。

用途: ヤコビアンが解析的に書けない / 計算コストが高い場合。 初期 $B_0$ は単位行列、 後は更新で改善。

auto r = broyden_method<double>(F, x0);  // J 不要、 F のみ与える

levenberg_marquardt

template<concepts::OrderedField T, typename V, typename M>
SystemRootFindingResult<V, T> levenberg_marquardt(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: $(J^\top J + \lambda I) \Delta\mathbf{x} = -J^\top F$ を解く。 ダンピングパラメータ $\lambda$ を成功時減らし、 失敗時増やす (Marquardt-Levenberg 更新規則)。 Newton と最急降下の補間。

用途: 初期値が真の根から遠い / ヤコビアンが特異近傍を通る / 最小二乗の感覚で根を探したい場合。 単なる方程式系解法というより、 残差最小化の混成手法。

powell_hybrid

template<concepts::OrderedField T, typename V, typename M>
SystemRootFindingResult<V, T> powell_hybrid(
    const std::function<V(const V&)>& F,
    const std::function<M(const V&)>& J,
    const V& x0,
    const ConvergenceCriteria<T>& criteria = ConvergenceCriteria<T>());

動作: Powell (1970) の dogleg-trust-region。 信頼領域内で Newton 方向と最急降下方向を組合せた dogleg 経路をたどる。 MINPACK hybrd ルーチン (SciPy fsolve の標準) と同系統。

用途: N-D 求根のデフォルト推奨。 局所収束性 (Newton 並み) と大域収束性 (信頼領域) を兼ね備える。

使用例

#include <math/roots/root_finding.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // 1D: f(x) = x^2 - 2 の根を Brent 法で求める
    auto result = brent_method<double>(
        [](double x) { return x * x - 2.0; },
        1.0, 2.0);
    if (result.converged && result.root)
        std::cout << "sqrt(2) = " << *result.root << std::endl;
    // 出力: sqrt(2) = 1.41421...

    // 多項式: x^3 - 6x^2 + 11x - 6 = (x-1)(x-2)(x-3) の全根
    Polynomial<double> p({-6.0, 11.0, -6.0, 1.0});  // 昇べき
    auto roots = solvePolynomial(p);
    for (auto& r : roots)
        std::cout << r.re << std::endl;
    // 出力: 1, 2, 3 (順序は実装依存)

    // 5次以上: solvePolynomial が自動的に jenkinsTraub を呼ぶ
    Polynomial<double> q({1.0, 0.0, 0.0, 0.0, 0.0, 1.0});  // x^5 + 1
    auto roots5 = solvePolynomial(q);
    for (auto& r : roots5)
        std::cout << r.re << " + " << r.im << "i" << std::endl;
}

関連する数学的背景

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

• 二分法 — 最も基本的な求根アルゴリズム
• ニュートン法 — 二次収束する反復法
• 割線法 — 導関数不要の超一次収束法
• 高度な多項式求根法 — Aberth-Ehrlich・companion matrix 法
• N 次根の精度倍増 Newton 法 — 多倍長 N 乗根
• Zimmermann 再帰平方根 — 高速多倍長平方根