LinAlg — 線形代数

概要

行列分解、固有値問題、線形方程式ソルバーを提供する。

• 分解 — LU, Cholesky, QR, SVD, LDL
• 固有値 — 対称/一般/一般化/エルミート、べき乗法
• ソルバー — 直接法 (LU/Cholesky/QR/SVD/LDL)、反復法 (CG, Jacobi, Gauss-Seidel, SOR)
• 最小二乗 — 通常/重み付き/正則化 (Ridge)
• 特殊行列 — 三重対角、帯行列、疎行列 (CSR)
• ベクトル空間 — Gram-Schmidt、基底完成、線形回帰

全関数は calx::algorithms 名前空間内のテンプレート関数 （ベクトル空間ユーティリティのみ calx）。

Complex モジュールについて: Complex モジュールの API ページは準備中だが、 Complex<T> クラス自体は利用可能であり、Matrix<Complex<double>> として 線形代数の各関数に渡すことができる。 std::complex<double> も同様に Matrix<std::complex<double>> として使用できる。

式テンプレート: 行列の要素単位演算（加減算、スカラー乗除）は式テンプレートで遅延評価される。 複数の演算を組み合わせても中間バッファは生成されず、代入時にまとめて評価される。 ただし、行列乗算 (operator*) は遅延評価されない（即座に計算される）。

次元チェック: 動的 Matrix<T> の次元不一致は実行時に DimensionError 例外で報告される。 Debug ビルド (NDEBUG 未定義) では例外に加えて assert が発動するため、デバッガで即座に停止できる。 StaticMatrix<T, Rows, Cols> では次元がテンプレート引数に含まれるため、一部の不一致はコンパイル時に検出される。 LinAlg 関数の StaticMatrix 用オーバーロード（static_assert による完全なコンパイル時チェック）は将来対応予定。

ヘッダ

#include <math/linalg/decomposition.hpp>          // 分解 + ユーティリティ
#include <math/linalg/eigenvalues.hpp>            // 固有値
#include <math/linalg/solvers.hpp>                // ソルバー
#include <math/linalg/vector_space_utilities.hpp> // ベクトル空間

ヘッダオンリー（リンク不要）。

cl /std:c++latest /EHsc /O2 /utf-8 /I<calx>/include your_code.cpp

Intel MKL の利用

Intel oneAPI MKL がインストールされている場合、LAPACK ルーチンを利用してパフォーマンスを向上させることができる。

• CMake オプション: -DCALX_HAS_MKL=ON
• 必要環境: Intel oneAPI MKL (oneAPI Base Toolkit に含まれる)
• リンクライブラリ: mkl_intel_lp64, mkl_sequential, mkl_core
• 有効になる機能: eigen, eigen_symmetric, eigen_hermitian, eigen_generalized 等の固有値ソルバーが LAPACK ルーチン (DSYEV, ZHEEV, DGEEV 等) を使用

cmake -B build -DCALX_HAS_MKL=ON -DCMAKE_BUILD_TYPE=Release
cmake --build build --config Release

MKL がない場合でも、calx 自前の実装 (Jacobi 法、QR 反復等) にフォールバックするため、全機能が利用可能である。

LU 分解

部分ピボット付き LU 分解。$PA = LU$。

• 計算量: $O(n^3)$
• 用途: 連立方程式、行列式、逆行列の汎用的な計算
• 制約: 正方行列

lu_decomposition

template<typename T>
std::pair<Matrix<T>, std::vector<size_type>>
    lu_decomposition(const Matrix<T>& a,
                     PivotStrategy pivot = PivotStrategy::Partial)

返り値: std::pair — 第 1 要素は L と U を 1 つの行列に格納した LU 行列 （対角より下が L、対角以上が U）、第 2 要素は行交換の記録（ピボット配列）。

例外: DimensionError — 行列が正方でない場合。

using namespace calx::algorithms;

Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
auto [LU, perm] = lu_decomposition(A);
// LU: L と U が 1 つの行列に格納されている
// perm: 行の交換順序

lu_solve

template<typename T>
Vector<T> lu_solve(const Matrix<T>& a, const Vector<T>& b,
                   PivotStrategy pivot = PivotStrategy::Partial)

返り値: 解ベクトル $\mathbf{x}$（$A\mathbf{x} = \mathbf{b}$ を満たす）。

例外: DimensionError — 次元不一致。

Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
Vector<double> b = {8, -11, -3};
auto x = lu_solve(A, b);
// x = {2, 3, -1}

lu_inverse

template<typename T>
Matrix<T> lu_inverse(const Matrix<T>& a)

返り値: 逆行列 $A^{-1}$。

Matrix<double> A = {{2, 1}, {5, 3}};
auto Ainv = lu_inverse(A);
// Ainv = {{3, -1}, {-5, 2}}
// A * Ainv = I

lu_determinant

template<typename T>
T lu_determinant(const Matrix<T>& a)

返り値: 行列式 $\det(A)$。LU 分解の対角要素の積にピボット交換の符号を乗じて計算。

Matrix<double> A = {{2, 1}, {5, 3}};
auto det = lu_determinant(A);
// det = 1.0  (2*3 - 1*5)

bareiss_determinant — Bareiss アルゴリズム

整数型の行列式を除算なしで計算する。 Sylvester の恒等式により、消去過程の除算が常に割り切れることを利用し、 中間結果も整数のまま $O(n^3)$ で行列式を得る。

• 計算量: $O(n^3)$
• アルゴリズム: 各ステップで $M_{ij}^{(k)} = \frac{M_{ij}^{(k-1)} \cdot M_{kk}^{(k-1)} - M_{ik}^{(k-1)} \cdot M_{kj}^{(k-1)}}{M_{k-1,k-1}^{(k-2)}}$ を計算。 この除算は Sylvester の恒等式により常に割り切れる
• 用途: Matrix<Int> や Matrix<int> など、除算が定義されない・または避けたい整数型
• 自動選択: Matrix::determinant() メンバ関数は、整数型 (numeric_traits<T>::is_integer) の場合に自動的に Bareiss を使用する

template<typename T>
T bareiss_determinant(const Matrix<T>& a)

返り値: 行列式 $\det(A)$。中間結果は整数のまま保持される。

// Int (多倍長整数) の行列式
Matrix<Int> A({{Int(100), Int(200), Int(300)},
               {Int(400), Int(500), Int(600)},
               {Int(700), Int(800), Int(1000)}});
auto det = bareiss_determinant(A);  // -3000000 (厳密)

// int でも使える
Matrix<int> B({{0, 1, 2}, {1, 0, 3}, {4, 5, 6}});
auto det2 = bareiss_determinant(B);  // 16

