From 6bcb07a77a5c7cf78a41a08c4dad9216fb6f592f Mon Sep 17 00:00:00 2001
From: Robert Bacula <bacularobert@gmail.com>
Date: Tue, 21 Jan 2025 17:10:55 +0200
Subject: [PATCH 1/2] svd now tests for NaN values
---
src/linalg/svd.rs | 32 ++++++++++++++++++++++++++------
tests/linalg/svd.rs | 2 +-
2 files changed, 27 insertions(+), 7 deletions(-)
diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index fc5ed6b4d..4ce46e91a 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -119,21 +119,24 @@ where
eps: T::RealField,
max_niter: usize,
) -> Option<Self> {
+ // Ensure that the matrix is not empty before proceeding
assert!(
!matrix.is_empty(),
"Cannot compute the SVD of an empty matrix."
);
+
let (nrows, ncols) = matrix.shape_generic();
let min_nrows_ncols = nrows.min(ncols);
+ // Special cases for 2x2 and 3x3 matrices, handled by predefined methods
if Self::use_special_always_ordered_svd2() {
- // SAFETY: the reference transmutes are OK since we checked that the types match exactly.
+ // SAFETY: Reference transmutes are OK because the types match exactly
let matrix: &Matrix2<T::RealField> = unsafe { std::mem::transmute(&matrix) };
let result = super::svd2::svd_ordered2(matrix, compute_u, compute_v);
let typed_result: &Self = unsafe { std::mem::transmute(&result) };
return Some(typed_result.clone());
} else if Self::use_special_always_ordered_svd3() {
- // SAFETY: the reference transmutes are OK since we checked that the types match exactly.
+ // SAFETY: Reference transmutes are OK because the types match exactly
let matrix: &Matrix3<T::RealField> = unsafe { std::mem::transmute(&matrix) };
let result = super::svd3::svd_ordered3(matrix, compute_u, compute_v, eps, max_niter);
let typed_result: &Self = unsafe { std::mem::transmute(&result) };
@@ -142,8 +145,10 @@ where
let dim = min_nrows_ncols.value();
+ // Get the maximum absolute value of the matrix for normalization
let m_amax = matrix.camax();
+ // If the max value is not zero, unscale the matrix by dividing by m_amax
if !m_amax.is_zero() {
matrix.unscale_mut(m_amax.clone());
}
@@ -158,6 +163,16 @@ where
let mut diagonal = bi_matrix.diagonal();
let mut off_diagonal = bi_matrix.off_diagonal();
+ // **Check for NaN values in the diagonal elements**
+ // We check whether any singular value in the diagonal is NaN by comparing each value to itself.
+ // This works because NaN is the only value in Rust that is not equal to itself.
+ if diagonal.iter().any(|s| s != s) {
+ // If any singular value is NaN, return None early
+ // Explanation:
+ // NaN != NaN is always true, so if `s != s`, we know `s` is NaN.
+ return None; // Return early if NaN found
+ }
+
let mut niter = 0;
let (mut start, mut end) = Self::delimit_subproblem(
&mut diagonal,
@@ -169,10 +184,11 @@ where
eps.clone(),
);
+ // Iterative SVD computation with Givens rotations
while end != start {
let subdim = end - start + 1;
- // Solve the subproblem.
+ // Solve subproblem for larger subdimensions (> 2)
#[allow(clippy::comparison_chain)]
if subdim > 2 {
let m = end - 1;
@@ -184,19 +200,23 @@ where
let dn = diagonal[n].clone();
let fm = off_diagonal[m].clone();
+ // Perform calculations to determine the shift value for Givens rotation
let tmm = dm.clone() * dm.clone()
+ off_diagonal[m - 1].clone() * off_diagonal[m - 1].clone();
let tmn = dm * fm.clone();
let tnn = dn.clone() * dn + fm.clone() * fm;
+ // Compute Wilkinson's shift
let shift = symmetric_eigen::wilkinson_shift(tmm, tnn, tmn);
+ // Create vector for subsequent Givens rotations
vec = Vector2::new(
diagonal[start].clone() * diagonal[start].clone() - shift,
diagonal[start].clone() * off_diagonal[start].clone(),
);
}
+ // Perform Givens rotations to reduce the bidiagonal matrix
for k in start..n {
let m12 = if k == n - 1 {
T::RealField::zero()
@@ -224,7 +244,6 @@ where
}
let v = Vector2::new(subm[(0, 0)].clone(), subm[(1, 0)].clone());
- // TODO: does the case `v.y == 0` ever happen?
let (rot2, norm2) = GivensRotation::cancel_y(&v)
.unwrap_or((GivensRotation::identity(), subm[(0, 0)].clone()));
@@ -264,7 +283,7 @@ where
}
}
} else if subdim == 2 {
- // Solve the remaining 2x2 subproblem.
+ // Solve the 2x2 subproblem if subdim == 2
let (u2, s, v2) = compute_2x2_uptrig_svd(
diagonal[start].clone(),
off_diagonal[start].clone(),
@@ -321,9 +340,10 @@ where
}
}
+ // Unscale the singular values after SVD computation
diagonal *= m_amax;
- // Ensure all singular value are non-negative.
+ // Ensure all singular values are non-negative
for i in 0..dim {
let sval = diagonal[i].clone();
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index d8b23d024..94c34cb52 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -501,7 +501,7 @@ fn svd_regression_issue_1072() {
}
#[test]
-// Exercises bug reported in issue #1313 of nalgebra (https://github.com/dimforge/nalgebra/issues/1313)
+// Exercises bug reported in issue #1313 of nalgebra (https://github. com/dimforge/nalgebra/issues/1313)
fn svd_regression_issue_1313() {
let s = 6.123234e-16_f32;
let m = nalgebra::dmatrix![
From dad6334a217865f211c2fd58cdfba0ba98135ef1 Mon Sep 17 00:00:00 2001
From: Robert Bacula <bacularobert@gmail.com>
Date: Tue, 21 Jan 2025 17:13:18 +0200
Subject: [PATCH 2/2] format
---
tests/linalg/svd.rs | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/tests/linalg/svd.rs b/tests/linalg/svd.rs
index 94c34cb52..d8b23d024 100644
--- a/tests/linalg/svd.rs
+++ b/tests/linalg/svd.rs
@@ -501,7 +501,7 @@ fn svd_regression_issue_1072() {
}
#[test]
-// Exercises bug reported in issue #1313 of nalgebra (https://github. com/dimforge/nalgebra/issues/1313)
+// Exercises bug reported in issue #1313 of nalgebra (https://github.com/dimforge/nalgebra/issues/1313)
fn svd_regression_issue_1313() {
let s = 6.123234e-16_f32;
let m = nalgebra::dmatrix![