Optimization Without Derivatives: Prima Fortran Version and Inclusion in SciPy

I cherish the following comment from you:

> if SciPy accepts this modern Fortran codebase, then I think the tide for Fortran will change.

I believe that you and the LFortran project are truly turning the tide on Fortran. I wish to see the production release of LFortran as soon as possible.

If PRIMA can also contribute to “turning the tide on Fortran”, I will feel honored. All the pain, depression, and frustration that I have experienced in the past three years will then be partially compensated.

Let us the community work together to get PRIMA included in SciPy and see whether the tide on Fortran will improve.

Many thanks.

https://github.com/uqfoundation/mystic

That's usually a bad assumption for wrapped math libraries outside of the foundational ones like numpy. They tend to have a pretty low ratio of "mileage" to complexity in the first place, it's easy for the project to become a zombie without users realizing, and the layers of challenges to fixing issues leave many of them unresolved. If you need reliability, you'll have to deliberately choose software that's mature and actively maintained.

I kinda tried this a few years back, but found the documentation/lack of modern tutorials appalling ... but there's been a nice explosion of activity over the last few years (e.g. fpm, lFortran), which makes me think this may have changed considerably.

Also, any tips on nice ways/resources that helped you get (back) into it?

https://arxiv.org/pdf/2302.13246.pdf

Upon examining the original FORTRAN code and accompanying papers, I realized that the complexity of the task was far from trivial. Working alone for three years, especially on such a complex codebase, requires a significant amount of dedication and perseverance.

“There do exist "translations" of Powell's Fortran 77 code into other languages. For example, NLopt contains a C version of COBYLA, NEWUOA, and BOBYQA, but the C code in NLopt is translated from the Fortran 77 code straightforwardly, if not automatically by f2c, and hence inherits the style, structure, and probably bugs of the original Fortran 77 implementation.”

See also https://github.com/stevengj/nlopt/issues/501 , where the author of NLopt talks.

Shameless plug: I happen to have recently written a package for the opposite limit! It finds roots when you can only measure the derivative.

https://github.com/EFavDB/inchwormrf

Sometimes you really don't know how to calculate the derivatives.

But perhaps more often, even if in principle you could do the additional work, formulate how to calculate derivatives, and add the required extra code, you just don't want to invest the additional resources, when you can just call a derivative-free optimization library. The optimization will run slower, but perhaps you have some extra seconds so spare.

Many common functions have non-derivative friendly discontinuities. Logic functions, rounding, conditional clauses, sawtooth functions, etc.)

Other non-derivative functions have flat constant zero-derivative areas. Like logic functions & rounding (again), etc.

Some functions have derivatives but have so many minima (or maxima) that non-optimal or non-acceptable minima are likely to trap learning. Sine and cosine waves, for instance.

Some functions even have repeating discontinuities going through +/-infinity. Tangent and secant.

Or, in some cases a model may exist of the system (or part of the system) to be optimized, which can be used to calculate the system response, but it contains relationships that are not easily accessible, so accurate derivatives are difficult or time costly to estimate. Since there is no symbolic expression to allow for efficient calculation of the derivatives.

Trying to use a derivative based optimization where a smooth derivative doesn't exist will likely result in getting stuck at a cliff where there's no change to either the function value or to the gradient.

And as a tangential gripe, I wish js had operator overloading. In python, I could just call my existing function with a symbolic argument and get a symbolic output to differentiate, because my `x+y*z` code would work on both kinds of arguments.

In my experience, it's because the client is trying to reuse an existing simulation code for something that requires optimization. This could be something like imaging (parameter estimation), calibration, or control. For example, say someone has something like a PDE solver that can simulate sound waves going through a media. If you can match the output of this simulation to some collected data, the parameters for the simulation will tell you what material is where.

