Struct Matrix 

pub struct Matrix {
    pub rows: usize,
    pub cols: usize,
    pub data: Vec<f64>,
}

A dense matrix stored in row-major order.

Storage

Elements are stored in a single flat vector in row-major order: data[row * cols + col]

Fields

rows: usize

Number of rows in the matrix

cols: usize

Number of columns in the matrix

data: Vec<f64>

Flat vector storing matrix elements in row-major order

Implementations

impl Matrix

pub fn new(rows: usize, cols: usize, data: Vec<f64>) -> Self

Creates a new matrix from the given dimensions and data.

Panics

Panics if data.len() != rows * cols.

Arguments

• rows - Number of rows
• cols - Number of columns
• data - Flat vector of elements in row-major order

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

Creates a matrix filled with zeros.

Arguments

• rows - Number of rows
• cols - Number of columns

pub fn get(&self, row: usize, col: usize) -> f64

Gets the element at the specified row and column.

Arguments

• row - Row index (0-based)
• col - Column index (0-based)

pub fn set(&mut self, row: usize, col: usize, val: f64)

Sets the element at the specified row and column.

Arguments

• row - Row index (0-based)
• col - Column index (0-based)
• val - Value to set

pub fn transpose(&self) -> Matrix

Returns the transpose of this matrix.

Swaps rows with columns: result\[col\]\[row\] = self\[row\]\[col\].

pub fn matmul(&self, other: &Matrix) -> Matrix

Performs matrix multiplication: self * other.

Panics

Panics if self.cols != other.rows.

pub fn mul_vec(&self, vec: &[f64]) -> Vec<f64>

Multiplies this matrix by a vector (treating vector as column matrix).

Computes self * vec where vec is treated as an n×1 column matrix.

Panics

Panics if self.cols != vec.len().

Arguments

• vec - Vector to multiply by

impl Matrix

pub fn qr(&self) -> (Matrix, Matrix)

Computes the QR decomposition using Householder reflections.

Factorizes the matrix as A = QR where Q is orthogonal and R is upper triangular.

Requirements

This implementation requires rows >= cols (tall matrix). For OLS regression, we always have more observations than predictors, so this requirement is satisfied.

Returns

A tuple (Q, R) where:

• Q is an orthogonal matrix (QᵀQ = I) of size m×m
• R is an upper triangular matrix of size m×n

pub fn identity(size: usize) -> Self

Creates an identity matrix of the given size.

Arguments

• size - Number of rows and columns (square matrix)

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

Inverts an upper triangular matrix (such as R from QR decomposition).

Uses back-substitution to compute the inverse. This is efficient for triangular matrices compared to general matrix inversion.

Panics

Panics if the matrix is not square.

Returns

None if the matrix is singular (has a zero or near-zero diagonal element). A matrix is considered singular if any diagonal element is below the internal tolerance (1e-10), which indicates the matrix does not have full rank.

Note

For upper triangular matrices, singularity is equivalent to having a zero (or near-zero) diagonal element. This is much simpler to check than for general matrices, which would require computing the condition number.

pub fn invert_upper_triangular_with_tolerance( &self, tolerance_mult: f64, ) -> Option<Matrix>

Inverts an upper triangular matrix with a custom tolerance multiplier.

The tolerance is computed as max_diag * epsilon * n * tolerance_mult. A higher tolerance_mult allows more tolerance for near-singular matrices.

Arguments

• tolerance_mult - Multiplier for the tolerance (1.0 = standard, higher = more tolerant)

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

Computes the inverse of a square matrix using QR decomposition.

For an invertible matrix A, computes A⁻¹ such that A * A⁻¹ = I. Uses QR decomposition for numerical stability.

Panics

Panics if the matrix is not square (i.e., self.rows != self.cols). Check dimensions before calling if the matrix shape is not guaranteed.

Returns

Returns Some(inverse) if the matrix is invertible, or None if the matrix is singular (non-invertible).

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

Computes the inverse of X’X given the QR decomposition of X (R’s chol2inv).

This is equivalent to computing (X'X)^(-1) using the QR decomposition of X. R’s chol2inv function is used for numerical stability in recursive residuals.

Arguments

• x - Input matrix (must have rows >= cols)

Returns

Some((X'X)^(-1)) if X has full rank, None otherwise.

Algorithm

Given QR decomposition X = QR where R is upper triangular:

1. Extract the upper p×p portion of R (denoted R₁)
2. Invert R₁ (upper triangular inverse)
3. Compute (X’X)^(-1) = R₁^(-1) × R₁^(-T)

This works because X’X = R’Q’QR = R’R, and R₁ contains the Cholesky factor.

pub fn chol2inv_from_qr_with_tolerance( &self, tolerance_mult: f64, ) -> Option<Matrix>

Computes the inverse of X’X given the QR decomposition with custom tolerance.

Similar to chol2inv_from_qr but allows specifying a tolerance multiplier for handling near-singular matrices.

Arguments

• tolerance_mult - Multiplier for the tolerance (higher = more tolerant)

impl Matrix

pub fn qr_linpack(&self, tol: Option<f64>) -> QRLinpack

Computes QR decomposition using R’s LINPACK dqrdc2 algorithm with column pivoting.

This is a port of R’s dqrdc2.f, which is a modification of LINPACK’s DQRDC. The algorithm:

1. Uses Householder transformations for QR factorization
2. Implements limited column pivoting based on column 2-norms
3. Moves columns with near-zero norm to the right-hand edge
4. Computes the rank (number of linearly independent columns)

Arguments

• tol - Tolerance for determining linear independence. Default is 1e-7 (R’s default). Columns with norm < tol * original_norm are considered negligible.

Returns

A QRLinpack struct containing the QR factorization, auxiliary information, pivot vector, and computed rank.

Algorithm Details

The decomposition is A * P = Q * R where:

• P is the permutation matrix coded by pivot
• Q is orthogonal (m × m)
• R is upper triangular in the first rank rows

The qr matrix contains:

• Upper triangle: R matrix (if pivoting was performed, this is R of the permuted matrix)
• Below diagonal: Householder vector information

Reference

• R source: src/appl/dqrdc2.f
• LINPACK documentation: https://www.netlib.org/linpack/dqrdc.f

pub fn qr_solve_linpack( &self, qr_result: &QRLinpack, y: &[f64], ) -> Option<Vec<f64>>

Solves a linear system using the QR decomposition with column pivoting.

This implements a least squares solver using the pivoted QR decomposition. For rank-deficient cases, coefficients corresponding to linearly dependent columns are set to f64::NAN.

Arguments

• qr_result - QR decomposition from Matrix::qr_linpack
• y - Right-hand side vector

Returns

A vector of coefficients, or None if the system is exactly singular.

Algorithm

This solver uses the standard QR decomposition approach:

1. Compute the QR decomposition of the permuted matrix
2. Extract R matrix (upper triangular with positive diagonal)
3. Compute qty = Q^T * y
4. Solve R * coef = qty using back substitution
5. Apply the pivot permutation to restore original column order

Note

The LINPACK QR algorithm stores R with mixed signs on the diagonal. This solver corrects for that by taking the absolute value of R’s diagonal.

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 Debug for Matrix

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

Formats the value using the given formatter.