Automatic differentiation in Fortran

[ivan-pi]

Recently, there has been a growth of interest in automatic differentiation tools used in adjoint modelling. Some popular projects in the Python and Julia worlds are JAX and
JuliaDiff, respectively.

Of course, plenty of such tools exist for other languages including (older) Fortran ones such as adjac, OpenAD, Tapenade, AUTO_DERIV, and ADIFOR just to name a few. At NAG they also support some algorithmic differentiation tools in the form of overloaded C++ classes.

The growth of interest in algorithmic differentiation appears to be fueled by the needs of the machine learning community. The Julia language developers are trying to push it so far to even have a differentiable programming language. A recent blog post touches on this topic from the perspective of a Fortran programmer (the idea for submitting opening this issue/proposal comes from reading this blog post).

I have no experience in this area (just plans to have a closer look in the future), but I am wondering whether some of the older tools (for example ADIFOR dates back to 1992) have achieved a level of maturity, where they might be implemented directly into the language. I am guessing the compiler developers already employ some experts in this field. It would be interesting to see whether an experimental branch of the LFortran compiler could make this idea take off. I can imagine that making these tools easy to use (while keeping the performance) would generate some new interest in Fortran and also send the message that Fortran is a living language (otherwise, we might end up in the situation where those rare Fortran programmers still left start calling ODE solvers through the Julia C API or simply abandoning ship).

Given that Fortran was named and developed for the purpose of "formula translation", one could imagine that in the current decade (or perhaps century) we also continue down this path of adding further mathematical abstractions such as differentiation. A related idea would be to introduce a new type for dual numbers, which can also be used for the purpose of automatic differentiation (see here and here).

I look forward to reading the comments from those of you who have used/authored such automatic differentiation tools in the past.

[jacobwilliams]

1000 thumbs up for this!

I see this as a major use case for some kind of generic programming/template feature. But it has to be well-designed. There are things we just cannot do right now.

For example: say we have a function that depends on an intrinsic real variable. This function could be using hundreds of modules, third-party codes, etc. I want to be able to call this function to get the result, but I also want to call the function and get the derivative, by replacing this real variable with some kind of autodiff class. I want both of these calls to be as efficient as they can be. Maybe this class is not even known about by some of the third-party code. I think the Julia multiple dispatch thing is the closest thing I've seen that can do something like this, but I've not actually used it before. Fortran needs that desperately. I don't want to have to clutter potentially hundreds of thousands of routines to put select type everywhere, or put all the logic in include files and duplicate all the interfaces. That's just not practical for anything complex.

[certik]

I have an issue open for this at LFortran: https://gitlab.com/lfortran/lfortran/issues/97. In the coming months I'll write a prototype, so that people can play with it, then we can go from there.

[LKedward]

This would be amazing!
If Fortran became a language with AD built in, it would be a game changer for scientific and numerical computing where derivatives are so important.
As you point out there are tools for performing AD in Fortran, however they are all old (don't support modern Fortran) and/or not open source.
As far as language modifications go, it's very complex, though there is a lot of research in the field.
I believe it can already be implemented simply (and inefficiently) in Fortran using operator overloading, but it really needs a source-code-transformation approach or equivalent if you want an efficient derivative code.
I don't know how this can be achieved at a language level, but it would set Fortran far apart if implemented.

[jacobwilliams]

I dare not dream that we'll ever get anything like this out of the Fortran committee. Which is why I think the best we could hope for is a really good generic programming capability so users could implement this themselves with (hopefully) a minimum of verbosity.

Having this baked into LFortran would also be amazing! That would really be a game changer.

[difference-scheme]

Yes, this can already be done today (to some extent) with the features already present in the language, by declaring a type for dual numbers and using operator overloading to declare arithmetic operations on these numbers (akin to how one would implement complex numbers if the language didn't provide them). The problem is it will be slow, not least due to the arithmetic overhead inherent in, e.g., dual number division (six times more expensive than division of reals).

[jacobwilliams]

@difference-scheme The main reason this is unworkable for anything complex is that you can't use any third-party library for any calculation unless it is also using your dual number class. There are a dozen different sort of classes you might devise to get derivatives besides dual numbers. All third party libs (say, an ODE solver or whatever) can't be expected to support all of them. The only way to make this workable is:

• Autodifferentiation is somehow baked into the language, so we don't have to worry about writing any custom classes.
• The language provides some amazing way to implement generics so that all third-party libs will be able to be written in such a way that they can be used for auto-differentiation without having to know about the particular implementation a user might use.

The Fortran committee needs to be looking at what Swift and Julia are doing. In this field, as in many other ways, Fortran is falling behind...

[ivan-pi]

Perhaps of interest: Algorithm 1005: Fortran Subroutines for Reverse Mode Algorithmic Differentiation of BLAS Matrix Operations. I could not find the actual Fortran routines. But it shows at least there is still large interest in having the possibility to automatically differentiate Fortran programs that rely on BLAS. Hopefully, LAPACK follows (in the paper they already show discuss algorithmic differentiation of the Cholesky decomposition).

This blog post by Christopher Rackauckas, The Essential Tools of Scientific Machine Learning (Scientific ML), also offers a summary of the automatic differentiation field with a focus on bringing together machine learning and ODE/PDE solvers.

[ivan-pi]

The ISO C++ committee is looking at differentiable programming in one of their study groups: https://docs.google.com/document/d/1_5TJCBvI6fZSdyuK7_Cpo5XwdoIS24DDbEPh2JjqQsg/edit?ts=5cf823e2#heading=h.t1fxx74w88nq

It seems as if the Fortran community is starting to fall behind on this issue compared to Julia and C++.

[ivan-pi]

Perhaps of interest: Algorithm 1005: Fortran Subroutines for Reverse Mode Algorithmic Differentiation of BLAS Matrix Operations. I could not find the actual Fortran routines. But it shows at least there is still large interest in having the possibility to automatically differentiate Fortran programs that rely on BLAS. Hopefully, LAPACK follows (in the paper they already show discuss algorithmic differentiation of the Cholesky decomposition).

The routines are now available on the ACM page: https://dl.acm.org/doi/10.1145/3382191

[dalon-work]

CppCon 2021 recently posted a talk about implementing AD in LLVM as a plugin, using the compiler directly, instead of any language-level constructs. I assume this would affect Flang as well.

https://youtu.be/1QQj1mAV-eY?t=2240

At the end he describes an effort to make this part of the C++ standard somehow.

[JuanDavidNavarro93]

We have implemented dual and hyper-dual numbers along with other kind of hypercomplex numbers in a library to perform AD in FORTRAN and PYTHON called HYPAD. Here is the link of the webpage: https://ceid.utsa.edu/HYPAD/

Stay tuned for further developments.