Merge pull request #85 from jturner314/add-determinant

Add determinant methods in solve module

Author: termoshtt

src/lapack_traits/solve.rs

@@ -10,6 +10,14 @@ use super::{Pivot, Transpose, into_result};
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`
 pub trait Solve_: Sized {
+    /// Computes the LU factorization of a general `m x n` matrix `a` using
+    /// partial pivoting with row interchanges.
+    ///
+    /// If the result matches `Err(LinalgError::Lapack(LapackError {
+    /// return_code )) if return_code > 0`, then `U[(return_code-1,
+    /// return_code-1)]` 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.
     unsafe fn lu(MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
     unsafe fn inv(MatrixLayout, a: &mut [Self], &Pivot) -> Result<()>;
     unsafe fn solve(MatrixLayout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;

src/layout.rs

@@ -75,6 +75,8 @@ pub trait AllocatedArray {
     type Elem;
     fn layout(&self) -> Result<MatrixLayout>;
     fn square_layout(&self) -> Result<MatrixLayout>;
+    /// Returns Ok iff the matrix is square (without computing the layout).
+    fn ensure_square(&self) -> Result<()>;
     fn as_allocated(&self) -> Result<&[Self::Elem]>;
 }
 
@@ -110,6 +112,14 @@ where
         }
     }
 
+    fn ensure_square(&self) -> Result<()> {
+        if self.is_square() {
+            Ok(())
+        } else {
+            Err(NotSquareError::new(self.rows() as i32, self.cols() as i32).into())
+        }
+    }
+
     fn as_allocated(&self) -> Result<&[A]> {
         Ok(self.as_slice_memory_order().ok_or(MemoryContError::new())?)
     }

src/solve.rs

@@ -330,3 +330,92 @@ where
         f.inv_into()
     }
 }
+
+/// An interface for calculating determinants of matrix refs.
+pub trait Determinant<A: Scalar> {
+    /// Computes the determinant of the matrix.
+    fn det(&self) -> Result<A>;
+}
+
+/// An interface for calculating determinants of matrices.
+pub trait DeterminantInto<A: Scalar> {
+    /// Computes the determinant of the matrix.
+    fn det_into(self) -> Result<A>;
+}
+
+fn lu_det<'a, A, P, U>(ipiv_iter: P, u_diag_iter: U) -> A
+where
+    A: Scalar,
+    P: Iterator<Item = i32>,
+    U: Iterator<Item = &'a A>,
+{
+    let pivot_sign = if ipiv_iter
+        .enumerate()
+        .filter(|&(i, pivot)| pivot != i as i32 + 1)
+        .count() % 2 == 0
+    {
+        A::one()
+    } 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()
+}
+
+impl<A, S> Determinant<A> for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn det(&self) -> Result<A> {
+        self.a.ensure_square()?;
+        Ok(lu_det(self.ipiv.iter().cloned(), self.a.diag().iter()))
+    }
+}
+
+impl<A, S> DeterminantInto<A> for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn det_into(self) -> Result<A> {
+        self.a.ensure_square()?;
+        Ok(lu_det(self.ipiv.into_iter(), self.a.into_diag().iter()))
+    }
+}
+
+impl<A, S> Determinant<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn det(&self) -> Result<A> {
+        self.ensure_square()?;
+        match self.factorize() {
+            Ok(fac) => fac.det(),
+            Err(LinalgError::Lapack(LapackError { return_code })) if return_code > 0 => Ok(A::zero()),
+            Err(err) => Err(err),
+        }
+    }
+}
+
+impl<A, S> DeterminantInto<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    fn det_into(self) -> Result<A> {
+        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()),
+            Err(err) => Err(err),
+        }
+    }
+}

src/types.rs

@@ -22,6 +22,7 @@ pub use num_complex::Complex64 as c64;
 /// - [abs_sqr](trait.Absolute.html#tymethod.abs_sqr)
 /// - [sqrt](trait.SquareRoot.html#tymethod.sqrt)
 /// - [exp](trait.Exponential.html#tymethod.exp)
+/// - [ln](trait.NaturalLogarithm.html#tymethod.ln)
 /// - [conj](trait.Conjugate.html#tymethod.conj)
 /// - [randn](trait.RandNormal.html#tymethod.randn)
 ///
@@ -33,6 +34,7 @@ pub trait Scalar
     + Absolute
     + SquareRoot
     + Exponential
+    + NaturalLogarithm
     + Conjugate
     + RandNormal
     + Neg<Output = Self>
@@ -118,6 +120,11 @@ pub trait Exponential {
     fn exp(&self) -> Self;
 }
 
+/// Define `ln()` more generally
+pub trait NaturalLogarithm {
+    fn ln(&self) -> Self;
+}
+
 /// Complex conjugate value
 pub trait Conjugate: Copy {
     fn conj(self) -> Self;
@@ -207,6 +214,18 @@ impl Exponential for $complex {
     }
 }
 
