Matrix Solver Integration Plan

Topic: internal.planning.matrix-solver-integration

Connect existing matrix decompositions (LU, QR, Cholesky) with solver infrastructure. Eliminates code duplication where SystemSolver.solve_nxn_system() reimplements Gaussian elimination despite Matrix.lu_decomposition() already existing.

Matrix Solver Integration Plan

Date: 2025-12-06T1530 Status: Ready for Execution Estimated Effort: 3.5 days

Executive Summary

Connect existing matrix decompositions (LU, QR, Cholesky) with solver infrastructure. Currently SystemSolver.solve_nxn_system() reimplements Gaussian elimination despite Matrix.lu_decomposition() already existing. This plan eliminates that duplication.

Phase 0: Forward/Backward Substitution (Foundation)

Priority: P0 - Required for all subsequent phases File: crates/mathhook-core/src/matrices/unified/operations.rs

Methods to Add

#![allow(unused)]
fn main() {
impl Matrix {
    /// Solve Lx = b for lower triangular L
    pub fn forward_substitution(&self, b: &[Expression]) -> Result<Vec<Expression>, MathError>;

    /// Solve Ux = b for upper triangular U
    pub fn backward_substitution(&self, b: &[Expression]) -> Result<Vec<Expression>, MathError>;
}
}

Algorithms

Forward substitution (Ly = b):

For i = 0 to n-1:
    sum = Σ(L[i][j] * x[j]) for j = 0 to i-1
    x[i] = (b[i] - sum) / L[i][i]
    if L[i][i] == 0: return Err(DivisionByZero)

Backward substitution (Ux = y):

For i = n-1 down to 0:
    sum = Σ(U[i][j] * x[j]) for j = i+1 to n-1
    x[i] = (b[i] - sum) / U[i][i]
    if U[i][i] == 0: return Err(DivisionByZero)

Tests Required

• Identity matrix (trivial)
• 2x2, 3x3 triangular matrices
• Zero diagonal element → DivisionByZero
• Symbolic coefficients

Phase 1: Matrix.solve(b) Using LU

Priority: P0 - Core integration File: crates/mathhook-core/src/matrices/unified/operations.rs

Method to Add

#![allow(unused)]
fn main() {
impl Matrix {
    /// Solve Ax = b using LU decomposition
    ///
    /// Algorithm:
    /// 1. Compute PA = LU
    /// 2. Solve Ly = Pb (forward substitution)
    /// 3. Solve Ux = y (backward substitution)
    pub fn solve(&self, b: &[Expression]) -> Result<Vec<Expression>, MathError>;
}
}

Helper Needed

#![allow(unused)]
fn main() {
fn apply_permutation(p: &Option<Matrix>, b: &[Expression]) -> Vec<Expression>
}

Error Cases

• Non-square matrix → DomainError
• Dimension mismatch → DomainError
• Singular matrix → DivisionByZero

Tests Required

• 2x2 integer solution
• 3x3 rational solution
• Identity matrix (trivial)
• Singular matrix → error
• Dimension mismatch → error

Phase 2: Replace SystemSolver.solve_nxn_system()

Priority: P0 - Eliminate duplication File: crates/mathhook-core/src/algebra/solvers/systems.rs

Current State (lines 407-518)

• Manually extracts coefficients into augmented matrix
• Reimplements Gaussian elimination with partial pivoting
• Performs back substitution inline

Refactored Approach

#![allow(unused)]
fn main() {
fn solve_nxn_system(&self, equations: &[Expression], variables: &[Symbol]) -> SolverResult {
    // 1. Build coefficient matrix A and RHS vector b
    let (a_matrix, b_vec) = self.build_system_matrix(equations, variables);

    // 2. Use Matrix.solve() which uses LU decomposition
    match a_matrix.solve(&b_vec) {
        Ok(solution) => SolverResult::Multiple(solution),
        Err(MathError::DivisionByZero) => self.check_singularity_type(&a_matrix, &b_vec),
        Err(_) => SolverResult::NoSolution,
    }
}
}