Anyway, the issue is that these existing codes are rarely instrumented or derived in a way to make it easy to get derivatives. Management doesn't want pay for a rewrite, so they demand that existing equity (the simulator) be used. Then, a derivative free optimizer (DFO) gets slapped on as a first attempt. To be clear, this doesn't work particularly well, but it tends to be the place where I most often see these algorithms used.

As to why it doesn't work very well, derivative based optimizers can scale to the order of hundreds of millions of variables without too many issues (modulo the kind of problem, but it works very well for a lot of problems of value). DFOs tend to have issues past a few dozen variables.

Now, for anyone not wanting to end up in this conundrum, when you write your original simulation code parameterize on the floating point type. Basically, don't write the code with double, use some kind of generic, template parameter, or whatever on something we'll call MyFloat. The purpose here is to allow the code to be more easily compatible with automatic differentiation (AD) tools, which will give a real and machine precision exact derivative. Don't do finite differences. Don't think reinstrumenting later is easy. It's not. It's a pain. It's expensive. And, to be clear, AD can be done poorly and it can be slow. However, I've never seen a code instrumented with AD and a good gradient based optimizer do more poorly than a derivative free optimization code. They've always done dramatically better. Beyond just giving the derivative, I'll also mention that even if the client wants to recode the simulator to give an analytic derivative for speed, it is vastly easier to do so if the results can be checked with a simulator instrumented with AD. There's also other tricks of value that can be used in a properly instrumented code such as interval arithmetic for some kinds of sensitivity analysis.

Anyway, there are other reasons to use derivative free optimization. Some functions just flat out don't have derivatives. However, in the world of mechanics, most do and DFOs gets abused due to poor prior planning. In either case, please, for the love of math, just parameterize on your floating point type. It would make my job much easier.

Unfortunately for many problems, they are either unavailable or expensive to compute. (Imagine a case when your “function” is a full CFD simulation or a 2 hr job — every time an optimizer iterates it has to run an entire simulation. You can’t see the function and you can’t really calculate the derivatives — except crudely through finite differences).

Derivative free optimizers choose to forgo derivative information and instead try to efficiently search the space to optimize the objective. Their convergence is often much slower than derivative based optimizers but they are also much simpler to implement and tend to not get stuck in valleys. All they really are are a way to cleverly evaluate f(X) to find the optimum, without needing f’(X) or f’’(X).

Caveat: derivative free optimizers are terrible for constrained optimization however (especially if your constraints are anything but simple bounds). You can use barrier methods but still they’re not great.

Also there’s an element of luck involved. It’s hard to benchmark derivative free methods because the methods themselves don’t matter too much (Nelder Mead Simplex method outperforms Powell’s method in some cases but not others. Powell’s method is well known but isn’t categorically better than any other algorithm) because when you don’t use derivatives, other higher order effects like initial point and shape of terrain dominate over the the choice of algorithm. I wouldn’t be too hung up on picking the right algorithm — instead I would focus my efforts on heuristics to get to a good initial point cheaply.

I'm trying to come up with a simple example. I think that for bunch of flying hard spheres in a box, any measurement that requires sufficiently time spaced samples of their positions will have trouble being differentiated.

[5] D. M. Bates, M. Mächler, B. Bolker, and S. Walker. Fitting linear mixed-effects models using lme4. J. Stat. Softw., 67:1–48, 2015.

[25] F. Gallard et al. GEMS: a Python library for automation of multidisciplinary design optimization process generation. In 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Kissimmee, FL, USA, 2018. AIAA.

[33] H. Izadinia, Q. Shan, and S. M. Seitz. IM2CAD. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5134–5143, San Juan, PR, USA, 2017. IEEE.

[40] K. Mombaur, A. Truong, and J. P. Laumond. From human to humanoid locomotion—an inverse optimal control approach. Auton. Robot., 28:369–383, 2010.

[39] G. A. Mamon, A. Biviano, and G. Boué. MAMPOSSt: modelling anisotropy and mass profiles of observed spherical systems I. Gaussian 3D velocities. Mon. Not. R. Astron. Soc., 429:3079–3098, 2013.