// determinant() メンバ関数は整数型で自動的に Bareiss を選択
auto det3 = A.determinant();  // bareiss_determinant と同じ結果

Cholesky 分解

正定値対称行列 $A = LL^\top$。LU 分解の約 2 倍高速。

• 計算量: $O(n^3/3)$（LU の半分）
• 用途: 最小二乗、カルマンフィルタ、ガウス過程、正規方程式
• 制約: 正定値対称行列のみ

cholesky_decomposition

template<typename T>
Matrix<T> cholesky_decomposition(const Matrix<T>& a)

返り値: 下三角行列 $L$（$A = LL^\top$ を満たす）。

例外: 正定値でない場合、分解に失敗する。

Matrix<double> A = {{4, 2}, {2, 3}};
auto L = cholesky_decomposition(A);
// L = {{2, 0}, {1, sqrt(2)}}
// L * L^T = A

cholesky_solve

template<typename T>
Vector<T> cholesky_solve(const Matrix<T>& a, const Vector<T>& b)

返り値: 解ベクトル $\mathbf{x}$。内部で Cholesky 分解 → 前進代入 → 後退代入を実行。

Matrix<double> A = {{4, 2}, {2, 3}};
Vector<double> b = {1, 2};
auto x = cholesky_solve(A, b);
// x ≈ {-0.125, 0.75}

QR 分解

Householder 反射を用いた QR 分解。$A = QR$。

• 計算量: $O(mn^2)$（$m \times n$ 行列）
• アルゴリズム: Householder 反射（Givens 回転より効率的）
• 用途: 最小二乗問題、固有値計算の前処理、数値的安定な連立方程式
• 特徴: LU より数値的に安定。長方形行列にも対応

qr_decomposition

template<typename T>
std::pair<Matrix<T>, Matrix<T>> qr_decomposition(const Matrix<T>& a)

返り値: std::pair — 第 1 要素は直交行列 $Q$（$m \times m$）、第 2 要素は上三角行列 $R$（$m \times n$）。

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
auto [Q, R] = qr_decomposition(A);
// Q: 3x3 直交行列 (Q^T Q = I)
// R: 3x2 上三角行列
// Q * R = A

qr_solve

template<typename T>
Vector<T> qr_solve(const Matrix<T>& a, const Vector<T>& b)

返り値: 解ベクトル $\mathbf{x}$。過決定系（$m > n$）では最小二乗解 $\min\|A\mathbf{x} - \mathbf{b}\|_2$ を返す。

// 過決定系 (3x2): 最小二乗解
Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
Vector<double> b = {1, 2, 2};
auto x = qr_solve(A, b);
// x ≈ {0.667, 0.5}

SVD（特異値分解）

Golub-Kahan アルゴリズムによる特異値分解。$A = U\Sigma V^\top$。

• 計算量: $O(mn \cdot \min(m,n))$
• アルゴリズム: Householder 二重対角化 → 陰的シフト QR 反復 （参考: Golub & Van Loan "Matrix Computations" §8.6）
• 用途: ランク決定、条件数計算、最小二乗、低ランク近似、主成分分析
• 特徴: 最も安定な分解。悪条件行列にも対応

svd_decomposition

template<typename T>
std::tuple<Matrix<T>, Vector<T>, Matrix<T>>
    svd_decomposition(const Matrix<T>& A, bool transV = true)

返り値: std::tuple<U, σ, V> — $U$ は左特異ベクトル行列（$m \times k$）、$\sigma$ は特異値ベクトル（$k$ 要素、降順）、 $V^\top$ は右特異ベクトル行列（$k \times n$）。$k = \min(m, n)$。

Matrix<double> A = {{1, 2}, {3, 4}, {5, 6}};
auto [U, sigma, Vt] = svd_decomposition(A);

// 再構成: A = U * diag(sigma) * Vt
auto A_rec = reconstruct_matrix_from_svd(U, sigma, Vt);

svd_solve

template<typename T>
Vector<T> svd_solve(const Matrix<T>& a, const Vector<T>& b,
    T rcond = std::numeric_limits<T>::epsilon() * T(100))

返り値: 解ベクトル $\mathbf{x}$。悪条件行列でもダンピングにより数値的に安定な解を返す。

使い分け: 条件数が良好なら lu_solve が最速。悪条件や特異に近い行列では svd_solve を使用。rcond を大きくするほどダンピングが強くなり安定性が増すが、精度は低下する。

Matrix<double> A = {{1, 2}, {3, 4}, {5, 6}};
Vector<double> b = {1, 2, 3};
auto x = svd_solve(A, b);
// x ≈ {0, 0.5}

reconstruct_matrix_from_svd

template<typename T>
Matrix<T> reconstruct_matrix_from_svd(
    const Matrix<T>& U, const Vector<T>& s, const Matrix<T>& Vt)

返り値: 再構成された行列 $A = U \cdot \mathrm{diag}(s) \cdot V^\top$（$m \times n$）。

Matrix<double> A = {{3, 0}, {0, 5}};
auto [U, s, Vt] = svd_decomposition(A);
auto A2 = reconstruct_matrix_from_svd(U, s, Vt);
// A2 ≈ A (再構成誤差 ≈ 0)

svd_determinant

template<typename T>
T svd_determinant(const Matrix<T>& a)

SVD を用いて行列式を計算する。特異値の積 $\prod \sigma_i$ に符号を乗じて返す。 悪条件行列で lu_determinant よりも数値的に安定。

Matrix<double> A = {{1, 2}, {3, 4}};
auto det = svd_determinant(A);
// det ≈ -2.0

LDL 分解

対称行列 $A = LDL^\top$。不定値行列にも対応。

• 計算量: $O(n^3/3)$
• 用途: 対称だが正定値とは限らない行列。最適化問題のヘッシアンなど
• 特徴: Cholesky（正定値のみ）の一般化

ldl_decomposition

template<typename T>
std::pair<Matrix<T>, Vector<T>> ldl_decomposition(const Matrix<T>& a)

返り値: std::pair — 第 1 要素は下三角行列 $L$（対角は 1）、第 2 要素は対角要素ベクトル $D$。

Matrix<double> A = {{4, 2, -1}, {2, 5, 3}, {-1, 3, 6}};
auto [L, D] = ldl_decomposition(A);
// A = L * diag(D) * L^T

ldl_solve

template<typename T>
Vector<T> ldl_solve(const Matrix<T>& a, const Vector<T>& b)

返り値: 解ベクトル $\mathbf{x}$。

Matrix<double> A = {{4, 2}, {2, 3}};
Vector<double> b = {1, 2};
auto x = ldl_solve(A, b);
// x ≈ {-0.125, 0.75}

