Short Codes Can Be Long on Insight
By Nicholas Higham
In 2013, MIT Press published a remarkable book entitled 10 PRINT CHR$(205.5+RND(1)); : GOTO 10 [8]. The title is a program created in 1982 in Commodore 64 Basic that prints a random maze to the screen. The book, written by a team of 10 authors with backgrounds in digital media, art, literature, and computer science, analyzes this program from every conceivable angle. It touches on the nature of patterns, repetition, and randomness; structured programming (including the polemic against the goto statement by Edsger Dijkstra); the Basic programming language; the history of computer art; and many other topics. For more on the program and a MATLAB version, see my blog post.
Cartoon created by mathematician John de Pillis.
The maze oneliner stands in stark contrast to much of today’s software, which can be thousands of lines long, multilayered, and far removed from the basic algorithm being employed. The classic \(Ax=b\) problem of \(n\) linear equations in \(n\) unknowns illustrates this point perfectly. In one sense, the problem is trivial: one can quickly code up the three nested loops of Gaussian elimination, and go on to implement partial pivoting with just a little more effort. Yet a modern \(Ax=b\) solver must be vastly more complicated than the mathematics suggests if it is to exploit modern computer architectures. Depending on the solver, this might mean employing a panel factorization with an appropriate blocking and use of level 3 BLAS operations, or using a tilebased factorization with dynamic task scheduling [1]. And exploiting sparsity introduces further complexities. Davis et al. [5] estimate that behind the MATLAB backslash with a sparse matrix lie about 120,000 lines of code implementing sparse direct methods.
Until the 1970s, it was possible to write practical codes for solving the standard numerical analysis problems within two to three pages; this is what George Forsythe and Cleve Moler did in their book Computer Solution of Linear Algebraic Systems (1967) and later, with Michael Malcolm, in Computer Methods for Mathematical Computations (1977). Moler has said, “One of the biggest reasons these two books were as successful as they were was because the programs in them were not only useful and correct, they were short and readable.”
In the 1990s, there was an effort to provide templates: descriptions of general algorithms that could be customized by the user. SIAM published two books of templates: one on eigenvalue problems [2] and the other on iterative methods for linear equations [3], the latter being accompanied by barebones Fortran and MATLAB codes. Both books have been very successful.
More recently, a number of SIAM book authors, especially those writing for the Fundamentals of Algorithms series, have built short, readable codes into their books. Nick Trefethen’s Spectral Methods in MATLAB (2000) contains 40 short programs, and he notes that “you can do an astonishing amount of serious computing in a few inches of computer code!” In a similar vein, his “Ten Digit Algorithms” essay [9] presents algorithms with three constraints: “The program can be at most one page long, and it has to solve your problem to at least ten digits of accuracy on your machine in less than five seconds.”
Providing a short code to help the reader understand the essence of an algorithm in research articles can also be beneficial. Indeed, a good test for an author of a survey paper is to produce simple implementations of the algorithms that are treated. My colleagues Stefan Güttel and Françoise Tisseur have recently written a survey of nonlinear eigenvalue problems that includes 12 very short MATLAB programs implementing the main methods discussed [6]. The reader can gain much insight by downloading and playing with the codes. It is also worth noting that some of SIAM Review’s most highly cited articles, such as “An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations” [7], include short codes.
What is my favorite short program? I offer funm_randomized in Figure 1. This MATLAB function computes \(f(A)\) for a square matrix \(A\) and a given scalarvalued function \(f\). The builtin MATLAB function funm does the same thing, but is a few hundred lines long and has some restrictions on \(f\). How can funm_randomized be so short? Adding a random perturbation to \(A\) produces a matrix that can be “safely” diagonalized. The error in the computed \(f(A)\) is typically at or below the level of the square root of the unit roundoff, so the accuracy will generally be quite a bit worse than for funm. Brian Davies proposed and analyzed this “approximate diagonalization” method [4]. Open questions about its behavior remain, and I hope it will be further studied.
Figure 1. MATLAB function funm_randomized.m. This code can be downloaded from https://gist.github.com/higham.
The 35year old “10 print” maze program may seem irrelevant, written in a moribund language. But according to the TIOBE Index for March 2017, Basic is still the sixth most popular programming language. As stateoftheart numerical software grows in length, in response to more complicated computer architectures, “10 print” serves as a reminder that we need to keep providing simple demonstrators that make the key ideas underlying our algorithms and software accessible. Not every algorithm can be as concisely described as funm_randomized.m, but simplifying down to the core of an idea is a great exercise to improve one’s understanding.
http://blogs.mathworks.com/cleve/2013/01/07/georgeforsythe/
References
[1] Abdelfattah, A., Anzt, H., Dongarra, J., Gates, M., Haidar, A., Kurzak, J.,…YarKhan, A. (2016). Linear Algebra Software for LargeScale Accelerated Multicore Computing. Acta Numerica, 25, 1160.
[2] Bai, Z., Demmel, J., Dongarra, J., Ruhe, A., & Van der Vorst, H. (Eds.). (2000). Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide. Philadelphia, PA: Society for Industrial and Applied Mathematics.
[3] Barrett, R., Berry, M., Chan, T., Demmel, J., Donato, J., Dongarra, J.,…Van der Vorst, H. (1994). Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Philadelphia, PA: Society for Industrial and Applied Mathematics.
[4] Davies, E.B. (2007). Approximate diagonalization. SIAM J. Matrix Anal. Appl., 29(4), 10511064.
[5] Davis, T., Rajamanickam, S., & SidLakhdar, W. (2016). A Survey of Direct Methods for Sparse Linear Systems. Acta Numerica, 25, 383566.
[6] Güttel, S., & Tisseur, F. (2017). The nonlinear eigenvalue problem. To appear in Acta Numerica. MIMS EPrint 2017.7, Manchester Institute for Mathematical Sciences, The University of Manchester.
[7] Higham, D.J. (2001). An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations. SIAM Review, 43, 525546.
[8] Montfort, N., Baudoin, P., Bell, J., Bogost, I., Douglass, J., Marino, M.,…Vawter, N. (2012). 10 PRINT CHR$(205.5+RND(1)); : GOTO 10. Cambridge, MA: MIT Press. Freely available as a PDF from http://nickm.com/trope_tank/10_PRINT_121114.pdf.
[9] Trefethen, L.N. (2005). Ten Digit Algorithms (Report Number 05/13). Oxford University Computing Laboratory. Retrieved from https://people.maths.ox.ac.uk/trefethen/tda.html.

Nicholas Higham is the Richardson Professor of Applied Mathematics at the University of Manchester. He is the current president of SIAM. 