Merge pull request #104 from jturner314/add-ln-det

Add methods to get natural log of determinant directly

Author: termoshtt

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())
     }
 }

src/solve.rs

@@ -47,6 +47,7 @@
 //! ```
 
 use ndarray::*;
+use num_traits::{Float, Zero};
 
 use super::convert::*;
 use super::error::*;
@@ -336,16 +337,54 @@ where
 /// An interface for calculating determinants of matrix refs.
 pub trait Determinant<A: Scalar> {
     /// Computes the determinant of the matrix.
-    fn det(&self) -> Result<A>;
+    fn det(&self) -> Result<A> {
+        let (sign, ln_det) = self.sln_det()?;
+        Ok(sign.mul_real(ln_det.exp()))
+    }
+
+    /// Computes the `(sign, natural_log)` of the determinant of the matrix.
+    ///
+    /// For real matrices, `sign` is `1`, `0`, or `-1`. For complex matrices,
+    /// `sign` is `0` or a complex number with absolute value 1. The
+    /// `natural_log` is the natural logarithm of the absolute value of the
+    /// determinant. If the determinant is zero, `sign` is 0 and `natural_log`
+    /// is negative infinity.
+    ///
+    /// To obtain the determinant, you can compute `sign * natural_log.exp()`
+    /// or just call `.det()` instead.
+    ///
+    /// This method is more robust than `.det()` to very small or very large
+    /// determinants since it returns the natural logarithm of the determinant
+    /// rather than the determinant itself.
+    fn sln_det(&self) -> Result<(A, A::Real)>;
 }
 
 /// An interface for calculating determinants of matrices.
-pub trait DeterminantInto<A: Scalar> {
+pub trait DeterminantInto<A: Scalar>: Sized {
     /// Computes the determinant of the matrix.
-    fn det_into(self) -> Result<A>;
+    fn det_into(self) -> Result<A> {
+        let (sign, ln_det) = self.sln_det_into()?;
+        Ok(sign.mul_real(ln_det.exp()))
+    }
+
+    /// Computes the `(sign, natural_log)` of the determinant of the matrix.
+    ///
+    /// For real matrices, `sign` is `1`, `0`, or `-1`. For complex matrices,
+    /// `sign` is `0` or a complex number with absolute value 1. The
+    /// `natural_log` is the natural logarithm of the absolute value of the
+    /// determinant. If the determinant is zero, `sign` is 0 and `natural_log`
+    /// is negative infinity.
+    ///
+    /// To obtain the determinant, you can compute `sign * natural_log.exp()`
+    /// or just call `.det_into()` instead.
+    ///
+    /// This method is more robust than `.det()` to very small or very large
+    /// determinants since it returns the natural logarithm of the determinant
+    /// rather than the determinant itself.
+    fn sln_det_into(self) -> Result<(A, A::Real)>;
 }
 
-fn lu_det<'a, A, P, U>(ipiv_iter: P, u_diag_iter: U) -> A
+fn lu_sln_det<'a, A, P, U>(ipiv_iter: P, u_diag_iter: U) -> (A, A::Real)
 where
     A: Scalar,
     P: Iterator<Item = i32>,
@@ -360,24 +399,27 @@ where
     } else {
         -A::one()
     };