分解の使い分け

ガウス消去法

拡大係数行列 $[A|\mathbf{b}]$ に対して前進消去と後退代入を行い、連立方程式 $A\mathbf{x} = \mathbf{b}$ を直接解く。 LU 分解と数学的に等価だが、拡大係数行列上で直接操作するため、分解結果を再利用しない場合に簡潔に使える。

gaussian_elimination

template<typename T>
Vector<T> gaussian_elimination(const Matrix<T>& a, const Vector<T>& b,
                               PivotStrategy pivot = PivotStrategy::Partial)

返り値: 解ベクトル $\mathbf{x}$。

例外:

• DimensionError — 行列が正方でない、または次元不一致
• MathError — 特異行列が検出された場合

Rational 型での特殊動作: T = Rational の場合、ピボット選択は絶対値最大ではなく 高さ最小（分子と分母の最大値が最も小さい要素）のピボットを選ぶ。 これにより Rational 演算での分母子の膨張を抑え、計算効率を維持する。

using namespace calx::algorithms;

// double で部分ピボット（既定）
Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
Vector<double> b = {8, -11, -3};
auto x = gaussian_elimination(A, b);
// x = {2, 3, -1}

using namespace calx::algorithms;

// Rational で厳密解
Matrix<Rational> A = {
    {Rational(2), Rational(1), Rational(-1)},
    {Rational(-3), Rational(-1), Rational(2)},
    {Rational(-2), Rational(1), Rational(2)}
};
Vector<Rational> b = {Rational(8), Rational(-11), Rational(-3)};
auto x = gaussian_elimination(A, b);
// x = {2, 3, -1} (厳密な有理数解)

ピボット戦略

PivotStrategy は LU 分解とガウス消去法の両方で使用される列挙型。

enum class PivotStrategy {
    None,       // ピボット選択なし (対角要素をそのまま使用)
    Partial,    // 部分ピボット選択 (列内で絶対値最大, Rational は高さ最小)
    Full        // 完全ピボット選択 (残り部分行列全体から選択)
};

Rational 型の場合: Partial と Full の選択基準が浮動小数点型と異なる。 絶対値最大ではなく、高さ（height）最小のピボットを選択する。 高さとは $\max(|\text{分子}|, |\text{分母}|)$ であり、小さいほど後続の演算で分母子の膨張が抑えられる。

ピボット戦略による結果の違い

悪条件行列では、ピボットなしだと数値精度が大幅に低下する。以下は (1,1) 成分が極めて小さい場合の例。

using namespace calx::algorithms;

// 悪条件行列: (1,1) 成分が極めて小さい
Matrix<float> A({{1e-8f, 1.0f}, {1.0f, 1.0f}});
Vector<float> b({1.0f, 2.0f});
// 真の解: x = {1.00000001..., 0.99999999...} ≈ {1.0, 1.0}

// ピボットなし → 1e-8 で割るため丸め誤差が爆発
auto x_none = gaussian_elimination(A, b, PivotStrategy::None);
// x_none = {0.000000, 1.000000}  ← x_none[0] が完全に壊れている

// 部分ピボット → 行交換で大きい成分を先に処理
auto x_partial = gaussian_elimination(A, b, PivotStrategy::Partial);
// x_partial = {1.000000, 1.000000}  ← 正確

ピボットなしでは第 1 行を $10^{-8}$ で割ることで $10^8$ 倍のスケーリングが発生し、 float の仮数部 23 bit (約 7 桁) では第 2 行の消去時に $1.0$ が $10^8$ に対して丸め落とされる。 部分ピボットでは行を交換して $|a_{11}| = 1$ とするため、この桁落ちを回避できる。

固有値

eigen — 一般固有値分解

一般（非対称を含む）行列の全固有値・固有ベクトルを計算。

• アルゴリズム: 内部で対称性を自動検出し、対称なら eigen_symmetric に委譲。 非対称の場合は Hessenberg 縮約 + Givens 回転ベースの QR 反復。 MKL が利用可能な場合は LAPACK ルーチンを使用
• 計算量: $O(n^3)$

template<typename T>
std::pair<Vector<std::complex<T>>, Matrix<std::complex<T>>>
    eigen(const Matrix<T>& a)

返り値: std::pair — 第 1 要素は固有値ベクトル（Vector<complex<T>>）、 第 2 要素は固有ベクトル行列（Matrix<complex<T>>、列ベクトルが各固有ベクトル）。

Matrix<double> A = {{0, -1}, {1, 0}};  // 回転行列
auto [eigenvalues, eigenvectors] = eigen(A);
// eigenvalues = {i, -i} (純虚数)

eigen_symmetric — 対称行列の固有値

実対称行列の全固有値・固有ベクトルを計算。固有値は実数。

• アルゴリズム: Jacobi 法（回転による対角化）。MKL 利用可能時は LAPACK
• 計算量: $O(n^3)$
• 特徴: 全固有値が実数保証。固有ベクトルは直交

template<typename T>
std::pair<Vector<T>, Matrix<T>>
    eigen_symmetric(const Matrix<T>& a)

Matrix<double> A = {{2, 1}, {1, 2}};
auto [vals, vecs] = eigen_symmetric(A);
// vals = {1, 3}
// vecs の各列は対応する固有ベクトル

eigen_hermitian — エルミート行列の固有値

エルミート行列（$A^H = A$）の全固有値・固有ベクトルを計算。固有値は実数。

• アルゴリズム: MKL 利用可能時は LAPACKE_zheev / LAPACKE_cheev。 MKL がない場合は、$n \times n$ エルミート行列を $2n \times 2n$ の実対称行列に変換し、 eigen_symmetric で解く代替実装を使用する
• MKL なしでも利用可能: 代替実装により全環境で動作する

template<typename T>
std::pair<Vector<T>, Matrix<std::complex<T>>>
    eigen_hermitian(const Matrix<std::complex<T>>& a)

// Hermitian 行列の固有値は実数
Matrix<std::complex<double>> H = {{{2,0}, {1,-1}}, {{1,1}, {3,0}}};
auto [vals, vecs] = eigen_hermitian(H);
// vals ≈ {1.268, 3.732} (実数)

eigen_generalized — 一般化固有値問題

$A\mathbf{x} = \lambda B\mathbf{x}$ を解く。$A$, $B$ は一般の正方行列。

• 用途: 振動解析、座屈問題、一般化主成分分析

template<typename T>
std::pair<Vector<std::complex<T>>, Matrix<std::complex<T>>>
    eigen_generalized(const Matrix<T>& a, const Matrix<T>& b)

