From cbb70f8a5ed4593514c1d5260e3783fa2834849c Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 14:39:54 +0000 Subject: [PATCH 01/15] figured out Cholesky with ndarray --- Cargo.toml | 3 +++ src/constraints/tests.rs | 27 +++++++++++++++++++++++++++ 2 files changed, 30 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index e130711f..7906fe43 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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. diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index a09b3c16..312c5d64 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1,4 +1,5 @@ use super::*; +use modcholesky::ModCholeskyGMW81; use rand; #[test] @@ -872,3 +873,29 @@ fn t_sphere2_center_projection_of_center() { fn t_ball1_alpha_negative() { let _ = Ball1::new(None, -1.); } + +#[test] +fn t_affine_space() { + // let m = 4; + // let n = 4; + // let mut a = Array2::::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::::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 = ndarray::arr2(&a_data); + + // Perform modified Cholesky decomposition + // The `Decomposition` struct holds L, E and P + let res = a.mod_cholesky_gmw81(); + + println!("{}", res.l); +} From bdd1f7d50f843ff6dcd0a72927fcfd7c66666a19 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 17:10:46 +0000 Subject: [PATCH 02/15] AffineSpace: computation of Cholesky --- src/constraints/affine_space.rs | 68 +++++++++++++++++++++++++++++++++ src/constraints/mod.rs | 2 + src/constraints/tests.rs | 50 ++++++++++++++---------- 3 files changed, 99 insertions(+), 21 deletions(-) create mode 100644 src/constraints/affine_space.rs diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs new file mode 100644 index 00000000..d884125b --- /dev/null +++ b/src/constraints/affine_space.rs @@ -0,0 +1,68 @@ +use super::Constraint; +use modcholesky::ModCholeskySE99; + +type OpenMat = ndarray::ArrayBase, ndarray::Dim<[usize; 2]>>; +type OpenVec = ndarray::ArrayBase, 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, + b_vec: OpenVec, + a_times_a_t: OpenMat, + l: OpenMat, + p: OpenVec, + e: OpenVec, +} + +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, b_data: Vec) -> 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 + } +} diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index ea8895da..2a02046d 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; @@ -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; diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 312c5d64..f6a911f2 100644 --- a/src/constraints/tests.rs +++ b/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::::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::::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 = ndarray::arr2(&a_data); + let data = [10., 2., 2., 15.].to_vec(); + let datb = vec![1., 2.]; + // let arr: ndarray::ArrayBase, 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::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 = 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); } From e27f87c9073e66d4535933d18758f5fd59af4290 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 18:13:13 +0000 Subject: [PATCH 03/15] AffineSpace: computation of y Solution of Ly = b(P) with backward substitution (tested) --- src/constraints/affine_space.rs | 27 +++++++++++++++++++++++++-- src/constraints/tests.rs | 33 +++++---------------------------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index d884125b..921db9d9 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -11,10 +11,11 @@ type OpenVec = ndarray::ArrayBase, ndarray::Dim<[usize; pub struct AffineSpace { a_mat: OpenMat, b_vec: OpenVec, - a_times_a_t: OpenMat, l: OpenMat, p: OpenVec, e: OpenVec, + n_rows: usize, + n_cols: usize, } impl AffineSpace { @@ -41,25 +42,47 @@ impl AffineSpace { let l = res.l; let p = res.p; let e = res.e; + // Print stuff + println!("A = {:?}", a_mat); + println!("AA' = {:?}", a_times_a_t); + println!("L = {:?}", l); + println!("P = {:?}", p); + println!("E = {:?}", e); // Construct and return new AffineSpace structure AffineSpace { a_mat, b_vec, - a_times_a_t, l, p, e, + n_rows, + n_cols, } } } impl Constraint for AffineSpace { fn project(&self, x: &mut [f64]) { + assert!(x.len() == self.n_cols, "x has wrong dimension"); 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); + // 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.; self.n_rows]; + y[0] = err[self.p[0]] / self.l[(0, 0)]; + for m in 1..self.n_rows { + let mut sum = 0.; + for i in 0..m { + sum += self.l[(m, i)] * y[i]; + } + y[m] = (err[self.p[m]] - sum) / self.l[(m, m)]; + } + println!("y = {:?}", y); + // Step 2: Solve L'x(P) = y } fn is_convex(&self) -> bool { diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index f6a911f2..13fbf1a9 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -876,34 +876,11 @@ fn t_ball1_alpha_negative() { #[test] fn t_affine_space() { - let data = [10., 2., 2., 15.].to_vec(); - let datb = vec![1., 2.]; - // let arr: ndarray::ArrayBase, 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::Dim<[usize; 1]>> = ndarray::Array1::from_shape_vec((nrows,), datb).unwrap(); - + let data = 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 datb = vec![1., 2., -0.5]; let aff = AffineSpace::new(data, datb); - let mut xx = [1., 1.]; + let mut xx = [1., -2., -0.3, 0.5]; aff.project(&mut xx); - - // let a: ndarray::Array2 = 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; - - // println!("{}", error); - - // println!("{}", res.e); } From 0d0209e4acb858d0c2d0eb59a50b05103a2f3d0a Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 18:30:27 +0000 Subject: [PATCH 04/15] AffineSpace: almost there backward substitution: WIP --- src/constraints/affine_space.rs | 14 +++++++++++++- src/constraints/tests.rs | 1 + 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 921db9d9..63be6125 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -82,7 +82,19 @@ impl Constraint for AffineSpace { y[m] = (err[self.p[m]] - sum) / self.l[(m, m)]; } println!("y = {:?}", y); - // Step 2: Solve L'x(P) = y + // Step 2: Solve L'z(P) = y + let mut z = vec![0.; self.n_rows]; + z[self.p[self.n_rows - 1]] = + y[self.n_rows - 1] / self.l[(self.n_rows - 1, self.n_rows - 1)]; + for m in (0..self.n_rows - 1).rev() { + // TODO! (WIP) + } + println!("z = {:?}", z); + // Step 3: Determine A' * z + let z_arr = ndarray::Array1::from_shape_vec((self.n_rows,), z).unwrap(); + let w = self.a_mat.t().dot(&z_arr); + println!("w = {:?}", w); + x.iter_mut().zip(w.iter()).for_each(|(xi, wi)| *xi -= wi); } fn is_convex(&self) -> bool { diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 13fbf1a9..513a72e2 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -883,4 +883,5 @@ fn t_affine_space() { let aff = AffineSpace::new(data, datb); let mut xx = [1., -2., -0.3, 0.5]; aff.project(&mut xx); + println!("x = {:?}", xx); } From 93c424b0fef8b74e86cc2f4f12bbaf7827af9734 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 19:23:59 +0000 Subject: [PATCH 05/15] backward substitution: wip --- src/constraints/affine_space.rs | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 63be6125..b2a003bd 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -63,33 +63,40 @@ impl AffineSpace { impl Constraint for AffineSpace { fn project(&self, x: &mut [f64]) { + let m = self.n_rows; + let n = self.n_cols; + let L = &self.l; + let P = &self.p; + assert!(x.len() == self.n_cols, "x has wrong dimension"); let x_vec = x.to_vec(); - let x_arr = ndarray::Array1::from_shape_vec((x_vec.len(),), x_vec).unwrap(); + let x_arr = ndarray::Array1::from_shape_vec((n,), x_vec).unwrap(); let ax = self.a_mat.dot(&x_arr); let err = ax - &self.b_vec; println!("err = {:?}", err); + // 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.; self.n_rows]; - y[0] = err[self.p[0]] / self.l[(0, 0)]; - for m in 1..self.n_rows { + let mut y = vec![0.; m]; + for i in 0..m { let mut sum = 0.; - for i in 0..m { - sum += self.l[(m, i)] * y[i]; + for j in 0..i { + sum += L[(i, j)] * y[j]; } - y[m] = (err[self.p[m]] - sum) / self.l[(m, m)]; + y[i] = (err[P[i]] - sum) / L[(i, i)]; } println!("y = {:?}", y); + // Step 2: Solve L'z(P) = y - let mut z = vec![0.; self.n_rows]; - z[self.p[self.n_rows - 1]] = - y[self.n_rows - 1] / self.l[(self.n_rows - 1, self.n_rows - 1)]; - for m in (0..self.n_rows - 1).rev() { - // TODO! (WIP) - } + let mut z = vec![0.; m]; + z[P[m - 1]] = y[m - 1] / L[(m - 1, m - 1)]; + z[P[m - 2]] = (y[m - 2] - L[(m - 1, m - 2)] * z[P[m - 1]]) / L[(m - 2, m - 2)]; + z[P[m - 3]] = + (y[m - 3] - L[(m - 2, m - 3)] * z[P[m - 2]] - L[(m - 1, m - 3)] * z[P[m - 1]]) + / L[(m - 3, m - 3)]; println!("z = {:?}", z); + // Step 3: Determine A' * z let z_arr = ndarray::Array1::from_shape_vec((self.n_rows,), z).unwrap(); let w = self.a_mat.t().dot(&z_arr); From e5bbf77f148c65f404d119ed365b625d650551c6 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 22:33:56 +0000 Subject: [PATCH 06/15] AffineSpace: first implementation Not fully tested --- src/constraints/affine_space.rs | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index b2a003bd..1050dfbb 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -65,15 +65,14 @@ impl Constraint for AffineSpace { fn project(&self, x: &mut [f64]) { let m = self.n_rows; let n = self.n_cols; - let L = &self.l; - let P = &self.p; + let chol: &ndarray::ArrayBase, ndarray::Dim<[usize; 2]>> = &self.l; + let perm = &self.p; assert!(x.len() == self.n_cols, "x has wrong dimension"); let x_vec = x.to_vec(); let x_arr = ndarray::Array1::from_shape_vec((n,), x_vec).unwrap(); let ax = self.a_mat.dot(&x_arr); let err = ax - &self.b_vec; - println!("err = {:?}", err); // Step 1: Solve Ly = b(P) // TODO: Make `y` into an attribute; however, to do this, we need to change @@ -82,25 +81,26 @@ impl Constraint for AffineSpace { for i in 0..m { let mut sum = 0.; for j in 0..i { - sum += L[(i, j)] * y[j]; + sum += chol[(i, j)] * y[j]; } - y[i] = (err[P[i]] - sum) / L[(i, i)]; + y[i] = (err[perm[i]] - sum) / chol[(i, i)]; } - println!("y = {:?}", y); // Step 2: Solve L'z(P) = y let mut z = vec![0.; m]; - z[P[m - 1]] = y[m - 1] / L[(m - 1, m - 1)]; - z[P[m - 2]] = (y[m - 2] - L[(m - 1, m - 2)] * z[P[m - 1]]) / L[(m - 2, m - 2)]; - z[P[m - 3]] = - (y[m - 3] - L[(m - 2, m - 3)] * z[P[m - 2]] - L[(m - 1, m - 3)] * z[P[m - 1]]) - / L[(m - 3, m - 3)]; - println!("z = {:?}", z); + 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 = ndarray::Array1::from_shape_vec((self.n_rows,), z).unwrap(); let w = self.a_mat.t().dot(&z_arr); - println!("w = {:?}", w); + + // Step 4: x <-- x - A'(AA')\(Ax-b) x.iter_mut().zip(w.iter()).for_each(|(xi, wi)| *xi -= wi); } From 6fb5b8bb7af827c680d6ecf1539d74d4cb8af28a Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 23:12:46 +0000 Subject: [PATCH 07/15] AffineSpace: more testing --- src/constraints/affine_space.rs | 57 +++++++++++++++++++++------------ src/constraints/tests.rs | 32 ++++++++++++++---- 2 files changed, 62 insertions(+), 27 deletions(-) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 1050dfbb..96f16729 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -1,19 +1,19 @@ use super::Constraint; use modcholesky::ModCholeskySE99; +use ndarray::{Array1, Array2, ArrayBase, Dim, OwnedRepr}; -type OpenMat = ndarray::ArrayBase, ndarray::Dim<[usize; 2]>>; -type OpenVec = ndarray::ArrayBase, ndarray::Dim<[usize; 1]>>; +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, $\{x\in\mathbb{R}^n: Ax = b\}$, which is an affine space. It is assumed that +/// 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, - e: OpenVec, n_rows: usize, n_cols: usize, } @@ -21,18 +21,27 @@ pub struct AffineSpace { 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, b_data: Vec) -> Self { + /// + /// ## 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_data.len(); - let n_elements_a = a_data.len(); + 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 = 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(); + 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' @@ -41,20 +50,13 @@ impl AffineSpace { let res = a_times_a_t.mod_cholesky_se99(); let l = res.l; let p = res.p; - let e = res.e; - // Print stuff - println!("A = {:?}", a_mat); - println!("AA' = {:?}", a_times_a_t); - println!("L = {:?}", l); - println!("P = {:?}", p); - println!("E = {:?}", e); + // Construct and return new AffineSpace structure AffineSpace { a_mat, b_vec, l, p, - e, n_rows, n_cols, } @@ -62,15 +64,27 @@ impl AffineSpace { } 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 + /// fn project(&self, x: &mut [f64]) { let m = self.n_rows; let n = self.n_cols; - let chol: &ndarray::ArrayBase, ndarray::Dim<[usize; 2]>> = &self.l; + let chol = &self.l; let perm = &self.p; - assert!(x.len() == self.n_cols, "x has wrong dimension"); + assert!(x.len() == n, "x has wrong dimension"); let x_vec = x.to_vec(); - let x_arr = ndarray::Array1::from_shape_vec((n,), x_vec).unwrap(); + 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; @@ -97,13 +111,14 @@ impl Constraint for AffineSpace { } // Step 3: Determine A' * z - let z_arr = ndarray::Array1::from_shape_vec((self.n_rows,), z).unwrap(); + 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/tests.rs b/src/constraints/tests.rs index 513a72e2..cada87fa 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -876,12 +876,32 @@ fn t_ball1_alpha_negative() { #[test] fn t_affine_space() { - let data = vec![ + 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 datb = vec![1., 2., -0.5]; - let aff = AffineSpace::new(data, datb); - let mut xx = [1., -2., -0.3, 0.5]; - aff.project(&mut xx); - println!("x = {:?}", xx); + 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] +#[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); } From 16e6d6b738ab37b033b51134a971a5d491b66ea8 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 23:15:57 +0000 Subject: [PATCH 08/15] AffineSpace: extern crate added --- CHANGELOG.md | 4 ++++ src/constraints/affine_space.rs | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 54b22078..0ea7803d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,10 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo --------------------- --> ## Unreleased +### Added + +- Implementation of `AffineSpace` + ### Fixed - Clippy fixes diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 96f16729..fc650d30 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -1,4 +1,8 @@ use super::Constraint; + +extern crate modcholesky; +extern crate ndarray; + use modcholesky::ModCholeskySE99; use ndarray::{Array1, Array2, ArrayBase, Dim, OwnedRepr}; From aa891e627f0ad3e6b807c68d41d4b74c36b86699 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Sun, 29 Oct 2023 23:29:05 +0000 Subject: [PATCH 09/15] update appveyor cargo add ndarray and modcholesky --- appveyor.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/appveyor.yml b/appveyor.yml index 614b3b7e..182dfc5f 100644 --- a/appveyor.yml +++ b/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% From 2234ab22b45a06c4a94f0225537d961038360df7 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Mon, 30 Oct 2023 00:30:43 +0000 Subject: [PATCH 10/15] AffineSpace: more testing --- src/constraints/affine_space.rs | 6 +++--- src/constraints/tests.rs | 25 +++++++++++++++++++++++++ 2 files changed, 28 insertions(+), 3 deletions(-) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index fc650d30..85e13da5 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -97,11 +97,11 @@ impl Constraint for AffineSpace { // &self to &mut self, which will require a mild refactoring let mut y = vec![0.; m]; for i in 0..m { - let mut sum = 0.; + y[i] = err[perm[i]]; for j in 0..i { - sum += chol[(i, j)] * y[j]; + y[i] -= chol[(i, j)] * y[j]; } - y[i] = (err[perm[i]] - sum) / chol[(i, i)]; + y[i] /= chol[(i, i)]; } // Step 2: Solve L'z(P) = y diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index cada87fa..72fa083d 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -898,6 +898,31 @@ fn t_affine_space() { ); } +#[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] #[should_panic] fn t_affine_space_wrong_dimensions() { From 44c41204fbb12dd4505db350cf560c30e7074cf8 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Mon, 30 Oct 2023 00:37:03 +0000 Subject: [PATCH 11/15] AffineSpace: another test A has single row --- src/constraints/tests.rs | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 72fa083d..eaab661c 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -923,6 +923,17 @@ fn t_affine_space_larger() { ); } +#[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() { From 580df713ff17d3c4689ac7df83aa4d5864c5bd3e Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Mon, 30 Oct 2023 11:56:14 +0000 Subject: [PATCH 12/15] Update CHANGELOG New version: 0.9.0 (unreleased) --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0ea7803d..08f21fb0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,7 +10,7 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo -## Unreleased +## [v0.9.0] - Unreleased ### Added @@ -279,6 +279,7 @@ This is a breaking API change. --------------------- --> +[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 From adc8f93236be1de940255bb80268ff66d8a8fe08 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Tue, 31 Oct 2023 01:29:08 +0000 Subject: [PATCH 13/15] affine space: example in docs --- src/constraints/affine_space.rs | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 85e13da5..0d17066e 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -80,6 +80,25 @@ impl Constraint for AffineSpace { /// /// - `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; From 9ce37f154384cb7ff09c276363b7a21fd8de4059 Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Tue, 31 Oct 2023 01:39:34 +0000 Subject: [PATCH 14/15] update website docs (rust constraints) --- docs/openrust-basic.md | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) 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 From 9b297e45be4007d66989cdb1c8e1e663d742750e Mon Sep 17 00:00:00 2001 From: Pantelis Sopasakis Date: Tue, 31 Oct 2023 01:40:46 +0000 Subject: [PATCH 15/15] [ci ckip] bump version in Cargo.toml --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 7906fe43..4e4ea073 100644 --- a/Cargo.toml +++ b/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"