From e9e0e055c8a97e575f5d4f8fb73c597abd6434fa Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:05:28 -0400 Subject: [PATCH 01/18] gram schmidt --- hoomd-linear-algebra/src/matrix.rs | 3 +++ hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 1 + 2 files changed, 4 insertions(+) create mode 100644 hoomd-linear-algebra/src/matrix/gram_schmidt.rs diff --git a/hoomd-linear-algebra/src/matrix.rs b/hoomd-linear-algebra/src/matrix.rs index e5ecccae3..27b372bb5 100644 --- a/hoomd-linear-algebra/src/matrix.rs +++ b/hoomd-linear-algebra/src/matrix.rs @@ -12,6 +12,9 @@ pub mod ops; /// ``qr`` decomposition for [`Matrix`] types. pub mod qr; +/// ``Gram-Schmidt`` orthonormalization for [`Matrix`] types. +pub mod gram_schmidt; + pub use crate::diagonal::DiagonalMatrix; /// A lightweight representation of a diagonal matrix. diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -0,0 +1 @@ + From 42ab16f68e890964573cad1c0c219c6f0cdcba25 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:21:57 -0400 Subject: [PATCH 02/18] WIP on implementing, needs 2 references --- .../src/matrix/gram_schmidt.rs | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 8b1378917..301ead4b1 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -1 +1,38 @@ +use crate::matrix::Matrix; +/// +/// +/// Implementation based on https://www.sfu.ca/~jtmulhol/py4math/linalg/np-gramschmidt/ +#[inline] +pub fn gram_schmidt(a: &Matrix) -> Matrix { + let mut a = a.clone(); + for j in 0..a.n_columns() { + // For the vector in column k, find the perpendicular of the projection onto + // the previous orthogonal vectors. + for k in (0..j) { + let j_dot_k = a + .get_col_slice_iter(k, 0..N) + .zip(a.get_col_slice_iter(j, 0..N)) + .fold(0.0, |acc, (l, r)| acc + (l * r)); + let proj_j_onto_k = a.get_col_slice_iter(k, 0..N).map(|x| x * j_dot_k); + a.get_col_slice_iter_mut(j, 0..N) + .zip(proj_j_onto_k) + .map(|(a_ji, proj_i)| *a_ji -= proj_i); + } + } + a +} + +// def gram_schmidt(A): +// for j in range(n): +// # For the vector in column j, find the perpendicular +// # of the projection onto the previous orthogonal vectors. +// for k in range(j): +// A[:, j] -= np.dot(A[:, k], A[:, j]) * A[:, k] +// # If original vectors aren't lin indep then we can check for this: +// # +// if np.isclose(np.linalg.norm(A[:, j]), 0, rtol=1e-15, atol=1e-14, equal_nan=False): +// A[:, j] = np.zeros(A.shape[0]) +// else: +// A[:, j] = A[:, j] / np.linalg.norm(A[:, j]) +// return A From e76c466a670c0059c20e46d9902822ad4dd7658f Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:25:51 -0400 Subject: [PATCH 03/18] I just need two references --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 301ead4b1..ff1cfabf3 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -14,10 +14,10 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< .get_col_slice_iter(k, 0..N) .zip(a.get_col_slice_iter(j, 0..N)) .fold(0.0, |acc, (l, r)| acc + (l * r)); - let proj_j_onto_k = a.get_col_slice_iter(k, 0..N).map(|x| x * j_dot_k); - a.get_col_slice_iter_mut(j, 0..N) - .zip(proj_j_onto_k) - .map(|(a_ji, proj_i)| *a_ji -= proj_i); + let mut proj_j_onto_k = a.get_col_slice_iter(k, 0..N).map(|x| x * j_dot_k); + for i in 0..N { + a[(i, j)] -= proj_j_onto_k.next().unwrap(); + } } } a From 28ddf05e4fd2d68333f21ee298a7d3d84c766e2a Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:31:50 -0400 Subject: [PATCH 04/18] wip --- .../src/matrix/gram_schmidt.rs | 21 +++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index ff1cfabf3..5ebb057cb 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -2,7 +2,8 @@ use crate::matrix::Matrix; /// /// -/// Implementation based on https://www.sfu.ca/~jtmulhol/py4math/linalg/np-gramschmidt/ +/// Implementation based on +#[must_use] #[inline] pub fn gram_schmidt(a: &Matrix) -> Matrix { let mut a = a.clone(); @@ -14,11 +15,23 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< .get_col_slice_iter(k, 0..N) .zip(a.get_col_slice_iter(j, 0..N)) .fold(0.0, |acc, (l, r)| acc + (l * r)); - let mut proj_j_onto_k = a.get_col_slice_iter(k, 0..N).map(|x| x * j_dot_k); for i in 0..N { - a[(i, j)] -= proj_j_onto_k.next().unwrap(); + a[(i, j)] -= a[(k, i)] * j_dot_k; } - } +// # If original vectors aren't lin indep then we can check for this: +// # +// if np.isclose(np.linalg.norm(A[:, j]), 0, rtol=1e-15, atol=1e-14, equal_nan=False): +// A[:, j] = np.zeros(A.shape[0]) +// else: +// A[:, j] = A[:, j] / np.linalg.norm(A[:, j]) + let column_j_norm = a + .get_col(j).iter_elements().fold(0.0, |acc, x| acc + x*x).sqrt(); + if column_j_norm.is_finite() { + a.get_col_slice_iter_mut(j, 0..N).for_each(|x|*x /= column_j_norm); + } + else { + + } } a } From 37110a5c05c6c158d59a143c5d1f80ddadf7124e Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:33:57 -0400 Subject: [PATCH 05/18] fix typo and bug --- .../src/matrix/gram_schmidt.rs | 27 ++++++++++--------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 5ebb057cb..71aa56338 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -15,23 +15,26 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< .get_col_slice_iter(k, 0..N) .zip(a.get_col_slice_iter(j, 0..N)) .fold(0.0, |acc, (l, r)| acc + (l * r)); + // Apply the projection for i in 0..N { - a[(i, j)] -= a[(k, i)] * j_dot_k; + a[(i, j)] -= a[(i, k)] * j_dot_k; } -// # If original vectors aren't lin indep then we can check for this: -// # -// if np.isclose(np.linalg.norm(A[:, j]), 0, rtol=1e-15, atol=1e-14, equal_nan=False): -// A[:, j] = np.zeros(A.shape[0]) -// else: -// A[:, j] = A[:, j] / np.linalg.norm(A[:, j]) + // # If original vectors aren't lin indep then we can check for this: + // # + // if np.isclose(np.linalg.norm(A[:, j]), 0, rtol=1e-15, atol=1e-14, equal_nan=False): + // A[:, j] = np.zeros(A.shape[0]) + // else: + // A[:, j] = A[:, j] / np.linalg.norm(A[:, j]) let column_j_norm = a - .get_col(j).iter_elements().fold(0.0, |acc, x| acc + x*x).sqrt(); + .get_col(j) + .iter_elements() + .fold(0.0, |acc, x| acc + x * x) + .sqrt(); if column_j_norm.is_finite() { - a.get_col_slice_iter_mut(j, 0..N).for_each(|x|*x /= column_j_norm); + a.get_col_slice_iter_mut(j, 0..N) + .for_each(|x| *x /= column_j_norm); } - else { - - } + } } a } From bb0244da989f97182e59507e02b4dab8c4305722 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:37:15 -0400 Subject: [PATCH 06/18] Working! Correct? --- .../src/matrix/gram_schmidt.rs | 31 +++++-------------- 1 file changed, 8 insertions(+), 23 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 71aa56338..cc6df0966 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -1,6 +1,6 @@ use crate::matrix::Matrix; -/// +/// Construct an orthonormal basis from the vectors in a [`Matrix`]. /// /// Implementation based on #[must_use] @@ -12,43 +12,28 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< // the previous orthogonal vectors. for k in (0..j) { let j_dot_k = a - .get_col_slice_iter(k, 0..N) - .zip(a.get_col_slice_iter(j, 0..N)) + .get_col(k) + .iter_elements() + .zip(a.get_col(j).iter_elements()) .fold(0.0, |acc, (l, r)| acc + (l * r)); // Apply the projection for i in 0..N { a[(i, j)] -= a[(i, k)] * j_dot_k; } - // # If original vectors aren't lin indep then we can check for this: - // # - // if np.isclose(np.linalg.norm(A[:, j]), 0, rtol=1e-15, atol=1e-14, equal_nan=False): - // A[:, j] = np.zeros(A.shape[0]) - // else: - // A[:, j] = A[:, j] / np.linalg.norm(A[:, j]) let column_j_norm = a .get_col(j) .iter_elements() .fold(0.0, |acc, x| acc + x * x) .sqrt(); + + // If the initial vectors are not linearly independent, zero out the col. if column_j_norm.is_finite() { a.get_col_slice_iter_mut(j, 0..N) .for_each(|x| *x /= column_j_norm); + } else { + a.get_col_slice_iter_mut(j, 0..N).for_each(|x| *x = 0.0); } } } a } - -// def gram_schmidt(A): -// for j in range(n): -// # For the vector in column j, find the perpendicular -// # of the projection onto the previous orthogonal vectors. -// for k in range(j): -// A[:, j] -= np.dot(A[:, k], A[:, j]) * A[:, k] -// # If original vectors aren't lin indep then we can check for this: -// # -// if np.isclose(np.linalg.norm(A[:, j]), 0, rtol=1e-15, atol=1e-14, equal_nan=False): -// A[:, j] = np.zeros(A.shape[0]) -// else: -// A[:, j] = A[:, j] / np.linalg.norm(A[:, j]) -// return A From 9c605db589179ca01eb3f9f04a13b2d337601410 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:37:33 -0400 Subject: [PATCH 07/18] lint --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index cc6df0966..b8b81f384 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -10,7 +10,7 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< for j in 0..a.n_columns() { // For the vector in column k, find the perpendicular of the projection onto // the previous orthogonal vectors. - for k in (0..j) { + for k in 0..j { let j_dot_k = a .get_col(k) .iter_elements() From a468a62db5c680f72f0221167f81194bf5d2db4b Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:39:25 -0400 Subject: [PATCH 08/18] move normalization out of loop --- .../src/matrix/gram_schmidt.rs | 27 +++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index b8b81f384..1afe90a60 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -20,20 +20,19 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< for i in 0..N { a[(i, j)] -= a[(i, k)] * j_dot_k; } - let column_j_norm = a - .get_col(j) - .iter_elements() - .fold(0.0, |acc, x| acc + x * x) - .sqrt(); - - // If the initial vectors are not linearly independent, zero out the col. - if column_j_norm.is_finite() { - a.get_col_slice_iter_mut(j, 0..N) - .for_each(|x| *x /= column_j_norm); - } else { - a.get_col_slice_iter_mut(j, 0..N).for_each(|x| *x = 0.0); - } + } // end loop over k + let column_j_norm = a + .get_col(j) + .iter_elements() + .fold(0.0, |acc, x| acc + x * x) + .sqrt(); + // If the initial vectors are not linearly independent, zero out the col. + if column_j_norm.is_finite() { + a.get_col_slice_iter_mut(j, 0..N) + .for_each(|x| *x /= column_j_norm); + } else { + a.get_col_slice_iter_mut(j, 0..N).for_each(|x| *x = 0.0); } - } + } // end loop over j a } From 3e2e214aa4f22f86f9abf925c25f82f2e18be879 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:41:25 -0400 Subject: [PATCH 09/18] sum instead of fold --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 1afe90a60..52148c88d 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -15,7 +15,8 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< .get_col(k) .iter_elements() .zip(a.get_col(j).iter_elements()) - .fold(0.0, |acc, (l, r)| acc + (l * r)); + .map(|(l, r)| l * r) + .sum::(); // Apply the projection for i in 0..N { a[(i, j)] -= a[(i, k)] * j_dot_k; @@ -24,7 +25,8 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< let column_j_norm = a .get_col(j) .iter_elements() - .fold(0.0, |acc, x| acc + x * x) + .map(|x| x * x) + .sum::() .sqrt(); // If the initial vectors are not linearly independent, zero out the col. if column_j_norm.is_finite() { From 482eb963f16ff3869893e69f19fb308e822863a1 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:42:08 -0400 Subject: [PATCH 10/18] docs --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 52148c88d..b3998fe3f 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -28,7 +28,7 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< .map(|x| x * x) .sum::() .sqrt(); - // If the initial vectors are not linearly independent, zero out the col. + // If the initial vectors are not linearly independent, zero out the column if column_j_norm.is_finite() { a.get_col_slice_iter_mut(j, 0..N) .for_each(|x| *x /= column_j_norm); From a839c1f9fcf58e060234e8089db8d363afcf787b Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:45:14 -0400 Subject: [PATCH 11/18] handle error --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index b3998fe3f..3ff42a371 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -1,5 +1,7 @@ use crate::matrix::Matrix; +pub const GRAM_SCHMIDT_EPSILON: f64 = 1e-12; + /// Construct an orthonormal basis from the vectors in a [`Matrix`]. /// /// Implementation based on @@ -29,7 +31,7 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< .sum::() .sqrt(); // If the initial vectors are not linearly independent, zero out the column - if column_j_norm.is_finite() { + if column_j_norm > GRAM_SCHMIDT_EPSILON { a.get_col_slice_iter_mut(j, 0..N) .for_each(|x| *x /= column_j_norm); } else { From 8f1f8dd32c85d91b45fb87615856ba4162a6b8e2 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:47:09 -0400 Subject: [PATCH 12/18] Final implementation --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 3ff42a371..3421d6b7f 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -1,8 +1,9 @@ use crate::matrix::Matrix; +/// Tolerance for norms below which a vector is considered to be the zero vector in Gram-Schmidt. pub const GRAM_SCHMIDT_EPSILON: f64 = 1e-12; -/// Construct an orthonormal basis from the vectors in a [`Matrix`]. +/// Construct an orthonormal basis from the vectors in a [`Matrix`] using the modified Gram-Schmidt procedure. /// /// Implementation based on #[must_use] From a02fcaad954d64ddf290c16228f37f11a7c0e39f Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:55:14 -0400 Subject: [PATCH 13/18] fix diags ulpseq to not use Diagonal --- hoomd-linear-algebra/src/matrix.rs | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix.rs b/hoomd-linear-algebra/src/matrix.rs index 27b372bb5..3058dbba9 100644 --- a/hoomd-linear-algebra/src/matrix.rs +++ b/hoomd-linear-algebra/src/matrix.rs @@ -1019,10 +1019,13 @@ pub mod test_utils { } } - pub(crate) fn assert_diags_ulps_eq( - m0: &T, - m1: &impl std::ops::Index, - ) { + pub(crate) fn assert_diags_ulps_eq( + m0: &T0, + m1: &T1, + ) where + T0: Index + ?Sized, + T1: Index + ?Sized, + { for i in 0..N { assert_ulps_eq!(m0[i], m1[i], epsilon = EPS); } @@ -1192,7 +1195,7 @@ mod tests { } assert_matrixes_ulps_eq::<2, 2, _, _>(&u, &faeru); - assert_diags_ulps_eq(&s, &faers); + assert_diags_ulps_eq::<2, _, _>(&s, &faers); // Note that faer returns V, not Vt assert_matrixes_ulps_eq::<2, 2, _, _>(&vt, &faerv.transpose()); } @@ -1226,7 +1229,7 @@ mod tests { let (nau, nas, navt) = (nasvd.u.unwrap(), nasvd.singular_values, nasvd.v_t.unwrap()); assert_matrixes_ulps_eq::<2, 2, _, _>(&u, &nau); - assert_diags_ulps_eq::<2>(&s, &nas); + assert_diags_ulps_eq::<2, _, _>(&s, &nas); assert_matrixes_ulps_eq::<2, 2, _, _>(&vt, &navt); } @@ -1259,7 +1262,7 @@ mod tests { let faers = faersvd.S(); // Our implementation allows negative singular value - assert_diags_ulps_eq( + assert_diags_ulps_eq::<3, _, _>( &DiagonalMatrix { elements: s.elements.map(f64::abs), }, @@ -1335,7 +1338,7 @@ mod tests { let expected_diag = DiagonalMatrix { elements: [1.0, 5.0, 9.0], }; - assert_diags_ulps_eq(&diag, &expected_diag); + assert_diags_ulps_eq::<3, _, _>(&diag, &expected_diag); let from_diag = diag.to_dense(); let expected_from_diag = Matrix { From 659e55f6d6dde29de8c7d2539e65b4343cf4fd9e Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 19:59:53 -0400 Subject: [PATCH 14/18] add tests --- .../src/matrix/gram_schmidt.rs | 82 +++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 3421d6b7f..b27ebc4ff 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -41,3 +41,85 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< } // end loop over j a } + +#[cfg(test)] +mod tests { + use super::*; + use crate::matrix::Matrix; + use crate::matrix::test_utils::assert_matrixes_ulps_eq; + + #[test] + fn test_gram_schmidt_3x3() { + // Example 1 from https://www.sfu.ca/~jtmulhol/py4math/linalg/np-gramschmidt/ + let a = Matrix::<3, 3> { + rows: [[1.0, -1.0, 0.0], [1.0, 2.0, 1.0], [0.0, 1.0, 1.0]], + }; + let q = gram_schmidt(&a); + + let sqrt2 = 2.0f64.sqrt(); + let sqrt11 = 11.0f64.sqrt(); + let sqrt22 = 22.0f64.sqrt(); + + let expected = Matrix::<3, 3> { + rows: [ + [1.0 / sqrt2, -3.0 / sqrt22, 1.0 / sqrt11], + [1.0 / sqrt2, 3.0 / sqrt22, -1.0 / sqrt11], + [0.0, 2.0 / sqrt22, 3.0 / sqrt11], + ], + }; + + assert_matrixes_ulps_eq::<3, 3, _, _>(&q, &expected); + } + + #[test] + fn test_gram_schmidt_linearly_dependent() { + // Example 2 from https://www.sfu.ca/~jtmulhol/py4math/linalg/np-gramschmidt/ + // A = [[1, 1, 2, 1], [1, 0, 1, 0], [0, 1, 1, 0], [1, 1, 2, 1]] (columns) + let a = Matrix::<4, 4> { + rows: [ + [1.0, 1.0, 2.0, 1.0], + [1.0, 0.0, 1.0, 0.0], + [0.0, 1.0, 1.0, 0.0], + [1.0, 1.0, 2.0, 1.0], + ], + }; + let q = gram_schmidt(&a); + + let sqrt3 = 3.0f64.sqrt(); + let sqrt15 = 15.0f64.sqrt(); + let sqrt10 = 10.0f64.sqrt(); + + let expected = Matrix::<4, 4> { + rows: [ + [1.0 / sqrt3, 1.0 / sqrt15, 0.0, 1.0 / sqrt10], + [1.0 / sqrt3, -2.0 / sqrt15, 0.0, -2.0 / sqrt10], + [0.0, 3.0 / sqrt15, 0.0, -2.0 / sqrt10], + [1.0 / sqrt3, 1.0 / sqrt15, 0.0, 1.0 / sqrt10], + ], + }; + + assert_matrixes_ulps_eq::<4, 4, _, _>(&q, &expected); + } + + #[test] + fn test_gram_schmidt_subspace() { + // Example 3 from https://www.sfu.ca/~jtmulhol/py4math/linalg/np-gramschmidt/ + let a = Matrix::<3, 2> { + rows: [[1.0, -1.0], [1.0, 2.0], [0.0, 1.0]], + }; + let q = gram_schmidt(&a); + + let sqrt2 = 2.0f64.sqrt(); + let sqrt22 = 22.0f64.sqrt(); + + let expected = Matrix::<3, 2> { + rows: [ + [1.0 / sqrt2, -3.0 / sqrt22], + [1.0 / sqrt2, 3.0 / sqrt22], + [0.0, 2.0 / sqrt22], + ], + }; + + assert_matrixes_ulps_eq::<3, 2, _, _>(&q, &expected); + } +} From fadb9d14ef2af40799dc98ffa2cbcae30b92ee3e Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 20:04:01 -0400 Subject: [PATCH 15/18] add bench --- hoomd-linear-algebra/benches/linalg.rs | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/hoomd-linear-algebra/benches/linalg.rs b/hoomd-linear-algebra/benches/linalg.rs index 7d3c6a53e..d7e164c3c 100644 --- a/hoomd-linear-algebra/benches/linalg.rs +++ b/hoomd-linear-algebra/benches/linalg.rs @@ -11,7 +11,10 @@ use divan::{self, Bencher, black_box, counter::ItemsCount}; use rand::{Rng, RngExt, SeedableRng, rngs::StdRng}; -use hoomd_linear_algebra::{Invertible, MatMul, matrix::Matrix}; +use hoomd_linear_algebra::{ + Invertible, MatMul, + matrix::{Matrix, gram_schmidt::gram_schmidt}, +}; fn main() { divan::main(); @@ -53,6 +56,16 @@ fn det_matn(bencher: Bencher) { .bench_local_values(|a| black_box(a.determinant())); } +#[divan::bench(consts = SQUARE_DIMENSIONS)] +fn gram_schmidt_matn(bencher: Bencher) { + let mut rng = StdRng::seed_from_u64(42); + + bencher + .counter(ItemsCount::from(1_u32)) + .with_inputs(|| create_random_matrix::(&mut rng)) + .bench_local_values(|a| black_box(gram_schmidt(&a))); +} + /// This benchmark is included as a reference implementation for comparison with ``SquareMatrix::<3, 3>::determinant``. /// We want to ensure the performance of the latter is comparable to this optimal implementation. #[divan::bench] From 471f91a90b839fc5c82980249acf57d767b3b7d7 Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 20:10:45 -0400 Subject: [PATCH 16/18] [WIP] Use raw loops, seems faster for small N --- .../src/matrix/gram_schmidt.rs | 34 ++++++++++--------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index b27ebc4ff..6b5440164 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -10,33 +10,35 @@ pub const GRAM_SCHMIDT_EPSILON: f64 = 1e-12; #[inline] pub fn gram_schmidt(a: &Matrix) -> Matrix { let mut a = a.clone(); - for j in 0..a.n_columns() { + for j in 0..M { // For the vector in column k, find the perpendicular of the projection onto // the previous orthogonal vectors. for k in 0..j { - let j_dot_k = a - .get_col(k) - .iter_elements() - .zip(a.get_col(j).iter_elements()) - .map(|(l, r)| l * r) - .sum::(); + let mut j_dot_k = 0.0; + for i in 0..N { + j_dot_k += a[(i, k)] * a[(i, j)]; + } // Apply the projection for i in 0..N { a[(i, j)] -= a[(i, k)] * j_dot_k; } } // end loop over k - let column_j_norm = a - .get_col(j) - .iter_elements() - .map(|x| x * x) - .sum::() - .sqrt(); + + let mut column_j_norm_sq = 0.0; + for i in 0..N { + column_j_norm_sq += a[(i, j)] * a[(i, j)]; + } + let column_j_norm = column_j_norm_sq.sqrt(); + // If the initial vectors are not linearly independent, zero out the column if column_j_norm > GRAM_SCHMIDT_EPSILON { - a.get_col_slice_iter_mut(j, 0..N) - .for_each(|x| *x /= column_j_norm); + for i in 0..N { + a[(i, j)] /= column_j_norm; + } } else { - a.get_col_slice_iter_mut(j, 0..N).for_each(|x| *x = 0.0); + for i in 0..N { + a[(i, j)] = 0.0; + } } } // end loop over j a From 7b19a4e25005022ace067c7995e0d7271535021b Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 20:38:55 -0400 Subject: [PATCH 17/18] lint --- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 6b5440164..370fd1497 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -1,3 +1,6 @@ +// Copyright (c) 2024-2026 The Regents of the University of Michigan. +// Part of hoomd-rs, released under the BSD 3-Clause License. + use crate::matrix::Matrix; /// Tolerance for norms below which a vector is considered to be the zero vector in Gram-Schmidt. From 7b8288732f67390edc7a990cc2e802a1dad5764e Mon Sep 17 00:00:00 2001 From: janbridley Date: Sun, 3 May 2026 20:39:09 -0400 Subject: [PATCH 18/18] +nightly format --- hoomd-linear-algebra/src/matrix.rs | 6 ++---- hoomd-linear-algebra/src/matrix/gram_schmidt.rs | 3 +-- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/hoomd-linear-algebra/src/matrix.rs b/hoomd-linear-algebra/src/matrix.rs index 3058dbba9..e04d4a657 100644 --- a/hoomd-linear-algebra/src/matrix.rs +++ b/hoomd-linear-algebra/src/matrix.rs @@ -1019,10 +1019,8 @@ pub mod test_utils { } } - pub(crate) fn assert_diags_ulps_eq( - m0: &T0, - m1: &T1, - ) where + pub(crate) fn assert_diags_ulps_eq(m0: &T0, m1: &T1) + where T0: Index + ?Sized, T1: Index + ?Sized, { diff --git a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs index 370fd1497..c99e4d5c1 100644 --- a/hoomd-linear-algebra/src/matrix/gram_schmidt.rs +++ b/hoomd-linear-algebra/src/matrix/gram_schmidt.rs @@ -50,8 +50,7 @@ pub fn gram_schmidt(a: &Matrix) -> Matrix< #[cfg(test)] mod tests { use super::*; - use crate::matrix::Matrix; - use crate::matrix::test_utils::assert_matrixes_ulps_eq; + use crate::matrix::{Matrix, test_utils::assert_matrixes_ulps_eq}; #[test] fn test_gram_schmidt_3x3() {