-    let (upper_sign, ln_det) = u_diag_iter.fold((A::one(), A::zero()), |(upper_sign, ln_det), &elem| {
-        let abs_elem = elem.abs();
-        (
-            upper_sign * elem.div_real(abs_elem),
-            ln_det.add_real(abs_elem.ln()),
-        )
-    });
-    pivot_sign * upper_sign * ln_det.exp()
+    let (upper_sign, ln_det) = u_diag_iter.fold(
+        (A::one(), A::Real::zero()),
+        |(upper_sign, ln_det), &elem| {
+            let abs_elem: A::Real = elem.abs();
+            (upper_sign * elem.div_real(abs_elem), ln_det + abs_elem.ln())
+        },
+    );
+    (pivot_sign * upper_sign, ln_det)
 }
 
 impl<A, S> Determinant<A> for LUFactorized<S>
 where
     A: Scalar,
     S: Data<Elem = A>,
 {
-    fn det(&self) -> Result<A> {
+    fn sln_det(&self) -> Result<(A, A::Real)> {
         self.a.ensure_square()?;
-        Ok(lu_det(self.ipiv.iter().cloned(), self.a.diag().iter()))
+        Ok(lu_sln_det(
+            self.ipiv.iter().cloned(),
+            self.a.diag().iter(),
+        ))
     }
 }
 
@@ -386,9 +428,12 @@ where
     A: Scalar,
     S: Data<Elem = A>,
 {
-    fn det_into(self) -> Result<A> {
+    fn sln_det_into(self) -> Result<(A, A::Real)> {
         self.a.ensure_square()?;
-        Ok(lu_det(self.ipiv.into_iter(), self.a.into_diag().iter()))
+        Ok(lu_sln_det(
+            self.ipiv.into_iter(),
+            self.a.into_diag().iter(),
+        ))
     }
 }
 
@@ -397,11 +442,14 @@ where
     A: Scalar,
     S: Data<Elem = A>,
 {
-    fn det(&self) -> Result<A> {
+    fn sln_det(&self) -> Result<(A, A::Real)> {
         self.ensure_square()?;
         match self.factorize() {
-            Ok(fac) => fac.det(),
-            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::zero()),
+            Ok(fac) => fac.sln_det(),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => {
+                // The determinant is zero.
+                Ok((A::zero(), A::Real::neg_infinity()))
+            }
             Err(err) => Err(err),
         }
     }
@@ -412,11 +460,14 @@ where
     A: Scalar,
     S: DataMut<Elem = A>,
 {
-    fn det_into(self) -> Result<A> {
+    fn sln_det_into(self) -> Result<(A, A::Real)> {
         self.ensure_square()?;
         match self.factorize_into() {
-            Ok(fac) => fac.det_into(),
-            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::zero()),
+            Ok(fac) => fac.sln_det_into(),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => {
+                // The determinant is zero.
+                Ok((A::zero(), A::Real::neg_infinity()))
+            }
             Err(err) => Err(err),
         }
     }

src/solveh.rs

@@ -254,20 +254,53 @@ where
 /// An interface for calculating determinants of Hermitian (or real symmetric) matrix refs.
 pub trait DeterminantH {
     type Output;
+    type SignLnOutput;
 
     /// Computes the determinant of the Hermitian (or real symmetric) matrix.
     fn deth(&self) -> Self::Output;
+
+    /// Computes the `(sign, natural_log)` of the determinant of the Hermitian
+    /// (or real symmetric) matrix.
+    ///
+    /// The `natural_log` is the natural logarithm of the absolute value of the
+    /// determinant. If the determinant is zero, `sign` is 0 and `natural_log`
+    /// is negative infinity.
+    ///
+    /// To obtain the determinant, you can compute `sign * natural_log.exp()`
+    /// or just call `.deth()` instead.
+    ///
+    /// This method is more robust than `.deth()` to very small or very large
+    /// determinants since it returns the natural logarithm of the determinant
+    /// rather than the determinant itself.
+    fn sln_deth(&self) -> Self::SignLnOutput;
 }
 
 /// An interface for calculating determinants of Hermitian (or real symmetric) matrices.
 pub trait DeterminantHInto {
     type Output;
+    type SignLnOutput;
 
     /// Computes the determinant of the Hermitian (or real symmetric) matrix.
     fn deth_into(self) -> Self::Output;
+
+    /// Computes the `(sign, natural_log)` of the determinant of the Hermitian
+    /// (or real symmetric) matrix.
+    ///
+    /// The `natural_log` is the natural logarithm of the absolute value of the
+    /// determinant. If the determinant is zero, `sign` is 0 and `natural_log`
+    /// is negative infinity.
+    ///
+    /// To obtain the determinant, you can compute `sign * natural_log.exp()`
+    /// or just call `.deth_into()` instead.
+    ///
+    /// This method is more robust than `.deth_into()` to very small or very
+    /// large determinants since it returns the natural logarithm of the
+    /// determinant rather than the determinant itself.
+    fn sln_deth_into(self) -> Self::SignLnOutput;
 }
 
-fn bk_det<P, S, A>(uplo: UPLO, ipiv_iter: P, a: &ArrayBase<S, Ix2>) -> A::Real
+/// Returns the sign and natural log of the determinant.
+fn bk_sln_det<P, S, A>(uplo: UPLO, ipiv_iter: P, a: &ArrayBase<S, Ix2>) -> (A::Real, A::Real)
 where
     P: Iterator<Item = i32>,
     S: Data<Elem = A>,
@@ -310,7 +343,7 @@ where
             ipiv_enum.next();
         }
     }
-    sign * ln_det.exp()
+    (sign, ln_det)
 }
 
 impl<A, S> DeterminantH for BKFactorized<S>
@@ -319,9 +352,15 @@ where
     S: Data<Elem = A>,
 {
     type Output = A::Real;
+    type SignLnOutput = (A::Real, A::Real);
 
     fn deth(&self) -> A::Real {
-        bk_det(UPLO::Upper, self.ipiv.iter().cloned(), &self.a)
+        let (sign, ln_det) = self.sln_deth();
+        sign * ln_det.exp()
+    }
+
+    fn sln_deth(&self) -> (A::Real, A::Real) {
+        bk_sln_det(UPLO::Upper, self.ipiv.iter().cloned(), &self.a)
     }
 }
 
@@ -331,9 +370,15 @@ where
     S: Data<Elem = A>,
 {
     type Output = A::Real;
+    type SignLnOutput = (A::Real, A::Real);
 
     fn deth_into(self) -> A::Real {
-        bk_det(UPLO::Upper, self.ipiv.into_iter(), &self.a)
+        let (sign, ln_det) = self.sln_deth_into();
+        sign * ln_det.exp()
+    }
+
+    fn sln_deth_into(self) -> (A::Real, A::Real) {
+        bk_sln_det(UPLO::Upper, self.ipiv.into_iter(), &self.a)
     }
 }
 
@@ -343,11 +388,20 @@ where
     S: Data<Elem = A>,
 {
     type Output = Result<A::Real>;
+    type SignLnOutput = Result<(A::Real, A::Real)>;
 
     fn deth(&self) -> Result<A::Real> {
+        let (sign, ln_det) = self.sln_deth()?;
+        Ok(sign * ln_det.exp())
+    }
+
+    fn sln_deth(&self) -> Result<(A::Real, A::Real)> {
         match self.factorizeh() {
-            Ok(fac) => Ok(fac.deth()),
-            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::Real::zero()),
+            Ok(fac) => Ok(fac.sln_deth()),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => {
+                // Determinant is zero.
+                Ok((A::Real::zero(), A::Real::neg_infinity()))
+            }
             Err(err) => Err(err),
         }
     }
@@ -359,11 +413,20 @@ where
     S: DataMut<Elem = A>,
 {
     type Output = Result<A::Real>;
+    type SignLnOutput = Result<(A::Real, A::Real)>;
 
     fn deth_into(self) -> Result<A::Real> {
+        let (sign, ln_det) = self.sln_deth_into()?;
+        Ok(sign * ln_det.exp())
+    }
+
+    fn sln_deth_into(self) -> Result<(A::Real, A::Real)> {
         match self.factorizeh_into() {
-            Ok(fac) => Ok(fac.deth_into()),
-            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::Real::zero()),
+            Ok(fac) => Ok(fac.sln_deth_into()),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => {
+                // Determinant is zero.
+                Ok((A::Real::zero(), A::Real::neg_infinity()))
+            }
             Err(err) => Err(err),
         }
     }