# Linear Algebra

numeris provides six matrix decompositions, each available as:

- **Free functions** operating in-place on `&mut impl MatrixMut<T>` — work on both `Matrix` and `DynMatrix`
- **Wrapper structs** with `solve()`, `inverse()`, `det()` convenience methods
- **Matrix convenience methods** — `a.lu()`, `a.cholesky()`, etc.

## Overview

| Decomposition | Struct | Dynamic wrapper | Use case |
|---|---|---|---|
| LU | `LuDecomposition` | `DynLu` | General square systems, det, inverse |
| Cholesky | `CholeskyDecomposition` | `DynCholesky` | SPD / Hermitian PD systems |
| QR | `QrDecomposition` | `DynQr` | Least-squares, orthogonal basis |
| SVD | `SvdDecomposition` | `DynSvd` | Rank, condition number, pseudoinverse |
| Symmetric Eigen | `SymmetricEigen` | `DynSymmetricEigen` | Real eigenvalues, eigenvectors |
| Schur | `SchurDecomposition` | `DynSchur` | General eigenvalues (complex pairs) |

## LU Decomposition

Partial pivoting LU:

$$
P A = L U,
$$

with `P` a row-permutation matrix, `L` unit lower-triangular, and `U` upper-triangular.

```rust
use numeris::{Matrix, Vector};

let a = Matrix::new([
    [2.0_f64, 1.0, -1.0],
    [-3.0,   -1.0,  2.0],
    [-2.0,    1.0,  2.0],
]);
let b = Vector::from_array([8.0, -11.0, -3.0]);

let lu = a.lu().unwrap();

// Solve Ax = b
let x   = lu.solve(&b);       // x = [2, 3, -1]

// Other operations
let inv = lu.inverse();        // A^{-1}
let det = lu.det();            // determinant

// Direct solve on Matrix (uses LU internally)
let x2 = a.solve(&b).unwrap();
```

!!! note "Small matrices"
    For N ≤ 4, `inverse()` and `det()` use direct closed-form formulas (adjugate method), bypassing LU entirely. This is 2–3x faster than the general path.

## Cholesky Decomposition

For symmetric positive-definite (SPD) matrices:

$$
A = L L^{\mathsf{H}},
$$

with `L` lower-triangular ($L^{\mathsf{H}} = L^{\top}$ in the real case). Works
with both real and complex (Hermitian positive-definite) matrices.

```rust
use numeris::Matrix;

let a = Matrix::new([
    [4.0_f64, 2.0],
    [2.0,     3.0],
]);

let chol = a.cholesky().unwrap();

// Solve Ax = b (forward + back substitution, no LU needed)
let b = numeris::Vector::from_array([1.0_f64, 2.0]);
let x = chol.solve(&b);

// Other operations
let inv    = chol.inverse();
let det    = chol.det();
let ln_det = chol.ln_det();    // log-determinant (stable for large matrices)
```

## QR Decomposition

Householder QR:

$$
A = Q R,
$$

with `Q` orthogonal ($Q^{\mathsf{H}} Q = I$) and `R` upper-triangular.

```rust
use numeris::{Matrix, Vector};

// Square system — exact solve
let a = Matrix::new([[2.0_f64, 1.0], [1.0, 3.0]]);
let b = Vector::from_array([5.0_f64, 10.0]);
let qr = a.qr().unwrap();
let x  = qr.solve(&b);

// Overdetermined system — least-squares
let a_rect = Matrix::new([
    [1.0_f64, 0.0],
    [1.0,     1.0],
    [1.0,     2.0],
]);
let b_rect = Vector::from_array([1.0_f64, 2.0, 4.0]);
let x_ls = a_rect.qr().unwrap().solve(&b_rect);

// Determinant
let det = a.qr().unwrap().det();
```

## SVD

Golub-Kahan implicit-shift QR on a bidiagonal form.

$$
A = U \Sigma V^{\top},
$$

where `A` is $M \times N$ ($M \ge N$), `U` and `V` are orthogonal, and
$\Sigma = \operatorname{diag}(\sigma_0 \ge \sigma_1 \ge \cdots \ge 0)$.

```rust
use numeris::Matrix;

let a = Matrix::new([
    [1.0_f64, 2.0],
    [3.0,     4.0],
    [5.0,     6.0],
]);

let svd = a.svd().unwrap();

// Components
let sigma = svd.singular_values();   // [σ₀, σ₁] sorted descending
let u = svd.u();                      // 3×3 orthogonal (or 3×2 thin)
let vt = svd.vt();                    // 2×2 orthogonal (Vᵀ)

// Rank and condition number
let rank = svd.rank(1e-10);          // rank with tolerance
let cond = svd.condition_number();   // σ_max / σ_min

// Shortcut: only singular values (faster, no U or V)
let sigma_only = a.singular_values_only().unwrap();
```

### SVD Solve (Least-Squares)

SVD-based solve computes the minimum-norm least-squares solution

$$
x = V \Sigma^{+} U^{\mathsf{H}} b,
$$

automatically handling rank-deficient systems ($\Sigma^{+}$ inverts only the
nonzero singular values):

```rust
use numeris::{Matrix, Vector};

// Overdetermined system (3 equations, 2 unknowns)
let a = Matrix::new([
    [1.0_f64, 0.0],
    [0.0,     1.0],
    [1.0,     1.0],
]);
let b = Vector::from_array([1.0_f64, 2.0, 3.0]);

// Via SVD decomposition
let svd = a.svd().unwrap();
let x = svd.solve(&b);

// Or directly on the matrix
let x2 = a.solve_svd(&b).unwrap();
```

### Pseudo-Inverse

