Is it possible for user written algorithms to outperform libraries' built-in optimized functions?

Question

I've always had this question in mind (even if it may sound vague), but in my numerical analysis courses we've always learned how to analyze and optimize code. However, since most linear algebra libraries (i.e. LAPACK, BLAS, etc.) have been continuously optimized by professional developers since the 1980s, I'm questioning what the purpose of having us write down algorithms and run them is. For example, we've had to write the QR iteration algorithm in MATLAB when there's already a built-in function (eig). Are there special cases where these built-in functions are vunerable for numerical instability where custom algorithms would outperform them?

Update : I want to thank you all for your interesting answers and comments I have read them all.

I also realized the beauty of this very narrow field of mathematics that has profound application in our work. The amount of effort that these mathematicians and computer scientist have invested in making these codes more powerful have made it beautiful for us (the ones who are learning about their theories) so that we can admire them while studying and knowing when to use these optimized built-in functions. While we write down a 15 line of codes to compute the Householder QR decomposition of some matrix $A$, I can only wonder how many lines of codes the built in function qr() has (probably in the 100s) so that it can serve us best in our computations in real life applications.

Answer 1

Although LAPACK has some incredibly optimized code, it can still be worth it to write your own version in a few cases.

1. The most important reason (and the reason they make you do it in your course). It's a great learning experience. You will never fully understand something like the QR algorithm (which is actually incredibly complicated to get the details right) unless you spend some time implementing it. A lot of the skills you pick up will translate when you write other numerical codes.

2. A lot of the code is contributed by the researchers who develop the algorithms and it takes a lot of effort. Researchers are usually evaluated on how many papers they publish, not the code they write for libraries, so there can be a lack of incentive for researchers to contribute these codes. As an example, it was only last year that an optimized version of the QZ algorithm was contributed, although a paper was published on it in 2006.

3. LAPACK code needs to be robust and general. This is mostly important for the small scale routines. I'll talk about DLARTG as an example. That routine calculates a Givens rotation. It has a lot of scaling, safety checks and exceptions to make sure it is always accurate for all the input you could ever give it. If you know that for your application, these safety checks will not be necessary or you don't need the strong accuracy guarantees, you can achieve faster results.

Answer 2

@ThiysSteel covers a lot, here is another perspective I find important:

Even if you have available excellent implementations of any algorithm you might need, you still need to understand some of the internals in order to make good use of your library. Some examples:

• Any linear-algebra library can compute an inverse matrix really fast and as accuratly as mathematics allows. But you have to know a little theory in order to know that you should still avoid using that if at all possible. (Matrix decomposition and solving without explicitly computing the inverse is better, generally speaking).
• Your library can probably compute eigenvalues of any arbitrary matrix, so it might be tempting to just use the general-purpose function and be done with it. But with some background knowledge you can appreciate how much(!) better the special-purpose algorithms for symmetric/hermitian matrices work. So it is reasonable to make a considerable effort to transform you problem to that case.
• In order to achieve optimal performance you have to think about data-layout. Row-Major? Column-Major? AoS or SoA (in case of complex numbers)? A good library will be able to work with any of these, but you need to know a little about the inner workings in order to make a informed choice (though in practice: try it out an measure performance. )

The situation is somewhat similar to writing assembly code. You should pretty much never do so by hand. But knowing (some) assembly will make you a better programmer in any high-level language.

Answer 3

To clarify @ThiysSteel's good answer:

The point is not to attempt to surpass the optimization of code written by very experienced people, who'd wrangled with it for decades.

The point is to acquaint yourself with the issue that were addressed in that code, to appreciate the "niceties", and not be naive about related matters in the future.

(It's wise/fair to learn from other people. You don't have to reinvent bad wheels yourself. "They sell them!" Both bad wheels and good... :)

Answer 4

I used to try to optimize code via using assembly language (as opposed to C). I had some clear success, where a real-time microphone array worked with assembly language, but it would not work at all in the C language.

However, since that time, more than 20 years have passed. It's very difficult to achieve this sort of improvement nowadays. Compilers have become better, and the same would be true for libraries. People still use packages based on LINPACK (e.g. linear algebra for eigenvalue/eigenvector routines), the core of which was written in Fortran in the 1960's! There have been new wrappers for C#, etc.

Any time one writes in a lower-level language, there are increased risks of errors and unreadable code, and of course a problem for anyone who has to maintain the code besides the original author.

