Enum Matrix 

pub enum Matrix {
    Dense(MatrixData),
    Identity(IdentityMatrixData),
    Zero(ZeroMatrixData),
    Diagonal(DiagonalMatrixData),
    Scalar(ScalarMatrixData),
    UpperTriangular(UpperTriangularMatrixData),
    LowerTriangular(LowerTriangularMatrixData),
    Symmetric(SymmetricMatrixData),
    Permutation(PermutationMatrixData),
}

Unified matrix type that can represent any matrix efficiently

This enum uses zero-cost abstractions to provide a single interface for all matrix types while maintaining optimal memory usage.

Variants

Dense(MatrixData)

Regular dense matrix: O(n²) memory

Identity(IdentityMatrixData)

Identity matrix: O(1) memory

Zero(ZeroMatrixData)

Zero matrix: O(1) memory

Diagonal(DiagonalMatrixData)

Diagonal matrix: O(n) memory

Scalar(ScalarMatrixData)

Scalar matrix: O(1) memory

UpperTriangular(UpperTriangularMatrixData)

Upper triangular: O(n²/2) memory

LowerTriangular(LowerTriangularMatrixData)

Lower triangular: O(n²/2) memory

Symmetric(SymmetricMatrixData)

Symmetric matrix: O(n²/2) memory

Permutation(PermutationMatrixData)

Permutation matrix: O(n) memory

Implementations

impl Matrix

Cholesky decomposition implementation

pub fn cholesky_decomposition(&self) -> Option<CholeskyDecomposition>

Perform Cholesky decomposition for positive definite matrices

Decomposes symmetric positive definite matrix A into A = LL^T where:

• L is lower triangular with positive diagonal elements

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::from_arrays([
    [4, 2],
    [2, 3]
]);

if let Some(chol) = matrix.cholesky_decomposition() {
    let (l_rows, l_cols) = chol.l.dimensions();
    assert_eq!(l_rows, 2);
    assert_eq!(l_cols, 2);
}

pub fn is_positive_definite_cholesky(&self) -> bool

Check if matrix is positive definite using Cholesky test

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let identity = Matrix::identity(3);
assert!(identity.is_positive_definite_cholesky());

let scalar = Matrix::scalar(2, Expression::integer(5));
assert!(scalar.is_positive_definite_cholesky());

impl Matrix

LU decomposition implementation

pub fn lu_decomposition(&self) -> Option<LUDecomposition>

Perform LU decomposition with partial pivoting

Decomposes matrix A into PA = LU where:

• P is a permutation matrix
• L is lower triangular with 1s on diagonal
• U is upper triangular

Examples

use mathhook_core::matrices::Matrix;

let matrix = Matrix::from_arrays([
    [2, 1],
    [4, 3]
]);

let lu = matrix.lu_decomposition().unwrap();
assert!(lu.p.is_some());

impl Matrix

QR decomposition implementation

pub fn qr_decomposition(&self) -> Option<QRDecomposition>

Perform QR decomposition using Gram-Schmidt process

Decomposes matrix A into A = QR where:

• Q is orthogonal (Q^T * Q = I)
• R is upper triangular

Examples

use mathhook_core::matrices::Matrix;

let matrix = Matrix::from_arrays([
    [1, 1],
    [0, 1]
]);

let qr = matrix.qr_decomposition().unwrap();
let (q_rows, q_cols) = qr.q.dimensions();
assert_eq!(q_rows, 2);
assert_eq!(q_cols, 2);

impl Matrix

SVD implementation

pub fn svd_decomposition(&self) -> Option<SVDDecomposition>

Perform Singular Value Decomposition

Decomposes matrix A into A = UΣV^T where:

• U contains left singular vectors (orthogonal)
• Σ contains singular values (diagonal, non-negative)
• V^T contains right singular vectors (orthogonal)

Examples

use mathhook_core::matrices::Matrix;

let matrix = Matrix::from_arrays([
    [1, 2],
    [3, 4]
]);

let svd = matrix.svd_decomposition().unwrap();
let (u_rows, u_cols) = svd.u.dimensions();
assert_eq!(u_rows, 2);

pub fn rank_via_svd(&self) -> usize

Get matrix rank using SVD

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let identity = Matrix::identity(3);
assert_eq!(identity.rank_via_svd(), 3);

let zero = Matrix::zero(3, 3);
assert_eq!(zero.rank_via_svd(), 0);

pub fn condition_number_via_svd(&self) -> Expression

