# Basic software libraries for numerical analysis using modern programming languages?

I'm looking for a software library with a scope similar to "numerical recipes", but implemented in a modern programming language. "Modern" in this context means to me: object oriented (not C or Fortran), running in a virtual machine with garbage collection (not C++), and with support for functional programming. Examples therefore would be: F#, C# from the .NET framework, Scala and to a lesser extend Java (no functional programming yet).

"Scope" means that it should cover all algorithms that are taught in, say, a general, one or two year(s) long introduction to numerical analysis, starting with linear algebra, up to partial differential equations and stochastic processes.

In addition I am interested if any research group using high performance computing uses a language that fits the description, and if not, why not. It is "common knowledge" in the software industry that one should no longer worry about performance problems of languages running in a virtual machine compared to languages compiled directly to machine code, e.g. Java versus C++. Any experience from high performance computing that affirms or contradicts this would be interesting to me, too.

Addendum after reading the first three answers (rather than commenting them individually): My motivation for asking about these "modern" programming language is this: It takes some effort to learn a programming language, much commercial software is written in Java, often not because Java is designed for the kind of problem people have to solve, but because there are many programmers out there that know Java, there is a big open source community, and, finally, "nobody gets fired for choosing Java" (paraphrasing "nobody gets fired for buying IBM"). This outweighs the fact that Java is often clearly not the best choice from a pure language point of view. A software library in a "modern" programming language could attract more people (maybe more contributors if it is open source), for this reason, than a library in a much better suited, but less known, programming language.

Edit: To a certain extend the Sage project and software is what I was looking for, although it employs Python as a higher level language which I don't think is a good choice. (http://www.sagemath.org/)

-
I agree with Alex's comment that programming languages are best if they're used for what they were intended for. From your description of "modern" languages, how about Matlab? I contend that any numerical algorithm within the scope of Recipes that is not already implemented as a function/extension in Matlab or Octave can be found, with a non-restrictive license, by spending a little time with google. –  Jerry Gagelman Oct 13 '10 at 9:42
Have a look at numerics.mathdotnet.com . –  Bruce Arnold Nov 30 '10 at 23:07
Thanks, will do! –  Tim van Beek Dec 1 '10 at 10:02
Supercomputers don't have Java compilers. Of course, one may argue about whether this is a cause or an effect. –  David Ketcheson Jul 4 '12 at 18:32

## 13 Answers

Try http://numpy.scipy.org/ (assuming you're fine with an interpreted language -- otherwise, there is a python compiler, but I know very little about it).

Why languages like Java are not widespread? Well, my view on the subject is the following. You have two kinds of people who do numerical mathematics:

• Some need to squeeze every instruction of performance out of their mission-critical software, and naturally they are worried about aliasing, garbage processing, cache misses. Thus they use plain C or even Fortran, and use all the voodoo they know to make it run faster.
• Some don't, as they develop algorithms for testing/research only, or where the sheer computing time is not the bottleneck. They generally use Matlab, since it is much easier to work with, and as long as you avoid (or recode in C) tight for loops has a reasonable performance.
-
Seconded. I use numpy to do numerical functional integration. It's not the fastest thing in the world, but it's good enough for nontrivial tasks. –  userN Oct 13 '10 at 12:01
Note 1: NumPy is part of the larger SciPy library. I think SciPy is a better fit for the question than NumPy per se. –  John D. Cook Nov 30 '10 at 20:54
Note 2: NumPy and SciPy are implemented in C under the hood. You have a nice higher-level language interface in Python, but you also have the efficiency of C. –  John D. Cook Nov 30 '10 at 20:55

I think, your idea of "modern" is not very practical, when it comes to computational mathematics. To people, who do research on algorithms, the one crucial measurement of their success is running time. This is also very important to most mathematicians, because hypotheses are tested for very large numbers/objects. Your "modern" languages are mainly geared at speeding up the design and coding process and not the running time of the programs. Implementing mathematical algorithms in C# would have more of an educational purpose than a practical one, and if someone tries to sell you such a thing, you shouldn't buy it.

I don't know much about the "common knowledge" that you are referring to, or in fact about the software industry, but fact is that when people write computer algebra packages, they still sometimes implement the crucial and most used parts in Assembler! Your comparison "Java versus C++" misses the point. Already C++ produces so much more overhead than C, say, that it would rarely be the language of choice for computer algebra packages, where efficiency is crucial. Don't get me wrong, I recognise the huge advantages of object oriented programming languages, but they are at their best if they are used for the things they were design for.

