datamole-ai/gomez

Framework and implementation for mathematical optimization and solving non-linear systems of equations.

gomez's Issues

`TrustRegion` solver has issues when using low-precision (`f32`) type as `Problem::Scalar`

This solver works:

struct SolveForOrigin;

impl Problem for SolveForOrigin {
    type Scalar = f64;
    type Dim = na::U2;

    fn dim(&self) -> Self::Dim {
        na::U2::name()
    }
}

impl System for SolveForOrigin {
    fn eval<Sx, Sfx>(
        &self,
        x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
        fx: &mut na::Vector<Self::Scalar, Self::Dim, Sfx>,
    ) -> Result<(), ProblemError>
    where
        Sx: na::storage::Storage<Self::Scalar, Self::Dim> + IsContiguous,
        Sfx: na::storage::StorageMut<Self::Scalar, Self::Dim>,
    {
        fx[0] = x[0];
        fx[1] = x[1];
        Ok(())
    }
}

iter = 1	|| fx || = 0.0000000004734010799665374	x = [3.3331393201052606e-10, 3.361702027859792e-10]
solved

Whereas this solver has issues:

struct SolveForOrigin;

impl Problem for SolveForOrigin {
    type Scalar = f32; // CHANGE is here
    type Dim = na::U2;

    fn dim(&self) -> Self::Dim {
        na::U2::name()
    }
}

impl System for SolveForOrigin {
    fn eval<Sx, Sfx>(
        &self,
        x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
        fx: &mut na::Vector<Self::Scalar, Self::Dim, Sfx>,
    ) -> Result<(), ProblemError>
    where
        Sx: na::storage::Storage<Self::Scalar, Self::Dim> + IsContiguous,
        Sfx: na::storage::StorageMut<Self::Scalar, Self::Dim>,
    {
        fx[0] = x[0];
        fx[1] = x[1];
        Ok(())
    }
}

Error: "neither newton nor steepest descent step can be taken from the point"

Add benchmarks to the documentation

Once the library will be rewritten to faer (#9), the benchmarks results should be added to the documentation comparing gomez with GSL and perhaps argmin (regardless if gomez is better or not).

Using gomez for solving geometric constraints / intecepts

Hi!

I'm wondering if gomez is appropriate to use for implementing a geometric solver (like this)?

Based on the description in the docs:

F(x) = 0,

where F(x) = { f1(x), ..., fn(x) }
and x = { x1, ..., xn }

I'm not sure if each Fn(x) can only be dependent on x, can it be dependent on other variables?

For instance if I have a system of equations that solve for the intersection of a line and a circle, I would need to reference the parameters of the line and circle across the same equation.

NB: Apologies if I'm totally missing the mark of understanding here, I'm learning :)

Tools for debugging

Analyzing why the solver can't find the solution or optimizer can't find the minimum is quite difficult process (at least for not so experienced people, like me) and often feels like black magic. gomez should provide ways how to help with this task.

