diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs index 158f899c7..c3081b6f7 100644 --- a/src/linalg/schur.rs +++ b/src/linalg/schur.rs @@ -135,6 +135,16 @@ where // Implicit double-shift QR method. let mut niter = 0; + // Number of QR steps since the active window last shrank due to a deflation. + let mut kdefl: usize = 0; + // Apply Wilkinson's exceptional shift every kexsh stalled QR steps. + let kexsh: usize = 10; + // This is used by the exceptional shift algorithm as mentioned in + // Wilkinson, J. H., & Reinsch, C. (1971). + // Handbook for automatic computation (Vol. 2: Linear algebra) + // P.373 + let exceptional_shift_a: T::RealField = crate::convert(0.75); + let exceptional_shift_b: T::RealField = crate::convert(-0.4375); let (mut start, mut end) = Self::delimit_subproblem(&mut t, eps.clone(), dim.value() - 1); while end != start { @@ -150,10 +160,28 @@ where let h22 = t[(start + 1, start + 1)].clone(); let h32 = t[(start + 2, start + 1)].clone(); - let hnn = t[(n, n)].clone(); - let hmm = t[(m, m)].clone(); - let hnm = t[(n, m)].clone(); - let hmn = t[(m, n)].clone(); + let (hnn, hmm, hnm, hmn) = if kdefl != 0 && kdefl % kexsh == 0 { + let s = t[(start + 1, start)].clone().norm1() + + t[(start + 2, start + 1)].clone().norm1(); + let h = t[(start, start)].clone() + + T::from_real(exceptional_shift_a.clone() * s.clone()); + // When the implicit QR step stalls, replace the trailing 2x2 Wilkinson shift + // with the exceptional-shift coefficients from Wilkinson & Reinsch. This + // perturbs the bulge start enough to escape the non-convergent cycle. + ( + h.clone(), + h, + T::from_real(exceptional_shift_b.clone() * s.clone()), + T::from_real(s), + ) + } else { + ( + t[(n, n)].clone(), + t[(m, m)].clone(), + t[(n, m)].clone(), + t[(m, n)].clone(), + ) + }; let tra = hnn.clone() + hmm.clone(); let det = hnn * hmm - hnm * hmn; @@ -256,10 +284,19 @@ where } } + let prev_start = start; + let prev_end = end; let sub = Self::delimit_subproblem(&mut t, eps.clone(), end); start = sub.0; end = sub.1; + if end < prev_end || start > prev_start { + // A split or deflation changed the active window, so restart the stall counter. + kdefl = 0; + } else { + // No deflation occurred: count another QR step toward an exceptional shift. + kdefl += 1; + } niter += 1; if niter == max_niter { diff --git a/tests/linalg/eigen.rs b/tests/linalg/eigen.rs index 2de280cd5..fc8385d06 100644 --- a/tests/linalg/eigen.rs +++ b/tests/linalg/eigen.rs @@ -170,6 +170,31 @@ fn eigenvalues_search_should_not_hang_issue_1528() { ); } +// Regression test for #611 +#[test] +#[rustfmt::skip] +fn eigenvalues_search_should_not_hang_issue_611() { + let m = nalgebra::Matrix4::::new( + 0.0, 0.0, 0.0, -0.8286, + 1.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 1.7094, + 0.0, 0.0, 1.0, 0.0, + ); + let complex_eigenvals = m.complex_eigenvalues(); + + assert_relative_eq!( + complex_eigenvals.iter().sum::>().re, + m.trace(), + epsilon = 1e-10 + ); + + assert_relative_eq!( + complex_eigenvals.iter().product::>().re, + m.determinant(), + epsilon = 1e-10 + ); +} + // #[cfg(feature = "arbitrary")] // quickcheck! { // TODO: full eigendecomposition is not implemented yet because of its complexity when some