AffineSpace: computation of Cholesky

src/constraints/affine_space.rs

@@ -0,0 +1,68 @@
+use super::Constraint;
+use modcholesky::ModCholeskySE99;
+
+type OpenMat<T> = ndarray::ArrayBase<ndarray::OwnedRepr<T>, ndarray::Dim<[usize; 2]>>;
+type OpenVec<T> = ndarray::ArrayBase<ndarray::OwnedRepr<T>, ndarray::Dim<[usize; 1]>>;
+
+#[derive(Clone)]
+/// An affine space here is defined as the set of solutions of a linear equation, $Ax = b$,
+/// that is, $\{x\in\mathbb{R}^n: Ax = b\}$, which is an affine space. It is assumed that
+/// the matrix $AA^\intercal$ is full-rank.
+pub struct AffineSpace {
+    a_mat: OpenMat<f64>,
+    b_vec: OpenVec<f64>,
+    a_times_a_t: OpenMat<f64>,
+    l: OpenMat<f64>,
+    p: OpenVec<usize>,
+    e: OpenVec<f64>,
+}
+
+impl AffineSpace {
+    /// Construct a new affine space given the matrix $A\in\mathbb{R}^{m\times n}$ and
+    /// the vector $b\in\mathbb{R}^m$
+    pub fn new(a_data: Vec<f64>, b_data: Vec<f64>) -> Self {
+        // Infer dimensions of A and b
+        let n_rows = b_data.len();
+        let n_elements_a = a_data.len();
+        assert!(
+            n_elements_a % n_rows == 0,
+            "A and b have incompatible dimensions"
+        );
+        let n_cols = n_elements_a / n_rows;
+        // Cast A and b as ndarray structures
+        let a_mat = ndarray::Array2::from_shape_vec((n_rows, n_cols), a_data).unwrap();
+        let b_vec = ndarray::Array1::from_shape_vec((n_rows,), b_data).unwrap();
+        // Compute a permuted Cholesky factorisation of AA'; in particular, we are looking for a
+        // minimum-norm matrix E, a permulation matrix P and a lower-trianular L, such that
+        //  P(AA' + E)P' = LL'
+        // and E should be 0 if A is full rank.
+        let a_times_a_t = a_mat.dot(&a_mat.t());
+        let res = a_times_a_t.mod_cholesky_se99();
+        let l = res.l;
+        let p = res.p;
+        let e = res.e;
+        // Construct and return new AffineSpace structure
+        AffineSpace {
+            a_mat,
+            b_vec,
+            a_times_a_t,
+            l,
+            p,
+            e,
+        }
+    }
+}
+
+impl Constraint for AffineSpace {
+    fn project(&self, x: &mut [f64]) {
+        let x_vec = x.to_vec();
+        let x_arr = ndarray::Array1::from_shape_vec((x_vec.len(),), x_vec).unwrap();
+        let ax = self.a_mat.dot(&x_arr);
+        let err = ax - &self.b_vec;
+        println!("err = {:?}", err);
+    }
+
+    fn is_convex(&self) -> bool {
+        true
+    }
+}

src/constraints/mod.rs

@@ -8,6 +8,7 @@
 //!
 //! [`Constraint`]: trait.Constraint.html
 
+mod affine_space;
 mod ball1;
 mod ball2;
 mod ballinf;
@@ -22,6 +23,7 @@ mod soc;
 mod sphere2;
 mod zero;
 
+pub use affine_space::AffineSpace;
 pub use ball1::Ball1;
 pub use ball2::Ball2;
 pub use ballinf::BallInf;

src/constraints/tests.rs

@@ -1,5 +1,5 @@
 use super::*;
-use modcholesky::ModCholeskyGMW81;
+use modcholesky::ModCholeskySE99;
 use rand;
 
 #[test]
@@ -876,26 +876,34 @@ fn t_ball1_alpha_negative() {
 
 #[test]
 fn t_affine_space() {
-    // let m = 4;
-    // let n = 4;
-    // let mut a = Array2::<f64>::zeros((m, n));
-    // a[(0, 0)] = 1.0;
-    // a[(1, 1)] = 1.0;
-    // a[(2, 2)] = 1.0;
-    // a[(3, 3)] = 1.0;
-    // let b = Array1::<f64>::zeros((n,));
-    // let aa = a.t().dot(&a);
-    let a_data = [
-        [1890.3, -1705.6, -315.8, 3000.3],
-        [-1705.6, 1538.3, 284.9, -2706.6],
-        [-315.8, 284.9, 52.5, -501.2],
-        [3000.3, -2706.6, -501.2, 4760.8],
-    ];
-    let a: ndarray::Array2<f64> = ndarray::arr2(&a_data);
+    let data = [10., 2., 2., 15.].to_vec();
+    let datb = vec![1., 2.];
+    // let arr: ndarray::ArrayBase<ndarray::OwnedRepr<f64>, ndarray::Dim<[usize; 2]>> =
+    //     ndarray::Array2::from_shape_vec((nrows, ncols), data).unwrap();
+
+    // let asq = arr.dot(&arr.tr());
+    // println!("A2 = {:?}", asq);
+
+    // let arrb: ndarray::ArrayBase<ndarray::OwnedRepr<f64>, ndarray::Dim<[usize; 1]>> = ndarray::Array1::from_shape_vec((nrows,), datb).unwrap();
+
+    let aff = AffineSpace::new(data, datb);
+    let mut xx = [1., 1.];
+    aff.project(&mut xx);
+
+    // let a: ndarray::Array2<f64> = ndarray::arr2(&a_data);
+    // let a_cp = a.clone();
+    // let mut a_sq = a.t().dot(&a_cp);
+
+    // println!("A'A = {:?}", a_sq);
+
+    // // Perform modified Cholesky decomposition
+    // // The `Decomposition` struct holds L, E and P
+    // let res = a_sq.mod_cholesky_se99();
+    // let x = a_sq[(res.p, 0)];
+    // let ll = res.l.dot(&res.l.t());
+    // let error = ll;
 
-    // Perform modified Cholesky decomposition
-    // The `Decomposition` struct holds L, E and P
-    let res = a.mod_cholesky_gmw81();
+    // println!("{}", error);
 
-    println!("{}", res.l);
+    // println!("{}", res.e);
 }