Currently, there are some log statements which can be enabled by injecting a global logger, but this is not perfect for number of reasons (it's not possible to switch different types of logs (log level is not a good approximation) and it gets very messy in a multi-threaded environment).

In some cases, it is especially hard to understand what's happening. An evidence is some non-standard things (1, 2, 3) inside gomez, which were often discovered by accident/bug, but proved to work much better than correct code in my cases. I would like to understand these things and come up with the real solutions or explain why they work.

Nonetheless, I don't really have a perfect idea what should be done in this regard. I will use this issue to collect ideas and would also appreciate any help/advice/thoughts.

Derivatives methods on the `System` and `Function` traits

Add jacobian method to System trait and gradient and hessian methods to Function trait, with default implementations that use finite differences. Change all solvers to use that instead of finite differences manually.

This will allow the users to specify analytical derivatives if they are available.

Next version

Ever since I discovered the faer linear algebra library I wanted to rewrite gomez into it. For two reasons:

1. Performance. The benchmarks for faer are very impressive. Moreover, the low level API allows to do operations, for which I am now forced to use cloning, in place.
2. API. nalgebra is very flexible with the matrix storage, but it also means that it's very trait-heavy (our public API is unnecessarily complex due to that, and the implementations of solvers are even worse). Interface when using faer will be much simpler.

Half a year ago I actually spent some time on implementing the algorithms using faer, and after initial struggles I managed to do most of them (and some more, and in a much more modular structure). The code that I have is for an old version of faer and a lot of changes were made since that, but they should be mostly cosmetic.

The plan is as follows:

1. Make some improvements to the API (more on that below).
2. Release 0.5.0.
3. Migrate everything to faer.
4. Release 0.6.0.
5. For a limited time period, make patches to 0.5.x branch if needed. This would be for people who want to continue using nalgebra.

Before fully migrating to faer, I would like to do some unrelated (or partly related) improvements to the API. Now I can think of these:

• Make function optimization a first-class citizen. Rename the next methods in Solver and Optimizer traits to solve_next and opt_next, respectively, so that they don't clash (see dd6524a).
• Remove apply_eval from the Function trait. A problem is either a system or function, not both. I guess this was a work around for being able to "solve" some systems by algorithms that only implemented the Optimizer trait. Instead of apply_eval, all Optimizers should also implement the Solver trait (not automatically, but as a guideline).
• Remove Dim associated type from the Problem trait. faer supports only dynamically sized matrices. Removing the type now should be an incremental change towards to final migration. Move the dim method from the Problem trait to the Domain type and make its return type usize.
• Add jacobian method to the System trait and gradient and hessian methods to the Function trait and provide default implementations using finite difference approximations. EDIT: With the current implementation of finite differences in gomez, this would allocate memory on every call. So I will postpone this to after faer rewrite, where the core API will be more functional and low-level.
• Remove ProblemError and don't return Result in System::eval and Function::apply. Encountering invalid dimensionality is a bug in the solver/optimizer not respecting the value in Domain::dim and so panic! is appropriate. Invalid values (infinity, NaN) should be detected by the solver. This will simplify APIs a bit and might also result in faster & smaller code (less branches). Transitively, other error types and corresponding Result returning methods will also be removed (JacobianError, etc.).
• Rename Problem::Scalar associated type to Problem::Field (seems like a better name considering its usage in Rust ecosystem). Introduce custom RealField trait having nalgebra's (later faer's) RealField as a super trait, and make Problem::Field type be constrained by our RealField trait. Move EPSILON_CBRT and EPSILON_SQRT constants to this trait (and use proper values for f32).
• Remove Variable and VariableBuilder and add specific constructors to Domain for unconstrained and rectangular domain. This will make Domain more flexible (we might add support for other domain types, for example spherical domain). But user should still be able to somehow specify custom scale/magnitude for variables.
• Make Problem::domain return &Domain instead of Domain. This will shift the responsibility of keeping the domain instance from the call-site (before creating the solver and throughout the solving process) to the problem type itself. Thanks to that, many functions will not need separate parameter for the domain. EDIT: This would be less convenient on the user side. Instead, keep the current low-level API and focus on high-level API which will hide these details.
• Remove cuckoo search algorithm. I don't think it provides much value, it was a relic from my experimentation in the past. I believe that LIPO or a more widely-used meta-heuristics algorithm (e.g., differential evolution) should be our focus.
• Remove the population module. As cuckoo search will be removed and it is currently the only population-based algorithm in gomez, there is no point in having this module in gomez (now). We may reintroduce it later, potentially in a different shape.
• Remove FunctionResultExt and VectorDomainExt. These are poor API.
• Hide the testing module behind a feature flag.
• Add TestFunction counterpart to the TestSystem in the testing module.
• Rename SolveError in testing module to something more generic (e.g., TestError) so the naming makes more sense in the context of optimization.
• Make Solver and Optimizer the only things exported from the prelude. Since Problem and System/Function are traits to be implemented by the user, it is probably better to explicitly import them. EDIT: Such a small prelude is probably not worthy of having a prelude. Instead, everything should be imported explicitly. Once there is high-level API, it will be encouraged and the Solver/Optimizer traits should not be necessary to import to use the high-level API.
• Introduce SolveIter<'a> and OptimizeIter<'a> which implement Iterator<Item = Result<('a x, 'a fx), S::Error>> (no lifetime 'a for fx in optimization case). Add solve(solver, problem, x, fx) and optimize(optimizer, problem, x) methods to the Solver and Optimizer traits, respectively, with default implementation that will simply construct the corresponding iterator and return it. This will be the high-level API, which still allows quite flexible control over the process (e.g., iter.take(max_iters).find(|result| result.map(|(x, fx)| fx.norm() <= 1e-6).unwrap_or(true))). Nevermind, I forgot that for Iterator::Item be a reference, I would need a lending iterator. So instead I implemented the drivers that to some extent mimic the iterator API.
• Move everything (what remains) from the core module to the root of the library.
• Remove RepulsiveSystem type. It is still a good idea, but it is currently a noop and I am not planning to spend time on it now. We may reintroduce it in the future.
• Migrate from rand to fastrand. It is designed to be fast (whereas in rand choosing fast generator requires conscious "effort") and has a bit simpler API. Nevertheless, still require the generator as an explicit parameter to support determinism (user sending generator initialized with a fixed seed).

After all these changes are done and 0.5.0 is released, then I will migrate the whole library from nalgebra to faer.

• Migrate to faer.

Residuals method on `System` trait

After return-position-impl-trait-in-trait feature gets into the stable, I would like to change the System trait to have the following API (some details omitted):

trait System: Problem {
    fn residuals<'a>(&self, x: &'a Vector) -> impl Iterator<Item = Self::Field> + 'a;

    fn eval(&self, x: &Vector, fx: &mut Vector) {
        // default impl, fill fx using self.residuals(x)
    }

    fn norm(&self, x: &Vector) -> Self::Field {
        // default impl, sqrt of sum of self.residuals(x)
    }
}

The main advantage is that by having residuals iterator the main interface, default implementation of norm can use that instead of eval and thus avoiding the need of allocating a temporary vector for fx. I believe that creating an iterator of residuals will not be more struggle than evaluating the system directly to fx. The iterator approach is already successfully used in the testing module.