Optimization Without Using Derivatives: the PRIMA Package, its Fortran Implementation, and Its Inclusion in SciPySkip to main content

As mentioned two weeks ago under the thread of LFortran 27, I have been developing a package named PRIMA 84 for solving general nonlinear optimization problems without using derivatives. @certik suggested that I should write about PRIMA, especially its inclusion in SciPy.

I think if you succeed to get your modern Fortran code into SciPy, that might turn the tide on Fortran. So I think this is a big deal. What needs to happen to get it done? Ask for help on this forum, just start a thread, create a plan and ask people to join your efforts. With enough people it will go quick.

What is PRIMA?

PRIMA is a package for solving general nonlinear optimization problems without using derivatives. It provides the reference implementation of Powell’s renowned derivative-free optimization methods, i.e., COBYLA, UOBYQA, NEWUOA, BOBYQA, and LINCOA. The “P” in the name stands for Powell 10, and “RIMA” is an acronym for “Reference Implementation with Modernization and Amelioration”.

PRIMA is a project that I have been working intensively on for the past three years. Almost all questions I posted on this forum come from PRIMA.

Who was Powell?

Michael James David Powell FRS 17 was “a British numerical analyst who was among the pioneers
of computational mathematics” 2. He was the inventor/early contributor of quasi-Newton method, trust region method, augmented Lagrangian method, and SQP method. Each of them is a pillar of modern numerical optimization. He also made significant contributions to approximation theory and methods 3, but I hesitate to mention more details as it is not my area.

Why PRIMA?

Professor Powell carefully implemented his derivative-free optimization methods into publicly available solvers. They are widely used by engineers and scientists. For instance, see Section 1 of a recent paper on Powell’s solvers 37 as well as the Google searches of COBYLA 30 and BOBYQA 10.

However, Professor Powell’s implementation was done in Fortran 77 5. The code is nontrivial to understand or maintain, let alone extend. For many practitioners, this has become an obstacle to exploiting these solvers in their applications. More seriously, bugs 25 keep emerging, but few of them get really fixed — it is highly challenging/frustrating to debug a maze of 244 GOTOs in 7939 lines of classic-style Fortran 77 code, especially if you do not know the sophisticated algorithms behind the code, which are nontrivial to understand by themselves.

Before he passed, Professor Powell had asked me and Professor Nick Gould 13 to maintain his solvers. This is an honorable mission. To make the solvers more accessible, I started PRIMA. It is a project somehow similar to the translation, interpretation, and annotation of Euclid’s Elements. It will make Powell’s solvers easily understandable to everyone, not only the experts. Few people remember who translated Elements 31, but it is a job that must be done.

Objective

PRIMA aims to provide the reference implementation of Powell’s methods in modern languages, first in modern Fortran 10 (F2008 or newer), and then in MATLAB, Python, C++, Julia, and R. It will be a faithful implementation, in the sense that the code will be mathematically equivalent to Powell’s, except for the bug fixes 25 and improvements 4 made intentionally.

The focus is to implement these methods in a structured and modularized way so that they are understandable, maintainable, extendable, fault tolerant, and future proof. The code will have no GOTO (of course) and will use matrix-vector procedures instead of loops whenever possible. In doing so, PRIMA codes the algorithms in a way that we would present them on a blackboard.

(N.B.: There are plenty of discussions on whether we should use matrix-vector procedures or loops in Fortran. Here, the former is the correct one by all means, as long as we note one fact — in derivative-free optimization, the dominating cost comes from the function evaluations implemented by the users rather than the numerical linear algebra of the solver. Each function evaluation takes minutes or months to finish.)

Current status

Modern Fortran

After almost three years of intensive coding, the modern Fortran version 67 of PRIMA has been finished by December 2022.

MATLAB

• An interface 1 is provided for using the modern Fortran implementation in MATLAB 3.
• A pure MATLAB version of NEWUOA 8 is implemented. It was generated straightforwardly (indeed, automatically) from an earlier version of the modern Fortran code.

Python

• The inclusion of PRIMA into SciPy is under discussion 66. It will replace the buggy 25 and unmaintained Fortran 77 version of COBYLA underlying scipy.optimize.minimize 27, and make the other four solvers available to all SciPy users.
• My ultimate objective is to have a native Python implementation of PRIMA independent of Fortran, similar to what we have done with NEWUOA in MATLAB as mentioned above.

Other languages

• Interfaces for using the modern Fortran implementation in other languages will be available later.
• Given the modern Fortran version, native implementations in other languages become much easier, because we now have a structured and modularized implementation as a reference. My team will implement the methods in other languages in this way. I am exploring how LLMs like ChatGPT can accelerate our progress in this process.

What is needed now?

Fortran

• Make PRIMA available with FPM.
• Generating documentation using FORD or other tools.

Python

How to include PRIMA in SciPy as soon as possible? This is the question. The major Scipy maintainers are positive about the inclusion of PRIMA solvers in SciPy. See the discussions on GitHub 66 for details. However, I do not know Python at all. So community efforts are greatly needed here.

The first step should be replacing the Fortran 77 implementation of COBYLA in SciPy 27 with the PRIMA version, respecting the current Python signature of COBYLA. The second step will be getting BOBYQA into SciPy, as this solver has been long requested by SciPy users 8. The third step is to get the other three solvers included.

Previously, my former student Tom @ragonneau and I developed Python interfaces 3 for the Fortran 77 implementation of Powell’s solvers under the PDFO project 2. They may provide a good starting point to begin with.

@zaikunzhang excellent, thanks for your efforts! You should try to get PRIMA into SciPy. Is that the first “modern” Fortran in SciPy? I don’t think SciPy uses any Fortran “modules” yet, correct?

I do not think there is any modern Fortran code with modules 5 in the SciPy code base, or even any .f90 code 2.