Matrix<double> A = {{1, 2}, {3, 4}};
Matrix<double> B = {{5, 6}, {7, 8}};
auto [vals, vecs] = eigen_generalized(A, B);
// A * vecs = B * vecs * diag(vals)

eigen_generalized_symmetric — 対称一般化固有値問題

$A\mathbf{x} = \lambda B\mathbf{x}$ を解く。$A$ は対称、$B$ は正定値対称。

• アルゴリズム: Cholesky 分解で $B = LL^\top$ とし、$C = L^{-1}AL^{-\top}$ の対称固有値問題に変換

template<typename T>
std::pair<Vector<T>, Matrix<T>>
    eigen_generalized_symmetric(const Matrix<T>& a, const Matrix<T>& b)

Matrix<double> A = {{4, 1}, {1, 3}};
Matrix<double> B = {{2, 0}, {0, 1}};
auto [vals, vecs] = eigen_generalized_symmetric(A, B);
// A * v = λ * B * v

eigen_complex — 複素行列の固有値

template<typename T>
std::pair<Vector<Complex<T>>, Matrix<Complex<T>>>
    eigen_complex(const Matrix<T>& a)

一般の実行列の固有値・固有ベクトルを複素数として返す。 eigen と同等だが、固有ベクトルも Complex<T> 型で返す。

Matrix<double> A = {{0, -1}, {1, 0}};  // 90° 回転
auto [vals, vecs] = eigen_complex(A);
// vals = {i, -i}

固有値問題の使い分け

べき乗法

特定の固有値を反復的に計算する手法。大行列で全固有値が不要な場合に有効。

power_method

最大固有値（絶対値最大）とその固有ベクトルを計算。

• 収束速度: $O(|\lambda_1/\lambda_2|^k)$ — 最大と 2 番目の固有値の比に依存
• 計算量: 反復あたり $O(n^2)$（行列-ベクトル積）

template<typename T>
std::pair<T, Vector<T>> power_method(const Matrix<T>& a,
    T tolerance = T(1e-10),
    std::size_t max_iterations = 1000)

Matrix<double> A = {{2, 1}, {1, 3}};
auto [lambda, v] = power_method(A);
// lambda ~ 3.618

inverse_power_method

指定したシフト $\sigma$ に最も近い固有値を計算。

template<typename T>
std::pair<T, Vector<T>> inverse_power_method(const Matrix<T>& a,
    T sigma = T(0),
    T tolerance = T(1e-10),
    std::size_t max_iterations = 1000)

Matrix<double> A = {{2, 1}, {1, 3}};
auto [lambda, v] = inverse_power_method(A);
// lambda ≈ 1.382 (最小固有値)

rayleigh_quotient_iteration

Rayleigh 商反復法。シフトを反復ごとに更新し、3 次収束を達成。

template<typename T>
std::pair<T, Vector<T>> rayleigh_quotient_iteration(const Matrix<T>& a,
    const Vector<T>& initial_guess,
    T tolerance = T(1e-10),
    std::size_t max_iterations = 100)

Matrix<double> A = {{2, 1}, {1, 3}};
Vector<double> v0 = {1, 1};
auto [lambda, v] = rayleigh_quotient_iteration(A, v0);
// lambda ≈ 3.618 (v0 に近い固有値に高速収束)

schur_decomposition — Schur 分解

template<typename T>
std::pair<Matrix<T>, Matrix<T>> schur_decomposition(const Matrix<T>& a)

実 Schur 分解 $A = QTQ^\top$ を計算する。$T$ は準上三角行列（実固有値は対角上、複素固有値は $2 \times 2$ ブロック）、 $Q$ は直交行列。行列関数 (expm, logm, sqrtm) の内部で使用される。

Matrix<double> A = {{1, 2}, {0, 3}};
auto [T, Q] = schur_decomposition(A);
// T ≈ {{1, *}, {0, 3}}  (上三角)
// Q * T * Q^T ≈ A

ソルバー

SolverType

enum class SolverType {
    Auto,       // 行列の性質から最適なソルバーを自動選択
    LU,         // LU 分解（汎用）
    Cholesky,   // Cholesky 分解（正定値対称）
    QR,         // QR 分解（数値安定性重視、過決定系対応）
    SVD,        // SVD（悪条件行列、ランク落ち対応）
    LDL,        // LDL 分解（対称不定値）
    Iterative   // 反復法（CG）
};

Auto の選択ロジック:

1. 対称かつ正定値 → Cholesky
2. 対称 → LDL
3. 正方 → LU
4. 過決定系 → QR

solve

template<typename T>
Vector<T> solve(const Matrix<T>& a, const Vector<T>& b,
    SolverType solver_type = SolverType::Auto)

返り値: 解ベクトル $\mathbf{x}$。

using namespace calx::algorithms;

Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
Vector<double> b = {8, -11, -3};

auto x1 = solve(A, b);                      // 自動選択
auto x2 = solve(A, b, SolverType::LU);      // LU 指定
auto x3 = solve(A, b, "cholesky");           // 文字列でも指定可能

文字列版オーバーロードも提供。使用可能な文字列: "auto", "lu", "cholesky", "qr", "svd", "ldl", "iterative"。

solve_multi — 複数右辺

template<typename T>
Matrix<T> solve_multi(const Matrix<T>& a, const Matrix<T>& b,
    SolverType solver_type = SolverType::Auto)

返り値: 解行列 $X$（$n \times m$）。$AX = B$ の各列を個別に solve で解く。

Matrix<double> A = {{2, 1}, {5, 3}};
Matrix<double> B = {{1, 0}, {0, 1}};  // 単位行列 → 逆行列を求める
auto X = solve_multi(A, B);
// X = {{3, -1}, {-5, 2}} = A^{-1}

反復法ソルバー

大規模行列で直接法が重い場合、または疎行列で有効な反復解法。

conjugate_gradient — 共役勾配法 (CG)

• 収束速度: $O(\sqrt{\kappa})$ 回の反復（$\kappa$: 条件数）
• 制約: 正定値対称行列のみ
• 計算量: 反復あたり $O(n^2)$（密行列）または $O(\mathrm{nnz})$（疎行列）

template<typename T>
Vector<T> conjugate_gradient(const Matrix<T>& a, const Vector<T>& b,
    T tolerance = std::numeric_limits<T>::epsilon() * 10,
    typename Matrix<T>::size_type max_iterations = 1000)

Matrix<double> A = {{4, 1}, {1, 3}};
Vector<double> b = {1, 2};
auto x = conjugate_gradient(A, b, 1e-10, 1000);

jacobi_method — Jacobi 反復法

• 収束条件: 対角優位（$|a_{ii}| > \sum_{j \ne i} |a_{ij}|$）で保証
• 収束速度: CG より遅いが実装が単純。並列化に適する