Lines to Delete

• Lines 426-479 (Gaussian elimination reimplementation)
• Lines 480-518 (back substitution reimplementation)

Import to Add

#![allow(unused)]
fn main() {
use crate::matrices::Matrix;
}

Phase 3: SPD Detection → Cholesky Routing

Priority: P1 - Performance optimization File: crates/mathhook-core/src/matrices/unified/operations.rs

Enhanced solve() Method

#![allow(unused)]
fn main() {
pub fn solve(&self, b: &[Expression]) -> Result<Vec<Expression>, MathError> {
    // Try Cholesky first for SPD matrices (2x faster)
    if self.is_symmetric() {
        if let Some(chol) = self.cholesky_decomposition() {
            let y = chol.l.forward_substitution(b)?;
            let lt = chol.l.transpose();
            return lt.backward_substitution(&y);
        }
    }

    // Fall back to LU
    self.solve_via_lu(b)
}
}

Expected Performance

• 2x speedup for symmetric positive definite matrices

Phase 4: QR-Based Least Squares

Priority: P1 - Overdetermined systems File: crates/mathhook-core/src/matrices/unified/operations.rs

Method to Add

#![allow(unused)]
fn main() {
impl Matrix {
    /// Solve min ||Ax - b||_2 using QR decomposition
    ///
    /// For m×n matrix (m >= n):
    /// 1. Compute A = QR
    /// 2. Compute Q^T * b
    /// 3. Solve Rx = (Q^T * b)[0:n]
    pub fn solve_least_squares(&self, b: &[Expression]) -> Result<Vec<Expression>, MathError>;
}
}

Helper Needed

#![allow(unused)]
fn main() {
fn matrix_vector_multiply(m: &Matrix, v: &[Expression]) -> Vec<Expression>
}

Phase 5: Matrix Inversion Using LU

Priority: P2 - Optimization File: crates/mathhook-core/src/matrices/unified/operations.rs

Enhanced inverse() Method

#![allow(unused)]
fn main() {
fn inverse(&self) -> Matrix {
    match self {
        Matrix::Identity(_) | Matrix::Scalar(_) | Matrix::Diagonal(_) => {
            // Keep existing fast paths
        }
        _ => self.inverse_via_lu().unwrap_or_else(|| self.gauss_jordan_inverse()),
    }
}

fn inverse_via_lu(&self) -> Option<Matrix> {
    // Solve A * X = I column by column
    for j in 0..n {
        let e_j = unit_vector(j, n);
        let col = self.solve(&e_j)?;
        // Store in result matrix
    }
}
}

Dependency Order

Phase 0: Forward/Backward Substitution
    ↓
Phase 1: Matrix.solve(b)
    ↓
Phase 2: Replace solve_nxn_system ←── Phase 3: Cholesky Routing
    ↓
Phase 4: QR Least Squares
    ↓
Phase 5: LU-based Inverse

Verification Strategy

SymPy Validation

# Create validation script at scripts/validate_matrix_solver.py
from sympy import Matrix

A = Matrix([[2, 1], [1, 3]])
b = Matrix([5, 7])
x = A.solve(b)  # Compare with Rust output

Rust Tests

#![allow(unused)]
fn main() {
#[test]
fn test_solve_via_lu_correctness() {
    let a = Matrix::from_arrays([[2, 1], [1, 3]]);
    let b = vec![expr!(5), expr!(7)];
    let x = a.solve(&b).unwrap();

    // Verify: A * x == b
    let ax = matrix_vector_multiply(&a, &x);
    assert_eq!(ax[0].simplify(), b[0].simplify());
    assert_eq!(ax[1].simplify(), b[1].simplify());
}
}

Risk Mitigation

Files Summary

Net change: ~+80 lines, elimination of code duplication

Examples

API Reference

• Rust: ``
• Python: ``
• JavaScript: ``

See Also

• internal.planning.fast-path-implementation

• internal.planning.fast-path-opportunities