Merge pull request #331 from alphaville/feature/projection-affine-space

Projection on affine space

CHANGELOG.md

@@ -10,7 +10,11 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo
 <!-- ---------------------
       Not released
      --------------------- -->
-## Unreleased
+## [v0.9.0] - Unreleased
+
+### Added
+
+- Implementation of `AffineSpace`
 
 ### Fixed
 
@@ -275,6 +279,7 @@ This is a breaking API change.
      --------------------- -->
 
 <!-- Releases -->
+[v0.9.0]: https://github.com/alphaville/optimization-engine/compare/v0.8.1...v0.9.0
 [v0.8.1]: https://github.com/alphaville/optimization-engine/compare/v0.8.0...v0.8.1
 [v0.8.0]: https://github.com/alphaville/optimization-engine/compare/v0.7.7...v0.8.0
 [v0.7.7]: https://github.com/alphaville/optimization-engine/compare/v0.7.6...v0.7.7 

Cargo.toml

@@ -42,7 +42,7 @@ homepage = "https://alphaville.github.io/optimization-engine/"
 repository = "https://github.com/alphaville/optimization-engine"
 
 # Version of this crate (SemVer)
-version = "0.8.1"
+version = "0.9.0"
 
 edition = "2018"
 
@@ -99,6 +99,9 @@ rpmalloc = { version = "0.2.0", features = [
 [target.'cfg(not(target_env = "msvc"))'.dependencies]
 jemallocator = { version = "0.5.0", optional = true }
 
+# Least squares solver
+ndarray = { version = "0.15", features = ["approx"] }
+modcholesky = "0.1.3"
 
 # --------------------------------------------------------------------------
 # F.E.A.T.U.R.E.S.

appveyor.yml

@@ -68,4 +68,7 @@ build: false
 #directly or perform other testing commands. Rust will automatically be placed in the PATH
 # environment variable.
 test_script:
+  - cargo add ndarray --features approx
+  - cargo add modcholesky
+  - cargo build
   - cargo test --verbose %cargoflags%

docs/openrust-basic.md

@@ -62,11 +62,15 @@ Constraints implement the namesake trait, [`Constraint`]. Implementations of [`C
 
 | Constraint           | Explanation                                          |
 |----------------------|------------------------------------------------------|
-| [`Ball2`]            | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_2 \leq r\\}$         |
-| [`BallInf`]          | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_\infty \leq r\\}$         |
+| [`AffineSpace`]      | $U {}={} \\{u\in\mathbb{R}^n : Au = b\\}$         |
+| [`Ball1`]            | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_1 \leq r\\}$  |
+| [`Ball2`]            | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_2 \leq r\\}$  |
+| [`Sphere2`]          | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_2 = r\\}$     |
+| [`BallInf`]          | $U {}={} \\{u\in\mathbb{R}^n : \Vert u-u^0\Vert_\infty \leq r\\}$ |
 | [`Halfspace`]        | $U {}={} \\{u\in\mathbb{R}^n : \langle c, u\rangle \leq b\\}$ |
 | [`Hyperplane`]       | $U {}={} \\{u\in\mathbb{R}^n : \langle c, u\rangle {}={} b\\}$ |
 | [`Rectangle`]        | $U {}={} \\{u\in\mathbb{R}^n : u_{\min} \leq u \leq u_{\max}\\}$ |
+| [`Simplex`]          | $U {}={} \\{u \in \mathbb{R}^n {}:{} u \geq 0, \sum_i u_i = \alpha\\}$ |
 | [`NoConstraints`]    | $U {}={} \mathbb{R}^n$                                   |
 | [`FiniteSet`]        | $U {}={} \\{u^{(1)}, u^{(2)},\ldots,u^{(N)}\\}$          |
 | [`SecondOrderCone`]  | $U {}={} \\{u=(z,t), t\in\mathbb{R}, \Vert{}z{}\Vert \leq \alpha t\\}$ |
@@ -353,6 +357,10 @@ the imposition of a maximum allowed duration, the exit status will be
 <!-- Links -->
 
 [`Constraint`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/trait.Constraint.html
+[`Simplex`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Simplex.html
+[`AffineSpace`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.AffineSpace.html
+[`Sphere2`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Sphere2.html
+[`Ball1`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Ball1.html
 [`Ball2`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Ball2.html
 [`BallInf`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.BallInf.html
 [`Halfspace`]: https://docs.rs/optimization_engine/*/optimization_engine/constraints/struct.Halfspace.html

src/constraints/affine_space.rs

@@ -0,0 +1,148 @@
+use super::Constraint;
+
+extern crate modcholesky;
+extern crate ndarray;
+
+use modcholesky::ModCholeskySE99;
+use ndarray::{Array1, Array2, ArrayBase, Dim, OwnedRepr};
+
+type OpenMat<T> = ArrayBase<OwnedRepr<T>, Dim<[usize; 2]>>;
+type OpenVec<T> = ArrayBase<OwnedRepr<T>, Dim<[usize; 1]>>;
+
+#[derive(Clone)]
+/// An affine space here is defined as the set of solutions of a linear equation, $Ax = b$,
+/// that is, $E=\\{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>,
+    l: OpenMat<f64>,
+    p: OpenVec<usize>,
+    n_rows: usize,
+    n_cols: usize,
+}
+
+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$
+    ///
+    /// ## Arguments
+    ///
+    /// - `a`: matrix $A$, row-wise data
+    /// - `b`: vector $b$
+    ///
+    /// ## Returns
+    /// New Affine Space structure
+    ///
+    pub fn new(a: Vec<f64>, b: Vec<f64>) -> Self {
+        // Infer dimensions of A and b
+        let n_rows = b.len();
+        let n_elements_a = a.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 = Array2::from_shape_vec((n_rows, n_cols), a).unwrap();
+        let b_vec = Array1::from_shape_vec((n_rows,), b).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;
+
+        // Construct and return new AffineSpace structure
+        AffineSpace {
+            a_mat,
+            b_vec,
+            l,
+            p,
+            n_rows,
+            n_cols,
+        }
+    }
+}
+
+impl Constraint for AffineSpace {
+    /// Projection onto the set $E = \\{x: Ax = b\\}$, which is computed by
+    /// $$P_E(x) = x - A^\intercal z(x),$$
+    /// where $z$ is the solution of the linear system
+    /// $$(AA^\intercal)z = Ax - b,$$
+    /// which has a unique solution provided $A$ has full row rank. The linear system
+    /// is solved by computing the Cholesky factorisation of $AA^\intercal$, which is
+    /// done using [`modcholesky`](https://crates.io/crates/modcholesky).
+    ///
+    /// ## Arguments
+    ///
+    /// - `x`: The given vector $x$ is updated with the projection on the set
+    ///
+    /// ## Example
+    ///
+    /// Consider the set $X = \\{x \in \mathbb{R}^4 :Ax = b\\}$, with $A\in\mathbb{R}^{3\times 4}$
+    /// being the matrix
+    /// $$A = \begin{bmatrix}0.5 & 0.1& 0.2& -0.3\\\\ -0.6& 0.3& 0 & 0.5 \\\\ 1.0& 0.1& -1& -0.4\end{bmatrix},$$
+    /// and $b$ being the vector
+    /// $$b = \begin{bmatrix}1 \\\\ 2 \\\\ -0.5\end{bmatrix}.$$
+    ///
+    /// ```rust
+    /// use optimization_engine::constraints::*;
+    ///
+    /// let a = vec![0.5, 0.1, 0.2, -0.3, -0.6, 0.3, 0., 0.5, 1.0, 0.1, -1.0, -0.4,];
+    /// let b = vec![1., 2., -0.5];
+    /// let affine_set = AffineSpace::new(a, b);
+    /// let mut x = [1., -2., -0.3, 0.5];
+    /// affine_set.project(&mut x);
+    /// ```
+    ///
+    /// The result is stored in `x` and it can be verified that $Ax = b$.
+    fn project(&self, x: &mut [f64]) {
+        let m = self.n_rows;
+        let n = self.n_cols;
+        let chol = &self.l;
+        let perm = &self.p;
+
+        assert!(x.len() == n, "x has wrong dimension");
+        let x_vec = x.to_vec();
+        let x_arr = Array1::from_shape_vec((n,), x_vec).unwrap();
+        let ax = self.a_mat.dot(&x_arr);
+        let err = ax - &self.b_vec;
+
+        // Step 1: Solve Ly = b(P)
+        // TODO: Make `y` into an attribute; however, to do this, we need to change
+        // &self to &mut self, which will require a mild refactoring
+        let mut y = vec![0.; m];
+        for i in 0..m {
+            y[i] = err[perm[i]];
+            for j in 0..i {
+                y[i] -= chol[(i, j)] * y[j];
+            }
+            y[i] /= chol[(i, i)];
+        }
+
+        // Step 2: Solve L'z(P) = y
+        let mut z = vec![0.; m];
+        for i in 1..=m {
+            z[perm[m - i]] = y[m - i];
+            for j in 1..i {
+                z[perm[m - i]] -= chol[(m - j, m - i)] * z[perm[m - j]];
+            }
+            z[perm[m - i]] /= chol[(m - i, m - i)];
+        }
+
+        // Step 3: Determine A' * z
+        let z_arr = Array1::from_shape_vec((self.n_rows,), z).unwrap();
+        let w = self.a_mat.t().dot(&z_arr);
+
+        // Step 4: x <-- x - A'(AA')\(Ax-b)
+        x.iter_mut().zip(w.iter()).for_each(|(xi, wi)| *xi -= wi);
+    }
+
+    /// Affine sets are convex.
+    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,4 +1,5 @@
 use super::*;
+use modcholesky::ModCholeskySE99;
 use rand;
 
 #[test]
@@ -872,3 +873,71 @@ fn t_sphere2_center_projection_of_center() {
 fn t_ball1_alpha_negative() {
     let _ = Ball1::new(None, -1.);
 }
+
+#[test]
+fn t_affine_space() {
+    let a = vec![
+        0.5, 0.1, 0.2, -0.3, -0.6, 0.3, 0., 0.5, 1.0, 0.1, -1.0, -0.4,
+    ];
+    let b = vec![1., 2., -0.5];
+    let affine_set = AffineSpace::new(a, b);
+    let mut x = [1., -2., -0.3, 0.5];
+    affine_set.project(&mut x);
+    let x_correct = [
+        1.888564346697095,
+        5.629857182200888,
+        1.796204902230790,
+        2.888362906715977,
+    ];
+    unit_test_utils::assert_nearly_equal_array(
+        &x_correct,
+        &x,
+        1e-10,
+        1e-12,
+        "projection on affine set is wrong",
+    );
+}
+
+#[test]
+fn t_affine_space_larger() {
+    let a = vec![
+        1.0f64, 1., 1., 0., 0., 0., 1., 1., 1., 0., 0., 0., 1., 1., 1., -1., 4., -1., 0., 2.,
+    ];
+    let b = vec![1., -2., 3., 4.];
+    let affine_set = AffineSpace::new(a, b);
+    let mut x = [10., 11., -9., 4., 5.];
+    affine_set.project(&mut x);
+    let x_correct = [
+        9.238095238095237,
+        -0.714285714285714,
+        -7.523809523809524,
+        6.238095238095238,
+        4.285714285714288,
+    ];
+    unit_test_utils::assert_nearly_equal_array(
+        &x_correct,
+        &x,
+        1e-10,
+        1e-12,
+        "projection on affine set is wrong",
+    );
+}
+
+#[test]
+fn t_affine_space_single_row() {
+    let a = vec![1., 1., 1., 1.];
+    let b = vec![1.];
+    let affine_set = AffineSpace::new(a, b);
+    let mut x = [5., 6., 10., 25.];
+    affine_set.project(&mut x);
+    let s = x.iter().sum();
+    unit_test_utils::assert_nearly_equal(1., s, 1e-12, 1e-14, "wrong sum");
+}
+
+#[test]
+#[should_panic]
+fn t_affine_space_wrong_dimensions() {
+    let a = vec![0.5, 0.1, 0.2, -0.3, -0.6, 0.3, 0., 0.5, 1.0, 0.1, -1.0];
+    let b = vec![1., 2., -0.5];
+    let _ = AffineSpace::new(a, b);
+}