template<typename T>
Vector<T> jacobi_method(const Matrix<T>& a, const Vector<T>& b,
    T tolerance = std::numeric_limits<T>::epsilon() * 10,
    typename Matrix<T>::size_type max_iterations = 1000)

Matrix<double> A = {{10, -1}, {-1, 10}};
Vector<double> b = {9, 9};
auto x = jacobi_method(A, b);
// x ≈ {1, 1}

gauss_seidel_method — Gauss-Seidel 反復法

• 収束速度: Jacobi の約 2 倍高速（同じ条件下）
• 収束条件: 対角優位 or 正定値対称で保証

template<typename T>
Vector<T> gauss_seidel_method(const Matrix<T>& a, const Vector<T>& b,
    T tolerance = std::numeric_limits<T>::epsilon() * 10,
    typename Matrix<T>::size_type max_iterations = 1000)

Matrix<double> A = {{10, -1}, {-1, 10}};
Vector<double> b = {9, 9};
auto x = gauss_seidel_method(A, b);
// x ≈ {1, 1} (Jacobi より速く収束)

sor_method — SOR 法（逐次過緩和法）

• 収束条件: $0 < \omega < 2$ で対角優位なら収束保証。$\omega = 1$ で Gauss-Seidel に一致

template<typename T>
Vector<T> sor_method(const Matrix<T>& a, const Vector<T>& b,
    T omega = static_cast<T>(1.5),
    T tolerance = std::numeric_limits<T>::epsilon() * 10,
    typename Matrix<T>::size_type max_iterations = 1000)

Matrix<double> A = {{10, -1}, {-1, 10}};
Vector<double> b = {9, 9};
auto x = sor_method(A, b, 1.2);  // omega = 1.2
// x ≈ {1, 1}

反復法の使い分け

最小二乗法

least_squares

template<typename T>
Vector<T> least_squares(const Matrix<T>& a, const Vector<T>& b,
    SolverType solver_type = SolverType::Auto)

返り値: 最小二乗解 $\mathbf{x} = \arg\min \|A\mathbf{x} - \mathbf{b}\|_2$（$n$ 要素）。

Auto 選択ロジック: 過剰決定系（$m \ge n$）では QR 分解、劣決定系（$m < n$）では SVD を使用する。 正規方程式 $A^\top A \mathbf{x} = A^\top \mathbf{b}$ を直接解く方式は数値的に不安定なため採用していない。

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
Vector<double> b = {1, 2, 2};
auto x = least_squares(A, b);
// x ≈ {0.667, 0.5}  (最小二乗フィット y = 0.667 + 0.5*x)

weighted_least_squares

template<typename T>
Vector<T> weighted_least_squares(const Matrix<T>& a, const Vector<T>& b,
    const Vector<T>& weights)

返り値: 重み付き最小二乗解 $\mathbf{x} = \arg\min \sum w_i (A_i\mathbf{x} - b_i)^2$。 内部的には $\sqrt{w_i}$ を $A$ の各行と $b$ の各要素に乗じて通常の最小二乗に変換。

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
Vector<double> b = {1, 2, 2};
Vector<double> w = {1, 1, 2};  // 3番目の観測を重視
auto x = weighted_least_squares(A, b, w);

regularized_least_squares — Ridge 回帰

template<typename T>
Vector<T> regularized_least_squares(const Matrix<T>& a, const Vector<T>& b,
    T lambda = static_cast<T>(0.1))

返り値: 正則化最小二乗解 $\mathbf{x} = (A^\top A + \lambda I)^{-1} A^\top \mathbf{b}$。

$\lambda$ が大きいほど解のノルムが小さくなり安定するが、バイアスが増大する。 交差検証で最適な $\lambda$ を選択するのが一般的。

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
Vector<double> b = {1, 2, 2};
auto x = regularized_least_squares(A, b, 0.1);
// λ = 0.1: 解 x のノルムが抑制される

特殊行列ソルバー

tridiagonal_solve — Thomas アルゴリズム

三重対角行列の高速ソルバー。

• 計算量: $O(n)$（LU の $O(n^3)$ に対して劇的に高速）
• 用途: スプライン補間、有限差分法 (FDM)、陰的時間積分

template<typename T>
Vector<T> tridiagonal_solve(const Vector<T>& a, const Vector<T>& b,
    const Vector<T>& c, const Vector<T>& d)

//  [2 1 0] [x0]   [1]
//  [1 3 1] [x1] = [2]
//  [0 1 2] [x2]   [3]
Vector<double> a = {1, 1};      // 下対角
Vector<double> b = {2, 3, 2};   // 主対角
Vector<double> c = {1, 1};      // 上対角
Vector<double> d = {1, 2, 3};   // 右辺
auto x = tridiagonal_solve(a, b, c, d);

band_solve — 帯行列ソルバー

部分ピボット付き帯行列 LU ソルバー。

• 計算量: $O(n \cdot (k_l + k_u)^2)$
• 最適化: $k_l = k_u = 1$（三重対角）の場合は Thomas アルゴリズムに自動委譲

template<typename T>
Vector<T> band_solve(const Matrix<T>& a, const Vector<T>& b,
    typename Matrix<T>::size_type kl,
    typename Matrix<T>::size_type ku)

// 三重対角行列 (帯幅 1)
Matrix<double> A = {{2,-1, 0}, {-1, 2,-1}, {0,-1, 2}};
Vector<double> b = {1, 0, 1};
auto x = band_solve(A, b, 1, 1);  // lower=1, upper=1
// x = {1, 1, 1}

疎行列ソルバー（CSR 形式）

CSR (Compressed Sparse Row) 形式で格納された疎行列の解法。

sparse_solve

template<typename T>
Vector<T> sparse_solve(
    const std::vector<T>& values,
    const std::vector<size_type>& row_ptr,
    const std::vector<size_type>& col_idx,
    const Vector<T>& b,
    T tolerance = std::numeric_limits<T>::epsilon() * T(100),
    typename Matrix<T>::size_type max_iterations = 0)

内部で sparse_bicgstab に委譲（非対称行列にも対応）。

SparseMatrix<double> A(3, 3);
A.insert(0,0, 4); A.insert(0,1,-1);
A.insert(1,0,-1); A.insert(1,1, 4); A.insert(1,2,-1);
A.insert(2,1,-1); A.insert(2,2, 4);
Vector<double> b = {1, 2, 1};
auto x = sparse_solve(A, b);

sparse_cg — 疎行列 CG

template<typename T>
Vector<T> sparse_cg(
    const std::vector<T>& values,
    const std::vector<size_type>& row_ptr,
    const std::vector<size_type>& col_idx,
    const Vector<T>& b,
    T tolerance = std::numeric_limits<T>::epsilon() * 10,
    typename Matrix<T>::size_type max_iterations = 0)

