Add random_hpd

Author: termoshtt

src/generate.rs

@@ -3,9 +3,23 @@ use ndarray::*;
 use rand::*;
 use std::ops::*;
 
+use super::layout::*;
 use super::types::*;
 use super::error::*;
 
+pub fn conjugate<A, Si, So>(a: &ArrayBase<Si, Ix2>) -> ArrayBase<So, Ix2>
+    where A: Conjugate,
+          Si: Data<Elem = A>,
+          So: DataOwned<Elem = A> + DataMut
+{
+    let mut a = replicate(&a.t());
+    for val in a.iter_mut() {
+        *val = Conjugate::conj(*val);
+    }
+    a
+}
+
+/// Random square matrix
 pub fn random_square<A, S>(n: usize) -> ArrayBase<S, Ix2>
     where A: Rand,
           S: DataOwned<Elem = A>
@@ -15,20 +29,31 @@ pub fn random_square<A, S>(n: usize) -> ArrayBase<S, Ix2>
     ArrayBase::from_shape_vec((n, n), v).unwrap()
 }
 
+/// Random Hermite matrix
 pub fn random_hermite<A, S>(n: usize) -> ArrayBase<S, Ix2>
     where A: Rand + Conjugate + Add<Output = A>,
           S: DataOwned<Elem = A> + DataMut
 {
     let mut a = random_square(n);
     for i in 0..n {
         a[(i, i)] = a[(i, i)] + Conjugate::conj(a[(i, i)]);
-        for j in i..n {
-            a[(i, j)] = a[(j, i)];
+        for j in (i + 1)..n {
+            a[(i, j)] = Conjugate::conj(a[(j, i)])
         }
     }
     a
 }
 
+/// Random Hermite Positive-definite matrix
+pub fn random_hpd<A, S>(n: usize) -> ArrayBase<S, Ix2>
+    where A: Rand + Conjugate + LinalgScalar,
+          S: DataOwned<Elem = A> + DataMut
+{
+    let a: Array2<A> = random_square(n);
+    let ah: Array2<A> = conjugate(&a);
+    replicate(&ah.dot(&a))
+}
+
 /// construct matrix from diag
 pub fn from_diag<A>(d: &[A]) -> Array2<A>
     where A: LinalgScalar

src/prelude.rs

@@ -3,6 +3,8 @@ pub use matrix::Matrix;
 pub use square::SquareMatrix;
 pub use triangular::*;
 pub use util::*;
+pub use types::*;
+pub use generate::*;
 pub use assert::*;
 
 pub use qr::*;

tests/cholesky.rs

@@ -1,24 +1,14 @@
 
-extern crate rand_extra;
 extern crate ndarray;
-extern crate ndarray_rand;
 #[macro_use]
 extern crate ndarray_linalg;
 
-use rand_extra::*;
 use ndarray::*;
-use ndarray_rand::RandomExt;
 use ndarray_linalg::prelude::*;
 
-pub fn random_hermite(n: usize) -> Array<f64, Ix2> {
-    let r_dist = RealNormal::new(0., 1.);
-    let a = Array::<f64, _>::random((n, n), r_dist);
-    a.dot(&a.t())
-}
-
 #[test]
 fn cholesky() {
-    let a = random_hermite(3);
+    let a: Array2<f64> = random_hpd(3);
     println!("a = \n{:?}", a);
     let c: Array2<_> = (&a).cholesky(UPLO::Upper).unwrap();
     println!("c = \n{:?}", c);
@@ -28,7 +18,7 @@ fn cholesky() {
 
 #[test]
 fn cholesky_t() {
-    let a = random_hermite(3);
+    let a: Array2<f64> = random_hpd(3);
     println!("a = \n{:?}", a);
     let c: Array2<_> = (&a).cholesky(UPLO::Upper).unwrap();
     println!("c = \n{:?}", c);