ForTimize: Trying to standardize and package different optimization procedures

Help
1

Hello! After having to deal with different optimization methods to solve one of my problems I became really tired of:

  • FInding source code / fpm libraries to use them
  • Restructuring the code to fit the specific optimization routine
  • Learning the special variables that some optimizers that don’t provide sane defaults have

During all that work (my research is not based on the optimization procedures, I just use them as a tool) I couldn’t stop thinking how nice it would be to have them all packaged in the same place, with a standard API to avoid all the extra trouble. Something like what scipy.optimize.minimize does, an easy-to-start API but with the possibility to specify details if the user knows about how the algorithms are implemented.

That’s why I’ve started drafting ForTimize (GitHub - fedebenelli/ForTimize). ForTimize intend to ease the usage of optimizers for the end-user. Right now it has a simple functionality like

program main
    use ForTimize, only: pr, minimize
    use my_objective, only: foo

    real(pr) :: x(3), F

    ! Initial guess
    x = [1, 2, 5]

    ! Minimize uses the Nelder-Mead algorithm as a default    
    call minimize(foo, x, F)
    
    ! Print results
    print *, x
    print *, F
end program

And it is also possible to use other algorithms just by adding the optional

type(MyOptimizer) :: optim
optim%some_setting = 24
call minimize(foo, x, F, optim=optim)

All the Optimizers must inherit from an abstract type, that ensures that the same API will be used. The base type is defined in the core.f90 file
Right now I’ve only implemented gradient descent and Nelder-Mead, as examples of two types of algorithms (with and without derivatives). Adding new Optimizers is easy, and now with fpm it is now possible to wrap all the other implemented algorithms and using them as dependencies

Since this is far from something I really need right now, just trying to help the FOSS Fortran community, I wanted to discuss your opinions on something like this. Any kind of code reviewing, now that the codebase is minimal and easy to follow and change, is welcomed too!

4 Likes
2

I’ve summarized some of my thinking on this topic in a few places:

I’ve gone down this path a bit too. Digging through my long list of stale private repositories there is one with the following README:

optimization-shed

Ivan’s optimization shed - a workspace for tinkering with optimization codes in Fortran, C, and C++

My goals for this repository include:

  • defining one or more canonical interfaces for (Fortran) callback routines used in mathematical optimization
  • experimenting with design patterns to develop robust, composable, flexible, and extensible solver libraries
  • develop adaptor patterns for mixed-language programming in Fortran/C/C++

Instead of re-implementing algorithms from scratch, we can start by integrating and evaluating existing libraries. The discovery process should feed information into achieving the goals above.

3 Likes
3

Thanks for your inputs! I’ve been using nlopt-f for some time on my projects and it was what inspired me to use a class(*) object to include the additional data for the function. Unfortunately, I can’t install nlopt on the system I’m working on now and that made me start writing ForTimize. While I think the approach that nlopt-f is one of the best of I’ve seen in this subject, I (personally) don’t like how the objective function must be defined as a nlopt_func type. I understand why it is needed but I think it sets a (little) barrier to a new user.

I know that my implementation needs (a lot) more polish, but I think we should all try to converge to an API that is easy to use and intuitive for a new user. Just as how I can now use solve(A, b) with stdlib instead of calling dgsev directly. I think that developing simple APIs for the things that Fortran works best (math problems) would help enormously to catch new users. If the Julia community could obtain this relatively quickly, why is it still “hard” for us?

4

I’ve played with writing a little wrapper on top of NLopt, which would make it a bit more generic:

program nlopt_alternative_demo

   use example_funcs
   use nlopt_alternative
   implicit none

   type(nlopt_prob(n = 2)) :: prob

   type :: constraint_data
      real(wp) :: a, b
   end type

   type(nlopt_func) :: cnstr(2)
   type(constraint_data), target :: d1, d2

   call create(prob, "LD_MMA")  ! initializes with default settings

   prob%lower_bounds(2) = 0.0_wp
   prob%xtol_rel = 1.e-4_wp

   prob%min_objective = nlopt_func(myfunc)

   d1 = constraint_data(2.0_wp, 0.0_wp)
   d2 = constraint_data(-1.0_wp, 1.0_wp)

   cnstr(1) = nlopt_func(myconstraint, d1)
   cnstr(2) = nlopt_func(myconstraint, d2)

   call add_inequality_constraints(opt, cnstr)

   x = [1.234_wp, 5.678_wp]

   call optimize(prob, x)

   print *, "Found minimum ", prob%optf, "at ", x
   print *, "Number of function evaluations ", prob%numevals

end program

The main thing to observe is the existence of an adapter that takes myfunc and myconstraint which have an interface different from what NLopt expects, but converts those into nlopt_funcs. Unfortunately, there are some caveats related to pointers, and the copying/ownership of data.

I fully agree.

1 Like
5

Thanks for the project. Looking at your module below, I wonder if

  1. The data argument needs to be intent(in out) rather than intent(in), since typically evaluation of the objective function should not modify data.
  2. The optimize subroutine should have optional arguments max_iter and tol specifying the maximum number of X values the optimizer can try and a convergence criterion based on changes on the objective function.
  3. The optimize subroutine could have optional arguments Xtry, Ftry, and dFtry containing arrays of the values of X that the optmizer tried and the values of F and dF at those tries. When an optimizer is not working, the user may want to see what values of X it tried. Alternatively, an iprint argument could be supplied to determine what is printed within the optimizer.
  4. The objective subroutine could have optional arguments Xmin and Xmax which if present give bound constraints on X.
  5. The objective_function subroutine could have an optional argument con(:) which if present has the values of constraints that must be non-negative.

I understand that all the suggestions other than the first may add clutter, but they could also increase generality.

module ForTimize__base
   use ForTimize__constants, only: pr
   implicit none

   type, abstract :: Optimizer
   contains
      procedure(optimize), deferred :: optimize
   end type Optimizer

   abstract interface
      subroutine objective_function(X, F, dF, data)
         import pr
         real(pr), intent(in) :: X(:)
         real(pr), intent(out) :: F
         real(pr), optional, intent(out) :: dF(:)
         class(*), optional, intent(in out) :: data
      end subroutine objective_function
   end interface

   abstract interface
      subroutine optimize(self, foo, X, F, data)
         import pr, objective_function, Optimizer
         class(Optimizer), intent(in out) :: self
         procedure(objective_function) :: foo
         real(pr), intent(in out) :: X(:)
         real(pr), intent(out) :: F
         class(*), optional, target, intent(in out) :: data
      end subroutine optimize
   end interface
end module ForTimize__base
6

Thanks for providing detailed feedback!

In 1 I agree that initially having it as intent(in) might be better, but the intent(in out) provides more flexibility to the user. I based this on personal cases when I might reuse most of the objective function logic but need some flexibility to modify the models I fit. This could also let the user save Xtry and Ftry in data at each function evaluation of they desire (but including it on the optimizer API could also be a good addition) without adding more clutter when implementing the optimizer. Maybe having two separate arguments mutable_data and immutable_data would be better?

Concerning the other points, I do think that all those kinds of parameters should be included, but I’m not sure if the best option would be to include them as arguments. I believe that the best option would be to have them as attributes of the base derived type. This would help keep everything encapsulated on the derived type and keep the general API simple.