LinAlg — Linear Algebra

Overview

Provides matrix decompositions, eigenvalue problems, and linear system solvers.

• Decompositions — LU, Cholesky, QR, SVD, LDL
• Eigenvalues — symmetric / general / generalized / Hermitian, power iteration
• Solvers — direct (LU/Cholesky/QR/SVD/LDL), iterative (CG, Jacobi, Gauss-Seidel, SOR)
• Least squares — ordinary / weighted / regularized (Ridge)
• Special matrices — tridiagonal, banded, sparse (CSR)
• Vector space — Gram-Schmidt, basis completion, linear regression

All functions reside in the calx::algorithms namespace (vector space utilities are in calx).

Complex module: While the Complex module API page is not yet available, the Complex<T> class itself is usable and can be passed to linear algebra functions as Matrix<Complex<double>>. std::complex<double> is also supported as Matrix<std::complex<double>>.

Expression templates: Element-wise matrix operations (addition, subtraction, scalar multiplication/division) are lazily evaluated using expression templates. No intermediate buffers are created when chaining multiple operations; evaluation occurs only at the point of assignment. However, matrix multiplication (operator*) is not lazily evaluated (it computes immediately).

Dimension checking: Dimension mismatches in dynamic Matrix<T> are reported at runtime via DimensionError exceptions. In Debug builds (NDEBUG not defined), assert fires in addition to the exception, allowing the debugger to stop immediately at the point of failure. With StaticMatrix<T, Rows, Cols>, some mismatches are caught at compile time since dimensions are template arguments. Full compile-time checking via static_assert in LinAlg function overloads for StaticMatrix is planned for a future release.

Header

#include <math/linalg/decomposition.hpp>          // Decompositions + utilities
#include <math/linalg/eigenvalues.hpp>            // Eigenvalues
#include <math/linalg/solvers.hpp>                // Solvers
#include <math/linalg/vector_space_utilities.hpp> // Vector space

Header-only (no linking required).

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

Using Intel MKL

If Intel oneAPI MKL is installed, calx can leverage LAPACK routines for improved performance.

• CMake option: -DCALX_HAS_MKL=ON
• Requirement: Intel oneAPI MKL (included in the oneAPI Base Toolkit)
• Link libraries: mkl_intel_lp64, mkl_sequential, mkl_core
• Enabled features: Eigenvalue solvers such as eigen, eigen_symmetric, eigen_hermitian, and eigen_generalized use LAPACK routines (DSYEV, ZHEEV, DGEEV, etc.)

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

Even without MKL, all features remain available through calx's own implementations (Jacobi method, QR iteration, etc.) as a fallback.

LU Decomposition

LU decomposition with partial pivoting. $PA = LU$.

• Complexity: $O(n^3)$
• Use cases: General-purpose computation of linear systems, determinants, and inverses
• Constraint: Square matrix

lu_decomposition

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

Return value: std::pair — the first element is the combined LU matrix (below the diagonal is L, the diagonal and above is U), the second element is the row permutation record (pivot array).

Exception: DimensionError — if the matrix is not square.

using namespace calx::algorithms;

Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
auto [LU, perm] = lu_decomposition(A);
// LU: L and U stored in a single matrix
// perm: row exchange order

lu_solve

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

Return value: Solution vector $\mathbf{x}$ satisfying $A\mathbf{x} = \mathbf{b}$.

Exception: DimensionError — dimension mismatch.

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)

Return value: Inverse matrix $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)

Return value: Determinant $\det(A)$. Computed as the product of diagonal elements of the LU decomposition, multiplied by the sign of the pivot permutation.

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

bareiss_determinant — Bareiss Algorithm

Computes the determinant of an integer matrix without division. By exploiting Sylvester's identity, the divisions arising during elimination are always exact, so all intermediate results remain integers while computing the determinant in $O(n^3)$.

• Complexity: $O(n^3)$
• Algorithm: At each step, computes $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)}}$. This division is always exact by Sylvester's identity
• Use cases: Integer types such as Matrix<Int> or Matrix<int> where division is undefined or should be avoided
• Auto-selection: The Matrix::determinant() member function automatically uses Bareiss for integer types (numeric_traits<T>::is_integer)

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