Get condition number using SVD

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let identity = Matrix::identity(2);
let cond = identity.condition_number_via_svd();
assert_eq!(cond, Expression::integer(1));

impl Matrix

pub fn eigen_decomposition(&self) -> Option<EigenDecomposition>

Compute eigenvalues and eigenvectors

Returns eigenvalues and corresponding eigenvectors for real matrices. For matrices with complex eigenvalues, use complex_eigen_decomposition.

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::diagonal(vec![
    Expression::integer(2),
    Expression::integer(3)
]);

let eigen = matrix.eigen_decomposition().unwrap();
assert_eq!(eigen.eigenvalues.len(), 2);
assert_eq!(eigen.eigenvalues[0], Expression::integer(2));
assert_eq!(eigen.eigenvalues[1], Expression::integer(3));

pub fn complex_eigen_decomposition(&self) -> Option<ComplexEigenDecomposition>

Compute complex eigenvalues and eigenvectors

Handles matrices that may have complex eigenvalues and eigenvectors.

Examples

use mathhook_core::matrices::Matrix;

let matrix = Matrix::from_arrays([
    [0, -1],
    [1, 0]
]);

// This matrix has complex eigenvalues ±i
let complex_eigen = matrix.complex_eigen_decomposition();
// Returns None as complex eigenvalue computation requires specialized algorithms
assert!(complex_eigen.is_none());

pub fn eigenvalues(&self) -> Vec<Expression>

Compute only eigenvalues (faster than full decomposition)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::scalar(3, Expression::integer(5));
let eigenvals = matrix.eigenvalues();
assert_eq!(eigenvals.len(), 3);
assert_eq!(eigenvals[0], Expression::integer(5));

pub fn is_diagonalizable(&self) -> bool

Check if matrix is diagonalizable

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let diagonal = Matrix::diagonal(vec![
    Expression::integer(1),
    Expression::integer(2),
    Expression::integer(3)
]);
assert!(diagonal.is_diagonalizable());

let identity = Matrix::identity(3);
assert!(identity.is_diagonalizable());

pub fn power_iteration_eigenvalues(&self) -> Option<EigenDecomposition>

Power iteration method for finding dominant eigenvalue

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::{Expression, expr};

let matrix = Matrix::diagonal(vec![expr!(2), expr!(3)]);
let eigen = matrix.power_iteration_eigenvalues().unwrap();
// Power iteration returns the dominant (largest) eigenvalue
// For symbolic computation, the result may be a complex expression
assert_eq!(eigen.eigenvalues.len(), 1);
// The eigenvalue exists but may not simplify to integer form symbolically
assert!(!eigen.eigenvalues[0].is_zero());

pub fn compute_vector_norm(&self, v: &[Expression]) -> Expression

Compute norm of a vector

pub fn is_small_value(&self, value: &Expression, tolerance: &Expression) -> bool

Check if a value is small (simplified convergence test)

impl Matrix

Matrix power computation using eigendecomposition

pub fn matrix_power_eigen(&self, n: i64) -> Option<Matrix>

Compute matrix power using eigendecomposition (A^n = P D^n P^(-1))

This method is particularly efficient for diagonal and diagonalizable matrices.

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::diagonal(vec![
    Expression::integer(2),
    Expression::integer(3)
]);

let power = matrix.matrix_power_eigen(3).unwrap();
let eigenvals = power.eigenvalues();
assert_eq!(eigenvals[0], Expression::integer(8)); // 2^3
assert_eq!(eigenvals[1], Expression::integer(27)); // 3^3

pub fn matrix_power_special(&self, n: i64) -> Option<Matrix>

Compute matrix power for special cases

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let identity = Matrix::identity(3);
let power = identity.matrix_power_special(5).unwrap();
assert!(matches!(power, Matrix::Identity(_)));

let scalar = Matrix::scalar(2, Expression::integer(3));
let power = scalar.matrix_power_special(2).unwrap();
// (3I)² = 9I

pub fn matrix_exponential_eigen(&self) -> Option<Matrix>

Compute matrix exponential using eigendecomposition exp(A) = P exp(D) P^(-1) where exp(D) = diag(exp(d_1), exp(d_2), …)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::zero(2, 2);
let exp_matrix = matrix.matrix_exponential_eigen().unwrap();
// exp(0) = 1, so result is diagonal(exp(0), exp(0))
let eigenvals = exp_matrix.eigenvalues();
assert_eq!(eigenvals.len(), 2);
// Eigenvalues are exp(0) in symbolic form
assert_eq!(eigenvals[0], Expression::function("exp", vec![Expression::integer(0)]));

