LinAlg デモ — 線形代数ソルバーの活用
LinAlg モジュールは直接法・反復法のソルバー、固有値分解、行列関数などを提供する。
Matrix<T> / Vector<T> と組み合わせて使う。
このページでは 4 つのデモを通じて、線形代数ソルバーの実用例を紹介する。 以下の出力はすべて calx で実際に計算したものである。 ソースコードはページ末尾に掲載しており、 calx をビルドすれば手元で同じ結果を確認できる。
Demo 1: 連立方程式 — LU / QR 統合インターフェース
連立方程式 $Ax = b$ を LU 分解と QR 分解の 2 つの方法で解く。
lu_solve() は直接呼び出し、solve() は統合インターフェースで
SolverType を指定して切り替えられる。
=== demo_linear_solve ===
x (LU) = [ 2 3 -1 ]
x (QR) = [ 2 3 -1 ]
LU 分解と QR 分解の両方で同じ解 $x = (2, 3, -1)$ が得られる。QR 分解は計算コストが高いが、悪条件系で数値安定性に優れる。
Demo 2: 対称行列の固有値分解
対称行列 $S$ の固有値と固有ベクトルを eigen_symmetric() で求める。
対称行列の固有値は常に実数であり、固有ベクトルは直交系を成す。
=== demo_eigen_symmetric ===
eigenvalues = [ 1.38197 3 4.61803 ]
eigenvectors (columns):
[ -0.371748 -0.707107 0.601501 ]
[ 0.601501 0 0.371748 ]
[ -0.371748 0.707107 0.601501 ]
固有値は昇順 $(\lambda_1 \approx 1.38,\; \lambda_2 = 3,\; \lambda_3 \approx 4.62)$ で返される。$\lambda_1 + \lambda_2 + \lambda_3 = 9 = \operatorname{tr}(S)$ が成り立つ。
Demo 3: 行列指数関数 — 回転行列の生成
交代行列 (skew-symmetric) の行列指数関数 $e^A$ は回転行列になる。
expm() で $2 \times 2$ の交代行列から回転行列を生成する。
=== demo_expm ===
expm(A) =
[ 0.540302 -0.841471 ]
[ 0.841471 0.540302 ]
expected: [[ 0.540302 -0.841471 ] [ 0.841471 0.540302 ]]
$e^A$ の結果が $(\cos 1, \sin 1) \approx (0.5403, 0.8415)$ を使った回転行列と一致する。Lie 群 $\mathrm{SO}(2)$ の指数写像そのものである。
Demo 4: 共役勾配法 — 反復ソルバー
対称正定値 (SPD) 行列に対して、conjugate_gradient() で反復的に解を求める。
直接法と異なり、行列を陽に分解せずに行列ベクトル積だけで解ける。
大規模疎行列に適する。
$5 \times 5$ の対角優位対称行列（三重対角）を構成する:
$$A = \begin{pmatrix} 10 & -1 & & & \\ -1 & 10 & -1 & & \\ & -1 & 10 & -1 & \\ & & -1 & 10 & -1 \\ & & & -1 & 10 \end{pmatrix}, \quad b = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{pmatrix}$$=== demo_iterative_solver ===
x (CG) = [ 0.109091 0.118182 0.12 0.118182 0.109091 ]
||Ax - b|| = 1.11022e-15
残差ノルム $\|Ax - b\| \approx 10^{-15}$ で機械精度まで収束している。解の対称性 $x_1 = x_5$, $x_2 = x_4$ は行列とベクトルの対称性を反映している。
ソースコードと実行方法
このページの全出力は以下の C++ プログラムから生成された。 calx をビルドすれば、手元で同じ結果を確認できる。
example_linalg.cpp (クリックで展開/折りたたみ)
// example_linalg.cpp
// Demonstrates linear algebra features from the calx API reference.
//
// Build:
// cd build && cmake --build . --config Release --target example-linalg
#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/linalg/solvers.hpp>
#include <math/linalg/eigenvalues.hpp>
#include <math/linalg/matrix_functions.hpp>
#include <iostream>
#include <cmath>
using namespace calx;
using namespace calx::algorithms;
// ---------------------------------------------------------------------------
// Snippet 1: Solve a linear system Ax = b (LU + unified interface)
// ---------------------------------------------------------------------------
void demo_linear_solve() {
std::cout << "=== demo_linear_solve ===" << std::endl;
Matrix<double> A({{2, 1, -1},
{-3, -1, 2},
{-2, 1, 2}});
Vector<double> b({8, -11, -3});
// LU solver
auto x = lu_solve(A, b);
std::cout << "x (LU) = [";
for (std::size_t i = 0; i < x.size(); ++i)
std::cout << " " << x[i];
std::cout << " ]" << std::endl;
// Unified interface with QR
auto x2 = solve(A, b, SolverType::QR);
std::cout << "x (QR) = [";
for (std::size_t i = 0; i < x2.size(); ++i)
std::cout << " " << x2[i];
std::cout << " ]" << std::endl;
std::cout << std::endl;
}
// ---------------------------------------------------------------------------
// Snippet 2: Eigenvalue decomposition of a symmetric matrix
// ---------------------------------------------------------------------------
void demo_eigen_symmetric() {
std::cout << "=== demo_eigen_symmetric ===" << std::endl;
Matrix<double> S({{4, 1, 1},
{1, 3, 0},
{1, 0, 2}});
auto [eigenvalues, eigenvectors] = eigen_symmetric(S);
std::cout << "eigenvalues = [";
for (std::size_t i = 0; i < eigenvalues.size(); ++i)
std::cout << " " << eigenvalues[i];
std::cout << " ]" << std::endl;
std::cout << "eigenvectors (columns):" << std::endl;
for (std::size_t i = 0; i < eigenvectors.rows(); ++i) {
std::cout << " [";
for (std::size_t j = 0; j < eigenvectors.cols(); ++j)
std::cout << " " << eigenvectors(i, j);
std::cout << " ]" << std::endl;
}
std::cout << std::endl;
}
// ---------------------------------------------------------------------------
// Snippet 3: Matrix exponential (rotation matrix from skew-symmetric)
// ---------------------------------------------------------------------------
void demo_expm() {
std::cout << "=== demo_expm ===" << std::endl;
Matrix<double> A({{0, -1},
{1, 0}});
auto E = expm(A);
std::cout << "expm(A) =" << std::endl;
for (std::size_t i = 0; i < E.rows(); ++i) {
std::cout << " [";
for (std::size_t j = 0; j < E.cols(); ++j)
std::cout << " " << E(i, j);
std::cout << " ]" << std::endl;
}
std::cout << "expected: [[ " << std::cos(1.0) << " " << -std::sin(1.0)
<< " ] [ " << std::sin(1.0) << " " << std::cos(1.0) << " ]]"
<< std::endl;
std::cout << std::endl;
}
// ---------------------------------------------------------------------------
// Snippet 4: Iterative solver (CG on a symmetric positive-definite matrix)
// ---------------------------------------------------------------------------
void demo_iterative_solver() {
std::cout << "=== demo_iterative_solver ===" << std::endl;
const int n = 5;
Matrix<double> A(n, n, 0.0);
for (int i = 0; i < n; ++i) {
A(i, i) = 10.0;
if (i > 0) { A(i, i - 1) = -1.0; A(i - 1, i) = -1.0; }
}
Vector<double> b(n, 1.0);
auto x = conjugate_gradient(A, b, 1e-12, 5000);
std::cout << "x (CG) = [";
for (std::size_t i = 0; i < x.size(); ++i)
std::cout << " " << x[i];
std::cout << " ]" << std::endl;
// Verify residual
auto r = A * x;
double res = 0.0;
for (std::size_t i = 0; i < b.size(); ++i)
res += (r[i] - b[i]) * (r[i] - b[i]);
std::cout << "||Ax - b|| = " << std::sqrt(res) << std::endl;
std::cout << std::endl;
}
// ---------------------------------------------------------------------------
int main() {
demo_linear_solve();
demo_eigen_symmetric();
demo_expm();
demo_iterative_solver();
return 0;
}API の詳細は LinAlg API リファレンス を参照のこと。
ビルドと実行
cd build && cmake --build . --config Release --target example-linalg