From f7498e29ccf16a260f40fd9b7913be97b9d1a6fa Mon Sep 17 00:00:00 2001 From: vmarcella Date: Thu, 12 Mar 2026 13:02:34 -0700 Subject: [PATCH 1/4] [update] smaller arrays to use gaussian elimination for computing the determinant --- crates/lambda-rs/src/math/matrix.rs | 166 +++++++++++++++++++++++----- 1 file changed, 139 insertions(+), 27 deletions(-) diff --git a/crates/lambda-rs/src/math/matrix.rs b/crates/lambda-rs/src/math/matrix.rs index b2cd7590..7263f366 100644 --- a/crates/lambda-rs/src/math/matrix.rs +++ b/crates/lambda-rs/src/math/matrix.rs @@ -259,6 +259,62 @@ pub fn identity_matrix< return result; } +/// Compute determinant via Gaussian elimination on fixed-size stack storage. +/// +/// This avoids heap allocations in hot fixed-size paths (for example 4x4 +/// transform matrices). +/// +/// # Arguments +/// +/// - `data`: Square matrix values stored in stack-allocated fixed-size arrays. +/// +/// # Returns +/// +/// - Determinant of `data`. +/// - `0.0` when elimination detects a near-zero pivot (singular matrix). +fn determinant_gaussian_stack(mut data: [[f32; N]; N]) -> f32 { + let mut sign = 1.0_f32; + for pivot in 0..N { + let mut pivot_row = pivot; + let mut pivot_abs = data[pivot][pivot].abs(); + for (candidate, row) in data.iter().enumerate().skip(pivot + 1) { + let candidate_abs = row[pivot].abs(); + if candidate_abs > pivot_abs { + pivot_abs = candidate_abs; + pivot_row = candidate; + } + } + + if pivot_abs <= f32::EPSILON { + return 0.0; + } + + if pivot_row != pivot { + data.swap(pivot, pivot_row); + sign = -sign; + } + + let pivot_value = data[pivot][pivot]; + let pivot_tail = data[pivot]; + for row in data.iter_mut().skip(pivot + 1) { + let factor = row[pivot] / pivot_value; + if factor == 0.0 { + continue; + } + for (dst, src) in row[pivot..].iter_mut().zip(pivot_tail[pivot..].iter()) + { + *dst -= factor * *src; + } + } + } + + let mut det = sign; + for (i, row) in data.iter().enumerate().take(N) { + det *= row[i]; + } + return det; +} + // -------------------------- ARRAY IMPLEMENTATION ----------------------------- /// Matrix implementations for arrays of f32 arrays. Including the trait Matrix into @@ -325,7 +381,15 @@ where todo!() } - /// Computes the determinant of any square matrix using Laplace expansion. + /// Computes the determinant of any square matrix using Gaussian elimination + /// with partial pivoting. + /// + /// # Behavior + /// + /// - Uses stack-allocated elimination for 3x3 and 4x4 matrices to avoid heap + /// allocation on common transform sizes. + /// - Uses a generic dense fallback for larger matrices. + /// - Overall asymptotic runtime is `O(n^3)` for square `n x n` matrices. fn determinant(&self) -> Result { let rows = self.as_ref().len(); if rows == 0 { @@ -341,34 +405,83 @@ where return Err(MathError::NonSquareMatrix { rows, cols }); } - return match rows { - 1 => Ok(self.as_ref()[0].as_ref()[0]), - 2 => { - let a = self.at(0, 0); - let b = self.at(0, 1); - let c = self.at(1, 0); - let d = self.at(1, 1); - return Ok(a * d - b * c); + if rows == 1 { + return Ok(self.as_ref()[0].as_ref()[0]); + } + + if rows == 2 { + let a = self.at(0, 0); + let b = self.at(0, 1); + let c = self.at(1, 0); + let d = self.at(1, 1); + return Ok(a * d - b * c); + } + + if rows == 3 { + let mut data = [[0.0_f32; 3]; 3]; + for (i, row) in data.iter_mut().enumerate() { + for (j, value) in row.iter_mut().enumerate() { + *value = self.at(i, j); + } + } + return Ok(determinant_gaussian_stack::<3>(data)); + } + + if rows == 4 { + let mut data = [[0.0_f32; 4]; 4]; + for (i, row) in data.iter_mut().enumerate() { + for (j, value) in row.iter_mut().enumerate() { + *value = self.at(i, j); + } + } + return Ok(determinant_gaussian_stack::<4>(data)); + } + + // Convert to a mutable dense matrix so we can perform in-place elimination. + let mut working = vec![vec![0.0_f32; rows]; rows]; + for (i, row) in working.iter_mut().enumerate() { + for (j, value) in row.iter_mut().enumerate() { + *value = self.at(i, j); } - _ => { - let mut result = 0.0; - for i in 0..rows { - let mut submatrix: Vec> = Vec::with_capacity(rows - 1); - for j in 1..rows { - let mut row = Vec::new(); - for k in 0..rows { - if k != i { - row.push(self.at(j, k)); - } - } - submatrix.push(row); - } - let sub_determinant = submatrix.determinant()?; - result += self.at(0, i) * sub_determinant * (-1.0_f32).powi(i as i32); + } + + let mut sign = 1.0_f32; + for pivot in 0..rows { + // Partial pivoting improves numerical stability and avoids division by 0. + let mut pivot_row = pivot; + let mut pivot_abs = working[pivot][pivot].abs(); + for (candidate, row) in working.iter().enumerate().skip(pivot + 1) { + let candidate_abs = row[pivot].abs(); + if candidate_abs > pivot_abs { + pivot_abs = candidate_abs; + pivot_row = candidate; } - return Ok(result); } - }; + + if pivot_abs <= f32::EPSILON { + return Ok(0.0); + } + + if pivot_row != pivot { + working.swap(pivot, pivot_row); + sign = -sign; + } + + let pivot_value = working[pivot][pivot]; + let pivot_tail: Vec = working[pivot][pivot..].to_vec(); + for row in working.iter_mut().skip(pivot + 1) { + let factor = row[pivot] / pivot_value; + if factor == 0.0 { + continue; + } + for (dst, src) in row[pivot..].iter_mut().zip(pivot_tail.iter()) { + *dst -= factor * *src; + } + } + } + + let diagonal_product = (0..rows).map(|i| working[i][i]).product::(); + return Ok(sign * diagonal_product); } /// Return the size as a (rows, columns). @@ -394,7 +507,6 @@ where #[cfg(test)] mod tests { - use super::{ filled_matrix, identity_matrix, From 1272f9933c153ceaa9bb06994a502c80f68b4e15 Mon Sep 17 00:00:00 2001 From: vmarcella Date: Thu, 12 Mar 2026 13:05:04 -0700 Subject: [PATCH 2/4] [update] implementation for larger matrices to use Vec instead of a Vec of Vecs --- crates/lambda-rs/src/math/matrix.rs | 37 ++++++++++++++++------------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/crates/lambda-rs/src/math/matrix.rs b/crates/lambda-rs/src/math/matrix.rs index 7263f366..788f34f9 100644 --- a/crates/lambda-rs/src/math/matrix.rs +++ b/crates/lambda-rs/src/math/matrix.rs @@ -437,11 +437,12 @@ where return Ok(determinant_gaussian_stack::<4>(data)); } - // Convert to a mutable dense matrix so we can perform in-place elimination. - let mut working = vec![vec![0.0_f32; rows]; rows]; - for (i, row) in working.iter_mut().enumerate() { - for (j, value) in row.iter_mut().enumerate() { - *value = self.at(i, j); + // Use a contiguous row-major buffer for larger matrices to keep the + // fallback cache-friendlier than `Vec>`. + let mut working = vec![0.0_f32; rows * rows]; + for i in 0..rows { + for j in 0..rows { + working[i * rows + j] = self.at(i, j); } } @@ -449,9 +450,9 @@ where for pivot in 0..rows { // Partial pivoting improves numerical stability and avoids division by 0. let mut pivot_row = pivot; - let mut pivot_abs = working[pivot][pivot].abs(); - for (candidate, row) in working.iter().enumerate().skip(pivot + 1) { - let candidate_abs = row[pivot].abs(); + let mut pivot_abs = working[pivot * rows + pivot].abs(); + for candidate in (pivot + 1)..rows { + let candidate_abs = working[candidate * rows + pivot].abs(); if candidate_abs > pivot_abs { pivot_abs = candidate_abs; pivot_row = candidate; @@ -463,24 +464,28 @@ where } if pivot_row != pivot { - working.swap(pivot, pivot_row); + for column in 0..rows { + working.swap(pivot * rows + column, pivot_row * rows + column); + } sign = -sign; } - let pivot_value = working[pivot][pivot]; - let pivot_tail: Vec = working[pivot][pivot..].to_vec(); - for row in working.iter_mut().skip(pivot + 1) { - let factor = row[pivot] / pivot_value; + let pivot_value = working[pivot * rows + pivot]; + for row in (pivot + 1)..rows { + let factor = working[row * rows + pivot] / pivot_value; if factor == 0.0 { continue; } - for (dst, src) in row[pivot..].iter_mut().zip(pivot_tail.iter()) { - *dst -= factor * *src; + for column in pivot..rows { + let row_idx = row * rows + column; + let pivot_idx = pivot * rows + column; + working[row_idx] -= factor * working[pivot_idx]; } } } - let diagonal_product = (0..rows).map(|i| working[i][i]).product::(); + let diagonal_product = + (0..rows).map(|i| working[i * rows + i]).product::(); return Ok(sign * diagonal_product); } From 02c5cd556ab1809a2290ea1ae0e4a6cbeeedf194 Mon Sep 17 00:00:00 2001 From: vmarcella Date: Thu, 12 Mar 2026 13:32:02 -0700 Subject: [PATCH 3/4] [fix] pivot checks so that non singular matrices aren't incorrectly marked as singular --- crates/lambda-rs/src/math/matrix.rs | 33 ++++++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/crates/lambda-rs/src/math/matrix.rs b/crates/lambda-rs/src/math/matrix.rs index 788f34f9..708b35d0 100644 --- a/crates/lambda-rs/src/math/matrix.rs +++ b/crates/lambda-rs/src/math/matrix.rs @@ -271,7 +271,8 @@ pub fn identity_matrix< /// # Returns /// /// - Determinant of `data`. -/// - `0.0` when elimination detects a near-zero pivot (singular matrix). +/// - `0.0` when elimination detects an exact zero pivot after pivoting +/// (singular matrix). fn determinant_gaussian_stack(mut data: [[f32; N]; N]) -> f32 { let mut sign = 1.0_f32; for pivot in 0..N { @@ -285,7 +286,7 @@ fn determinant_gaussian_stack(mut data: [[f32; N]; N]) -> f32 { } } - if pivot_abs <= f32::EPSILON { + if pivot_abs == 0.0 { return 0.0; } @@ -459,7 +460,7 @@ where } } - if pivot_abs <= f32::EPSILON { + if pivot_abs == 0.0 { return Ok(0.0); } @@ -571,6 +572,32 @@ mod tests { assert_eq!(m2.determinant(), Ok(-306.0)); } + #[test] + fn determinant_preserves_small_non_zero_pivot_in_stack_path() { + let m = [[1.0e-8, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]]; + crate::assert_approximately_equal!( + m.determinant().expect("determinant should succeed"), + 6.0e-8, + 1.0e-12 + ); + } + + #[test] + fn determinant_preserves_small_non_zero_pivot_in_dynamic_path() { + let m = [ + [1.0e-8, 0.0, 0.0, 0.0, 0.0], + [0.0, 2.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 3.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 4.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 5.0], + ]; + crate::assert_approximately_equal!( + m.determinant().expect("determinant should succeed"), + 1.2e-6, + 1.0e-10 + ); + } + #[test] fn non_square_matrix_determinant() { let m = [[3.0, 8.0], [4.0, 6.0], [0.0, 1.0]]; From 9bdf71e9829fe3e0793d9faec3c0a23a4a6313e2 Mon Sep 17 00:00:00 2001 From: vmarcella Date: Thu, 12 Mar 2026 13:34:37 -0700 Subject: [PATCH 4/4] [add] tests for 4x4 and 5x5 row swaps --- crates/lambda-rs/src/math/matrix.rs | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/crates/lambda-rs/src/math/matrix.rs b/crates/lambda-rs/src/math/matrix.rs index 708b35d0..683c7767 100644 --- a/crates/lambda-rs/src/math/matrix.rs +++ b/crates/lambda-rs/src/math/matrix.rs @@ -572,6 +572,29 @@ mod tests { assert_eq!(m2.determinant(), Ok(-306.0)); } + #[test] + fn determinant_4x4_stack_path_handles_row_swap_sign() { + let m = [ + [0.0, 2.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 3.0, 0.0], + [0.0, 0.0, 0.0, 4.0], + ]; + assert_eq!(m.determinant(), Ok(-24.0)); + } + + #[test] + fn determinant_5x5_dynamic_path_handles_row_swap_sign() { + let m = [ + [0.0, 2.0, 0.0, 0.0, 0.0], + [1.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 3.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 4.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 5.0], + ]; + assert_eq!(m.determinant(), Ok(-120.0)); + } + #[test] fn determinant_preserves_small_non_zero_pivot_in_stack_path() { let m = [[1.0e-8, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 3.0]];