Fixes for nalgebra-lapack (faulty test, bugs in implementations, compile errors) (#1541)

* fix bugs in lu and schur decomposition, update tests for singularity detection

* improve cholesky tests for better output

* make sure in cholesky tests to have strictly positive definite matrices

* relax numerical accuracy in cholesky tests

* make all cholesky tests use the strictly sym pos def matrices

* add nalgebra-lapack tests to CI

* fix copy paste error in ci yml

* fix serde-serialize bounds

* add serde-serialize to CI features for na-lapack

* add slow-tests feature

* fix unused imports

Author: geo-ant

.github/workflows/nalgebra-ci-build.yml

@@ -36,8 +36,8 @@ jobs:
         run: cargo build;
       - name: Build --features serde-serialize
         run: cargo build --features serde-serialize
-      - name: Build nalgebra-lapack
-        run: cd nalgebra-lapack; cargo build;
+      - name: Check nalgebra-lapack
+        run: cd nalgebra-lapack; cargo check --features serde-serialize;
       - name: Build nalgebra-sparse --no-default-features
         run: cd nalgebra-sparse; cargo build --no-default-features;
       - name: Build nalgebra-sparse (default features)
@@ -105,6 +105,12 @@ jobs:
       - uses: actions/checkout@v4
       - name: test nalgebra-macros
         run: cargo test -p nalgebra-macros
+  test-nalgebra-lapack:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: test nalgebra-lapack
+        run: cargo test -p nalgebra-lapack --features proptest-support;
   build-wasm:
     runs-on: ubuntu-latest
     #    env:

nalgebra-lapack/Cargo.toml

@@ -20,6 +20,7 @@ maintenance = { status = "actively-developed" }
 serde-serialize = ["serde", "nalgebra/serde-serialize"]
 proptest-support = ["nalgebra/proptest-support"]
 arbitrary = ["nalgebra/arbitrary"]
+slow-tests = []
 
 # For BLAS/LAPACK
 default = ["netlib"]

nalgebra-lapack/src/eigen.rs

@@ -23,7 +23,7 @@ use lapack;
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(deserialize = "DefaultAllocator: Allocator<D, D> + Allocator<D>,
-         OVector<T, D>: Serialize,
+         OVector<T, D>: Deserialize<'de>,
          OMatrix<T, D, D>: Deserialize<'de>"))
 )]
 #[derive(Clone, Debug)]

nalgebra-lapack/src/hessenberg.rs

@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde::{Deserialize, Serialize};
+
 use num::Zero;
 use num_complex::Complex;
 

nalgebra-lapack/src/lu.rs

@@ -1,3 +1,6 @@
+#[cfg(feature = "serde-serialize")]
+use serde::{Deserialize, Serialize};
+
 use num::{One, Zero};
 use num_complex::Complex;
 
@@ -23,14 +26,14 @@ use lapack;
     serde(bound(serialize = "DefaultAllocator: Allocator<R, C> +
                            Allocator<DimMinimum<R, C>>,
          OMatrix<T, R, C>: Serialize,
-         PermutationSequence<DimMinimum<R, C>>: Serialize"))
+         OVector<i32, DimMinimum<R, C>>: Serialize"))
 )]
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(deserialize = "DefaultAllocator: Allocator<R, C> +
                            Allocator<DimMinimum<R, C>>,
          OMatrix<T, R, C>: Deserialize<'de>,
-         PermutationSequence<DimMinimum<R, C>>: Deserialize<'de>"))
+         OVector<i32, DimMinimum<R, C>>: Deserialize<'de>"))
 )]
 #[derive(Clone, Debug)]
 pub struct LU<T: Scalar, R: DimMin<C>, C: Dim>
@@ -39,6 +42,9 @@ where
 {
     lu: OMatrix<T, R, C>,
     p: OVector<i32, DimMinimum<R, C>>,
+    /// the LU decomposition can be computed for singular matrices, but we
+    /// cannot solve linear systems for singular matrices.
+    singular: bool,
 }
 
 impl<T: Scalar + Copy, R: DimMin<C>, C: Dim> Copy for LU<T, R, C>
@@ -78,9 +84,30 @@ where
             ipiv.as_mut_slice(),
             &mut info,
         );