+impl NaturalLogarithm for $real {
+    fn ln(&self) -> Self {
+        Float::ln(*self)
+    }
+}
+
+impl NaturalLogarithm for $complex {
+    fn ln(&self) -> Self {
+        Complex::ln(self)
+    }
+}
+
 impl Conjugate for $real {
     fn conj(self) -> Self {
         self

tests/solve.rs

@@ -0,0 +1,141 @@
+extern crate ndarray;
+#[macro_use]
+extern crate ndarray_linalg;
+extern crate num_traits;
+
+use ndarray::*;
+use ndarray_linalg::*;
+use num_traits::{One, Zero};
+
+/// Returns the matrix with the specified `row` and `col` removed.
+fn matrix_minor<A, S>(a: ArrayBase<S, Ix2>, (row, col): (usize, usize)) -> Array2<A>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    let mut select_rows = (0..a.rows()).collect::<Vec<_>>();
+    select_rows.remove(row);
+    let mut select_cols = (0..a.cols()).collect::<Vec<_>>();
+    select_cols.remove(col);
+    a.select(Axis(0), &select_rows).select(
+        Axis(1),
+        &select_cols,
+    )
+}
+
+/// Computes the determinant of matrix `a`.
+///
+/// Note: This implementation is written to be clearly correct so that it's
+/// useful for verification, but it's very inefficient.
+fn det_naive<A, S>(a: ArrayBase<S, Ix2>) -> A
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    assert_eq!(a.rows(), a.cols());
+    match a.cols() {
+        0 => A::one(),
+        1 => a[(0, 0)],
+        cols => {
+            (0..cols)
+                .map(|col| {
+                    let sign = if col % 2 == 0 { A::one() } else { -A::one() };
+                    sign * a[(0, col)] * det_naive(matrix_minor(a.view(), (0, col)))
+                })
+                .fold(A::zero(), |sum, subdet| sum + subdet)
+        }
+    }
+}
+
+#[test]
+fn det_empty() {
+    macro_rules! det_empty {
+        ($elem:ty) => {
+            let a: Array2<$elem> = Array2::zeros((0, 0));
+            assert_eq!(a.factorize().unwrap().det().unwrap(), One::one());
+            assert_eq!(a.factorize().unwrap().det_into().unwrap(), One::one());
+            assert_eq!(a.det().unwrap(), One::one());
+            assert_eq!(a.det_into().unwrap(), One::one());
+        }
+    }
+    det_empty!(f64);
+    det_empty!(f32);
+    det_empty!(c64);
+    det_empty!(c32);
+}
+
+#[test]
+fn det_zero() {
+    macro_rules! det_zero {
+        ($elem:ty) => {
+            let a: Array2<$elem> = Array2::zeros((1, 1));
+            assert_eq!(a.det().unwrap(), Zero::zero());
+            assert_eq!(a.det_into().unwrap(), Zero::zero());
+        }
+    }
+    det_zero!(f64);
+    det_zero!(f32);
+    det_zero!(c64);
+    det_zero!(c32);
+}
+
+#[test]
+fn det_zero_nonsquare() {
+    macro_rules! det_zero_nonsquare {
+        ($elem:ty, $shape:expr) => {
+            let a: Array2<$elem> = Array2::zeros($shape);
+            assert!(a.det().is_err());
+            assert!(a.det_into().is_err());
+        }
+    }
+    for &shape in &[(1, 2).into_shape(), (1, 2).f()] {
+        det_zero_nonsquare!(f64, shape);
+        det_zero_nonsquare!(f32, shape);
+        det_zero_nonsquare!(c64, shape);
+        det_zero_nonsquare!(c32, shape);
+    }
+}
+
+#[test]
+fn det() {
+    macro_rules! det {
+        ($elem:ty, $shape:expr, $rtol:expr) => {
+            let a: Array2<$elem> = random($shape);
+            println!("a = \n{:?}", a);
+            let det = det_naive(a.view());
+            assert_rclose!(a.factorize().unwrap().det().unwrap(), det, $rtol);
+            assert_rclose!(a.factorize().unwrap().det_into().unwrap(), det, $rtol);
+            assert_rclose!(a.det().unwrap(), det, $rtol);
+            assert_rclose!(a.det_into().unwrap(), det, $rtol);
+        }
+    }
+    for rows in 1..5 {
+        for &shape in &[(rows, rows).into_shape(), (rows, rows).f()] {
+            det!(f64, shape, 1e-9);
+            det!(f32, shape, 1e-4);
+            det!(c64, shape, 1e-9);
+            det!(c32, shape, 1e-4);
+        }
+    }
+}
+
+#[test]
+fn det_nonsquare() {
+    macro_rules! det_nonsquare {
+        ($elem:ty, $shape:expr) => {
+            let a: Array2<$elem> = random($shape);
+            assert!(a.factorize().unwrap().det().is_err());
+            assert!(a.factorize().unwrap().det_into().is_err());
+            assert!(a.det().is_err());
+            assert!(a.det_into().is_err());
+        }
+    }
+    for &dims in &[(1, 0), (1, 2), (2, 1), (2, 3)] {
+        for &shape in &[dims.clone().into_shape(), dims.clone().f()] {
+            det_nonsquare!(f64, shape);
+            det_nonsquare!(f32, shape);
+            det_nonsquare!(c64, shape);
+            det_nonsquare!(c32, shape);
+        }
+    }
+}