Add .ln_detc*() methods to DeterminantC*

Author: jturner314

src/cholesky.rs

@@ -102,12 +102,15 @@ where
     type Output = <A as AssociatedReal>::Real;
 
     fn detc(&self) -> Self::Output {
+        self.ln_detc().exp()
+    }
+
+    fn ln_detc(&self) -> Self::Output {
         self.factor
             .diag()
             .iter()
             .map(|elem| elem.abs_sqr().ln())
             .sum::<Self::Output>()
-            .exp()
     }
 }
 
@@ -121,6 +124,10 @@ where
     fn detc_into(self) -> Self::Output {
         self.detc()
     }
+
+    fn ln_detc_into(self) -> Self::Output {
+        self.ln_detc()
+    }
 }
 
 impl<A, S> InverseC for CholeskyFactorized<S>
@@ -391,6 +398,14 @@ pub trait DeterminantC {
     /// Computes the determinant of the Hermitian (or real symmetric) positive
     /// definite matrix.
     fn detc(&self) -> Self::Output;
+
+    /// Computes the natural log of the determinant of the Hermitian (or real
+    /// symmetric) positive definite matrix.
+    ///
+    /// This method is more robust than `.detc()` to very small or very large
+    /// determinants since it returns the natural logarithm of the determinant
+    /// rather than the determinant itself.
+    fn ln_detc(&self) -> Self::Output;
 }
 
 
@@ -401,6 +416,14 @@ pub trait DeterminantCInto {
     /// Computes the determinant of the Hermitian (or real symmetric) positive
     /// definite matrix.
     fn detc_into(self) -> Self::Output;
+
+    /// Computes the natural log of the determinant of the Hermitian (or real
+    /// symmetric) positive definite matrix.
+    ///
+    /// This method is more robust than `.detc_into()` to very small or very
+    /// large determinants since it returns the natural logarithm of the
+    /// determinant rather than the determinant itself.
+    fn ln_detc_into(self) -> Self::Output;
 }
 
 impl<A, S> DeterminantC for ArrayBase<S, Ix2>
@@ -411,7 +434,11 @@ where
     type Output = Result<<A as AssociatedReal>::Real>;
 
     fn detc(&self) -> Self::Output {
-        Ok(self.factorizec(UPLO::Upper)?.detc())
+        Ok(self.ln_detc()?.exp())
+    }
+
+    fn ln_detc(&self) -> Self::Output {
+        Ok(self.factorizec(UPLO::Upper)?.ln_detc())
     }
 }
 
@@ -423,6 +450,10 @@ where
     type Output = Result<<A as AssociatedReal>::Real>;
 
     fn detc_into(self) -> Self::Output {
-        Ok(self.factorizec_into(UPLO::Upper)?.detc_into())
+        Ok(self.ln_detc_into()?.exp())
+    }
+
+    fn ln_detc_into(self) -> Self::Output {
+        Ok(self.factorizec_into(UPLO::Upper)?.ln_detc_into())
     }
 }

tests/cholesky.rs

@@ -104,10 +104,16 @@ fn cholesky_det() {
         ($elem:ty, $atol:expr) => {
             let a: Array2<$elem> = random_hpd(3);
             println!("a = \n{:?}", a);
-            let det = a.eigvalsh(UPLO::Upper).unwrap().mapv(|elem| elem.ln()).scalar_sum().exp();
-            assert_aclose!(a.detc().unwrap(), det, $atol);
+            let ln_det = a.eigvalsh(UPLO::Upper).unwrap().mapv(|elem| elem.ln()).scalar_sum();
+            let det = ln_det.exp();
             assert_aclose!(a.factorizec(UPLO::Upper).unwrap().detc(), det, $atol);
-            assert_aclose!(a.factorizec(UPLO::Lower).unwrap().detc(), det, $atol);
+            assert_aclose!(a.factorizec(UPLO::Upper).unwrap().ln_detc(), ln_det, $atol);
+            assert_aclose!(a.factorizec(UPLO::Lower).unwrap().detc_into(), det, $atol);
+            assert_aclose!(a.factorizec(UPLO::Lower).unwrap().ln_detc_into(), ln_det, $atol);
+            assert_aclose!(a.detc().unwrap(), det, $atol);
+            assert_aclose!(a.ln_detc().unwrap(), ln_det, $atol);
+            assert_aclose!(a.clone().detc_into().unwrap(), det, $atol);
+            assert_aclose!(a.ln_detc_into().unwrap(), ln_det, $atol);
         }
     }
     cholesky_det!(f64, 1e-9);