-
Thanks for this answer! Of course, in the software industry, it is most of the time much cheaper to by better hardware than to pay people that optimize the running time of the software. I'm interested what mathematicians working in computational numerical analysis think about this, so your answer helpful. The problem with these "performance prejudice" that I referred to in my question is of course that no one specifies the kind of computations they refer to. Assembler? Really? –  Tim van Beek Oct 13 '10 at 10:52
@Tim: Example, consider basic functions such as sin, cos, log, exp, pow --- I would be surprised if they were not implemented in Assembler. –  Suvrit Oct 13 '10 at 12:37
Tim, I should have said that I am not working in numerical analysis. I regularly use computers to test hypotheses in number theory or group theory and the programs are usually much quicker to write than they will run. It is not uncommon to let a program run over the weekend or over a week, just to try lots of examples (and it often pays off). I don't know, whether numerical analysis people do that sort of thing, but I imagine so. –  Alex B. Oct 13 '10 at 13:49
@Alex: thanks for pointing out your research focus :-). Actually I'm interested in all kinds of uses of computer programs in mathematics including number theory, not only numerical analysis (but I thought it would be helpful to focus on numerical analysis for this question). –  Tim van Beek Oct 13 '10 at 15:17
Alex, you write "Already C++ produces so much more overhead than C [...]". This is not true in general, it clearly depends on how you write your programs. You can use a proper subset of C++ to write programs which are as fast as C programs (or in fact, generate the same machine code), but which are much more readable and both easier to write and maintain; for example, by using classes (but not inheritance!) and templates. (Of course, you shouldn't accidently pass objects by value.) Anyway, for me, C++ is the language of choice if I want to implement something which should better be efficient. –  felix Oct 13 '10 at 22:20

This is still something of an open problem. The reason is that achieving high performance on numerical codes takes fairly close attention to things like data locality and the details of the cache hierarchy -- failing to do so can easily cost a factor of 10 or 20 in performance. So compiling high-performance numerical code from high-level programs requires fairly sophisticated language and compiler technology.

• The state of the art at the moment for high-level numerical programming is probably something like SAC, which is a purely functional language which uses ideas from functional programming, linear types, and nested data parallelism to efficiently parallelize and SIMD-ize high level specifications of array processing programs. This builds on ideas from older research languages like Sisal and NESL.

• Another example is the ATS project, which is a functional language which uses dependent type theory and linear types to support generating very efficient code. I don't think they have explicitly focused on numerical codes, though they have investigated high-performance systems programming more generally.

• Another commonly-used approach is to write a high-level program which generates the low-level code, rather than writing it by hand. This is used by the FFTW package ("Fastest Fourier Transform in the West"), which is an Ocaml program (written in a very higher-order functional style) which generates the C code in the FFTW package.

-
Thanks! Never heard of those...F# is sometimes marketed as exactly the kind of thing as SAC, as you describe it, which is why I wanted to know that the professionals think about it :-) –  Tim van Beek Oct 13 '10 at 10:56
F# is a very nice language. However, it plays in the same space as Java or C#: it offers excellent bindings to platform linear algebra packages (such as BLAS) and the base language is fast enough that statistical programs that have significant non-matrix components will run acceptably fast. Actually writing the matrix code itself in it won't work well. (Compare this with R, Matlab, or Numpy, where the non-matrix part of the language is hopelessly slow. This sucks for (eg) robust statistics, where algorithms often do many small matrix ops, rather than a few big ones.) –  Neel Krishnaswami Oct 15 '10 at 7:58

A number of numerical libraries in Java can be found at JavaNumerics

"Modern" in this context means to me: object oriented (not C or Fortran)

Fortran has object-oriented features since 2003, see, for instance Fortran working group note 5 (IMHO object-orientation has been oversold, in general, and does not belong in Fortran, specifically, but the Fortran standards committee didn't ask me, alas)

Commercial libraries that have been ported to Fortran 2003 include IMSL and NAG. Opensource library LAPACK is F90, but if you are going to do numerical work, chances are, you'll have to use it.

A software library in a "modern" programming language could attract more people (maybe more contributors if it is open source), for this reason, than a library in a much better suited, but less known, programming language.

It is much harder to write robust numerical code, than to contribute to an average opensource project. I believe, that the choice of language is a tertiary concern, - understanding the mathematical subject matter and error anlysis takes more work than learning another (procedural) programming language.

-
Thanks! I learned about the object oriented features of FORTRAN several weeks after asking the question, but haven't looked into it yet. I agree with your overall assesment concerning open cource projects versus numerical code, but I'm nevertheless concerned about the growing language chasm of numerical code projects and others. –  Tim van Beek Dec 1 '10 at 10:01
There's an article "Matrix Multiply: A Case Study" which is highly illuminating wrt cost one may pay for using popular prescriptions, such as OO or immutability. They end up with 4 orders of magnitude speedup obtained by "playing dirty". –  Mikhail Kagalenko Dec 1 '10 at 12:33

An advantage of sage is that it wraps high performance libraries (asm/C/C++) in python.

I don't mean to start a language flamewar, but when I run some computation for more than a day, I definitely pick underlying C/C++ or something professionally coded in assembler and I suppose this is common practice.

-
Sure, this question is not about language wars, that is, I'm not interested in one either. However, I am still interested if there are software projects under way that would profit from a higher level of abstraction compared to C. Especially "glue code" or "housekeeping code" (io, testing, logging, exception handling etc.) is much easier to write and maintain in higher level languages. This becomes interesting for big and long living software projects like Sage (or, much more, for global climate models), compared to the "hack'n forget an algorithm for a single paper" - kind of coding. –  Tim van Beek Nov 7 '11 at 10:34
Weather models are coded in C++ or Matlab. Check out eigen.tuxfamily.org/index.php?title=Main_Page for a use of C++ template metaprogramming to speed things up. –  Vít Tuček Jul 14 '12 at 10:19

It is not entirely clear to me what are you looking for and why.

If you want the best performance you are almost surely bound to use compiled language such as Fortran1, C++ or C. Of course, you can always use almost any "higher level language" such as Python, Ruby or whatnot to glue together routines from libraries written in some "low level language" as C, Fortran, etc. Octave, Matlab and Sage come to mind. NumPy is quite good example of this approach, since almost all of its core functions are written in C (i.e. LAPACK).

If, on the other hand, you want to experiment with algorithms themselves, you can implement them in Haskell or (oca)ML or some other "mathematicians friendly" language. Also, succinct syntax and lack of side effects means it's much easier to prove correctness. Moreover there are area specific systems/languages as LiE, Macaulay2, Singular, GAP, ...

Writing basic routines (such as those from Numerical Recipes) in any other language than C, C++ or Fortran means that

• you are doing an exercise
• you are trying to solve a problem in a language where you can't use functions from libraries written in C/C++/Fortran
• you need numerical routines in a big (i.e. not feasible to rewrite in C/C++/Fortran) project where calling an external library causes unwanted overhead (this can be the case for example with Java).

1 Please note that modern Fortran is quite high level language with functional and object-oriented features.

edit:

I am not a Fortran programmer myself, so what follows are information relayed from a friend of mine (physicist and active Fortran programmer). I realize that this is not what you are asking for, but it may be helpful nevertheless.

books:

When starting, it's best to learn Fortran 95/90 first. For a gentle introduction see http://www.cs.mtu.edu/~shene/COURSES/cs201/NOTES/fortran.html

For an overview of the language see http://en.wikipedia.org/wiki/Fortran_language_features (mostly written by Michael Metcalf; see the book above).

As for OOP, it was possible in Fortran 95 using various hacks, but another revision of the language (Fortran 2003) included OOP features into the language norm itself. For technical account on Fortran 2003 features see ftp://ftp.nag.co.uk/sc22wg5/N1551-N1600/N1579.pdf For tutorial on OOP in Fortran 2003 see http://www.pgroup.com/lit/articles/insider/v3n1a3.htm and http://www.pgroup.com/lit/articles/insider/v3n2a2.htm

Fortran is still evolving, the latest standard being Fortran 2008 which introduced concurrent programming techniques:
http://www.training.prace-ri.eu/uploads/tx_pracetmo/coarrayvideo1.pdf
http://www.training.prace-ri.eu/uploads/tx_pracetmo/coarrayvideo2.pdf
http://www.training.prace-ri.eu/uploads/tx_pracetmo/coarrayvideo3.pdf

I should add that there probably isn't a compiler that supports all the features of the standard (2003).

-
Thanks, I infer from your answer that while it may not be "entirely clear" what I am looking for, it is clear enough to provide useful information :-) Could you expand your footnote about Fortran a bit? (Literature, development environments, online references, experiences?) –  Tim van Beek Jul 16 '12 at 9:29

There are several good answers already, but none of them so far deal with high performance computing, which nowadays means computing on tens or hundreds of thousands of cores. I am aware of only three Python codes that have scaled to this level:

Of course, you can't write a performant numerical code entirely in Python, so typically these use a compiled language for performance-critical parts (in PyClaw, we use both Fortran and C). One reason for choosing Python is that automated tools like f2py and Cython make it easy to incorporate compiled code.

I do not believe there are any high performance computing codes written in the languages you mention. I suspect modern supercomputers don't support those languages.

-
Thanks! Yes, I wondered if anyone would bother to write virtual machines (.NET or Java) for high performance architectures. –  Tim van Beek Jul 5 '12 at 7:29

You might also check out the GNU Scientific Library (GSL) at http://www.gnu.org/software/gsl/, for which there are bindings and wrappers for other languages. Or you could forgo OOP and write your programs in C, which does at least have the advantage of being extremely fast and efficient. And I believe you can compile GSL so it is optimized for your system. I believe there is (or was) a Haskell wrapper for GSL, but I can't find it.

Most high level systems (Matlab, Octave, Sage etc) allow you to write routines in C and link them in, thus getting the best of both worlds: the speed of C and the ease of use of a higher level computing environment. I guess it all comes down to what sort of algorithms you are using: big slow ones (such as Monte Carlo routines with 10^20 random numbers) or something quick and straightforward, such as quadrature.

Over the years I've used Maple, Matlab, MuPAD, Sage, Octave, C, Java, and others, and I've written (or at least blogged) about them all. I can't prescribe one over the other - you have to decide for yourself which best suits your usage, expertise, and requirements.

-
This is the only answer so far that makes the main point that needed to be made. It's basically irrelevant whether the language in which the library is implemented has fancy, stylish features. If the user has a favorite language and wants to use it, all that matters is whether the library has an API for that language. Similarly, I use linux, and the fact that the linux kernel is implemented in plain C has zero effect on how I use the OS. –  Ben Crowell Jun 5 at 18:51

I think jHepWork is very attractive and modern math software (written in Java and uses Jython as a macro language)

-
Thanks for your anwer. –  Tim van Beek Nov 2 '11 at 7:47

In my opinion, high performance computing focuses more on the efficiency of computing, especially in the area of scientific computing. So the C with high performance is undisputed a good choice.

Moreover, in parallel computing which is a branch of high performance computing, whose enviroment is Linux terminal, programmers should deal with much more details of the data transfering and allocation by using MPI library which is more convenient implemented by C or Fortran.

From my view point, especially in parallel programming, high performance computing is still at the status of developing and exploring. And a lot of details need to be dealt with including the standard, mechanism,abstract model of the parallel programming.(I think maybe the MPI is not the appropriate approach to describing parallel programming ultimately but is a good way nowadays.)

When people make them clear, I think the modern programming languages as you described will be adopted generally in high performance computing.

-
Thanks, but the focus of my question is more on languages and environments that are intuitive for mathematicians who intend to implement algorithms - languages that provide a better infrastructure than, for example, C. "Infrastructure" is meant is a couple of different ways here, both "code for housekeeping, GUI, networking", level of abstraction (e.g. pointer handling, OOP), easy of build and deployment etc. My question is less about what someone interested in high performance computing needs to focus on. –  Tim van Beek Jul 16 '12 at 9:34

I recommend the Enthought Python Distribution, available here. From the EPD homepage:

The Enthought Python Distribution provides scientists with a comprehensive set of tools to perform rigorous data analysis and visualization. Python, distinguished by its flexibility, coherence, and ease-of-use, is rapidly becoming the programming language of choice for researchers worldwide. EPD extends this capacity with a powerful collection of Python libraries to enable interactive technical computing and cross-platform rapid application development.

There are free academic versions (Windows, Mac OS, Linux), very easy to install.

DISCLAIMER: I am not affiliated directly or indirectly with Enthought.

-

Another contender that has appeared very recently is Julia. It is a Matlab-like language designed from scratch for scientific computing by compiler and programming language experts, with a special eye to speed potential, parallelizability and syntax consistency. The language is JIT-compiled, and seems blazing fast on simple CPU-intensive benchmarks (several factors faster than Python/Matlab/R, within a small factor from Fortran and C).

It is garbage-collected; its object-oriented features are based on multiple dispatch and supports lambdas and some functional programming,

Library support is still not as large as the contenders, but I'd say it already includes everything that you can study in a first course on numerical analysis.

-

You should use code from Numerical Recipes. If you do this, then when people ask you for your code you can tell them that you aren't allowed to give it to them because it would violate the licence. Then people will stop bugging you and you won't have to reveal how messy your code is.

-
Sarcasm is cheap, but I share your frustration. This is not the place to indulge in it, though. –  Tim van Beek Jul 16 '12 at 9:25