Return value: Determinant $\det(A)$. All intermediate results are kept as integers.

// Determinant of a multi-precision integer matrix
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 (exact)

// Also works with int
Matrix<int> B({{0, 1, 2}, {1, 0, 3}, {4, 5, 6}});
auto det2 = bareiss_determinant(B);  // 16

// The determinant() member function automatically selects Bareiss for integer types
auto det3 = A.determinant();  // same result as bareiss_determinant

Cholesky Decomposition

Symmetric positive definite matrix $A = LL^\top$. About twice as fast as LU decomposition.

• Complexity: $O(n^3/3)$ (half of LU)
• Use cases: Least squares, Kalman filters, Gaussian processes, normal equations
• Constraint: Symmetric positive definite matrices only

cholesky_decomposition

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

Return value: Lower triangular matrix $L$ satisfying $A = LL^\top$.

Exception: Decomposition fails if the matrix is not positive definite.

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)

Return value: Solution vector $\mathbf{x}$. Internally performs Cholesky decomposition followed by forward and back substitution.

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

QR Decomposition

QR decomposition using Householder reflections. $A = QR$.

• Complexity: $O(mn^2)$ (for an $m \times n$ matrix)
• Algorithm: Householder reflections (more efficient than Givens rotations)
• Use cases: Least squares problems, preprocessing for eigenvalue computation, numerically stable linear systems
• Properties: More numerically stable than LU. Also handles rectangular matrices

qr_decomposition

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

Return value: std::pair — the first element is the orthogonal matrix $Q$ ($m \times m$), the second element is the upper triangular matrix $R$ ($m \times n$).

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
auto [Q, R] = qr_decomposition(A);
// Q: 3x3 orthogonal matrix (Q^T Q = I)
// R: 3x2 upper triangular matrix
// Q * R = A

qr_solve

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

Return value: Solution vector $\mathbf{x}$. For overdetermined systems ($m > n$), returns the least squares solution $\min\|A\mathbf{x} - \mathbf{b}\|_2$.

// Overdetermined system (3x2): least squares solution
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 (Singular Value Decomposition)

Singular value decomposition using the Golub-Kahan algorithm. $A = U\Sigma V^\top$.

• Complexity: $O(mn \cdot \min(m,n))$
• Algorithm: Householder bidiagonalization followed by implicit-shift QR iteration (ref: Golub & Van Loan "Matrix Computations" §8.6)
• Use cases: Rank determination, condition number estimation, least squares, low-rank approximation, principal component analysis
• Properties: The most stable decomposition. Handles ill-conditioned matrices

svd_decomposition

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

Return value: std::tuple<U, σ, V> — $U$ is the left singular vector matrix ($m \times k$), $\sigma$ is the singular value vector ($k$ elements, in descending order), $V^\top$ is the right singular vector matrix ($k \times n$). $k = \min(m, n)$.

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

// Reconstruction: 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))

Return value: Solution vector $\mathbf{x}$. Returns a numerically stable solution via damping, even for ill-conditioned matrices.

When to use: If the condition number is good, lu_solve is the fastest. For ill-conditioned or near-singular matrices, use svd_solve. Increasing rcond strengthens the damping and improves stability, at the cost of accuracy.

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)

Return value: Reconstructed matrix $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 (reconstruction error ≈ 0)

svd_determinant

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

Computes the determinant using SVD. Returns the product of singular values $\prod \sigma_i$ multiplied by the sign. More numerically stable than lu_determinant for ill-conditioned matrices.

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

LDL Decomposition

Symmetric matrix $A = LDL^\top$. Also handles indefinite matrices.

• Complexity: $O(n^3/3)$
• Use cases: Symmetric but not necessarily positive definite matrices, such as Hessians in optimization problems
• Properties: Generalization of Cholesky (which requires positive definiteness)

ldl_decomposition

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

Return value: std::pair — the first element is the lower triangular matrix $L$ (with unit diagonal), the second element is the diagonal element vector $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)

