Add .deth*() for determinants of Hermitian matrices

Author: jturner314

src/solveh.rs

@@ -50,6 +50,7 @@
 //! ```
 
 use ndarray::*;
+use num_traits::{Float, One, Zero};
 
 use super::convert::*;
 use super::error::*;
@@ -258,3 +259,121 @@ where
         f.invh_into()
     }
 }
+
+/// An interface for calculating determinants of Hermitian (or real symmetric) matrix refs.
+pub trait DeterminantH {
+    type Output;
+
+    /// Computes the determinant of the Hermitian (or real symmetric) matrix.
+    fn deth(&self) -> Self::Output;
+}
+
+/// An interface for calculating determinants of Hermitian (or real symmetric) matrices.
+pub trait DeterminantHInto {
+    type Output;
+
+    /// Computes the determinant of the Hermitian (or real symmetric) matrix.
+    fn deth_into(self) -> Self::Output;
+}
+
+fn bk_det<P, S, A>(uplo: UPLO, ipiv_iter: P, a: &ArrayBase<S, Ix2>) -> A::Real
+where
+    P: Iterator<Item = i32>,
+    S: Data<Elem = A>,
+    A: Scalar,
+{
+    let mut sign = A::Real::one();
+    let mut ln_det = A::Real::zero();
+    let mut ipiv_enum = ipiv_iter.enumerate();
+    while let Some((k, ipiv_k)) = ipiv_enum.next() {
+        debug_assert!(k < a.rows() && k < a.cols());
+        if ipiv_k > 0 {
+            // 1x1 block at k, must be real.
+            let elem = unsafe { a.uget((k, k)) }.real();
+            debug_assert_eq!(elem.imag(), Zero::zero());
+            sign = sign * elem.signum();
+            ln_det = ln_det + elem.abs().ln();
+        } else {
+            // 2x2 block at k..k+2.
+
+            // Upper left diagonal elem, must be real.
+            let upper_diag = unsafe { a.uget((k, k)) }.real();
+            debug_assert_eq!(upper_diag.imag(), Zero::zero());
+
+            // Lower right diagonal elem, must be real.
+            let lower_diag = unsafe { a.uget((k + 1, k + 1)) }.real();
+            debug_assert_eq!(lower_diag.imag(), Zero::zero());
+
+            // Off-diagonal elements, can be complex.
+            let off_diag = match uplo {
+                UPLO::Upper => unsafe { a.uget((k, k + 1)) },
+                UPLO::Lower => unsafe { a.uget((k + 1, k)) },
+            };
+
+            // Determinant of 2x2 block.
+            let block_det = upper_diag * lower_diag - off_diag.abs_sqr();
+            sign = sign * block_det.signum();
+            ln_det = ln_det + block_det.abs().ln();
+
+            // Skip the k+1 ipiv value.
+            ipiv_enum.next();
+        }
+    }
+    sign * ln_det.exp()
+}
+
+impl<A, S> DeterminantH for BKFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    type Output = A::Real;
+
+    fn deth(&self) -> A::Real {
+        bk_det(UPLO::Upper, self.ipiv.iter().cloned(), &self.a)
+    }
+}
+
+impl<A, S> DeterminantHInto for BKFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    type Output = A::Real;
+
+    fn deth_into(self) -> A::Real {
+        bk_det(UPLO::Upper, self.ipiv.into_iter(), &self.a)
+    }
+}
+
+impl<A, S> DeterminantH for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    type Output = Result<A::Real>;
+
+    fn deth(&self) -> Result<A::Real> {
+        match self.factorizeh() {
+            Ok(fac) => Ok(fac.deth()),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::Real::zero()),
+            Err(err) => Err(err),
+        }
+    }
+}
+
+impl<A, S> DeterminantHInto for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    type Output = Result<A::Real>;
+
+    fn deth_into(self) -> Result<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()),
+            Err(err) => Err(err),
+        }
+    }
+}

tests/solveh.rs

@@ -0,0 +1,99 @@
+extern crate ndarray;
+#[macro_use]
+extern crate ndarray_linalg;
+extern crate num_traits;
+
+use ndarray::*;
+use ndarray_linalg::*;
+use num_traits::{One, Zero};
+
+#[test]
+fn deth_empty() {
+    macro_rules! deth_empty {
+        ($elem:ty) => {
+            let a: Array2<$elem> = Array2::zeros((0, 0));
+            assert_eq!(a.factorizeh().unwrap().deth(), One::one());
+            assert_eq!(a.factorizeh().unwrap().deth_into(), One::one());
+            assert_eq!(a.deth().unwrap(), One::one());
+            assert_eq!(a.deth_into().unwrap(), One::one());
+        }
+    }
+    deth_empty!(f64);
+    deth_empty!(f32);
+    deth_empty!(c64);
+    deth_empty!(c32);
+}
+
+#[test]
+fn deth_zero() {
+    macro_rules! deth_zero {
+        ($elem:ty) => {
+            let a: Array2<$elem> = Array2::zeros((1, 1));
+            assert_eq!(a.deth().unwrap(), Zero::zero());
+            assert_eq!(a.deth_into().unwrap(), Zero::zero());
+        }
+    }
+    deth_zero!(f64);
+    deth_zero!(f32);
+    deth_zero!(c64);
+    deth_zero!(c32);
+}
+
+#[test]
+fn deth_zero_nonsquare() {
+    macro_rules! deth_zero_nonsquare {
+        ($elem:ty, $shape:expr) => {
+            let a: Array2<$elem> = Array2::zeros($shape);
+            assert!(a.deth().is_err());
+            assert!(a.deth_into().is_err());
+        }
+    }
+    for &shape in &[(1, 2).into_shape(), (1, 2).f()] {
+        deth_zero_nonsquare!(f64, shape);
+        deth_zero_nonsquare!(f32, shape);
+        deth_zero_nonsquare!(c64, shape);
+        deth_zero_nonsquare!(c32, shape);
+    }
+}
+
+#[test]
+fn deth() {
+    macro_rules! deth {
+        ($elem:ty, $rows:expr, $atol:expr) => {
+            let a: Array2<$elem> = random_hermite($rows);
+            println!("a = \n{:?}", a);
+            let det = a.eigvalsh(UPLO::Upper).unwrap().iter().product();
+            assert_aclose!(a.factorizeh().unwrap().deth(), det, $atol);
+            assert_aclose!(a.factorizeh().unwrap().deth_into(), det, $atol);
+            assert_aclose!(a.deth().unwrap(), det, $atol);
+            assert_aclose!(a.deth_into().unwrap(), det, $atol);
+        }
+    }
+    for rows in 1..6 {
+        deth!(f64, rows, 1e-9);
+        deth!(f32, rows, 1e-3);
+        deth!(c64, rows, 1e-9);
+        deth!(c32, rows, 1e-3);
+    }
+}
+
+#[test]
+fn deth_nonsquare() {
+    macro_rules! deth_nonsquare {
+        ($elem:ty, $shape:expr) => {
+            let a: Array2<$elem> = Array2::zeros($shape);
+            assert!(a.factorizeh().is_err());
+            assert!(a.factorizeh().is_err());
+            assert!(a.deth().is_err());
+            assert!(a.deth_into().is_err());
+        }
+    }
+    for &dims in &[(1, 0), (1, 2), (2, 1), (2, 3)] {
+        for &shape in &[dims.clone().into_shape(), dims.clone().f()] {
+            deth_nonsquare!(f64, shape);
+            deth_nonsquare!(f32, shape);
+            deth_nonsquare!(c64, shape);
+            deth_nonsquare!(c32, shape);
+        }
+    }
+}