pub fn matrix_logarithm_eigen(&self) -> Option<Matrix>

Compute matrix logarithm using eigendecomposition log(A) = P log(D) P^(-1) where log(D) = diag(log(d_1), log(d_2), …)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let identity = Matrix::identity(2);
let log_matrix = identity.matrix_logarithm_eigen().unwrap();
// log(I) has eigenvalues log(1) = 0, so result is diagonal matrix with zeros
let eigenvals = log_matrix.eigenvalues();
assert_eq!(eigenvals.len(), 2);
assert_eq!(eigenvals[0], Expression::function("log", vec![Expression::integer(1)]));
assert_eq!(eigenvals[1], Expression::function("log", vec![Expression::integer(1)]));

pub fn matrix_sqrt_eigen(&self) -> Option<Matrix>

Compute matrix square root using eigendecomposition sqrt(A) = P sqrt(D) P^(-1) where sqrt(D) = diag(sqrt(d_1), sqrt(d_2), …)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::diagonal(vec![
    Expression::integer(4),
    Expression::integer(9)
]);
let sqrt_matrix = matrix.matrix_sqrt_eigen().unwrap();
let eigenvals = sqrt_matrix.eigenvalues();
// Eigenvalues are sqrt(4) and sqrt(9) in symbolic form
assert_eq!(eigenvals.len(), 2);
assert_eq!(eigenvals[0], Expression::pow(Expression::integer(4), Expression::rational(1, 2)));
assert_eq!(eigenvals[1], Expression::pow(Expression::integer(9), Expression::rational(1, 2)));

pub fn is_nilpotent(&self) -> bool

Check if matrix is nilpotent (A^k = 0 for some positive integer k)

Examples

use mathhook_core::matrices::Matrix;

let zero_matrix = Matrix::zero(3, 3);
assert!(zero_matrix.is_nilpotent());

let identity = Matrix::identity(3);
assert!(!identity.is_nilpotent());

pub fn minimal_polynomial(&self) -> CharacteristicPolynomial

Compute the minimal polynomial (smallest degree polynomial that annihilates the matrix)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::diagonal(vec![
    Expression::integer(2),
    Expression::integer(2),
    Expression::integer(3)
]);
let min_poly = matrix.minimal_polynomial();
// For this matrix, minimal polynomial is (λ-2)(λ-3)
assert!(!min_poly.coefficients.is_empty());

impl Matrix

pub fn dimensions(&self) -> (usize, usize)

Get matrix dimensions efficiently

This method provides O(1) dimension lookup for all matrix types.

pub fn get_element(&self, i: usize, j: usize) -> Expression

Get element at position (i, j) efficiently

This method provides optimized element access for each matrix type.

pub fn as_numeric(&self) -> Option<NumericMatrix>

Try to convert this matrix to a NumericMatrix for fast numeric operations.

Returns Some(NumericMatrix) if all elements can be converted to f64, None otherwise (e.g., if matrix contains symbolic expressions).

pub fn is_square(&self) -> bool

Check if this is a square matrix

pub fn is_zero(&self) -> bool

Check if this is a zero matrix

pub fn is_identity(&self) -> bool

Check if this is an identity matrix

pub fn is_diagonal(&self) -> bool

Check if this is a diagonal matrix

pub fn is_symmetric(&self) -> bool

Check if this is symmetric

pub fn optimize(self) -> Matrix

Convert to the most efficient representation

This method analyzes the matrix and converts it to the most memory-efficient representation possible.

pub fn dense(rows: Vec<Vec<Expression>>) -> Self

Create a dense matrix from rows

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::dense(vec![
    vec![Expression::integer(1), Expression::integer(2)],
    vec![Expression::integer(3), Expression::integer(4)]
]);

pub fn identity(size: usize) -> Self

Create an identity matrix of given size Memory efficient: O(1) storage vs O(n²) for dense matrix

Examples

use mathhook_core::matrices::Matrix;

let identity = Matrix::identity(3);
assert_eq!(identity.dimensions(), (3, 3));
assert!(identity.is_identity());

pub fn zero(rows: usize, cols: usize) -> Self

Create a zero matrix of given dimensions Memory efficient: O(1) storage vs O(n*m) for dense matrix

Examples

use mathhook_core::matrices::Matrix;

let zero = Matrix::zero(2, 3);
assert_eq!(zero.dimensions(), (2, 3));
assert!(zero.is_zero());