Return value: Solution vector $\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}

Choosing a Decomposition

Gaussian Elimination

Performs forward elimination and back substitution on the augmented matrix $[A|\mathbf{b}]$ to directly solve the linear system $A\mathbf{x} = \mathbf{b}$. Mathematically equivalent to LU decomposition, but operates directly on the augmented matrix, making it more concise when the decomposition will not be reused.

gaussian_elimination

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

Return value: Solution vector $\mathbf{x}$.

Exceptions:

• DimensionError — matrix is not square, or dimension mismatch
• MathError — singular matrix detected

Special behavior for Rational type: When T = Rational, pivot selection chooses the element with the minimum height (the one whose numerator and denominator are smallest in absolute value) rather than the element with maximum absolute value. This suppresses numerator/denominator growth in Rational arithmetic and maintains computational efficiency.

using namespace calx::algorithms;

// double with partial pivoting (default)
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 for exact solution
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} (exact rational solution)

Pivot Strategies

PivotStrategy is an enumeration used by both LU decomposition and Gaussian elimination.

enum class PivotStrategy {
    None,       // No pivoting (use diagonal element as-is)
    Partial,    // Partial pivoting (max absolute value in column; min height for Rational)
    Full        // Full pivoting (search entire remaining submatrix)
};

For Rational type: The selection criterion for Partial and Full differs from floating-point types. Instead of maximum absolute value, the pivot with minimum height is selected. The height is defined as $\max(|\text{numerator}|, |\text{denominator}|)$; a smaller height suppresses numerator/denominator growth in subsequent operations.

How Pivot Strategy Affects Results

For ill-conditioned matrices, omitting pivoting can severely degrade numerical accuracy. The following example shows a case where the (1,1) entry is extremely small.

using namespace calx::algorithms;

// Ill-conditioned: (1,1) entry is extremely small
Matrix<float> A({{1e-8f, 1.0f}, {1.0f, 1.0f}});
Vector<float> b({1.0f, 2.0f});
// True solution: x = {1.00000001..., 0.99999999...} ≈ {1.0, 1.0}

// No pivoting → divides by 1e-8, rounding error explodes
auto x_none = gaussian_elimination(A, b, PivotStrategy::None);
// x_none = {0.000000, 1.000000}  ← x_none[0] is completely wrong

// Partial pivoting → swaps rows to process the larger entry first
auto x_partial = gaussian_elimination(A, b, PivotStrategy::Partial);
// x_partial = {1.000000, 1.000000}  ← accurate

Without pivoting, dividing row 1 by $10^{-8}$ introduces a scaling factor of $10^8$. Since float has only 23 mantissa bits (~7 digits), $1.0$ is rounded away when subtracted from $10^8$ during elimination of row 2. Partial pivoting swaps the rows so that $|a_{11}| = 1$, avoiding this catastrophic cancellation.

Eigenvalues

eigen — General Eigenvalue Decomposition

Computes all eigenvalues and eigenvectors of a general (including non-symmetric) matrix.

• Algorithm: Automatically detects symmetry; if symmetric, delegates to eigen_symmetric. For non-symmetric matrices, uses Hessenberg reduction + Givens rotation-based QR iteration. Uses LAPACK routines when MKL is available
• Complexity: $O(n^3)$

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

Return value: std::pair — the first element is the eigenvalue vector (Vector<complex<T>>), the second element is the eigenvector matrix (Matrix<complex<T>>, where each column is an eigenvector).

Matrix<double> A = {{0, -1}, {1, 0}};  // rotation matrix
auto [eigenvalues, eigenvectors] = eigen(A);
// eigenvalues = {i, -i} (purely imaginary)

eigen_symmetric — Symmetric Matrix Eigenvalues

Computes all eigenvalues and eigenvectors of a real symmetric matrix. Eigenvalues are real.

• Algorithm: Jacobi method (diagonalization by rotations). Uses LAPACK when MKL is available
• Complexity: $O(n^3)$
• Properties: All eigenvalues are guaranteed real. Eigenvectors are orthogonal

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}
// each column of vecs is the corresponding eigenvector