Also, there may be special cases, edge cases, etc. in the programming problem to be solved. However, these would probably be best handled in the 'calling program' (i.e.:

if (such and such) 
     call the library 
  else 
     implement special rule for edge case 

... etc.....

For example, one could avoid divide by zero if the absolute value of a floating point number is less than some threshold, and then call a calculation involving a division.

Answer 5

This might be tangent but i feel it worth adding. A common way for user code to beat libraries is to exploit structure the libraries are not aware of. Block diagonal matrices are easier to work in structured form then dense matrix that contain the complete matrix and lots of zeros. Cyclic boundary conditions in PDEs can be solved more easily than working on the over the dense system matrix. In addition the rise of new datatypes in hardware as Float16/BFloat16 means old libraries need to catch up to support those. Similarly if you have some composed numerical type like complex numbers or dual numbers writing and optimizing your own functionality will beat emulating the math with real number matrices due to data locality.

Answer 6

You mention LAPACK and BLAS. One really sobering thing to realise is that they have been written in FORTRAN and that means code looking straightforward (and thus easy to write by people specialising in mathematics rather than programming) can be heavily optimised by the compiler.

The same does not work in C or C++. For starters, there is no standardised way to represent multidimensional matrices to subroutines (some recent C standard got variable-dimensioned arrays in, just to have the next standard declare support optional again). For another, arrays passed to FORTRAN subroutines are guaranteed to point to non-overlapping structures. That allows significant optimisations and reordering of load/store instructions. In connection with standard representation and passing conventions, strength reduction and pointer arithmetic are expected to be done by the compiler as a rule.

Recent C implementations have gotten a "restrict" keyword which can be employed to similar effect but conflicts with variable-dimension matrix syntax and leads to inscrutable code only writable and readable by experts.

C++ can get around some of those problems with wagonloads of template programming, but every linear algebra package implements its own matrix representations which are not compatible with one another.

As a result, creating good interoperable linear algebra libraries in C/C++ is still playing catch-up with FORTRAN code half a century old and many of the better C/C++ implementations of those libraries essentially call extern "FORTRAN" code routines that have been compiled using a good FORTRAN compiler.

To a good degree, that is the problem of C interchanging arrays and pointers particularly in function call elements and exposing the internals of array implementations in that manner, making it hard for the compiler to assume anything useful about storage layout and making more-dimensional arrays with sizes not fixed at compilation awkward to deal with.

So when trying to outperform "standard routines", you may not just be fighting against a lot of accumulated knowledge but you may also be hampered down by competing with code that has not even been forced through a compiler that got its hands tied behind its back due to language semantics unsuitable for numeric work.

Answer 7

The point of learning how to write "the QR iteration algorithm" is not because eig will be slower. It is because "the QR iteration algorithm" is easy.

When you learn things, you are given easy tasks.

We don't teach children how to add 1+1 because we need the answer. We don't teach them how to manually carry on paper when they add 7+17 because that is the optimal way to do it. We don't teach long division of 3 into 36 for the same reason.

We teach those because it is easy, and the skills you learn are useful for non-easy problems.

Not every problem you'll want to solve in MATLAB will have a pre-written solution for it. Some of the problems with pre-written solutions will be harder to find than (if you have some knowledge) just writing it yourself, and the version you write yourself will be good enough.

Optimization techniques often cross from one algorithm to another. So when you learn how to write an algorithm, and how to modify it to make it faster, and even learn the insane things the "professional" versions do, you gain skills that may apply the next time you need an algorithm that isn't pre-written for you.

Beyond that, it is possible that this isn't the end. You may be the next one who implements or optimizes those algorithms. You might invent a new mathematical algorithm, and have to write out a working sketch of it to see how it acts. You might leave scientific computing, and become a computer programmer (I hear they are hiring), and writing algorithms and optimizing them is a useful skill for that job.

Answer 8

In some special cases it may be if you have a special case application that's very different for where the libraries are intended.

I have written an optimization application that needs to do 80 million very simple linear programming solutions per second on 8 cores. That's 10 million per second per core. I tested a number of available linear programming libraries. Every single one of them was too slow by a factor of over 10x.

In the end I had to write a system that picks some of the inequalities, treats them like equations, solves them (I solved them algebraically to do as much work as I can in program writing stage leaving as little work as possible in execution stage), calculates the value of the optimized function, and then the same is done for a different set of picked inequalities, and the one with most optimal value wins.

If I had no knowledge about the nature of linear programming problems, namely that the optimal point has to lie on the boundary of the inequalities, this custom solver would have been impossible to write.

It seems that the available libraries were intended for cases where the linear programming problem is complex as opposed to being simple, and the sheer amount of CPU time used to solve the problem is needed to solve one complex problem rather than solving 10 million simple problems with slightly different parameters.