Is it common not to use libraries for standard numerical algorithms, and why?

Question

A lot of numerical algorithms (integration, differentiation, interpolation, special functions, etc.) are available in scientific computation libraries like GSL. But I often see code with "hand-rolled" implementations of these functions. For small programs which are not necessarily intended for public distribution, is it common practice among computational scientists to just implement numerical algorithms yourself (by which I mean copying or transcribing from a website, Numerical Recipes, or similar) when you need them? If so, is there a particular reason to avoid linking to something like GSL, or is it just more "tradition" than anything else?

I ask because I'm a big fan of code reuse, which would suggest that I should try to use existing implementations when possible. But I'm curious whether there are reasons that principle is less valuable in scientific computation than in general programming.

---

Forgot to mention: I'm specifically asking about C and C++, as opposed to languages like Python where there is a clear benefit (speed of execution) to using a library.

Answer 1

I used to implement everything myself, but lately have begun using libraries much more. I think there are several very important advantages of using a library, beyond just the issue of whether you have to write a routine yourself or not. If you use a library, you get

• Code that has been tested by hundreds/thousands/more users
• Code that will continue to be updated and improved in the future, without any work on your part
• Optimized code that is more efficient and perhaps more scalable than what you would write in a first attempt
• Depending on the library, by establishing an interface to it in your code you may get access to many algorithms that you currently don't use but may want to in the future

In the last bullet point above, I'm thinking of large libraries like Trilinos or PETSc. I can reinforce this with a couple of concrete personal examples in development of PyClaw. Although it would have been straightforward to parallelize Clawpack with MPI calls, we chose to use PETSc. This allowed us to limit the paralle code in the package to less than 300 lines of Python, but even better, by putting our data in PETSc's format we gained immediate access to PETSc's implicit solvers, enabling current work on an implicit solver in PyClaw. As a second example, PyClaw initially included hand-code fifth-order WENO reconstruction, but we eventually decided to rely on the PyWENO package for this. This was a huge gain, since PyWENO can automatically generate WENO routines of any order in several languages.

Finally, if you use libraries, you can contribute back by developing improvements or finding bugs, which will benefit many other people, whereas debugging or improving your own code only benefits you.

Answer 2

There is substantial programmer overhead involved in linking to a library function, especially if that library is new to the programmer. It is often simpler to just rewrite simple algorithms rather than figure out the specifics of a particular library. As the algorithms become more complex this behavior switches.

Python has excelled at reducing this overhead with tools like pip/easy_install and a uniform data structure interface (i.e. every library seems to take and produce a numpy array).

Answer 3

One of the projects I'm involved in right now is writing a flexible simulation and analysis package for a class of particle physics detectors. One of the goals of this project is to provide the code base to be used in these things for decades to come.

At this point point we already have two dozen dependencies, making the build process such a nightmare that it has spun off a separate project managed out of the Fermilab computing center just to provide a reliable tool-chain.

Now imagine that you encounter a need for some tool that isn't in that tool-chain (happened to me just last month). You have three choices

1. Roll you own. With all the risks and hassles that involves.
2. Scrape some code out of a library somewhere and include it directing in The Project. Meaning that you have to take on maintenance going forward, and that you'll have to understand someone else's code when that happens.
3. Go to the people who maintain the tool-chain, beg them for what you need and then wait for a release cycle to get it. These guys are pretty responsive, but you have to make the case for it without working code or after you've done (1) or (2).

It's very easy to choose (1). Maybe too easy.

Answer 4

I think it is quite common, with some algorithms more likely to be re-implemented than others.

There's a tricky trade-off between how annoying a library is to install, how hard it is to implement the algorithm yourself, how hard it is to optimize it, and how well the library fits your needs. Also, sometimes using a library is just overkill: I used the slow bisection algorithm in one of my programs because I called it only a few times and I didn't want to add a library just for that.

Is it easy for you to write a well-optimized version? If it is, you might be better off doing it. You're going to get exactly what you need and you won't depend on anyone's work. But of course you really need to know what you're doing: even simple algorithms can be tricky to implement.

I would be curious to see a study on this, but from my biased perspective, scientists often use libraries for linear algebra and random number generators, with most of the remaining code being homemade.

Answer 5

I think that implementing an algorithm instead of using a library can sometimes give a better understanding and control of the model. When I am coding some program for scientific computations, it's important for me to understand what I am doing. Implementing the important algorithms helps me to get a better knowledge of the problem and achieve better control of it.