eigen_hermitian — Hermitian Matrix Eigenvalues

Computes all eigenvalues and eigenvectors of a Hermitian matrix ($A^H = A$). Eigenvalues are real.

• Algorithm: Uses LAPACKE_zheev / LAPACKE_cheev when MKL is available. Without MKL, the $n \times n$ Hermitian matrix is converted to a $2n \times 2n$ real symmetric matrix and solved via eigen_symmetric
• Works without MKL: A fallback implementation is available in all environments

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

// Eigenvalues of a Hermitian matrix are real
Matrix<std::complex<double>> H = {{{2,0}, {1,-1}}, {{1,1}, {3,0}}};
auto [vals, vecs] = eigen_hermitian(H);
// vals ≈ {1.268, 3.732} (real)

eigen_generalized — Generalized Eigenvalue Problem

Solves $A\mathbf{x} = \lambda B\mathbf{x}$. $A$ and $B$ are general square matrices.

• Use cases: Vibration analysis, buckling problems, generalized PCA

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 — Symmetric Generalized Eigenvalue Problem

Solves $A\mathbf{x} = \lambda B\mathbf{x}$. $A$ is symmetric and $B$ is symmetric positive definite.

• Algorithm: Decomposes $B = LL^\top$ via Cholesky, then solves the symmetric eigenvalue problem for $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 — Complex Matrix Eigenvalues

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

Returns the eigenvalues and eigenvectors of a general real matrix as complex numbers. Equivalent to eigen, but returns eigenvectors as Complex<T> type.

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

Choosing an Eigenvalue Solver

Power Iteration

Iterative methods for computing specific eigenvalues. Useful when not all eigenvalues are needed for large matrices.

power_method

Computes the dominant eigenvalue (largest in absolute value) and its eigenvector.

• Convergence rate: $O(|\lambda_1/\lambda_2|^k)$ — depends on the ratio of the largest to second-largest eigenvalue
• Complexity: $O(n^2)$ per iteration (matrix-vector product)

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

Computes the eigenvalue closest to a given shift $\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 (smallest eigenvalue)

rayleigh_quotient_iteration

Rayleigh quotient iteration. Updates the shift at each iteration, achieving cubic convergence.

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 (converges rapidly to the eigenvalue nearest to v0)

schur_decomposition — Schur Decomposition

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

Computes the real Schur decomposition $A = QTQ^\top$. $T$ is a quasi-upper-triangular matrix (real eigenvalues on the diagonal, complex eigenvalues as $2 \times 2$ blocks), and $Q$ is an orthogonal matrix. Used internally by matrix functions (expm, logm, sqrtm).

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

Solvers

SolverType

enum class SolverType {
    Auto,       // Automatically selects the best solver based on matrix properties
    LU,         // LU decomposition (general-purpose)
    Cholesky,   // Cholesky decomposition (symmetric positive definite)
    QR,         // QR decomposition (numerical stability, overdetermined systems)
    SVD,        // SVD (ill-conditioned, rank-deficient)
    LDL,        // LDL decomposition (symmetric indefinite)
    Iterative   // Iterative method (CG)
};

Auto selection logic:

1. Symmetric and positive definite → Cholesky
2. Symmetric → LDL
3. Square → LU
4. Overdetermined → QR

solve

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

Return value: Solution vector $\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-selection
auto x2 = solve(A, b, SolverType::LU);      // explicit LU
auto x3 = solve(A, b, "cholesky");           // string specification also available

A string overload is also provided. Available strings: "auto", "lu", "cholesky", "qr", "svd", "ldl", "iterative".

solve_multi — Multiple Right-Hand Sides

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

Return value: Solution matrix $X$ ($n \times m$). Solves each column of $AX = B$ individually via solve.

Matrix<double> A = {{2, 1}, {5, 3}};
Matrix<double> B = {{1, 0}, {0, 1}};  // identity → computes the inverse
auto X = solve_multi(A, B);
// X = {{3, -1}, {-5, 2}} = A^{-1}

Iterative Solvers