pub fn diagonal(diagonal_elements: Vec<Expression>) -> Self

Create a diagonal matrix from diagonal elements Memory efficient: O(n) storage vs O(n²) for dense matrix

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let diag = Matrix::diagonal(vec![
    Expression::integer(1),
    Expression::integer(2),
    Expression::integer(3)
]);
assert_eq!(diag.dimensions(), (3, 3));
assert!(diag.is_diagonal());

pub fn scalar(size: usize, scalar_value: Expression) -> Self

Create a scalar matrix (c*I) Memory efficient: O(1) storage vs O(n²) for dense matrix

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let scalar = Matrix::scalar(3, Expression::integer(5));

pub fn upper_triangular(size: usize, elements: Vec<Expression>) -> Self

Create an upper triangular matrix Memory efficient: ~50% storage vs dense matrix

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let upper = Matrix::upper_triangular(3, vec![
    Expression::integer(1), Expression::integer(2), Expression::integer(3),
    Expression::integer(4), Expression::integer(5),
    Expression::integer(6)
]);

pub fn lower_triangular(size: usize, elements: Vec<Expression>) -> Self

Create a lower triangular matrix Memory efficient: ~50% storage vs dense matrix

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let lower = Matrix::lower_triangular(3, vec![
    Expression::integer(1),
    Expression::integer(2), Expression::integer(3),
    Expression::integer(4), Expression::integer(5), Expression::integer(6)
]);

pub fn symmetric(size: usize, elements: Vec<Expression>) -> Self

Create a symmetric matrix Memory efficient: ~50% storage vs dense matrix

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let symmetric = Matrix::symmetric(3, vec![
    Expression::integer(1), Expression::integer(2), Expression::integer(3),
    Expression::integer(4), Expression::integer(5),
    Expression::integer(6)
]);

pub fn permutation(permutation: Vec<usize>) -> Self

Create a permutation matrix Memory efficient: O(n) storage vs O(n²) for dense matrix

Examples

use mathhook_core::matrices::Matrix;

let perm = Matrix::permutation(vec![2, 0, 1]);

pub fn from_arrays<const R: usize, const C: usize>( arrays: [[i64; C]; R], ) -> Self

Create matrix from nested arrays (convenience method)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::from_arrays([
    [1, 2, 3],
    [4, 5, 6]
]);

pub fn from_flat(rows: usize, cols: usize, elements: &[Expression]) -> Self

Create matrix from flat vector (row-major order)

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::Expression;

let matrix = Matrix::from_flat(2, 3, &[
    Expression::integer(1), Expression::integer(2), Expression::integer(3),
    Expression::integer(4), Expression::integer(5), Expression::integer(6)
]);

impl Matrix

pub fn trace(&self) -> Expression

Get the trace (sum of diagonal elements) efficiently

pub fn determinant(&self) -> Result<Expression, MathError>

Get the determinant efficiently (for square matrices)

Returns

Result containing the determinant expression, or MathError for non-square matrices

Errors

Returns DomainError if matrix is not square

Algorithm

• Special matrices (Identity, Zero, Scalar, Diagonal): O(1) or O(n)
• Small matrices (1x1, 2x2): Direct formulas
• Numeric matrices: NumericMatrix fast-path with O(n³) LU decomposition
• Larger symbolic matrices (n≥3): LU decomposition O(n³)

pub fn scalar_multiply(&self, scalar: &Expression) -> Matrix

Scalar multiplication

impl Matrix

pub fn forward_substitution( &self, b: &[Expression], ) -> Result<Vec<Expression>, MathError>

Solve Lx = b for lower triangular L using forward substitution

Arguments

• b - Right-hand side vector

Returns

Solution vector x

Errors

• DivisionByZero if any diagonal element is zero
• DomainError if dimensions don’t match

Algorithm

For i = 0 to n-1: x[i] = (b[i] - Σ(L[i][j] * x[j]) for j < i) / L[i][i]

pub fn backward_substitution( &self, b: &[Expression], ) -> Result<Vec<Expression>, MathError>

Solve Ux = b for upper triangular U using backward substitution

Arguments

• b - Right-hand side vector

Returns

Solution vector x

Errors

• DivisionByZero if any diagonal element is zero
• DomainError if dimensions don’t match

For i = n-1 down to 0: x[i] = (b[i] - Σ(U[i][j] * x[j]) for j > i) / U[i][i]

pub fn solve(&self, b: &[Expression]) -> Result<Vec<Expression>, MathError>

Solve Ax = b using optimal decomposition