On the other hand, sometimes it's not a trivial task to select a library that is needed for getting a solution, so it's better to search for already-implemented algorithms when you are sure what are you trying to achieve and why do you want it.

If the algorithms are complex, then coding them by hand gives you the opportunity to improve performance/quality of solution using task-specific features. And sometimes it is necessary to change the algorithm a bit, which is easier if you know the code that you wrote and you can edit it in the way you want.

Answer 6

One answer is that there are so many slight variations to numerical code that it is really hard to encapsulate that in a library. Take this in comparison to web software, which is often easy to install and has a clear set of inputs and outputs. I think more common is people grabbing a framework, or big library that acts like a framework (Trilinos/PETSc), and using that ecosystem to get the benefits of using community codes.

Answer 7

Before deciding whether or not to use libraries, I think you'd also want to figure out how much the use of a library will help your code. If you're going to be using a well-optimized library for a key computational kernel, then it's probably a lot more efficient than trying to write your own.

However, if you're writing a specialized routine that's only going to be called once during a program's execution, it may not be worth it to adapt your code to fit the framework required by a library.

(Another thing to think about: how much re-architecting will you need to do to take advantage of the library? Unless the man-hours you spend to fix the code are compensated for by corresponding gains in efficiency or numerical accuracy, it may not be worth it in the long run. Ideally, though, this is something you plan for when initially designing data structures and algorithms, so that the library's "flow" is taken into account from the ground up.)

Answer 8

My 2 cents.

I think it is easier to write generally about this, rather than just about C/C++. First, libraries in languages like Python are not necessarily used to get a speed benefit, even if that is a consequence. I think @David covered the reasons pretty well.

Taking it from the top, the language implementation to some extent dictates what libraries you have access to. Commonly used languages in computational science include C, C++, Python, Perl, Java, Fortran, and R. Less common examples might be Ocaml and Common Lisp. Now, since most of these languages are written in C, they have a natural Foreign function interface to C. However, it is not so easy to call, say, a Perl library from Python or vice versa. So in practice people tend to either

1. Use a library written in their implementation language, usally something that is part of the standard libraries, or otherwise widely available, or

2. Call a C/C++ library through the languages FFI. This assumes that a wrapper does not already exist, since if it does, it is not easily distinguishable from (1).

(2) is usually harder, because you have to wrap the C/C++ function yourself. Also, you have to either bundle the library, or add an extra dependency. For that reason, people are more likely to use the builtin language libraries rather than use GSL for example, which is in C.

For very generic routines, say generating random samples from distributions, or basic numerical routines like quadrature of integrals, it is easy and common to reuse some library. As the functionality one is trying to implement becomes more complex, it becomes exponentially more unlikely that one is going to find the exact function one wants in another library, and even one does, one could spend lots of time searching and finally adapting the function as necessary (the code style/design could be a problem for example). And as discussed above, one has access to only a subset of the libraries out there. On the other hand, implementing an algorithm oneself if it is complex and not the main focus can be daunting, and of course one has to deal with those pesky speed issues.

So, this becomes an optimization problem in cost/benefit analysis. My experience is that even for comparatively standard techniques like MCMC, I usually wind up writing my own code, because it fits better with how I am designing the overall software.

Of course, even if you end up not using the code, it is possible to learn from other people's code. I don't know how often scientists actually bother to do this, though. My impression is that reading other people's code to learn is more a software engineer thing.

Answer 9

Thinking back to my second-year mechanics course, it occurs to me that part of reason I have implemented my own versions of well-known algorithms is that I was taught to do it that way. I can't think of a single example where I was taught how to interface to and link in a library in my undergrad physics education. I do have a fond memory of seeing for the first time a list of coordinates of a spinning golf ball, having computed the solution of the coupled Newton equations in FORTRAN myself. There is a certain thrill and satisfaction (even pride) that comes from calculating things from scratch.

Answer 10

I think one should use tested libraries as much as possible. Most people are not experts in numerical computing and probably will not be able to implement a solution as correctly and carefully as what's available in well-tested libraries. That said, however, it is sometimes the case that there are no available libraries that implement the combination of capabilities needed in a given application. I've seen this happen in the technical area in which I work; existing codes didn't solve all the cases, and someone ended up implementing a solver from scratch that did.

Answer 11

The fundamental problem is often with the interface between the application and the library. An application programmer has knowledge about the problem that is often lost when passing the problem (for example as a matrix) to a library. This knowledge is such that exploiting it more than offsets the benefits of using the highly optimized library. As a result, the application programmer "rolls" his/her own implementation that does exploit the knowledge.

Thus, a really good library needs for such knowledge to be passed from the application into the library, so that the library, too, can take advantage of it.

Answer 12

In addition to all the things already said above, I will repeat my answer from the "Fortran vs C++" question: The most valuable asset a programmer has is her time. Yes, external dependencies are often awkward. But spending time to re-implement, debug and test algorithms that others have already implemented is almost always stupid, and the result will also rarely be as good as code that was written by experts on a particular topic. Re-use what others have done!

Answer 13

The group I work with uses libraries as much as possible. I'm one of the few programmers, and the rest of the folks picked up programming on the job. They know enough of their own limitations to know where they shouldn't be dabbling. IMSL is the preferred library. Stuff like GSL would be forbidden due to licensing restrictions, even though this is a federal agency and we give our software away anyway.

Answer 14

"Re-use is primarily a social phenomenon. I can use someone else's software provided that

1. it works
2. it is comprehensible
3. it can co-exist
4. it is supported (or I'm willing to support it myself, mostly I'm not)
5. it is economical
6. I can find it.

" — B. Stroustrup, The C++ Programming Language 2 ed. (1991) p. 383.

Answer 15

Several good reasons have been given by others for using libraries, and also for rolling your own routines.

You can sometimes speed up computations for a specific application because you know in advance that you will never need the wide range of values the library routine covers, or the accuracy those routines deliver.

Of course, a lot depends on the particular application and how many times the library routine will be called. Why would you call a library routine for Bessel functions billions of times if you only need a few significant figures for a small range of x, and some simpler technique will be sufficient for your needs?

Answer 16

Is little to add, we have to reuse code, it is about code sustainability and contribution to society, but that is all above.

The reason why we do not reuse code is that if you are beginning programmer is hard to understand others code. It is particularly difficult with advanced C++, and you can do some tricks in pure C as well.

Very often at the beginning, one understand the method, but not as it is implemented in the library, or simply how to use the library with its generic interface, error control and conventions, very often documented for experienced programmers if at all. This gives the illusion that is better to implement a standard method, like LU factorization by yourself. Moreover, new programmers underestimate the value of code testing, validation and portability for different operating system. So at the end reason is laziness, writing own code looks like as a faster and easier solution.

The reality is that we can learn more by using and reading code than programming by yourself from scratch.

Laziness drives me for most of the time, I think also the majority of people. For the same reason, some write code from scratch and other use existing libraries.

Answer 17

Libraries algorithms provide in contrast to own implementations:

• They are generic and templated. You can later reparameterize your implementation without worrying changing your own code which is supposed to have a lot of constraints.
• There is savety versus degenerate cases of input data. A lot of computational geometry algorithms, e.g. convex hull ones, need to handle for example colinearity of three points. You might be able to ignore those cases if you never plan to distribute your code and also do not want to reuse your code in future often.
• They provide the minimal runtime complexity for either expected or worst case input configurations. Higher level algorithms have as building bricks often lower level algorithms,e.g. sorting algorithms, or special data types. Quick sort might be the most common choice for sorting data, but if your implementation of the algorithm has to garantee n(log(n)) you cannot use it.
• They are memory usage efficient
• They are further runtime optmized
• If supported, by far more closed to be "bug" free in general, especially if you work with the main branch. Nothing is more tested than a well distributed library. Not every bug crashes, not every bug produces unreasonable results. Your algorithm's implementation might still produce acceptable results, just not as good as it is designed for. The less visible a bug is, the less likely it is that you as a single person can even detect it.

I still consider it good when entering a new field to implement one version of a well understandable algorithm on your own. I takes a lot of time in total. I bought and read books, the named Press et al. I always read much theory before and during those implementations. And after understanding the general concepts of a field and experiencing the traps in practice for me it is time to jump to the in all aspects better library implementations. I think you will become a better user of the library if you wrote an "hello world" algorithm in the libraries field on your own.

If you work in a bigger team it might not be your own choice wether your team uses a specific library or not. The core team might do the decision. And there might be a person responsible for the library binding in your project with its own time plannings. Rewriting one algorithm you can do with your own time planning, not relying on other peoples decision.

If you are on your own and like to distribute there is another problem. I consider as well as many other source code as most useful resource. Many to all informaticians might agree here. In an applied field outside informatics it might be needed to provide a precompiled exectuable on windows. Under Linux you can set up things relatively easy on your own in case of open source usage libraries.

Rewriting an algorithm on your own gives lisence freedom. Your project might not support GPL lisence of GSL for example.

The lisence might be the only constraint which is independent from the reseachers point of view.