Iterative methods that are effective when direct solvers are too costly for large matrices, or when the matrix is sparse.

conjugate_gradient — Conjugate Gradient (CG)

• Convergence rate: $O(\sqrt{\kappa})$ iterations ($\kappa$: condition number)
• Constraint: Symmetric positive definite matrices only
• Complexity: $O(n^2)$ per iteration (dense matrix) or $O(\mathrm{nnz})$ (sparse matrix)

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 Iteration

• Convergence condition: Guaranteed for diagonally dominant matrices ($|a_{ii}| > \sum_{j \ne i} |a_{ij}|$)
• Convergence rate: Slower than CG but simple to implement. Well-suited for parallelization

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 Iteration

• Convergence rate: About twice as fast as Jacobi (under the same conditions)
• Convergence condition: Guaranteed for diagonally dominant or symmetric positive definite matrices

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} (converges faster than Jacobi)

sor_method — SOR (Successive Over-Relaxation)

• Convergence condition: Guaranteed for diagonally dominant matrices when $0 < \omega < 2$. Reduces to Gauss-Seidel when $\omega = 1$

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}

Choosing an Iterative Method

Least Squares

least_squares

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

Return value: Least squares solution $\mathbf{x} = \arg\min \|A\mathbf{x} - \mathbf{b}\|_2$ ($n$ elements).

Auto selection logic: Uses QR decomposition for overdetermined systems ($m \ge n$), and SVD for underdetermined systems ($m < n$). Solving the normal equations $A^\top A \mathbf{x} = A^\top \mathbf{b}$ directly is not used due to numerical instability.

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}  (least squares fit 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)

Return value: Weighted least squares solution $\mathbf{x} = \arg\min \sum w_i (A_i\mathbf{x} - b_i)^2$. Internally multiplies each row of $A$ and each element of $b$ by $\sqrt{w_i}$ and reduces to an ordinary least squares problem.

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
Vector<double> b = {1, 2, 2};
Vector<double> w = {1, 1, 2};  // emphasize the 3rd observation
auto x = weighted_least_squares(A, b, w);

regularized_least_squares — Ridge Regression

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

Return value: Regularized least squares solution $\mathbf{x} = (A^\top A + \lambda I)^{-1} A^\top \mathbf{b}$.

Larger $\lambda$ reduces the solution norm and improves stability, but increases bias. It is common practice to select the optimal $\lambda$ via cross-validation.

Matrix<double> A = {{1, 1}, {1, 2}, {1, 3}};
Vector<double> b = {1, 2, 2};
auto x = regularized_least_squares(A, b, 0.1);
// lambda = 0.1: the norm of solution x is suppressed

Special Matrix Solvers

tridiagonal_solve — Thomas Algorithm

Fast solver for tridiagonal matrices.

• Complexity: $O(n)$ (dramatically faster than LU's $O(n^3)$)
• Use cases: Spline interpolation, finite difference methods (FDM), implicit time integration

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};      // sub-diagonal
Vector<double> b = {2, 3, 2};   // main diagonal
Vector<double> c = {1, 1};      // super-diagonal
Vector<double> d = {1, 2, 3};   // right-hand side
auto x = tridiagonal_solve(a, b, c, d);

band_solve — Banded Matrix Solver

Banded matrix LU solver with partial pivoting.

• Complexity: $O(n \cdot (k_l + k_u)^2)$
• Optimization: Automatically delegates to the Thomas algorithm when $k_l = k_u = 1$ (tridiagonal)

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)

// Tridiagonal matrix (bandwidth 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}

Sparse Matrix Solvers (CSR Format)

Solvers for sparse matrices stored in CSR (Compressed Sparse Row) format.

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)

Internally delegates to sparse_bicgstab (supports non-symmetric matrices).

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 — Sparse 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)

Constraint: Symmetric positive definite matrices only. Use sparse_bicgstab for non-symmetric matrices.

// SPD sparse matrix
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 — Sparse 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)

Unlike CG, supports non-symmetric matrices. More memory-efficient than 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}

Matrix Functions

Analytic functions applied to matrices. Internally uses Schur decomposition and Parlett recurrence.

