Thursday, May 3, 2012

Derivative-free nonlinear optimization for .NET and Java

Optimization of non-linear objective functions and potentially also non-linear constraints is a common problem in for example radiation therapy.

More often than not, the objective and constraint functions are not easily differentiable, which limits the possibility of applying the more efficient derivative-based optimization algorithms (of which IPOPT is a prominent open-source example).

There are several more or less time-efficient ways of solving (constrained) non-linear optimization problems, including simulated annealing, genetic algorithms and direct search methods such as the Nelder-Mead method (a good overview of mathematical optimization with several links is provided here.)

One rather popular code for solving non-linear optimization problems involving (potentially nonlinear) constraints is the COBYLA2 method. Originally formulated and implemented in Fortran 77 by Michael Powell, this method has been ported to both Fortran 90 and C (through f2c conversion), and it has also been integrated into the Python scientific calculations library scipy. COBYLA2 is simple to use, but it exhibits a rather slow convergence rate; the number of function and constraint evaluations required to locate an optimum that meets the constraints is often very large. As a first-attempt solver when objective and constraint function gradients are complex or time-consuming to derive, COBYLA2 is however often a good choice.

Being primarily a .NET developer, it has until recently not been straightforward to incorporate COBYLA2 though. It has of course been possible to build a native DLL and call COBYLA2 via P/Invoke. This works relatively well in a full .NET Framework environment, but it is only barely applicable in Silverlight applications and completely inapplicable in Windows Phone applications. On the other hand, the Fortran code is relatively well structured and well documented, so I recently decided to port COBYLA2 to C#.

I based the porting on the Fortran 90 code, since this implementation already utilizes more C/C++/C# like constructs. The porting has been successful and I have released the end result as an open-source project cscobyla on Github. Included in this project are unit tests demonstrating the use of COBYLA2 in C# applications. The C# code is very straightforward and does not make use of any advanced functions, so it should be possible to incorporate the code in Silverlight and Windows Phone applications just as easy as it can be incorporated in regular .NET applications. When built into a C# class library, it should be straightforward to reference the COBYLA2 optimization from any other .NET language as well.

Encouraged by the successful porting to C#, I then embarked on a Java porting experience! I have not been able to find a pure Java implementation of COBYLA2, so I considered this a good Java development exercise. Java does not provide delegates, multidimensional arrays can only be represented as jagged arrays, and Java does not support the goto statement, so I was forced to make some design changes. Nevertheless, the Java porting effort also succeeded, and I have also made this result available as open source on Github, the jcobyla project.

And this is not all there is! As mentioned, COBYLA2 converges slowly. Michael Powell has made several efforts to overcome this issue for specialized cases by making use of a local quadratic approximation rather than a linear approximation that is being used in COBYLA2. These improvements have been made available in the NEWUOA code for unconstrained optimization of nonlinear objective functions, and more recently for bound constrained optimization of nonlinear objective functions in the BOBYQA (Bound Optimization BY Quadratic Approximation) code. These codes exhibits substantially improved convergence properties compared to COBYLA2, albeit for unconstrained or bound constrained problems only.

In particular, I consider the BOBYQA being a very interesting development, as many problems encountered in for example radiation therapy can be formulated as bound constrained problems. BOBYQA has not been available for the .NET platform other than through P/Invoke, so again I decided to port Fortran code to C#. This time I had to rely on the Fortran 77 implementation, since I have not been able to identify any Fortran 90 port of BOBYQA. It took some additional effort, but even BOBYQA is now available as open source in C#. I have denoted this project csbobyqa and made it available on Github. Also written using standard C# constructs, the code should be easily incorporated in any .NET, Silverlight or Windows Phone application.

In the case of BOBYQA, there actually already is an open-source Java implementation available as part of the Apache Commons Math project. The API for this Java class can be found here

To confirm that my C# implementation is sufficiently efficient, I have also identified the longest-running unit test (bound-constrained Rosen with varying number of interpolation points) in the Java implementation and implemented the corresponding unit test in the C# unit test library. On my laptop, the C# unit test takes around 15 seconds in each consecutive run, whereas the corresponding unit test on the Java implementation ranges from 15 to 40 seconds. I am not sure why the Java timing varies this much, but it is however encouraging that the C# implementation consistently performs equal to or better than the Java implementation in this case.