-        lapack_panic!(info);
 
-        Self { lu: m, p: ipiv }
+        // @note(geo-ant)
+        // the lapack documentation for xGETRF states:
+        //
+        // !>          INFO is INTEGER
+        // !>          = 0:  successful exit
+        // !>          < 0:  if INFO = -i, the i-th argument had an illegal value
+        // !>          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization
+        // !>                has been completed, but the factor U is exactly
+        // !>                singular, and division by zero will occur if it is used
+        // !>                to solve a system of equations.
+        // !>
+        //
+        // This is a bit unusual for lapack routines, so we have to check whether it is
+        // negative before panicking.
+        if info < 0 {
+            lapack_panic!(info);
+        }
+
+        Self {
+            lu: m,
+            p: ipiv,
+            singular: info > 0,
+        }
     }
 
     /// Gets the lower-triangular matrix part of the decomposition.
@@ -155,6 +182,10 @@ where
     where
         DefaultAllocator: Allocator<R2, C2> + Allocator<R2>,
     {
+        if self.singular {
+            return false;
+        }
+
         let dim = self.lu.nrows();
 
         assert!(
@@ -195,7 +226,7 @@ where
         DefaultAllocator: Allocator<R2, C2> + Allocator<R2>,
     {
         let mut res = b.clone_owned();
-        if self.generic_solve_mut(b'T', &mut res) {
+        if self.generic_solve_mut(b'N', &mut res) {
             Some(res)
         } else {
             None
@@ -231,7 +262,7 @@ where
         DefaultAllocator: Allocator<R2, C2> + Allocator<R2>,
     {
         let mut res = b.clone_owned();
-        if self.generic_solve_mut(b'T', &mut res) {
+        if self.generic_solve_mut(b'C', &mut res) {
             Some(res)
         } else {
             None
@@ -245,7 +276,7 @@ where
     where
         DefaultAllocator: Allocator<R2, C2> + Allocator<R2>,
     {
-        self.generic_solve_mut(b'T', b)
+        self.generic_solve_mut(b'N', b)
     }
 
     /// Solves in-place the linear system `self.transpose() * x = b`, where `x` is the unknown to be
@@ -267,7 +298,7 @@ where
     where
         DefaultAllocator: Allocator<R2, C2> + Allocator<R2>,
     {
-        self.generic_solve_mut(b'T', b)
+        self.generic_solve_mut(b'C', b)
     }
 }
 
@@ -279,6 +310,10 @@ where
 {
     /// Computes the inverse of the decomposed matrix.
     pub fn inverse(mut self) -> Option<OMatrix<T, D, D>> {
+        if self.singular {
+            return None;
+        }
+
         let dim = self.lu.nrows() as i32;
         let mut info = 0;
         let lwork = T::xgetri_work_size(

nalgebra-lapack/src/schur.rs

@@ -24,7 +24,7 @@ use lapack;
 #[cfg_attr(
     feature = "serde-serialize",
     serde(bound(deserialize = "DefaultAllocator: Allocator<D, D> + Allocator<D>,
-         OVector<T, D>: Serialize,
+         OVector<T, D>: Deserialize<'de>,
          OMatrix<T, D, D>: Deserialize<'de>"))
 )]
 #[derive(Clone, Debug)]
@@ -82,7 +82,11 @@ where
 
         let lwork = T::xgees_work_size(
             b'V',
-            b'T',
+            //@note(geo-ant): this interface is bad because we could pass
+            // 'S' here, but we cannot give a function to SELECT, which
+            // isn't exposed here.
+            // https://www.netlib.org/lapack/explore-html/d5/d38/group__gees.html
+            b'N',
             n as i32,
             m.as_mut_slice(),
             lda,
@@ -100,7 +104,7 @@ where
 
         T::xgees(
             b'V',
-            b'T',
+            b'N',
             n as i32,
             m.as_mut_slice(),
             lda,
@@ -169,6 +173,7 @@ pub trait SchurScalar: Scalar {
     fn xgees(
         jobvs: u8,
         sort: u8,
+        //@note(geo-ant) see other comments in this file on this method
         // select: ???
         n: i32,
         a: &mut [Self],
@@ -208,6 +213,10 @@ macro_rules! real_eigensystem_scalar_impl (
             #[inline]
             fn xgees(jobvs:  u8,
                      sort:   u8,
+                     //@note(geo-ant) see the documentation for xGEES here:
+                     // https://www.netlib.org/lapack/explore-html/d5/d38/group__gees.html
+                     // If we pass a sort argument of b'S', then we must offer a
+                     // SELECT implementation.
                      // select: ???
                      n:      i32,
                      a:      &mut [$N],

nalgebra-lapack/src/svd.rs

@@ -17,17 +17,17 @@ use lapack;
     serde(bound(serialize = "DefaultAllocator: Allocator<DimMinimum<R, C>> +
                            Allocator<R, R> +
                            Allocator<C, C>,
-         OMatrix<T, R>: Serialize,
-         OMatrix<T, C>: Serialize,
+         OMatrix<T, R, R>: Serialize,
+         OMatrix<T, C, C>: Serialize,
          OVector<T, DimMinimum<R, C>>: Serialize"))
 )]
 #[cfg_attr(
     feature = "serde-serialize",
-    serde(bound(serialize = "DefaultAllocator: Allocator<DimMinimum<R, C>> +
-                           Allocator<R, R> +
-                           Allocator<C, C>,
-         OMatrix<T, R>: Deserialize<'de>,
-         OMatrix<T, C>: Deserialize<'de>,
+    serde(bound(deserialize = "DefaultAllocator: Allocator<DimMinimum<R, C>> +
+                             Allocator<R, R> +
+                             Allocator<C, C>,
+         OMatrix<T, R, R>: Deserialize<'de>,
+         OMatrix<T, C, C>: Deserialize<'de>,
          OVector<T, DimMinimum<R, C>>: Deserialize<'de>"))
 )]
 #[derive(Clone, Debug)]

nalgebra-lapack/tests/linalg/cholesky.rs

@@ -1,15 +1,44 @@
-use std::cmp;
-
-use na::{DMatrix, DVector, Matrix4x3, Vector4};
+use crate::proptest::*;
+use core::f64;
+use na::{Const, DMatrix, Matrix4, Matrix4xX};
 use nl::Cholesky;
+use proptest::prelude::*;
+
+fn positive_definite_dmatrix() -> impl Strategy<Value = DMatrix<f64>> {
+    // @note(geo-ant) to get positive definite matrices we use M*M^T + alpha*I,
+    // where alpha is a constant that is chosen so that the eigenvales stay
+    // positive.
+    dmatrix().prop_map(|m| {
+        let alpha = f64::EPSILON.sqrt() * m.norm_squared();
+        let nrows = m.nrows();
+        &m * m.transpose() + alpha * DMatrix::identity(nrows, nrows)
+    })
+}
 
-use crate::proptest::*;
-use proptest::{prop_assert, proptest};
+fn positive_definite_matrix4() -> impl Strategy<Value = Matrix4<f64>> {
+    matrix4().prop_map(|m| {
+        let alpha = f64::EPSILON.sqrt() * m.norm_squared();
+        &m * m.transpose() + alpha * Matrix4::identity()
+    })
+}
+
+fn positive_definite_linear_system() -> impl Strategy<Value = (DMatrix<f64>, DMatrix<f64>)> {
+    positive_definite_dmatrix().prop_flat_map(|a| {
+        let b = matrix(PROPTEST_F64, a.nrows(), PROPTEST_MATRIX_DIM);
+        (Just(a), b)
+    })
+}
+
+fn positive_definite_linear_system_4() -> impl Strategy<Value = (Matrix4<f64>, Matrix4xX<f64>)> {
+    positive_definite_matrix4().prop_flat_map(|a| {
+        let b = matrix(PROPTEST_F64, Const::<4>, PROPTEST_MATRIX_DIM);
+        (Just(a), b)
+    })
+}
 
 proptest! {
     #[test]
-    fn cholesky(m in dmatrix()) {
-        let m = &m * m.transpose();
+    fn cholesky(m in positive_definite_dmatrix()) {
         if let Some(chol) = Cholesky::new(m.clone()) {
             let l = chol.unpack();
             let reconstructed_m = &l * l.transpose();
@@ -19,8 +48,7 @@ proptest! {
     }
 
     #[test]
-    fn cholesky_static(m in matrix3()) {
-        let m = &m * m.transpose();
+    fn cholesky_static(m in positive_definite_matrix4()) {
         if let Some(chol) = Cholesky::new(m) {
             let l = chol.unpack();
             let reconstructed_m = &l * l.transpose();
@@ -30,61 +58,37 @@ proptest! {
     }
 
     #[test]
-    fn cholesky_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
-        let n  = cmp::min(n, 15);  // To avoid slowing down the test too much.
-        let nb = cmp::min(nb, 15); // To avoid slowing down the test too much.
-        let m  = DMatrix::<f64>::new_random(n, n);
-        let m   = &m * m.transpose();
+    fn cholesky_solve((a,b) in positive_definite_linear_system()) {
 
-        if let Some(chol) = Cholesky::new(m.clone()) {
-            let b1 = DVector::new_random(n);
-            let b2 = DMatrix::new_random(n, nb);
-
-            let sol1 = chol.solve(&b1).unwrap();
-            let sol2 = chol.solve(&b2).unwrap();
-
-            prop_assert!(relative_eq!(&m * sol1, b1, epsilon = 1.0e-6));
-            prop_assert!(relative_eq!(&m * sol2, b2, epsilon = 1.0e-6));
+        if let Some(chol) = Cholesky::new(a.clone()) {
+            let sol = chol.solve(&b).unwrap();
+            prop_assert!(relative_eq!(&a * sol, b, epsilon = 1.0e-5));
         }
     }
 
     #[test]
-    fn cholesky_solve_static(m in matrix4()) {
-        let m = &m * m.transpose();
-        if let Some(chol) = Cholesky::new(m) {
-            let b1 = Vector4::new_random();
-            let b2 = Matrix4x3::new_random();
-
-            let sol1 = chol.solve(&b1).unwrap();
-            let sol2 = chol.solve(&b2).unwrap();
-
-            prop_assert!(relative_eq!(m * sol1, b1, epsilon = 1.0e-4));
-            prop_assert!(relative_eq!(m * sol2, b2, epsilon = 1.0e-4));
+    fn cholesky_solve_static((a,b) in positive_definite_linear_system_4()) {
+        if let Some(chol) = Cholesky::new(a) {
+            let sol = chol.solve(&b).unwrap();
+            prop_assert!(relative_eq!(a * sol, b, epsilon = 1.0e-5));
         }
     }
 
     #[test]
-    fn cholesky_inverse(n in PROPTEST_MATRIX_DIM) {
-        let n = cmp::min(n, 15); // To avoid slowing down the test too much.
-        let m = DMatrix::<f64>::new_random(n, n);
-        let m = &m * m.transpose();
-
-        if let Some(m1) = Cholesky::new(m.clone()).unwrap().inverse() {
-            let id1 = &m  * &m1;
-            let id2 = &m1 * &m;
+    fn cholesky_inverse(a in positive_definite_dmatrix()) {
+        let minv = Cholesky::new(a.clone()).unwrap().inverse().unwrap();
+        let id1 = &a  * &minv;
+        let id2 = &minv * &a;
 
-            prop_assert!(id1.is_identity(1.0e-6) && id2.is_identity(1.0e-6));
-        }
+        prop_assert!(id1.is_identity(1.0e-6) && id2.is_identity(1.0e-6));
     }
 
     #[test]
-    fn cholesky_inverse_static(m in matrix4()) {
-        let m = m * m.transpose();
-        if let Some(m1) = Cholesky::new(m.clone()).unwrap().inverse() {
-            let id1 = &m  * &m1;
-            let id2 = &m1 * &m;
+    fn cholesky_inverse_static(a in positive_definite_matrix4()) {
+        let minv = Cholesky::new(a.clone()).unwrap().inverse().unwrap();
+        let id1 = &a  * &minv;
+        let id2 = &minv * &a;
 
-            prop_assert!(id1.is_identity(1.0e-4) && id2.is_identity(1.0e-4))
-        }
+        prop_assert!(id1.is_identity(1.0e-6) && id2.is_identity(1.0e-6));
     }
 }