expm — Matrix Exponential

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

Computes $e^A$ using the scaling-and-squaring method with Padé approximation (Higham 2005). The Padé order is automatically selected from [3/3] to [13/13] based on the 1-norm of the matrix. For matrices with a small norm, a low-order approximation is sufficient and fast; for large norms, scaling (halving) combined with the maximum order [13/13] ensures high accuracy. There is currently no parameter to manually specify the order.

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

sqrtm — Matrix Square Root

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

Computes the principal square root $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 — Matrix Logarithm

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

Computes the principal logarithm $\log(A)$. $e^{\log(A)} = A$.

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

Gershgorin Discs

Estimates the region where eigenvalues lie. For each row, returns a disc with center $a_{ii}$ and radius $\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)

Returns the interval $[\lambda_{\min}, \lambda_{\max}]$ that contains all eigenvalues.

auto [lo, hi] = gershgorin_bounds(A);
// all eigenvalues lie within [lo, hi]

Matrices with Complex Eigenvalues

gershgorin_bounds uses std::abs() for radius computation, so it works correctly with matrices that have complex-valued elements. However, the returned bounds represent an interval on the real axis and do not reflect the imaginary parts of eigenvalues.

To locate complex eigenvalues precisely, inspect the individual Gershgorin discs via gershgorin_discs (center: $a_{ii}$, radius: $R_i = \sum_{j \ne i}|a_{ij}|$). Each eigenvalue lies within at least one disc.

Vector Space Utilities

Provides fundamental vector space operations: Gram-Schmidt orthonormalization to convert a set of vectors into an orthogonal basis, projection of one vector onto another, linear independence testing, and linear regression via least squares.

Typical use cases: constructing and orthogonalizing bases (post-processing eigenvectors, subspace projection), linear fitting of data, and computing angles between vectors.

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 is the orthogonalized set of vectors

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)

Linear regression via least squares: $\beta = (X^\top X)^{-1} X^\top y$.

// Fit 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)

Returns the angle (in radians) between two vectors.

BLAS (Basic Linear Algebra Subprograms)

BLAS (Basic Linear Algebra Subprograms) is a standard interface defining fundamental linear algebra operations. calx implements the key Level 1 through Level 3 BLAS functions in the calx::blas namespace.

These functions are used internally by decomposition algorithms and solvers, and can also be called directly when you need explicit control over matrix-vector operations (e.g., performing update operations like $y \leftarrow \alpha A x + \beta y$ without buffer allocation).

Level 1: Vector-Vector Operations

dot — Inner Product

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-Norm

Computes the Euclidean norm (2-norm) of a vector: $\|x\|_2 = \sqrt{\sum x_i^2}$.

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

asum — 1-Norm

Computes the sum of absolute values (1-norm) of a vector: $\|x\|_1 = \sum |x_i|$.

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

scal — Scaling

Scales a vector by a scalar: $x \leftarrow \alpha x$. In-place operation.

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

axpy — AXPY

$y \leftarrow \alpha x + y$. Performs vector scaling and addition in a single pass. One of the most frequently used BLAS functions.

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: Matrix-Vector Operations

gemv — General Matrix-Vector Product

General matrix-vector product: $y \leftarrow \alpha \cdot \mathrm{op}(A) \cdot x + \beta \cdot y$. Set trans=true for transpose. Buffer reuse eliminates temporary objects.

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 — Triangular Solve

Solves the triangular system $\mathrm{op}(A) x = b$ via forward/back substitution. The upper/trans/unit_diag flags specify the triangular matrix type.

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: Matrix-Matrix Operations

gemm — General Matrix Multiply

General matrix multiplication: $C \leftarrow \alpha \cdot \mathrm{op}(A) \cdot \mathrm{op}(B) + \beta \cdot C$. The core Level 3 BLAS routine. For large matrices, memory access patterns dominate performance.

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 — Triangular Solve (Matrix Version)

Solves a triangular matrix equation: $B \leftarrow \alpha \cdot \mathrm{op}(A)^{-1} \cdot B$ or $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)⁻¹