制約: 対称正定値行列のみ。非対称なら sparse_bicgstab を使用。

// SPD 疎行列
SparseMatrix<double> A(2, 2);
A.insert(0,0, 4); A.insert(0,1, 1);
A.insert(1,0, 1); A.insert(1,1, 3);
Vector<double> b = {1, 2};
auto x = sparse_cg(A, b);
// x ≈ {0.091, 0.636}

sparse_bicgstab — 疎行列 BiCGSTAB

template<typename T>
Vector<T> sparse_bicgstab(
    const std::vector<T>& values,
    const std::vector<size_type>& row_ptr,
    const std::vector<size_type>& col_idx,
    const Vector<T>& b,
    T tolerance = std::numeric_limits<T>::epsilon() * 10,
    typename Matrix<T>::size_type max_iterations = 0)

CG と異なり非対称行列にも対応。GMRES よりメモリ効率が良い。

SparseMatrix<double> A(2, 2);
A.insert(0,0, 3); A.insert(0,1, 1);
A.insert(1,0, 2); A.insert(1,1, 4);
Vector<double> b = {1, 2};
auto x = sparse_bicgstab(A, b);
// x ≈ {0.2, 0.4}

行列関数

行列に対する解析関数。Schur 分解と Parlett 再帰を内部で使用。

expm — 行列指数関数

template<typename T>
Matrix<T> expm(const Matrix<T>& A)

$e^A$ を Scaling-and-squaring 法 + Padé 近似で計算する (Higham 2005)。 行列の 1-ノルムに応じて Padé 次数 [3/3]〜[13/13] を自動選択する。 ノルムが小さい行列では低次近似で十分なため高速に、 大ノルムの行列ではスケーリング (二分割) と最大次数 [13/13] を組み合わせて高精度を確保する。 次数を手動で指定するパラメータは現在提供していない。

Matrix<double> A = {{0, -1}, {1, 0}};
auto E = expm(A);
// E ≈ {{cos(1), -sin(1)}, {sin(1), cos(1)}}

sqrtm — 行列平方根

template<typename T>
Matrix<T> sqrtm(const Matrix<T>& A)

主値平方根 $A^{1/2}$ を計算する。$\text{sqrtm}(A) \cdot \text{sqrtm}(A) = A$。

Matrix<double> A = {{4, 1}, {1, 3}};
auto S = sqrtm(A);
// S * S ≈ A

logm — 行列対数

template<typename T>
Matrix<T> logm(const Matrix<T>& A)

主値対数 $\log(A)$ を計算する。$e^{\log(A)} = A$。

Matrix<double> A = {{2, 1}, {0, 2}};
auto L = logm(A);
auto check = expm(L);
// check ≈ A

Gershgorin 円盤

固有値の存在範囲を推定する。各行について中心 $a_{ii}$ と半径 $\sum_{j \ne i}|a_{ij}|$ の円盤を返す。

gershgorin_discs

template<typename T>
std::vector<GershgorinDisc<T>> gershgorin_discs(const Matrix<T>& a)

// GershgorinDisc<T> { T center; T radius; }

Matrix<double> A = {{5, 0.5, 0.2}, {0.5, 3, 0.3}, {0.2, 0.3, 4}};
auto discs = gershgorin_discs(A);
// discs[0]: center=5.0, radius=0.7
// discs[1]: center=3.0, radius=0.8
// discs[2]: center=4.0, radius=0.5

gershgorin_bounds

template<typename T>
std::pair<T, T> gershgorin_bounds(const Matrix<T>& a)

全固有値が含まれる区間 $[\lambda_{\min}, \lambda_{\max}]$ を返す。

auto [lo, hi] = gershgorin_bounds(A);
// 全固有値は [lo, hi] に含まれる

複素固有値を持つ行列への適用

gershgorin_bounds は半径の計算に std::abs() を使用するため、 複素数の要素を持つ行列でも正常に動作する。 ただし、返される bounds は実軸上の区間であり、固有値の虚部は反映されない。

複素固有値の位置を正確に把握したい場合は、gershgorin_discs で 個別の Gershgorin 円盤（中心: $a_{ii}$、半径: $R_i = \sum_{j \ne i}|a_{ij}|$）を確認すること。 各固有値は少なくとも 1 つの円盤内に存在する。

ベクトル空間ユーティリティ

ベクトル空間の基本操作を提供する。与えられたベクトル集合を直交基底に変換する Gram-Schmidt 正規直交化、 あるベクトルを別のベクトル方向に射影する projection、ベクトル集合の線形独立性の判定、 そして最小二乗法による線形回帰などが含まれる。

典型的な使用場面: 基底の構成・直交化 (固有ベクトルの後処理、部分空間射影)、 データの線形フィッティング、ベクトル間の角度計算など。

gram_schmidt

template<typename T, typename V>
std::vector<V> gram_schmidt(std::vector<V> vectors)

std::vector<Vector<double>> basis = {{1,0,1}, {0,1,1}, {1,1,0}};
auto ortho = gram_schmidt<double, Vector<double>>(basis);
// ortho は直交化されたベクトル列

projection

template<typename T, typename V>
V projection(const V& v, const V& onto)

Vector<double> v = {3, 4};
Vector<double> u = {1, 0};
auto proj = projection<double>(v, u);
// proj = {3, 0}

normalize

template<typename T, typename V>
V normalize(const V& v)

linear_combination

template<typename T, typename V>
V linear_combination(const std::vector<V>& vectors, const std::vector<T>& coefficients)

is_linearly_independent

template<typename T, typename V, typename M>
bool is_linearly_independent(const std::vector<V>& vectors)

linear_regression

template<typename T, typename V, typename M>
V linear_regression(const M& X, const V& y)

最小二乗法による線形回帰: $\beta = (X^\top X)^{-1} X^\top y$。

// y = a + b*x のフィッティング
Matrix<double> X = {{1,1}, {1,2}, {1,3}, {1,4}};
Vector<double> y = {2.1, 3.9, 6.2, 7.8};
auto beta = linear_regression<double, Vector<double>, Matrix<double>>(X, y);
// beta ≈ {0.05, 1.97} → y ≈ 0.05 + 1.97*x

angle

template<typename T, typename V>
T angle(const V& u, const V& v)

2 ベクトル間の角度 (ラジアン) を返す。

BLAS (Basic Linear Algebra Subprograms)

BLAS (Basic Linear Algebra Subprograms) は線形代数の基本演算を定義した標準インターフェースである。 calx では calx::blas 名前空間に BLAS Level 1〜3 の主要関数を実装している。

