diff --git a/CHANGELOG.md b/CHANGELOG.md index 64ac85cf..c8fa985d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo ### Added - Rust implementation of epigraph of squared Euclidean norm (constraint) +- Implementation of `AffineSpace` ### Fixed diff --git a/Cargo.toml b/Cargo.toml index 963f797b..d4ca6231 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -99,10 +99,16 @@ rpmalloc = { version = "0.2.0", features = [ [target.'cfg(not(target_env = "msvc"))'.dependencies] jemallocator = { version = "0.5.0", optional = true } + # computation of roots of cubic equation needed for the projection on the # epigraph of the squared Eucliean norm roots = "0.0.8" +# Least squares solver +ndarray = { version = "0.15", features = ["approx"] } +modcholesky = "0.1.3" + + # -------------------------------------------------------------------------- # F.E.A.T.U.R.E.S. # -------------------------------------------------------------------------- diff --git a/appveyor.yml b/appveyor.yml index 20bd0b3d..e216a306 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -69,5 +69,7 @@ build: false # environment variable. test_script: - cargo add roots + - cargo add ndarray --features approx + - cargo add modcholesky - cargo build - cargo test --verbose %cargoflags% diff --git a/docs/openrust-basic.md b/docs/openrust-basic.md index 10b8cbd0..412f7bd0 100644 --- a/docs/openrust-basic.md +++ b/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 [`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 diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs new file mode 100644 index 00000000..0d17066e --- /dev/null +++ b/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 = ArrayBase, Dim<[usize; 2]>>; +type OpenVec = ArrayBase, 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, + b_vec: OpenVec, + l: OpenMat, + p: OpenVec, + 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, b: Vec) -> 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 + } +} diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index 66ac7b33..dd0b7038 100644 --- a/src/constraints/mod.rs +++ b/src/constraints/mod.rs @@ -8,6 +8,7 @@ //! //! [`Constraint`]: trait.Constraint.html +mod affine_space; mod ball1; mod ball2; mod ballinf; @@ -23,6 +24,7 @@ mod soc; mod sphere2; mod zero; +pub use affine_space::AffineSpace; pub use ball1::Ball1; pub use ball2::Ball2; pub use ballinf::BallInf; diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs index 2befc82b..8ea4d344 100644 --- a/src/constraints/sphere2.rs +++ b/src/constraints/sphere2.rs @@ -38,7 +38,7 @@ impl<'a> Constraint for Sphere2<'a> { if let Some(center) = &self.center { let norm_difference = crate::matrix_operations::norm2_squared_diff(x, center).sqrt(); if norm_difference <= epsilon { - x.copy_from_slice(¢er); + x.copy_from_slice(center); x[0] += self.radius; return; } diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index d0cf1d0b..1533ec55 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1,6 +1,7 @@ use crate::matrix_operations; use super::*; +use modcholesky::ModCholeskySE99; use rand; #[test] @@ -921,3 +922,69 @@ fn t_epigraph_squared_norm_correctness() { "wrong projection on epigraph of squared norm", ); } +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); +}