Randomized SVD

Halko-Martinsson-Tropp (2011) algorithm. Rapidly computes the top $k$ singular values and vectors of a large matrix.

randomizedSVD

Parameters: k = number of singular values to compute, oversampling = oversampling count (improves accuracy), powerIter = number of power iterations (effective when singular value decay is slow), seed = random 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);  // large matrix
// ... fill A ...
auto result = randomizedSVD(A, 10);  // top 10 singular values
// result.U: (1000×10), result.S: 10 values, result.Vt: (10×500)

Sparse Direct Solvers

Direct decomposition solvers for sparse matrices. #include <math/linalg/simplicial_cholesky.hpp> etc.

SimplicialLLT — Sparse 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 decomposition $A = LL^\top$ for symmetric positive definite sparse matrices.

SparseMatrix<double> A(100, 100);
// ... build symmetric positive definite matrix ...
SimplicialLLT<double> chol(A);
auto x = chol.solve(b);

SimplicialLDLT — Sparse 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$ decomposition for symmetric sparse matrices. No sqrt needed; numerically stable.

SparseLU — Sparse 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 decomposition $PA = LU$ for non-symmetric sparse matrices.

SparseQR — Sparse 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 decomposition for sparse matrices. Also usable for least squares problems.

Sparse Iterative Solvers

Namespace calx::sparse_algorithms. Iterative solvers for large-scale sparse matrices.

IterativeSolver — Builder Pattern

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, or MINRES
    IterativeSolver& preconditioner(PreconditionerType type); // None, Jacobi, ILU0, or IC0
    IterativeSolver& tolerance(T abs_tol, T rel_tol);
    IterativeSolver& maxIterations(std::size_t n);
    IterativeSolver& restart(std::size_t m);            // GMRES restart interval
    IterativeSolver& initialGuess(Vector<T> x0);
    IterativeSolver& enableLog();                       // Record residual at each iteration
    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

For symmetric indefinite matrices. Uses Lanczos tridiagonalization + 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

Allows different preconditioners at each iteration. For non-symmetric matrices.

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 (Sparse Version)

For non-symmetric matrices. More memory-efficient than 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 (Sparse Version)

For symmetric positive definite matrices.

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

Examples

1. Solving a Linear System with LU Decomposition

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

int main() {
    // 3x3 linear system Ax = b
    Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
    Vector<double> b = {8, -11, -3};

    // Method 1: auto-selection
    auto x = solve(A, b);
    std::cout << "x = [" << x[0] << ", " << x[1] << ", " << x[2] << "]" << std::endl;
    // x = [2, 3, -1]

    // Method 2: direct LU decomposition (when reusing the decomposition)
    auto [LU, perm] = lu_decomposition(A);
    auto det = lu_determinant(A);
    auto inv = lu_inverse(A);
    std::cout << "det(A) = " << det << std::endl;
}

2. Low-Rank Approximation with SVD

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

int main() {
    // 3x3 matrix
    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;

    // Rank-1 approximation: use only the largest singular value
    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() {
    // Full pivoting for maximum numerical stability
    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};

    // Partial pivoting (default)
    auto x_partial = gaussian_elimination(A, b, PivotStrategy::Partial);

    // Full pivoting (maximum stability)
    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. Exact Solution of a Hilbert Matrix with 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 Hilbert matrix: H(i,j) = 1/(i+j+1)
    // A classic example of an ill-conditioned matrix (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));

    // Right-hand side: b = H * [1, 1, 1, 1]^T (so the solution is all ones)
    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;
    }

    // Exact solution using Rational (minimum-height pivot selection)
    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)
}

Related Mathematical Background

The following articles explain the mathematical concepts underlying the linear algebra module.

• Matrix Norms — Definitions and properties of various norms
• Gaussian Elimination — Fundamental method for solving linear systems
• Gauss-Jordan Elimination — Matrix inversion
• Lanczos Algorithm — Eigenvalue computation for large sparse matrices
• Randomized SVD — Efficient singular value decomposition for large matrices
• Exact Linear Algebra — Linear algebra over rationals and integers