To re-iterate, here are the links to my most recent open-source optimization projects:

Good luck with the derivative-free optimization!


  1. Interesting work! Did you port the code manually, or is there a tool to assist you with porting Fortran code to ?

    The varying run time of the Java unit test could be explained by the differences between the JVM and the CLR, since the CLR always JITs each method exactly once, while the JVM might interpret the code the first time and only JIT it later when it finds that the method is called more than once. If you include a few warm up iterations before timing the code I'm sure the timings will stabilize.

    1. Thanks Bas for the explanation regarding the Java unit test run times!

      The code is manually ported, I have not been able to find any automated tool for Fortran-to-C# conversion. For Fortran-to-C conversion there is f2c ( Theoretically it might be possible to create an f2c# application from the f2c code base. In practice (and by looking at typical C code generated by f2c), that kind of development hardly seems to be worth the effort.

  2. I have downloaded the java implementation and tested it with the following codes.

    public static void main(String[] args)
    CobylaTest test = new CobylaTest();

    But I keep getting the following error messages. Can you tell me what's happening or better what actions I should take to get the results?

    Output from test problem 1 (Simple quadratic)
    Exception in thread "main" java.lang.VerifyError: Inconsistent stackmap frames at branch target 2864 in method com.cureos.numerics.Cobyla.trstlp(II[[D[DD[D)Z at offset 314
    at com.cureos.numerics.CobylaTest.test01FindMinimum(
    at com.cureos.numerics.Main.main(

  3. Strange, I have not seen that kind of error before.

    The test samples have been written to run within the context of a JUnit test runner, which is more or less automatically managed by Netbeans and Eclipse. I have not tried to run the samples directly from a Main method, but I will give it a try.

    Note that the test code relies on the presence of a JUnit JAR; have you added a reference to this library?

    1. Mr. Gustafsson, thank you for your quick reply. I have added a reference to JUnit4. I understand that without a reference to JUnit, your codes will not simply compile.

    2. Dear Gunsang,

      actually it is only necessary to reference the JUnit4 library if you run the unit test source code unmodified. jcobyla itself is not dependent on JUnit4.

      For example, if you extract one of the unit test methods into your own class, remove the "@Test" annotation and replace the last line "assertArrayEquals(..." with your own output line for example, you will see that you are not at all required to reference JUnit.

    3. Mr. Gustafsson, I found that the errors I reported above have something to do with the Java version 7 bug reports. I changed the runtime to version 6 and the programs work just fine. Thank you for your efforts in translating the Fortran version to Java, and for your kind reply to me. By the way, are there any "answer sheets" availabel on the web for the sample test problems that you wrote?

    4. Oh! It seems that the answers to the test problems are in the assertArrayEquals method of CobylaTest. Please just disregard the last question above. Thanks a lot.

  4. It would be interesting to look at the comparison between Fortran and C# effectiveness. Actually, I'm surprised you haven't done/shown it in your work.

  5. Many thanks for your feedback and sorry for not getting back sooner. The main purpose of this work was to provide a standalone managed implementation of the two optimization algorithms. If speed is the bottleneck, I am fairly convinced that a native implementation should always be the primary choice. However, in a managed environment a 100% managed implementation of the optimization algorithms may very well be a plausible choice, since it would relieve development from mixing managed and native code.

    In a language-agnostic sense, the original algorithms are re-implemented as close to the originals as possible. There are some additional steps of copying arrays back and forth to adapt to the base-1 index scenario in the original FORTRAN implementation. Thus, the number of operations is always slightly larger in the managed implementations, and it is most likely that a sufficiently speed optimized native compilation will always be faster than the managed analog.

  6. Could you give some explanation of the different error codes that are returned by the Bobyqa.FindMinimum method in C#? I see no documentation on them on the various sites where you posted this code.

    In particular, I'm getting "DenominatorCancellation" returned from my method, and I'm not sure what that means or how to correct it.

  7. Could you please specify the meaning of the int argument in Func calfun?