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]

Example

let m = Matrix::new(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
assert_eq!(m.rows, 2);
assert_eq!(m.cols, 3);
assert_eq!(m.get(0, 0), 1.0);
assert_eq!(m.get(1, 2), 6.0);

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

Example

let m = Matrix::new(2, 2, vec![1.0, 2.0, 3.0, 4.0]);
assert_eq!(m.get(0, 0), 1.0);
assert_eq!(m.get(0, 1), 2.0);
assert_eq!(m.get(1, 0), 3.0);
assert_eq!(m.get(1, 1), 4.0);

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

Creates a matrix filled with zeros.

Arguments

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

Example

let m = Matrix::zeros(3, 2);
assert_eq!(m.rows, 3);
assert_eq!(m.cols, 2);
assert_eq!(m.get(1, 1), 0.0);

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].

Example

let m = Matrix::new(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
let t = m.transpose();
assert_eq!(t.rows, 3);
assert_eq!(t.cols, 2);
assert_eq!(t.get(0, 1), 4.0);

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

Performs matrix multiplication: self * other.

Panics

Panics if self.cols != other.rows.

Example

let a = Matrix::new(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
let b = Matrix::new(3, 2, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
let c = a.matmul(&b);
assert_eq!(c.rows, 2);
assert_eq!(c.cols, 2);
assert_eq!(c.get(0, 0), 22.0); // 1*1 + 2*3 + 3*5

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

Example

let m = Matrix::new(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]);
let v = vec![1.0, 2.0, 3.0];
let result = m.mul_vec(&v);
assert_eq!(result.len(), 2);
assert_eq!(result[0], 14.0); // 1*1 + 2*2 + 3*3

pub fn col_dot(&self, col: usize, v: &[f64]) -> f64

Computes the dot product of a column with a vector: Σ(data[i * cols + col] * v[i]).

For a row-major matrix, this iterates through all rows at a fixed column.

Arguments

• col - Column index
• v - Vector to dot with (must have length equal to rows)

Panics

Panics if col >= cols or v.len() != rows.

pub fn col_axpy_inplace(&self, col: usize, alpha: f64, v: &mut [f64])

Performs the column-vector operation in place: v += alpha * column_col.

This is the AXPY operation where the column is treated as a vector. For row-major storage, we iterate through rows at a fixed column.

Arguments

• col - Column index
• alpha - Scaling factor for the column
• v - Vector to modify in place (must have length equal to rows)

Panics

Panics if col >= cols or v.len() != rows.

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

Computes the squared L2 norm of a column: Σ(data[i * cols + col]²).

Arguments

• col - Column index

Panics

Panics if col >= cols.

pub fn add_diagonal_in_place(&mut self, alpha: f64, start_index: usize)

Adds a value to diagonal elements starting from a given index.

This is useful for ridge regression where we add lambda * I to X^T X, but the intercept column should not be penalized.

Arguments

• alpha - Value to add to diagonal elements
• start_index - Starting diagonal index (0 = first diagonal element)

Panics

Panics if the matrix is not square.

Example

For a 3×3 identity matrix with intercept in first column (unpenalized):

add_diagonal_in_place(lambda, 1) on:
[1.0, 0.0, 0.0]       [1.0,   0.0,   0.0  ]
[0.0, 1.0, 0.0]  ->   [0.0,  1.0+λ, 0.0  ]
[0.0, 0.0, 1.0]       [0.0,   0.0,  1.0+λ]

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)

Example

let i = Matrix::identity(3);
assert_eq!(i.get(0, 0), 1.0);
assert_eq!(i.get(1, 1), 1.0);
assert_eq!(i.get(2, 2), 1.0);
assert_eq!(i.get(0, 1), 0.0);

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.