これらの関数は分解アルゴリズムやソルバーの内部で使用されるほか、 ユーザーが行列・ベクトル演算を明示的に制御したい場合（例: $y \leftarrow \alpha A x + \beta y$ のような 更新演算をバッファ割り当てなしで実行）にも直接利用できる。

Level 1: ベクトル-ベクトル演算

dot — 内積

template<typename T> T dot(const Vector<T>& x, const Vector<T>& y)

Vector<double> x = {1, 2, 3}, y = {4, 5, 6};
double d = blas::dot(x, y);  // 32.0

nrm2 — 2-ノルム

ベクトルのユークリッドノルム（2-ノルム）$\|x\|_2 = \sqrt{\sum x_i^2}$ を計算する。

template<typename T> T nrm2(const Vector<T>& x)

asum — 1-ノルム

ベクトルの絶対値和（1-ノルム）$\|x\|_1 = \sum |x_i|$ を計算する。

template<typename T> T asum(const Vector<T>& x)

scal — スケーリング

ベクトルをスカラー倍する。$x \leftarrow \alpha x$。in-place 演算。

template<typename T> void scal(T alpha, Vector<T>& x)  // x ← α·x

axpy — AXPY

$y \leftarrow \alpha x + y$。ベクトルのスケーリングと加算を1回の走査で実行する。BLAS で最も頻繁に使われる関数の一つ。

template<typename T> void axpy(T alpha, const Vector<T>& x, Vector<T>& y)  // y ← α·x + y

Vector<double> x = {1, 2, 3}, y = {4, 5, 6};
blas::axpy(2.5, x, y);  // y = {6.5, 10.0, 13.5}

Level 2: 行列-ベクトル演算

gemv — 一般行列-ベクトル積

一般行列とベクトルの積。$y \leftarrow \alpha \cdot \mathrm{op}(A) \cdot x + \beta \cdot y$。trans=true で転置。バッファ再利用により一時オブジェクト不要。

template<typename T>
void gemv(bool trans, T alpha, const Matrix<T>& A,
          const Vector<T>& x, T beta, Vector<T>& y)  // y ← α·op(A)·x + β·y

trsv — 三角行列ソルバー

三角連立方程式 $\mathrm{op}(A) x = b$ を前進/後退代入で解く。upper/trans/unit_diag で三角行列の種類を指定。

template<typename T>
void trsv(bool upper, bool trans, bool unit_diag,
          const Matrix<T>& A, Vector<T>& x)  // x ← op(A)⁻¹·x

Level 3: 行列-行列演算

gemm — 一般行列積

一般行列積。$C \leftarrow \alpha \cdot \mathrm{op}(A) \cdot \mathrm{op}(B) + \beta \cdot C$。Level 3 BLAS の中核。大規模行列ではメモリアクセスパターンが性能を支配する。

template<typename T>
void gemm(bool transA, bool transB, T alpha,
          const Matrix<T>& A, const Matrix<T>& B,
          T beta, Matrix<T>& C)  // C ← α·op(A)·op(B) + β·C

trsm — 三角行列ソルバー (行列版)

三角行列を係数とする行列方程式を解く。$B \leftarrow \alpha \cdot \mathrm{op}(A)^{-1} \cdot B$ または $B \cdot \mathrm{op}(A)^{-1}$。

template<typename T>
void trsm(bool side_left, bool upper, bool trans, bool unit_diag,
          T alpha, const Matrix<T>& A, Matrix<T>& B)  // B ← α·op(A)⁻¹·B or B·op(A)⁻¹

ランダム化 SVD

Halko-Martinsson-Tropp (2011) アルゴリズム。大規模行列の上位 $k$ 個の特異値・ベクトルを高速に計算する。

randomizedSVD

パラメータ: k=求める特異値の数、oversampling=オーバーサンプリング数（精度向上）、powerIter=パワー反復回数（特異値の減衰が遅い場合に有効）、seed=乱数シード。

template<typename T>
RandomizedSVDResult<T> randomizedSVD(const Matrix<T>& A, size_t k,
    size_t oversampling = 5, size_t powerIter = 1, unsigned seed = 42)

// RandomizedSVDResult<T> { Matrix<T> U; std::vector<T> S; Matrix<T> Vt; }

Matrix<double> A(1000, 500);  // 大規模行列
// ... fill A ...
auto result = randomizedSVD(A, 10);  // 上位 10 特異値
// result.U: (1000×10), result.S: 10個, result.Vt: (10×500)

疎行列直接法

疎行列の直接分解ソルバー。#include <math/linalg/simplicial_cholesky.hpp> 等。

SimplicialLLT — 疎行列 Cholesky

template<typename T, typename IndexType = std::size_t>
class SimplicialLLT {
public:
    explicit SimplicialLLT(const SparseMatrix<T, IndexType>& A);
    void compute(const SparseMatrix<T, IndexType>& A);
    Vector<T> solve(const Vector<T>& b) const;
    bool success() const;
};

対称正定値疎行列の Cholesky 分解 $A = LL^\top$。

SparseMatrix<double> A(100, 100);
// ... 対称正定値行列を構築 ...
SimplicialLLT<double> chol(A);
auto x = chol.solve(b);

SimplicialLDLT — 疎行列 LDL

template<typename T, typename IndexType = std::size_t>
class SimplicialLDLT {
public:
    explicit SimplicialLDLT(const SparseMatrix<T, IndexType>& A);
    void compute(const SparseMatrix<T, IndexType>& A);
    Vector<T> solve(const Vector<T>& b) const;
    bool success() const;
};

対称疎行列の $A = LDL^\top$ 分解。sqrt 不要で数値的に安定。

SparseLU — 疎行列 LU

template<typename T, typename IndexType = std::size_t>
class SparseLU {
public:
    void compute(const SparseMatrix<T, IndexType>& A);
    Vector<T> solve(const Vector<T>& b) const;
    bool success() const;
};

非対称疎行列の LU 分解 $PA = LU$。

SparseQR — 疎行列 QR

template<typename T, typename IndexType = std::size_t>
class SparseQR {
public:
    void compute(const SparseMatrix<T, IndexType>& A);
    Vector<T> solve(const Vector<T>& b) const;
    bool success() const;
};

疎行列の QR 分解。最小二乗問題にも使用可能。

疎行列反復法

名前空間 calx::sparse_algorithms。大規模疎行列の反復解法。

IterativeSolver — ビルダーパターン