The Moore-Penrose pseudo-inverse `A⁺ = V Σ⁺ Uᴴ`, where singular values below `ε · σ_max · max(M, N)` are treated as zero:

```rust
use numeris::Matrix;

let a = Matrix::new([
    [1.0_f64, 0.0],
    [0.0,     2.0],
    [0.0,     0.0],
]);

let a_pinv = a.pinv().unwrap();  // 2×3 pseudo-inverse
```

!!! info "M < N matrices"
    `DynSvd` handles `M < N` by transposing internally. The fixed `SvdDecomposition` requires `M ≥ N` at compile time.

## Matrix Exponential

Computes the matrix exponential

$$
e^{A} = \sum_{k=0}^{\infty} \frac{A^{k}}{k!}
$$

for square matrices using a [13/13] Padé approximation with scaling and squaring (Higham 2005).

```rust
use numeris::Matrix;

// Rotation matrix via expm of an antisymmetric matrix
let theta = std::f64::consts::FRAC_PI_4;
let a = Matrix::new([
    [0.0, -theta],
    [theta, 0.0],
]);
let r = a.expm().unwrap();  // 45° rotation matrix

// Diagonal: expm(diag(a,b)) = diag(e^a, e^b)
let d = Matrix::new([[2.0_f64, 0.0], [0.0, 3.0]]);
let ed = d.expm().unwrap();
assert!((ed[(0,0)] - 2.0_f64.exp()).abs() < 1e-12);
```

Also available as a free function: `numeris::linalg::expm(&a)`.

## Symmetric Eigendecomposition

Householder tridiagonalization + implicit QR with Wilkinson shift.

For real symmetric (or Hermitian) matrices only. Produces real eigenvalues:

$$
A = Q \Lambda Q^{\top},
\qquad \Lambda = \operatorname{diag}(\lambda_0 \le \lambda_1 \le \cdots),
$$

with `Q` orthogonal — its columns are the orthonormal eigenvectors.

```rust
use numeris::Matrix;

let a = Matrix::new([
    [4.0_f64, 1.0, 0.0],
    [1.0,     3.0, 1.0],
    [0.0,     1.0, 2.0],
]);

let eig = a.eig_symmetric().unwrap();

let vals = eig.eigenvalues();    // sorted ascending
let vecs = eig.eigenvectors();   // columns = orthonormal eigenvectors
// Reconstruction: A ≈ vecs * diag(vals) * vecs^T

// Shortcut: eigenvalues only (faster, no eigenvectors)
let vals_only = a.eigenvalues_symmetric().unwrap();
```

## Schur Decomposition

Francis double-shift QR on the upper Hessenberg form.

Produces quasi-upper-triangular `S` with 1×1 (real eigenvalue) and 2×2 (conjugate complex pair) diagonal blocks, and orthogonal `Q` such that

$$
A = Q S Q^{\top}.
$$

```rust
use numeris::Matrix;

// 90° rotation — eigenvalues ±i
let a = Matrix::new([[0.0_f64, -1.0], [1.0, 0.0]]);

let schur = a.schur().unwrap();
let s = schur.s();   // quasi-upper-triangular
let q = schur.q();   // orthogonal factor

// Extract eigenvalues (real and imaginary parts)
let (re, im) = a.eigenvalues().unwrap();
// re = [0, 0], im = [1, -1]  (conjugate pair ±i)
```

## DynMatrix Variants

All six decompositions have dynamic wrappers in `DynMatrix`:

```rust
use numeris::{DynMatrix, DynVector};

let a = DynMatrix::from_rows(3, 3, &[
    4.0_f64, 2.0, 0.0,
    2.0,     3.0, 1.0,
    0.0,     1.0, 2.0,
]);
let b = DynVector::from_slice(&[6.0_f64, 7.0, 3.0]);

let lu   = a.lu().unwrap();
let chol = a.cholesky().unwrap();
let qr   = a.qr().unwrap();
let svd  = a.svd().unwrap();
let eig  = a.eig_symmetric().unwrap();
let sch  = a.schur().unwrap();

let x   = a.solve(&b).unwrap();
let inv = a.inverse().unwrap();
```

## Complex Matrices

Enable the `complex` feature to use decompositions with complex elements:

```toml
numeris = { version = "0.5", features = ["complex"] }
```

```rust
use numeris::{Complex, Matrix, Vector};
type C = Complex<f64>;

let a = Matrix::new([
    [C::new(2.0, 1.0), C::new(1.0, -1.0)],
    [C::new(1.0, 0.0), C::new(3.0, 2.0)],
]);
let b = Vector::from_array([C::new(5.0, 3.0), C::new(7.0, 4.0)]);
let x = a.solve(&b).unwrap();

// Hermitian positive-definite Cholesky (A = L L^H)
let hpd = Matrix::new([
    [C::new(4.0, 0.0), C::new(2.0,  1.0)],
    [C::new(2.0,-1.0), C::new(5.0,  0.0)],
]);
let chol = hpd.cholesky().unwrap();
```

Complex overhead is zero for real code paths — `conj()` and `re()` on `f64` are identity functions, fully inlined and erased by the compiler.

## Error Handling

All fallible operations return `Result<_, LinalgError>`:

```rust
use numeris::linalg::LinalgError;

match a.cholesky() {
    Ok(chol) => { /* use chol */ }
    Err(LinalgError::NotPositiveDefinite) => { /* not SPD */ }
    Err(LinalgError::SingularMatrix)      => { /* singular */ }
    Err(LinalgError::ConvergenceFailure)  => { /* QR iteration didn't converge */ }
    Err(e) => { /* other */ }
}
```
