diff --git a/CHANGELOG.md b/CHANGELOG.md index a71ece5a..83c1cb8f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,13 +9,19 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo -## [v0.9.2] - Unreleased +## [v0.10.0] - 2026-03-10 + +### Added + +- Custom implementation of Cholesky factorisation (and solver); this is used in `AffineSpace` now. +- New function in `matrix_operations` to compute AA' given a matrix A ### Changed -- Update version of `rand`, `ndarray`, and `modcholesky` in `Cargo.toml` +- Update version of `ndarray`, in `Cargo.toml` +- Removed `modcholesky` because it was causing a bug (see issue #378) -[v0.9.2]: https://github.com/alphaville/optimization-engine/compare/v0.9.1...v0.9.2 +[v0.10.0]: https://github.com/alphaville/optimization-engine/compare/v0.9.1...v0.10.0 [v0.9.1]: https://github.com/alphaville/optimization-engine/compare/v0.9.0...v0.9.1 [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 diff --git a/Cargo.toml b/Cargo.toml index b967529b..70da22a4 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.9.1" +version = "0.10.0" edition = "2018" @@ -80,7 +80,7 @@ rustdoc-args = ["--html-in-header", "katex-header.html"] num = "0.4" # Our own stuff - L-BFGS: limited-memory BFGS directions -lbfgs = "0.2" +lbfgs = "0.3" # Instant is a generic timer that works on Wasm (with wasm-bindgen) instant = { version = "0.1" } @@ -98,10 +98,7 @@ rpmalloc = { version = "0.2", features = [ # epigraph of the squared Euclidean norm roots = "0.0.8" -# Least squares solver (NOTE: ndarray must be version 0.15 - not 0.16) -# Bug report: https://github.com/argmin-rs/modcholesky/issues/34 -ndarray = { version = "=0.16.1", features = ["approx"] } -modcholesky = "0.2" +ndarray = { version = "=0.17.2", features = ["approx"] } # jemallocator is an optional feature; it will only be loaded if the feature # `jem` is used (i.e., if we compile with `cargo build --features jem`) @@ -127,15 +124,15 @@ wasm = ["wasm-bindgen", "instant/wasm-bindgen", "instant/inaccurate"] # -------------------------------------------------------------------------- # These dependencies are only used for testing [dev-dependencies] -unit_test_utils = "0.1.3" +unit_test_utils = "0.1.4" # Note: path dependencies (even dev-dependencies) cannot be local in order # to publish a crate, so the following line can only be used when # testing locally # icasadi_test = { path = "test/icasadi_test/" } # instead, use: -icasadi_test = "0.0.2" +icasadi_test = "0.0.3" # Random number generators for unit tests: -rand = "0.9" +rand = "0.10" # -------------------------------------------------------------------------- diff --git a/appveyor.yml b/appveyor.yml index e216a306..e91d5115 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -70,6 +70,5 @@ build: false test_script: - cargo add roots - cargo add ndarray --features approx - - cargo add modcholesky - cargo build - cargo test --verbose %cargoflags% diff --git a/open-codegen/CHANGELOG.md b/open-codegen/CHANGELOG.md index b9e7ccde..a5be9058 100644 --- a/open-codegen/CHANGELOG.md +++ b/open-codegen/CHANGELOG.md @@ -11,6 +11,7 @@ Note: This is the Changelog file of `opengen` - the Python interface of OpEn ### Changed +- Removed dependence on `pkg_resources` (issue #380) - Additional unit tests: increased coverage to 92% ## [0.9.4] - 2025-05-08 diff --git a/open-codegen/VERSION b/open-codegen/VERSION index 2bd77c74..03834411 100644 --- a/open-codegen/VERSION +++ b/open-codegen/VERSION @@ -1 +1 @@ -0.9.4 \ No newline at end of file +0.9.5 \ No newline at end of file diff --git a/open-codegen/opengen/builder/optimizer_builder.py b/open-codegen/opengen/builder/optimizer_builder.py index 67956171..0ff153e7 100644 --- a/open-codegen/opengen/builder/optimizer_builder.py +++ b/open-codegen/opengen/builder/optimizer_builder.py @@ -9,9 +9,9 @@ import os import jinja2 import logging -import pkg_resources import sys +from importlib.metadata import version from .ros_builder import RosBuilder _AUTOGEN_COST_FNAME = 'auto_casadi_cost.c' @@ -750,7 +750,7 @@ def __generate_yaml_data_file(self): target_yaml_file_path = os.path.join(target_dir, "optimizer.yml") - opengen_version = pkg_resources.require("opengen")[0].version + opengen_version = version("opengen") tcp_details = None if tcp_config is None \ else {'ip': tcp_config.bind_ip, 'port': tcp_config.bind_port} diff --git a/open-codegen/opengen/definitions.py b/open-codegen/opengen/definitions.py index 40baed33..7949027b 100644 --- a/open-codegen/opengen/definitions.py +++ b/open-codegen/opengen/definitions.py @@ -1,21 +1,22 @@ -import pkg_resources +from importlib.resources import files, as_file def templates_dir(): """Directory where the templates are found (for internal use, mainly)""" - return pkg_resources.resource_filename('opengen', 'templates/') + return files("opengen") / "templates" def templates_subdir(subdir=None): """ Directory where the templates are found and subfolder relative - to that path(for internal use, mainly) + to that path (for internal use, mainly) """ - if subdir is None: - return templates_dir() - return pkg_resources.resource_filename('opengen', 'templates/%s/' % subdir) + resource = files("opengen") / "templates" + if subdir is not None: + resource = resource / subdir + return resource def original_icasadi_dir(): """Directory where the original icasadi files are found (for internal use)""" - return pkg_resources.resource_filename('opengen', 'icasadi/') + return files("opengen") / "icasadi" diff --git a/open-codegen/opengen/tcp/optimizer_tcp_manager.py b/open-codegen/opengen/tcp/optimizer_tcp_manager.py index 8345c6aa..68df889f 100644 --- a/open-codegen/opengen/tcp/optimizer_tcp_manager.py +++ b/open-codegen/opengen/tcp/optimizer_tcp_manager.py @@ -6,11 +6,10 @@ import logging import time import math -import pkg_resources from threading import Thread from retry import retry from .solver_response import SolverResponse - +from importlib.metadata import version class OptimizerTcpManager: """Client for TCP interface of parametric optimizers @@ -66,8 +65,7 @@ def __init__(self, optimizer_path=None, ip=None, port=None): # Check whether the optimizer was built with the current version of opengen # We can only check the optimizer version if the optimizer runs locally opengen_version = self.__optimizer_details['build']['opengen_version'] - current_opengen_version = pkg_resources.require("opengen")[ - 0].version + current_opengen_version = version("opengen") if current_opengen_version != opengen_version: logging.warn( 'the target optimizer was build with a different version of opengen (%s)' % opengen_version) diff --git a/open-codegen/test/test_raspberry_pi.py b/open-codegen/test/test_raspberry_pi.py index e6236241..f040f95f 100644 --- a/open-codegen/test/test_raspberry_pi.py +++ b/open-codegen/test/test_raspberry_pi.py @@ -1,3 +1,4 @@ +import os import opengen as og import unittest import casadi.casadi as cs @@ -7,6 +8,11 @@ class RaspberryPiTest(unittest.TestCase): + @staticmethod + def get_open_local_absolute_path(): + cwd = os.getcwd() + return cwd.split('open-codegen')[0] + # ----------------------------------------------------------------------- # Cross-compile to Raspberry Pi # ----------------------------------------------------------------------- @@ -25,6 +31,7 @@ def test_compile_rpi(self): meta = og.config.OptimizerMeta()\ .with_optimizer_name(optimizer_name) build_config = og.config.BuildConfiguration() \ + .with_open_version(local_path=RaspberryPiTest.get_open_local_absolute_path()) \ .with_build_directory(optimizers_dir) \ .with_build_mode(og.config.BuildConfiguration.DEBUG_MODE) \ .with_tcp_interface_config() \ diff --git a/src/alm/alm_cache.rs b/src/alm/alm_cache.rs index e014c3bc..2f5211bf 100644 --- a/src/alm/alm_cache.rs +++ b/src/alm/alm_cache.rs @@ -51,7 +51,7 @@ impl AlmCache { /// # Arguments /// /// - `panoc_cache`: an instance of `PANOCCache` that will be used by - /// the inner problem + /// the inner problem /// - `n1`, `n2`: range dimensions of mappings `F1` and `F2` respectively /// /// # Panics @@ -77,9 +77,9 @@ impl AlmCache { w_pm: if n2 > 0 { Some(vec![0.0; n2]) } else { None }, iteration: 0, delta_y_norm: 0.0, - delta_y_norm_plus: std::f64::INFINITY, + delta_y_norm_plus: f64::INFINITY, f2_norm: 0.0, - f2_norm_plus: std::f64::INFINITY, + f2_norm_plus: f64::INFINITY, inner_iteration_count: 0, last_inner_problem_norm_fpr: -1.0, available_time: None, diff --git a/src/alm/alm_factory.rs b/src/alm/alm_factory.rs index 735c2be5..21776453 100644 --- a/src/alm/alm_factory.rs +++ b/src/alm/alm_factory.rs @@ -13,7 +13,7 @@ use crate::{constraints::Constraint, matrix_operations, FunctionCallResult}; /// # Types /// /// - `Cost`: cost function $f:\mathbb{R}^{n_u} \to \mathbb{R}$ which is computed -/// by a function with signature: +/// by a function with signature: /// ///```rust,ignore ///fn f(u: &[f64], cost: &mut f64) -> FunctionCallResult @@ -22,7 +22,7 @@ use crate::{constraints::Constraint, matrix_operations, FunctionCallResult}; /// where `cost` is updated with the value $f(u)$, /// /// - `CostGradient`: gradient of the cost function, $\nabla f: \mathbb{R}^{n_u} \to \mathbb{R}^{n_u}$, -/// which is computed by a function with signature +/// which is computed by a function with signature /// /// ```rust,ignore /// fn df(u: &[f64], grad: &mut [f64]) -> FunctionCallResult @@ -31,19 +31,19 @@ use crate::{constraints::Constraint, matrix_operations, FunctionCallResult}; /// where on exit `grad` stores the /// /// - `MappingF1` and `MappingF2`: mappings $F_1:\mathbb{R}^n\to\mathbb{R}^{n_1}$ -/// and $F_2:\mathbb{R}^n\to\mathbb{R}^{n_2}$ which -/// are computed by functions with signature +/// and $F_2:\mathbb{R}^n\to\mathbb{R}^{n_2}$ which are computed by functions +/// with signature /// /// ```rust,ignore /// fn mapping(u: &[f64], fu: &mut [f64]) -> FunctionCallResult /// ``` /// /// - `JacobianMappingF1Trans` and `JacobianMappingF2Trans`: functions that compute -/// product of the form $JF_i(u)^\top{}d$ for given $d\in\mathbb{R}^{n_i}$ and -/// $u\in\mathbb{R}^{n_u}$ +/// product of the form $JF_i(u)^\top{}d$ for given $d\in\mathbb{R}^{n_i}$ and +/// $u\in\mathbb{R}^{n_u}$ /// /// - `SetC`: A set $C\subseteq \mathbb{R}^{n_1}$, which is used in the definition -/// of the constraints $F_1(u) \in C$ +/// of the constraints $F_1(u) \in C$ /// /// The above are used to compute $\psi:\mathbb{R}^{n_u}\to\mathbb{R}$ for given /// $u\in\mathbb{R}^{n_u}$ and $\xi=(c, y)\in\mathbb{R}^{n_1+1}$, where $c\in\mathbb{R}$ diff --git a/src/alm/alm_optimizer.rs b/src/alm/alm_optimizer.rs index 0fcb1fde..c7d9ec08 100644 --- a/src/alm/alm_optimizer.rs +++ b/src/alm/alm_optimizer.rs @@ -13,7 +13,7 @@ const DEFAULT_PENALTY_UPDATE_FACTOR: f64 = 5.0; const DEFAULT_EPSILON_UPDATE_FACTOR: f64 = 0.1; const DEFAULT_INFEAS_SUFFICIENT_DECREASE_FACTOR: f64 = 0.1; const DEFAULT_INITIAL_TOLERANCE: f64 = 0.1; -const SMALL_EPSILON: f64 = std::f64::EPSILON; +const SMALL_EPSILON: f64 = f64::EPSILON; /// Internal/private structure used by method AlmOptimizer.step /// to return some minimal information about the inner problem @@ -81,8 +81,8 @@ impl InnerProblemStatus { /// - For $\nu=0,\ldots, \nu_{\max}$ /// - $y \gets \Pi_Y(y)$ /// - $u \gets \arg\min_{u\in U} \psi(u, \xi)$, where $\psi(u, \xi)$ is a given function: this problem is -/// solved with tolerance $\bar\epsilon$ -/// (see [`AlmFactory`](./struct.AlmFactory.html) regarding how this is constructed) +/// solved with tolerance $\bar\epsilon$ +/// (see [`AlmFactory`](./struct.AlmFactory.html) regarding how this is constructed) /// - $y^+ \gets y + c(F_1(u) - \Pi_C(F_1(u) + y/c))$ /// - Define $z^+ \gets \Vert y^+ - y \Vert$ and $t^+ = \Vert F_2(u) \Vert$ /// - If $z^+ \leq c\delta$, $t^+ \leq \delta$ and $\epsilon_\nu \leq \epsilon$, return $(u, y^+)$ @@ -205,10 +205,10 @@ where /// # Arguments /// /// - `alm_cache`: a reuseable instance of [`AlmCache`](./struct.AlmCache.html), which is borrowed by - /// `AlmOptimizer` + /// `AlmOptimizer` /// - `alm_problem`: the problem specification (data for $\psi(u, \xi)$, - /// $\nabla_u \psi(u, \xi)$, $F_1(u)$ (if any), $F_2(u)$ (if any), and sets - /// $C$, $U$ and $Y$) + /// $\nabla_u \psi(u, \xi)$, $F_1(u)$ (if any), $F_2(u)$ (if any), and sets + /// $C$, $U$ and $Y$) /// /// /// # Example @@ -467,7 +467,7 @@ where /// # Arguments /// /// - `initial_inner_tolerance`: the initial value of the inner tolerance, that is, - /// the value $\espilon_0$ + /// the value $\espilon_0$ /// /// # Returns /// @@ -530,7 +530,7 @@ where /// # Arguments /// /// - `y_init`: initial vector of Lagrange multipliers (type: `&[f64]`) of - /// length equal to `n1` + /// length equal to `n1` /// /// # Returns /// @@ -732,7 +732,7 @@ where .with_max_duration( alm_cache .available_time - .unwrap_or_else(|| std::time::Duration::from_secs(std::u64::MAX)), + .unwrap_or_else(|| std::time::Duration::from_secs(u64::MAX)), ) // Set the maximum number of inner iterations .with_max_iter(self.max_inner_iterations); @@ -1059,7 +1059,7 @@ mod tests { // Test: the initial value of the penalty parameter is positive if let Some(xi) = &alm_optimizer.alm_cache.xi { - assert!(xi[0] > std::f64::EPSILON); + assert!(xi[0] > f64::EPSILON); } // Test: with_initial_penalty diff --git a/src/alm/alm_optimizer_status.rs b/src/alm/alm_optimizer_status.rs index e6c9bbed..5c10e477 100644 --- a/src/alm/alm_optimizer_status.rs +++ b/src/alm/alm_optimizer_status.rs @@ -19,7 +19,7 @@ pub struct AlmOptimizerStatus { num_inner_iterations: usize, /// Norm of the fixed-point residual of the the problem last_problem_norm_fpr: f64, - /// + /// Lagrange multipliers vector lagrange_multipliers: Option>, /// Total solve time solve_time: std::time::Duration, @@ -122,7 +122,7 @@ impl AlmOptimizerStatus { /// # Arguments /// /// - `lagrange_multipliers`: vector of Lagrange multipliers (which is copied - /// into an internal field of `AlmOptimizerStatus`) + /// into an internal field of `AlmOptimizerStatus`) /// /// # Panics /// diff --git a/src/alm/alm_problem.rs b/src/alm/alm_problem.rs index 7827069e..9f28bd2e 100644 --- a/src/alm/alm_problem.rs +++ b/src/alm/alm_problem.rs @@ -103,7 +103,7 @@ where /// - `constraints`: hard constraints, set $U$ /// - `alm_set_c`: Set $C$ of ALM-specific constraints (convex, closed) /// - `alm_set_y`: Compact, convex set $Y$ of Lagrange multipliers, which needs to be a - /// compact subset of $C^*$ (the convex conjugate of the convex set $C{}\subseteq{}\mathbb{R}^{n_1}$) + /// compact subset of $C^*$ (the convex conjugate of the convex set $C{}\subseteq{}\mathbb{R}^{n_1}$) /// - `parametric_cost`: Parametric cost function, $\psi(u, \xi)$, where $\xi = (c, y)$ /// - `parametric_gradient`: Gradient of cost function wrt $u$, that is $\nabla_x \psi(u, \xi)$ /// - `mapping_f1`: Mapping `F1` of ALM-specific constraints ($F1(u) \in C$) diff --git a/src/cholesky_factorizer.rs b/src/cholesky_factorizer.rs new file mode 100644 index 00000000..12cad81c --- /dev/null +++ b/src/cholesky_factorizer.rs @@ -0,0 +1,302 @@ +//! Cholesky factorization and linear solves for symmetric positive-definite matrices. +//! +//! This module provides a [`CholeskyFactorizer`] type for computing and storing the +//! Cholesky factorization of a square matrix +//! +//! $$A = L L^\intercal$$ +//! +//! where $L$ is lower triangular. +//! +//! The matrix is provided as a flat slice in **row-major** order. Internally, the +//! computed Cholesky factor `L` is also stored in row-major order in a dense +//! `Vec`. Only the lower-triangular part is meaningful after factorization. +//! +//! Once a factorization has been computed with [`CholeskyFactorizer::factorize`], +//! the stored factor can be reused to solve linear systems of the form +//! +//! $$A x = b$$ +//! +//! via [`CholeskyFactorizer::solve`]. +//! +//! # Requirements +//! +//! The input matrix must: +//! +//! - be square of dimension `n × n` +//! - be stored in row-major order +//! - be symmetric positive definite +//! +//! If the matrix is not positive definite, factorization fails with +//! [`CholeskyError::NotPositiveDefinite`]. +//! +//! # Errors +//! +//! The module uses [`CholeskyError`] to report the following conditions: +//! +//! - [`CholeskyError::DimensionMismatch`] when the input matrix or right-hand side +//! has an incompatible size +//! - [`CholeskyError::NotPositiveDefinite`] when factorization fails because the +//! matrix is not positive definite +//! - [`CholeskyError::NotFactorized`] when attempting to solve a system before a +//! valid factorization has been computed +//! +//! # Example +//! +//! ```rust +//! use optimization_engine::CholeskyFactorizer; +//! +//! let a = vec![ +//! 4.0_f64, 1.0, +//! 1.0, 3.0, +//! ]; +//! +//! let b = vec![1.0_f64, 2.0]; +//! +//! let mut factorizer = CholeskyFactorizer::new(2); +//! factorizer.factorize(&a).unwrap(); +//! +//! let x = factorizer.solve(&b).unwrap(); +//! +//! assert!((x[0] - 0.0909090909).abs() < 1e-10); +//! assert!((x[1] - 0.6363636364).abs() < 1e-10); +//! ``` +//! +//! # Notes +//! +//! - The implementation is generic over `T: Float`. +//! - Storage is preallocated when constructing the factorizer with +//! [`CholeskyFactorizer::new`]. +//! - The factorizer keeps the computed factor internally so that multiple right-hand +//! sides can be solved efficiently after a single factorization. + +use num::Float; + +#[derive(Debug, Clone)] +/// Cholesky factoriser +pub struct CholeskyFactorizer { + n: usize, + cholesky_factor: Vec, // row-major storage of the lower-triangular factor L + is_factorized: bool, // whether factorization has been completed +} + +#[derive(Debug, Clone, PartialEq, Eq)] +/// Cholesky errors +pub enum CholeskyError { + /// Not positive definite matrix + NotPositiveDefinite, + /// RHS dimension mismatch + DimensionMismatch, + /// Attempting to solve without having factorized + NotFactorized, +} + +impl CholeskyFactorizer { + /// Create a factorizer with storage preallocated for an n x n matrix. + pub fn new(n: usize) -> Self { + Self { + n, + cholesky_factor: vec![T::zero(); n * n], + is_factorized: false, + } + } + + /// Compute the Cholesky factorization $A = L L^\intercal$ + /// from a square matrix in row-major order. + /// + /// The input matrix must have dimension equal to the one used in `new(n)`. + pub fn factorize(&mut self, a: &[T]) -> Result<(), CholeskyError> { + self.is_factorized = false; + if a.len() != self.n * self.n { + return Err(CholeskyError::DimensionMismatch); + } + self.cholesky_factor.fill(T::zero()); + let n = self.n; + + for i in 0..n { + let row_i = i * n; + + for j in 0..=i { + let row_j = j * n; + let mut sum = a[row_i + j]; + + for k in 0..j { + sum = sum - self.cholesky_factor[row_i + k] * self.cholesky_factor[row_j + k]; + } + + if i == j { + if sum <= T::zero() { + return Err(CholeskyError::NotPositiveDefinite); + } + self.cholesky_factor[row_i + i] = sum.sqrt(); + } else { + self.cholesky_factor[row_i + j] = sum / self.cholesky_factor[row_j + j]; + } + } + } + self.is_factorized = true; + Ok(()) + } + + #[inline] + #[must_use] + /// dimension + pub fn dimension(&self) -> usize { + self.n + } + + #[inline] + #[must_use] + /// Cholesky factor, L + pub fn cholesky_factor(&self) -> &[T] { + &self.cholesky_factor + } + + /// Solves the linear system $A x = b$ using the stored Cholesky factorization + /// $A = L L^\intercal$. + /// + /// This method assumes that [`Self::factorize`] has already been called successfully, + /// so that the lower-triangular Cholesky factor `L` is available internally. + /// + /// The solution is computed in two steps: + /// + /// 1. Forward substitution to solve $L y = b$ + /// 2. Back substitution to solve $L^\intercal x = y$ + /// + /// # Arguments + /// + /// * `b` - Right-hand-side vector of length `n` + /// + /// # Returns + /// + /// Returns the solution vector `x` such that $A x = b$. + /// + /// # Errors + /// + /// Returns an error in the following cases: + /// + /// * `CholeskyError::NotFactorized` if no valid Cholesky factorization is currently + /// stored in the factorizer + /// * `CholeskyError::DimensionMismatch` if `b.len() != n`, where `n` is the matrix dimension + /// + /// # Notes + /// + /// * The matrix `A` itself is not accessed directly; only its Cholesky factor `L` is used. + /// * The Cholesky factor is stored internally in row-major order. + /// * This method allocates temporary vectors for the intermediate solution `y` and the + /// final solution `x`. + pub fn solve(&self, b: &[T]) -> Result, CholeskyError> { + if !self.is_factorized { + return Err(CholeskyError::NotFactorized); + } + if b.len() != self.n { + return Err(CholeskyError::DimensionMismatch); + } + + let n = self.n; + + // Forward substitution: L y = b + let mut y = vec![T::zero(); n]; + for i in 0..n { + let row_i = i * n; + let mut sum = b[i]; + for (&lij, &yj) in self.cholesky_factor[row_i..row_i + i] + .iter() + .zip(y[..i].iter()) + { + sum = sum - lij * yj; + } + y[i] = sum / self.cholesky_factor[row_i + i]; + } + + // Back substitution: L^T x = y + let mut x = vec![T::zero(); n]; + for i in (0..n).rev() { + let mut sum = y[i]; + for (row_j, &xj) in self + .cholesky_factor + .chunks_exact(n) + .skip(i + 1) + .zip(x.iter().skip(i + 1)) + { + sum = sum - row_j[i] * xj; + } + x[i] = sum / self.cholesky_factor[i * n + i]; + } + + Ok(x) + } +} + +/* ---------------------------------------------------------------------------- */ +/* TESTS */ +/* ---------------------------------------------------------------------------- */ +#[cfg(test)] +mod tests { + use crate::*; + + #[test] + fn t_cholesky_basic() { + let a = vec![4.0_f64, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 98.0]; + let mut factorizer = CholeskyFactorizer::new(3); + let _ = factorizer.factorize(&a); + assert!(3 == factorizer.dimension(), "wrong dimension"); + let expected_l = [2.0, 0.0, 0.0, 6.0, 1.0, 0.0, -8.0, 5.0, 3.0]; + unit_test_utils::nearly_equal_array( + &expected_l, + &factorizer.cholesky_factor(), + 1e-10, + 1e-12, + ); + } + + #[test] + fn t_cholesky_solve_linear_system() { + let a = vec![4.0_f64, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 98.0]; + let mut factorizer = CholeskyFactorizer::new(3); + let _ = factorizer.factorize(&a); + let rhs = vec![-5.0_f64, 2.0, -3.0]; + let x = factorizer.solve(&rhs).unwrap(); + let expected_sol = [-280.25_f64, 77., -12.]; + unit_test_utils::nearly_equal_array(&expected_sol, &x, 1e-10, 1e-12); + } + + #[test] + fn t_cholesky_not_square_matrix() { + let a = vec![1.0_f64, 2., 7., 5., 9.]; + let mut factorizer = CholeskyFactorizer::new(3); + let result = factorizer.factorize(&a); + assert_eq!(result, Err(CholeskyError::DimensionMismatch)); + } + + #[test] + fn t_cholesky_solve_wrong_dimension_rhs() { + let a = vec![4.0_f64, 12.0, -16.0, 12.0, 37.0, -43.0, -16.0, -43.0, 98.0]; + let mut factorizer = CholeskyFactorizer::new(3); + let _ = factorizer.factorize(&a); + let rhs = vec![-5.0_f64, 2.0]; + let result = factorizer.solve(&rhs); + assert_eq!(result, Err(CholeskyError::DimensionMismatch)); + } + + #[test] + fn t_cholesky_solve_not_factorized_1() { + let factorizer = CholeskyFactorizer::new(3); + let rhs = vec![-5.0_f64, 2.0]; + let result = factorizer.solve(&rhs); + assert_eq!(result, Err(CholeskyError::NotFactorized)); + } + + #[test] + fn t_cholesky_solve_not_factorized_2() { + let a = vec![1.0_f64, 1.0, 1.0, 1.0]; + let mut factorizer = CholeskyFactorizer::new(2); + let factorization_result = factorizer.factorize(&a); + assert_eq!( + factorization_result, + Err(CholeskyError::NotPositiveDefinite) + ); + let rhs = vec![-5.0_f64, 2.0]; + let result = factorizer.solve(&rhs); + assert_eq!(result, Err(CholeskyError::NotFactorized)); + } +} diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 0d17066e..2f25ceba 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -1,23 +1,17 @@ use super::Constraint; +use crate::matrix_operations; +use crate::CholeskyFactorizer; -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]>>; +use ndarray::{ArrayView1, ArrayView2}; #[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, + a_mat: Vec, + b_vec: Vec, + factorizer: CholeskyFactorizer, n_rows: usize, n_cols: usize, } @@ -35,32 +29,21 @@ impl AffineSpace { /// 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_rows > 0, "b must not be empty"); assert!( - n_elements_a % n_rows == 0, + n_elements_a.is_multiple_of(n_rows), "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 + let aat = matrix_operations::mul_a_at(&a, n_rows, n_cols).unwrap(); + let mut factorizer = CholeskyFactorizer::new(n_rows); + factorizer.factorize(&aat).unwrap(); AffineSpace { - a_mat, - b_vec, - l, - p, + a_mat: a, + b_vec: b, + factorizer, n_rows, n_cols, } @@ -73,8 +56,8 @@ impl Constraint for AffineSpace { /// 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). + /// is solved by computing the Cholesky factorization of $AA^\intercal$, which is + /// done using `CholeskyFactorizer` /// /// ## Arguments /// @@ -100,45 +83,25 @@ impl Constraint for AffineSpace { /// /// 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 1: Compute e = Ax - b + let a = ArrayView2::from_shape((self.n_rows, self.n_cols), &self.a_mat) + .expect("invalid A shape"); + let x_view = ArrayView1::from(&x[..]); + let b = ArrayView1::from(&self.b_vec[..]); + let e = a.dot(&x_view) - b; + let e_slice: &[f64] = e.as_slice().unwrap(); - // 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 2: Solve AA' z = e and compute z + let z = self.factorizer.solve(e_slice).unwrap(); - // Step 4: x <-- x - A'(AA')\(Ax-b) - x.iter_mut().zip(w.iter()).for_each(|(xi, wi)| *xi -= wi); + // Step 3: Compute x = x - A'z + let at_z = a.t().dot(&ArrayView1::from(&z[..])); + for (xi, corr) in x.iter_mut().zip(at_z.iter()) { + *xi -= *corr; + } } /// Affine sets are convex. diff --git a/src/constraints/halfspace.rs b/src/constraints/halfspace.rs index 5d9d4a38..5442ca54 100644 --- a/src/constraints/halfspace.rs +++ b/src/constraints/halfspace.rs @@ -72,7 +72,7 @@ impl<'a> Constraint for Halfspace<'a> { /// # Arguments /// /// - `x`: (in) vector to be projected on the current instance of a halfspace, - /// (out) projection on the second-order cone + /// (out) projection on the second-order cone /// /// # Panics /// diff --git a/src/constraints/hyperplane.rs b/src/constraints/hyperplane.rs index d3e6dd46..fa466c61 100644 --- a/src/constraints/hyperplane.rs +++ b/src/constraints/hyperplane.rs @@ -69,7 +69,7 @@ impl<'a> Constraint for Hyperplane<'a> { /// # Arguments /// /// - `x`: (in) vector to be projected on the current instance of a hyperplane, - /// (out) projection on the second-order cone + /// (out) projection on the second-order cone /// /// # Panics /// diff --git a/src/core/fbs/fbs_cache.rs b/src/core/fbs/fbs_cache.rs index de134799..83dd61bb 100644 --- a/src/core/fbs/fbs_cache.rs +++ b/src/core/fbs/fbs_cache.rs @@ -43,7 +43,7 @@ impl FBSCache { work_u_previous: vec![0.0; n.get()], gamma, tolerance, - norm_fpr: std::f64::INFINITY, + norm_fpr: f64::INFINITY, } } } diff --git a/src/core/panoc/panoc_cache.rs b/src/core/panoc/panoc_cache.rs index 3f3c0b13..9bd38f70 100644 --- a/src/core/panoc/panoc_cache.rs +++ b/src/core/panoc/panoc_cache.rs @@ -72,7 +72,7 @@ impl PANOCCache { u_plus: vec![0.0; problem_size], gamma: 0.0, tolerance, - norm_gamma_fpr: std::f64::INFINITY, + norm_gamma_fpr: f64::INFINITY, lbfgs: lbfgs::Lbfgs::new(problem_size, lbfgs_memory_size) .with_cbfgs_alpha(DEFAULT_CBFGS_ALPHA) .with_cbfgs_epsilon(DEFAULT_CBFGS_EPSILON) @@ -158,7 +158,7 @@ impl PANOCCache { /// It checks whether: /// - the FPR condition, `gamma*||fpr|| < epsilon` , /// - (if activated) the AKKT condition `||gamma*fpr + (df - df_prev)|| < eps_akkt` - /// are satisfied. + /// are satisfied. pub fn exit_condition(&self) -> bool { self.fpr_exit_condition() && self.akkt_exit_condition() } diff --git a/src/core/solver_status.rs b/src/core/solver_status.rs index 9098d4bf..e4f7a138 100644 --- a/src/core/solver_status.rs +++ b/src/core/solver_status.rs @@ -32,7 +32,7 @@ impl SolverStatus { /// the specified tolerance /// - `num_iter` number of iterations /// - `fpr_norm` norm of the fixed-point residual; a gauge of the solution - /// quality + /// quality /// - `cost_value` the value of the cost function at the solution /// pub fn new( diff --git a/src/lib.rs b/src/lib.rs index 2363cfc1..e0f61bb0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -53,11 +53,13 @@ pub enum SolverError { pub type FunctionCallResult = Result<(), SolverError>; pub mod alm; +pub mod cholesky_factorizer; pub mod constraints; pub mod core; pub mod lipschitz_estimator; pub mod matrix_operations; +pub use crate::cholesky_factorizer::{CholeskyError, CholeskyFactorizer}; pub use crate::core::fbs; pub use crate::core::panoc; pub use crate::core::{AlgorithmEngine, Optimizer, Problem}; diff --git a/src/lipschitz_estimator.rs b/src/lipschitz_estimator.rs index 22d45e43..d8657062 100644 --- a/src/lipschitz_estimator.rs +++ b/src/lipschitz_estimator.rs @@ -74,13 +74,13 @@ where /// # Arguments /// /// - `u_` On entry: point where the Lipschitz constant is estimated, - /// On exit: the provided slice is modified (this is why it is a mutable - /// reference). The value of `u_` at exit is slightly perturbed. If you need - /// to keep the original value of `u_`, you need to make a copy of the variable - /// before you provide it to this method. + /// On exit: the provided slice is modified (this is why it is a mutable + /// reference). The value of `u_` at exit is slightly perturbed. If you need + /// to keep the original value of `u_`, you need to make a copy of the variable + /// before you provide it to this method. /// - `f_` given closure /// - `function_value_` externally allocated memory which on exit stores the - /// value of the given function at `u_`, that is `f_(u_)` + /// value of the given function at `u_`, that is `f_(u_)` /// /// # Returns /// diff --git a/src/matrix_operations.rs b/src/matrix_operations.rs index f1eafe08..5420eb5b 100644 --- a/src/matrix_operations.rs +++ b/src/matrix_operations.rs @@ -36,6 +36,7 @@ //! ``` //! +use ndarray::ArrayView2; use num::{Float, Zero}; use std::iter::Sum; use std::ops::Mul; @@ -140,6 +141,55 @@ where !a.iter().any(|&xi| !xi.is_finite()) } +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +/// Matrix error +pub enum MatrixError { + /// Matrix dimension mismatch + DimensionMismatch, +} + +/// Computes the matrix product `A * Aᵀ` for a matrix `A` stored in row-major order. +/// +/// The input slice `a` is interpreted as a matrix with shape `rows × cols`. +/// The returned matrix has shape `rows × rows` and is also stored in row-major order. +/// +/// # Arguments +/// +/// * `a` - Input matrix data in row-major order +/// * `rows` - Number of rows of `A` +/// * `cols` - Number of columns of `A` +/// +/// # Returns +/// +/// Returns a vector containing the product `A * Aᵀ` in row-major order. +/// +/// If `a.len() != rows * cols`, the function returns +/// `Err(MatrixError::DimensionMismatch)`. +/// +/// +/// +/// # Notes +/// +/// * The result is symmetric by construction. +/// * This implementation computes the full matrix explicitly. +/// * For better performance, a specialized version can compute only one triangle +/// and mirror it. +pub fn mul_a_at( + a: &[T], + rows: usize, + cols: usize, +) -> Result, MatrixError> { + if a.len() != rows * cols { + return Err(MatrixError::DimensionMismatch); + } + let a_mat = + ArrayView2::from_shape((rows, cols), a).map_err(|_| MatrixError::DimensionMismatch)?; + let out = a_mat.dot(&a_mat.t()); + let (vec, offset) = out.into_raw_vec_and_offset(); + debug_assert_eq!(offset, Some(0)); + Ok(vec) +} + /* ---------------------------------------------------------------------------- */ /* TESTS */ /* ---------------------------------------------------------------------------- */ @@ -241,4 +291,45 @@ mod tests { let norm2sq = matrix_operations::norm2_squared_diff(&x, &y); unit_test_utils::assert_nearly_equal(190., norm2sq, 1e-10, 1e-12, "norm sq diff"); } + + #[test] + fn t_matmul_a_at_tall() { + let a = vec![1.0_f64, 2.0, 3.0, 4.0, 5.0, 6.0]; + let aat = matrix_operations::mul_a_at(&a, 3, 2).unwrap(); + let expected = vec![5.0_f64, 11., 17., 11., 25., 39., 17., 39., 61.]; + unit_test_utils::nearly_equal_array(&expected, &aat, 1e-10, 1e-12); + } + + #[test] + fn t_matmul_a_at_fat() { + let a = vec![1.0_f64, 2.0, 3.0, 4.0, 5.0, 6.0]; + let aat = matrix_operations::mul_a_at(&a, 2, 3).unwrap(); + let expected = vec![14.0_f64, 32., 32., 77.]; + unit_test_utils::nearly_equal_array(&expected, &aat, 1e-10, 1e-12); + } + + #[test] + fn t_matmul_a_at_column_vec() { + let a = vec![1.0_f64, 2.0, 3.0]; + let aat = matrix_operations::mul_a_at(&a, 3, 1).unwrap(); + let expected = vec![1.0_f64, 2., 3., 2., 4., 6., 3., 6., 9.]; + unit_test_utils::nearly_equal_array(&expected, &aat, 1e-10, 1e-12); + } + + #[test] + fn t_matmul_a_at_column_rowvec() { + let a = vec![1.0_f64, 2.0, 3.0]; + let aat = matrix_operations::mul_a_at(&a, 1, 3).unwrap(); + unit_test_utils::nearly_equal(14.0, aat[0], 1e-10, 1e-12); + } + + #[test] + fn t_matmul_a_at_wrong_dimension() { + let a = vec![1.0_f64, 2.0, 3.0, 4.0, 5.0, 6.0]; + let result = matrix_operations::mul_a_at(&a, 2, 2); + assert_eq!( + result, + Err(matrix_operations::MatrixError::DimensionMismatch) + ); + } }