template<typename T, typename IndexType = std::size_t>
class IterativeSolver {
public:
    IterativeSolver(const SparseMatrix<T, IndexType>& A, const Vector<T>& b);
    IterativeSolver& solver(SolverType type);          // CG, BiCGSTAB, GMRES, MINRES から選択
    IterativeSolver& preconditioner(PreconditionerType type); // 前処理: None, Jacobi, ILU0, IC0
    IterativeSolver& tolerance(T abs_tol, T rel_tol);
    IterativeSolver& maxIterations(std::size_t n);
    IterativeSolver& restart(std::size_t m);            // GMRES の再起動間隔
    IterativeSolver& initialGuess(Vector<T> x0);
    IterativeSolver& enableLog();                       // 反復ごとの残差を記録
    SolverResult<T> solve();
};

auto result = IterativeSolver<double>(A, b)
    .solver(SolverType::BiCGSTAB)
    .preconditioner(PreconditionerType::ILU0)
    .tolerance(1e-10, 1e-8)
    .maxIterations(500)
    .solve();
// result.success, result.x, result.final_residual

minres — MINRES

対称不定行列用。Lanczos 三重対角化 + Givens QR。

template<typename T, typename IndexType>
SolverResult<T> minres(const SparseMatrix<T, IndexType>& A, const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria, ...)

fgmres — Flexible GMRES

反復ごとに異なる前処理を適用可能。非対称行列用。

template<typename T, typename IndexType>
SolverResult<T> fgmres(const SparseMatrix<T, IndexType>& A, const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria, std::size_t restart = 30, ...)

bicgstab — BiCGSTAB (疎行列版)

非対称行列用。GMRES よりメモリ効率が良い。

template<typename T, typename IndexType>
SolverResult<T> bicgstab(const SparseMatrix<T, IndexType>& A, const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria, ...)

conjugate_gradient — CG (疎行列版)

対称正定値行列用。

template<typename T, typename IndexType>
SolverResult<T> conjugate_gradient(const SparseMatrix<T, IndexType>& A, const Vector<T>& b,
    const ConvergenceCriteria<T>& criteria, ...)

使用例

1. LU 分解で連立方程式を解く

#include <math/linalg/decomposition.hpp>
#include <math/linalg/solvers.hpp>
#include <iostream>
using namespace calx;
using namespace calx::algorithms;

int main() {
    // 3x3 連立方程式 Ax = b
    Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
    Vector<double> b = {8, -11, -3};

    // 方法 1: 自動選択
    auto x = solve(A, b);
    std::cout << "x = [" << x[0] << ", " << x[1] << ", " << x[2] << "]" << std::endl;
    // x = [2, 3, -1]

    // 方法 2: LU 分解を直接使用（分解結果を再利用する場合）
    auto [LU, perm] = lu_decomposition(A);
    auto det = lu_determinant(A);
    auto inv = lu_inverse(A);
    std::cout << "det(A) = " << det << std::endl;
}

2. SVD で低ランク近似

#include <math/linalg/decomposition.hpp>
#include <iostream>
using namespace calx;
using namespace calx::algorithms;

int main() {
    // 3x3 行列
    Matrix<double> A = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    auto [U, sigma, Vt] = svd_decomposition(A);

    std::cout << "Singular values: ";
    for (std::size_t i = 0; i < sigma.size(); ++i)
        std::cout << sigma[i] << " ";
    std::cout << std::endl;

    // ランク 1 近似: 最大特異値のみ使用
    Matrix<double> U1(U.rows(), 1), Vt1(1, Vt.cols());
    for (std::size_t i = 0; i < U.rows(); ++i) U1(i, 0) = U(i, 0);
    for (std::size_t j = 0; j < Vt.cols(); ++j) Vt1(0, j) = Vt(0, j);
    Vector<double> s1 = {sigma[0]};
    auto A1 = reconstruct_matrix_from_svd(U1, s1, Vt1);
    std::cout << "Rank-1 approximation:" << std::endl;
    for (std::size_t i = 0; i < A1.rows(); ++i) {
        for (std::size_t j = 0; j < A1.cols(); ++j)
            std::cout << A1(i, j) << " ";
        std::cout << std::endl;
    }
}

3. gaussian_elimination + PivotStrategy::Full

#include <math/linalg/decomposition.hpp>
#include <iostream>
using namespace calx;
using namespace calx::algorithms;

int main() {
    // 完全ピボットで数値安定性を最大化
    Matrix<double> A = {
        {1e-18, 1.0,  1.0},
        {1.0,   1.0,  2.0},
        {1.0,   2.0,  1.0}
    };
    Vector<double> b = {2.0, 4.0, 4.0};

    // 部分ピボット (既定)
    auto x_partial = gaussian_elimination(A, b, PivotStrategy::Partial);

    // 完全ピボット (最大安定性)
    auto x_full = gaussian_elimination(A, b, PivotStrategy::Full);

    std::cout << "Partial pivot: [" << x_partial[0] << ", "
              << x_partial[1] << ", " << x_partial[2] << "]" << std::endl;
    std::cout << "Full pivot:    [" << x_full[0] << ", "
              << x_full[1] << ", " << x_full[2] << "]" << std::endl;
}

4. Matrix<Rational> でヒルベルト行列の厳密解

#include <math/linalg/decomposition.hpp>
#include <math/core/mp/Rational.hpp>
#include <iostream>
using namespace calx;
using namespace calx::algorithms;

int main() {
    // 4x4 ヒルベルト行列: H(i,j) = 1/(i+j+1)
    // 悪条件行列の代表例 (cond ~ 15514)
    const std::size_t n = 4;
    Matrix<Rational> H(n, n);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            H(i, j) = Rational(1, static_cast<int64_t>(i + j + 1));

    // 右辺: b = H * [1, 1, 1, 1]^T (解が全て 1 になるように設定)
    Vector<Rational> b(n);
    for (std::size_t i = 0; i < n; ++i) {
        Rational sum(0);
        for (std::size_t j = 0; j < n; ++j)
            sum = sum + H(i, j);
        b[i] = sum;
    }

    // Rational で厳密解 (高さ最小ピボット選択)
    auto x = gaussian_elimination(H, b);

    std::cout << "Exact solution of Hilbert system:" << std::endl;
    for (std::size_t i = 0; i < n; ++i)
        std::cout << "  x[" << i << "] = " << x[i] << std::endl;
    // x = {1, 1, 1, 1} (exact)
}

関連する数学的背景

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

• 行列ノルム — 各種ノルムの定義と性質
• ガウスの消去法 — 連立方程式の基本解法
• ガウス＝ジョルダン消去法 — 逆行列の計算
• ランチョスのアルゴリズム — 大規模疎行列の固有値計算
• ランダム化 SVD — 大規模行列の効率的な特異値分解
• 厳密線形代数 — 有理数・整数上の線形代数