Arguments

• b - Right-hand side vector

Returns

Solution vector x

Errors

• DomainError if matrix is not square or dimensions don’t match
• DivisionByZero if matrix is singular

Algorithm Selection

• Symmetric positive definite matrices: Cholesky (LL^T), ~2x faster
• General square matrices: LU decomposition with partial pivoting

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::expr;

let a = Matrix::from_arrays([[2, 1], [1, 3]]);
let b = vec![expr!(5), expr!(7)];
let x = a.solve(&b).unwrap();

pub fn solve_least_squares( &self, b: &[Expression], ) -> Result<Vec<Expression>, MathError>

Solve least squares problem: min ||Ax - b||₂ using QR decomposition

Arguments

• b - Right-hand side vector

Returns

Solution vector x that minimizes ||Ax - b||₂

Errors

• DomainError if dimensions don’t match or m < n
• DivisionByZero if R has zero diagonal elements

Algorithm

For m×n matrix A (m >= n):

1. Compute A = QR (Q is m×n, R is n×n upper triangular)
2. Compute c = Q^T * b
3. Solve Rx = c[0:n] using backward substitution

Examples

use mathhook_core::matrices::Matrix;
use mathhook_core::expr;

// Overdetermined system: 3 equations, 2 unknowns
let a = Matrix::from_arrays([[1, 0], [0, 1], [1, 1]]);
let b = vec![expr!(1), expr!(2), expr!(2)];
let x = a.solve_least_squares(&b).unwrap();

Trait Implementations

impl Clone for Matrix

fn clone(&self) -> Matrix

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl CoreMatrixOps for Matrix

fn add(&self, other: &Matrix) -> Result<Matrix, MathError>

fn multiply(&self, other: &Matrix) -> Result<Matrix, MathError>

fn transpose(&self) -> Matrix

fn inverse(&self) -> Matrix

impl Debug for Matrix

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<'de> Deserialize<'de> for Matrix

fn deserialize<__D>(__deserializer: __D) -> Result<Self, __D::Error>

where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl EigenOperations for Matrix

Implementation of EigenOperations trait for Matrix

fn eigen_decomposition(&self) -> Option<EigenDecomposition>

Compute eigenvalues and eigenvectors

fn complex_eigen_decomposition(&self) -> Option<ComplexEigenDecomposition>

Compute complex eigenvalues and eigenvectors

fn eigenvalues(&self) -> Vec<Expression>

Compute only eigenvalues (faster than full decomposition)

fn characteristic_polynomial(&self) -> CharacteristicPolynomial

Compute characteristic polynomial det(A - λI)

fn trace(&self) -> Expression

Get the trace (sum of eigenvalues)

fn determinant_via_eigenvalues(&self) -> Expression

Get the determinant (product of eigenvalues)

fn is_diagonalizable(&self) -> bool

Check if matrix is diagonalizable

fn matrix_power_eigen(&self, n: i64) -> Option<Matrix>

Compute matrix power using eigendecomposition

fn matrix_exponential(&self) -> Option<Matrix>

Compute matrix exponential using eigendecomposition

fn matrix_logarithm(&self) -> Option<Matrix>

Compute matrix logarithm using eigendecomposition

fn matrix_sqrt(&self) -> Option<Matrix>

Compute matrix square root using eigendecomposition

fn is_nilpotent(&self) -> bool

Check if matrix is nilpotent

impl MatrixDecomposition for Matrix

Implementation of MatrixDecomposition trait for Matrix

fn lu_decomposition(&self) -> Option<LUDecomposition>

Perform LU decomposition with partial pivoting

fn qr_decomposition(&self) -> Option<QRDecomposition>

Perform QR decomposition using Gram-Schmidt process

fn cholesky_decomposition(&self) -> Option<CholeskyDecomposition>

Perform Cholesky decomposition for positive definite matrices

fn svd_decomposition(&self) -> Option<SVDDecomposition>

Perform Singular Value Decomposition

fn rank(&self) -> usize

Get matrix rank using SVD

fn is_positive_definite(&self) -> bool

Check if matrix is positive definite

fn condition_number(&self) -> Expression

Get condition number (ratio of largest to smallest singular value)

impl PartialEq for Matrix

fn eq(&self, other: &Matrix) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl Serialize for Matrix

fn serialize<__S>(&self, __serializer: __S) -> Result<__S::Ok, __S::Error>

where __S: Serializer,

Serialize this value into the given Serde serializer.

impl StructuralPartialEq for Matrix