I got the idea for these test matrices from Gene Golub years ago, but the mathematical foundation comes from a paper by Divakar Viswanath and Nick Trefethen.... read more >>

]]>I got the idea for these test matrices from Gene Golub years ago, but the mathematical foundation comes from a paper by Divakar Viswanath and Nick Trefethen.

This post is about two MATLAB programs included in the collection from *Numerical Computing with MATLAB*. The program `golub` generates badly conditioned integer test matrices. The program `lugui` is an interactive graphical interface that allows you to experiment with pivot strategies in Gaussian elimination.

Set the random number generator to the state it has when MATLAB starts and generate a 6-by-6 example.

```
rng('default')
A = golub(6)
```

A = 1 3 11 0 -11 -15 18 55 209 15 -198 -277 -23 -33 144 532 259 82 9 55 405 437 -100 -285 3 -4 -111 -180 39 219 -13 -9 202 346 401 253

It's not obvious that this matrix is badly conditioned, but it is.

cond(A)

ans = 3.0115e+12

In case you needed another demonstration that the determinant does not indicate near singularity.

```
format short
det(A)
```

ans = 1.0000

The elements of the inverse vary over nine orders of magnitude and show perverse scaling that is a direct consequence of the way this matrix was created.

format short e X = inv(A)

X = 2.8508e+09 -1.5513e+08 3.1110e+06 1.4854e+06 -1.6559e+05 -2.0380e+04 -1.3362e+09 7.2711e+07 -1.4582e+06 -6.9624e+05 7.7615e+04 9.5520e+03 1.0422e+08 -5.6713e+06 1.1373e+05 5.4306e+04 -6.0540e+03 -7.4500e+02 1.2586e+07 -6.8489e+05 1.3735e+04 6.5580e+03 -7.3100e+02 -9.0000e+01 -8.4225e+05 4.5833e+04 -9.1900e+02 -4.3900e+02 4.9000e+01 6.0000e+00 -1.3830e+05 7.5260e+03 -1.5100e+02 -7.2000e+01 8.0000e+00 1.0000e+00

In a 1998 paper Divakar Viswanath and Nick Trefethen demonstrate that the 2-norm condition number of an *n* -by- *n* triangular matrix with normally distributed entries below the diagonal and ones on the diagonal grows exponentially almost surely at a rate *c^n* with *c* = 1.30568...

I want my matrices to have integer entries so I use randomly distributed integers, but the exponential behavior must be nearly the same.

Here is the code for `golub`.

```
type golub
```

function A = golub(n) %GOLUB Badly conditioned integer test matrices. % GOLUB(n) is the product of two random integer n-by-n matrices, % one of them unit lower triangular and one unit upper triangular. % LU factorization without pivoting fails to reveal that such % matrices are badly conditioned. % See also LUGUI. % Copyright 2014 Cleve Moler % Copyright 2014 The MathWorks, Inc. s = 10; L = tril(round(s*randn(n)),-1)+eye(n); U = triu(round(s*randn(n)),1)+eye(n); A = L*U;

`Lugui` lets you use the mouse to select the pivot elements for Gaussian elimination. Or, you can choose diagonal pivoting, partial pivoting or complete pivoting. Here are animated gifs of these three choices.

Diagonal pivoting actually defies the conventional wisdom and picks the *smallest* possible element at each stage as the pivot, namely 1.0. But this choice means that division by the pivot involves no roundoff. Moreover the choice leads to the reconstruction of the original triangular factors.

```
[L,U] = lugui_noshow(A,'diagonal')
cond_L = cond(L)
cond_U = cond(U)
```

L = 1 0 0 0 0 0 18 1 0 0 0 0 -23 36 1 0 0 0 9 28 -2 1 0 0 3 -13 -1 7 1 0 -13 30 15 16 -8 1 U = 1 3 11 0 -11 -15 0 1 11 15 0 -7 0 0 1 -8 6 -11 0 0 0 1 11 24 0 0 0 0 1 -6 0 0 0 0 0 1 cond_L = 8.4806e+06 cond_U = 7.9889e+05

Partial pivoting gives only a weak indication that this matrix is badly conditioned.

```
[L,U,p] = lugui_noshow(A,'partial');
condition_estimate_partial = max(abs(diag(U)))/min(abs(diag(U)))
```

condition_estimate_partial = 5.8208e+06

Complete pivoting does a much better job.

```
[L,U,p,q] = lugui_noshow(A,'complete');
condition_estimate_complete = max(abs(diag(U)))/min(abs(diag(U)))
```

condition_estimate_complete = 1.5166e+12

The QR decomposition needs column pivoting to reveal the ill conditioning.

[Q,R] = qr(A); condition_estimate_without = max(abs(diag(R)))/min(abs(diag(R))) [Q,R,E] = qr(A); condition_estimate_with = max(abs(diag(R)))/min(abs(diag(R)))

condition_estimate_without = 6.6062e+06 condition_estimate_with = 2.2595e+12

Divakar Viswanath and Nick Trefethen, "Condition Numbers of Random Triangular Matrices," *SIAM J. Matrix Anal. Appl.* 19, 564-581, 1998, <http://epubs.siam.org/doi/pdf/10.1137/S0895479896312869>. Also available at https://people.maths.ox.ac.uk/trefethen/publication/PDF/1998_77.pdf

Cleve Moler, *Numerical Computing with MATLAB*, 2004. <http://www.mathworks.com/moler/ncmfilelist.html>.

Get
the MATLAB code

Published with MATLAB® R2015a

A rectangular box, such as a book or a cell phone, thrown in the air can tumble stably about its longest axis, or about its shortest axis, but not about its middle axis.... read more >>

]]>A rectangular box, such as a book or a cell phone, thrown in the air can tumble stably about its longest axis, or about its shortest axis, but not about its middle axis.

In his latest book, *Differential Equations and Linear Algebra*, my colleague Gilbert Strang describes an interesting model problem made famous among MIT students by another MIT professor, Alar Toomre. If you throw a rectangular box in the air with a twist, you can make it tumble stably about its longest axis, or about its shortest axis. But, if the lengths of the three sides of the box are different, you cannot make it tumble about its middle-sized axis.

To describe the dynamics of the box, discount the effects of gravity by using a coordinate system centered in the box and that moves with it. Let $x(t)$, $y(t)$, and $z(t)$ be the angular momenta about the three principal axes, and let $I_1$, $I_2$, and $I_3$ be the moments of inertia about these axes. Then Euler's equations are

$$ { d \over dt} \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{r} ayz \\ bxz\\ cxy \end{array} \right) $$

where

$$ a = 1/I_3 - 1/I_2, \ \ b = 1/I_1 - 1/I_3, \ \ c = 1/I_2 - 1/I_1 $$

If we have

$$ I_1 = 1, \ \ I_2 = 1/2, \ \ I_3 = 1/3 $$

the equations become simply

$$ { d \over dt} \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{r} yz \\ -2xz\\ xy \end{array} \right) $$

There are six critical points. Any solution to the Euler equations that starts exactly at one of these points remains there.

$$ X = \pm \left( \begin{array}{c} 1 \\ 0 \\ 0 \end{array} \right), \ \ Y = \pm \left( \begin{array}{c} 0 \\ 1 \\ 0 \end{array} \right), \ \ Z = \pm \left( \begin{array}{c} 0 \\ 0 \\ 1 \end{array} \right) $$.

It turns out that any solution that starts out with $x^2+y^2+z^2=1$ retains that property, so solutions travel on the surface of the unit sphere. If you think of the sphere as the earth with $+Z$ at the North Pole, then $+X$ is the point where the Greenwich meridian crosses the equator. This is in the eastern Atlantic, off the coast of West Africa. $+Y$ is the point where the $90^\circ$ meridian east crosses the equator. This is in Indian Ocean, west of Sumatra.

To see what happens near a critical point, we need the Jacobian.

$$ J = \left[ \begin{array}{rrr} 0 & z & y \\ -2z & 0 & -2x \\ y & x & 0 \end{array} \right] $$

At $+X$

$$ J = \left[ \begin{array}{rrr} 0 & 0 & 0 \\ 0 & 0 & -2 \\ 0 & 1 & 0 \end{array} \right] $$

The eigenvalues of $J$ are pure imaginary, namely $\lambda = \pm \sqrt{2}i$, along with $\lambda = 0$. So near $+X$ the solutions behave locally like $\cos{(\sqrt{2}t)}$ and $\sin{(\sqrt{2}t)}$. This critical point, which corresponds to rotation about the long axis, is a stable center. The critical point $-X$, which corresponds to rotation about the long axis in the opposite direction, is also a stable center.

At $+Z$

$$ J = \left[ \begin{array}{rrr} 0 & 1 & 0 \\ -2 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right] $$

We have the same eigenvalues and so rotation about the short axis is also a stable center.

The middle axis, $+Y$, is a different story. Here

$$ J = \left[ \begin{array}{rrr} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{array} \right] $$

The eigenvalues of $J$ are real, namely $\lambda = 1, 0,$ and $-1$. Near $+Y$ and $-Y$ the solutions have components that behave locally like $e^t$. These are unstable saddle points.

The surface of the sphere can be divided into four quadrants, with one of the stable critical points, $\pm X$ or $\pm Z$, at the center of each quadrant. Any solution to the differential equation is periodic and describes an orbit that remains within one of the quadrants, circulating around its center.

If the initial point is near the center of a quadrant, the orbit is nearly elliptical and the orbit has period close to $\sqrt{2}\pi$. If the initial point is near the boundary of a quadrant, the orbit sticks near that boundary, making sharp turns at the two saddle points $\pm Y$. The period increases as the initial point gets further from the center.

I think of these as giant atmospheric currents that circulate around one-fourth of a planet.

Here is a link to a program, `tumbling_box.m`, that uses `ode45` to solve these equations and produce the following graphic. I invite you to download the program and experiment yourself. You are presented with a sphere that is empty except for dots at three critical points. You can click repeatedly to provide the initial conditions for the orbits.

You can also click on icons in the figure window toolbar to zoom in and out, and rotate and pan over the sphere. If you zoom in to the stable centers points, you can see the elliptical orbits associated with stability. If you zoom in to the saddle points, you can see the hyperbolic orbits associated with instability. If you rotate the sphere, you can follow one of the unstable orbits around the sphere, back to its starting point.

`Tumbling_box` includes `ode45` event handling code to determine the length of a period of the solution. It also has mouse button down code to interpret the coordinates of a mouse click as a point on the surface of the unit sphere.

Gilbert Strang, *Differential Equations and Linear Algebra*, Wellesley-Cambridge Press, 2014, 502pp. <http://www.wellesleycambridge.com>

Get
the MATLAB code

Published with MATLAB® R2015a

An interactive graphical experiment lets you discover the value of one of the most important numerical quantities in mathematics.... read more >>

]]>An interactive graphical experiment lets you discover the value of one of the most important numerical quantities in mathematics.

One of my favorite graphical experiments allows you to find the numerical value of $e$. The program is described in Chapter 8 of Experiments with MATLAB and is included in the toolbox available with that book. It was originally called `expgui`; today it is known as `expshow`. Here is a link to the code. I hope you are able to start up MATLAB and actually run `expshow` while you are reading this blog.

Let's say you've forgotten how to find the derivative of $a^t$ with respect to $t$. Be careful. Don't blindly follow your memory that the derivative of $t^n$ is $n t^{n-1}$ and claim that the derivative of $a^t$ is $t a^{t-1}$.

For graphical purposes, we can use an approximate derivative, replacing the tangent by a secant with a step size of $10^{-4}$. Here is the core of `expshow` and the resulting screen shot with $a = 2$

```
t = 0:1/64:2;
h = .0001;
% Compute y = a^t and its approximate derivative.
y = a.^t;
yp = (a.^(t+h) - a.^t)/h;
```

The blue line is the graph of $2^t$. At $t = 1$ it passes through $y = 2$ and at $t = 2$ it hits $y = 4$. The sienna line is the graph of the derivative of $2^t$. The important facts are that it has the same shape as the blue line and lies entirely below it.

Now take your mouse and move the blue line. It's very satisfying to be actually interact with this experiment. Push the blue line until the sienna line moves above it, to somewhere around $a = 3$.

With $a = 3$, at $t = 1$ the graph passes through $y = 3$ and at $t = 2$ it hits $y = 9$. The sienna line is now the graph of the derivative of $3^t$. The important facts are that it still has the same shape as the blue line and now lies entirely above it.

Now you know what to do.

Using the mouse, move the blue line until the two lines lie on top of each other. You have found the only function in the world that is equal to its own derivative and, in the process, discovered that, to three decimal places, the crucial value of the base is 2.718. And you did this without touching the keyboard or typing in any numbers.

Get
the MATLAB code

Published with MATLAB® R2015a

I have just returned from two meetings in Europe, the 26th Biennial Conference on Numerical Analysis at the University of Strathclyde in Glasgow, Scotland, and Sparse Days III in Saint-Girons, France.... read more >>

]]>I have just returned from two meetings in Europe, the 26th Biennial Conference on Numerical Analysis at the University of Strathclyde in Glasgow, Scotland, and Sparse Days III in Saint-Girons, France.

The Biennial Conference on Numerical Analysis in Scotland has a long history. The first two meetings were held at the University of St. Andrews in 1965 and 1967. The meetings moved to the University of Dundee in 1969. They were held there every two years, under the leadership of Ron Mitchell, until his retirement in 2007. They moved to the University of Strathclyde in 2009, under the leadership of Alison Ramage. This meeting was the 26th and the 50th anniversary. Although this is a premier conference series in numerical analysis, it was the first time I participated.

The University of Strathclyde is a public research university located in Glasgow, Scotland. The university's name comes from the "Valley of the River Clyde". The school is Scotland's foremost institution of science and engineering.

The conference was held June 23rd through 26th. I was one of a dozen invited speakers -- five from the US, four from continental Europe, and three from the UK.

Almost all of the participants gave half-hour talks. There were 168 talks organized in seven parallel sessions. Many of the talks were in one of twelve minisymposia. Sample minisymposia titles include "Stable and accurate discretisations for convection-dominated problems", "City analytics", "Chebfun: new developments, cool applications and on the horizon", and "Numerical linear algebra for optimsation and data assimilation".

Two of the invited lectures are given special emphasis to honor British numerical analysts. The A. R. Mitchell Lecture honors the University of Dundee's Ron Mitchell who was the dominant force in this conference for forty years.

This year the lecture was given by Mike Giles of the University of Oxford on "Multilevel Monte Carlo Methods". I knew nothing about this subject before his talk and I learned a great deal. His abstract provides a link to a web page featuring a survey paper and MATLAB codes, <http://people.maths.ox.ac.uk/gilesm/acta>

The Fletcher-Powell Lecture honors two mathematicians known for their work in optimization algorithms, including the Davidon-Fletcher-Powell, DFP, formula. Mike Powell, from Cambridge University, had passed away in April. Roger Fletcher, from the University of Dundee, attended the lecture.

This year the lecture was given by my long-time friend Mike Saunders from Stanford. He talked about "Experiments with linear and nonlinear optimization using Quad precision." He has used the GFortran compiler to build a version of his MINOS optimization software with the REAL*16 floating point datatype. He tackled flux balance analysis models of metabolic networks that involve coefficients ranging over 15 or 16 orders of magnitude. Ordinary double precision, that is Fortran's REAL*8, cannot do the job.

The history of the Sparse Days in St Girons conference is itself sparse. There have been just two previous conferences, one in 1994 and one in 2003. The conference is part of a program on "High Performance Linear and Nonlinear Methods for Large Scale Applications" organized by CIMI, the French Centre International de Mathematique et d'Informatique in Toulouse.

Organizers of Sparse Days in St Girons III included my buddies Iain Duff, who splits his time between Rutherford Appleton Laboratory in the UK and the Parallel Algorithms Group at CERFACS in Toulouse; Jack Dongarra, from the University of Tennessee; and Jim Demmel, from the University of California, Berkeley.

St Girons is a picturesque town in the French Pyrenees, not far from the border with Spain. The Tour de France sometimes passes through here. We met in the town's only cinema, which is air conditioned. This turned out to be fortunate planning because a heat wave hit that week, June 28th through July 2nd, with temperatures in the mid 30s C, which is mid 90s F.

The format for this conference was quite different from the one the previous week in Scotland. There were only about half as many attendees, around 100. There were no invited talks and no parallel sessions. Anyone who wanted to give a talk gave one. There were 42 talks, most of them half an hour.

Tim Davis, MathWorks consultant on sparse matrices, who recently moved to Texas A&M, talked about "Sparse SVD, and a GPU-accelerated sparse QR".

Bora Ucar of Ecoles Normales Superieures in Lyon talked about "Two approximation algorithms for bipartite matching on multicore architectures." The analysis of his algorithm happened to involve my old friend, the Lambert W function.

Joost Rommes, from Mentor Graphics in the Netherlands, talked about "Challenges in numerical simulation of electrical networks."

Alex Pothen, from Purdue, talked about "The Virtual Scalpel: real-time matrix computations and finite elements for surgery."

Dan Sorensen, from Rice, talked about "A DEIM induced CUR factorization." This title means that C is some columns of a matrix A, R is some rows of A, and DEIM is an algorithm that finds U so that the product CUR is a good low approximation to A.

Martin Gander, from Universite de Geneve, talked about "Five decades of time parallel integration." I was particularly interested in his talk because I knew something about the subject many years ago, but was not familiar with recent developments. At first, it might seem impossible to parallelize a computation that involves stepping along in time. But Martin gave an overview of four different approaches to how this can be done.

There is another Cleve in this business. Cleve Ashcraft of LSTC, Livermore Software Technology Corporation, is actually a "grandstudent". His PhD advisor at Yale was Stan Eisenstat, who was my PhD advisee when I was a visiting professor at Stanford. At Sparse Days, Cleve talked about "Separability, partitions and coverings.". His talk didn't mention the applications he is deeply involved in at LSTC, which are the developers, among other things, of the crash analysis software LS-DYNA. If you want to see a hint of some serious use of dynamic finite element calculations, take at the LSTC web pages.

This is the group photo from Sparse Days in St Girons III. There are two Cleves. Cleve Ashcraft is the guy with the green shirt and grey beard in the front row. Can you find the other one? (Thanks to Pierre-Henri Cros for the photo, and for handling local arrangements.)

Get
the MATLAB code

Published with MATLAB® R2015a

Augustin (Austin) Dubrulle deserves to be better known in the numerical linear algebra community. His version of the implicit QR algorithm for computing the eigenvalues of a symmetric tridiagonal matrix that was published in a half-page paper in *Numerische Mathematik* in 1970 is faster than Wilkinson's version published earlier. It is still a core algorithm in MATLAB today.... read more >>

Augustin (Austin) Dubrulle deserves to be better known in the numerical linear algebra community. His version of the implicit QR algorithm for computing the eigenvalues of a symmetric tridiagonal matrix that was published in a half-page paper in *Numerische Mathematik* in 1970 is faster than Wilkinson's version published earlier. It is still a core algorithm in MATLAB today.

"QR" and "QL" are right- and left-handed, or forward and backward, versions of the same algorithm. We are used to thinking of factoring a matrix into an orthogonal factor, Q, and an upper or right triangular factor, R. This leads to QR algorithms. But for reasons having to do with starting a loop at $1$ rather than $n-1$, the authors of the Handbook Series on Linear Algebra decided to use left triangular and QL algorithms.

Augustin (Austin) Dubrulle is the only guy that I know who improved on an algorithm of Jim Wilkinson. Austin was born in France. He began his career in the early 1960's with IBM at their *Institut de Calcul Scientifique* in Paris. In 1968 he transferred to an IBM center in Houston, and began work on SSP, IBM's Scientific Subroutine Package.

After studying papers and books by Wilkinson and Parlett, Austin derived his own version of the implicit QR algorithm for computing the eigenvalues of a symmetric tridiagonal matrix. He was careful about the use of common subexpressions to save operations. This is especially important when computing eigenvectors because the transformations are applied to an accompanying full matrix.

Austin told me that he wrote up his result in the style of the Linear Algebra series in *Numerische Mathematik* and asked IBM for permission to submit the paper for publication. He was told by company lawyers that, in view of the Justice Department antitrust suit, IBM could not give away software, and that the ALGOL included in the paper was software. Austin tried to make the case that ALGOL was primarily a means of human communication and secondarily a programming language of little practical use, but that argument did not fly.

When Martin and Wilkinson published their implicit algorithm, Austin was encouraged to realize that his version had fewer operations in the inner loop and so would be faster. He arranged to meet Wilkinson at the 1969 University of Michigan numerical analysis summer school. Jim encouraged Austin to submit a note with just the algorithm, no provocative software, to the series.

Here is a recreation of Dubrulle's half-page 1970 paper in *Numerische Mathematik*. This is the entire paper. You can get the original here or purchase it here.

(You can get the Martin and Wilkinson contribution here or purchase it here.)

'********************************************************************************************'

The following is a formulation of the *QL* algorithm with implicit shift which requires fewer operations than the explicit and implicit algorithms described in [1] and [2].

Let $d_i^{(s)}$ $(i=1,...,n)$ and $e_i^{(s)}$ $(i=1,...,n-1)$ be the diagonal and codiagonal elements of the matrix at the $s$ -th iteration and $k_s$ the shift. Then the $(s+1)$ -st iteration can be expressed as follows.

*Acknowledgements*. The author wishes to thank Dr. J.H. Wilkinson for his helpful comments

References

1. Bowdler, H., Martin, R.S., Reinsch, C., Wilkinson, J.H.: The QR and QL algorithms for symmetric matrices. Numerische Mathematik 11, 293-306 (1968).

2. Martin, R. S., Wilkinson, J.H.: The implicit QL algorithm. Numerische Mathematik 12, 377-383 (1968).

A. Dubrulle IBM Corporation Industry Development-Scientific Applications 6900 Fannin Houston, Texas 77025, U.S.A.

'********************************************************************************************'

Dubrulle's algorithm still needs a little work. The computation of the next $e_{i+1}$ as the square root of the sum of squares is dangerous. There is potential for unnecessary underflow or overflow. I discussed this in a post almost two years ago on Pythagorean Addition. I don't suggest actually using the *pythag* algorithm. Instead, do something like this.

if abs(p) >= abs(q) c = q/p; t = sqrt(1+c^2); e(i+1) = t*p; s = 1/t; c = c*s; else s = p/q; t = sqrt(1+s^2); e(i+1) = t*q; c = 1/t; s = s*c; end

When the papers from *Numerische Mathematik* were collected to form the 1971 *Handbook for Automatic Computation, Volume II, Linear Algebra*, almost all of them where reprinted without change. But, despite what its footnote says, Contribution II/4, *The Implicit QL Algorithm*, never appeared in the journal. The paper is the merger of Dubrulle's note and the Martin-Wilkinson contribution that it references. The ALGOL procedures, *imtql1* and *imtql2*, are new.

The authors of Contribution II/4 are A. Dubrulle, R.S.Martin, and J.H.Wilkinson, although the three of them never worked together. Proper credit is given, but I'm afraid that an interesting little bit of history has been lost.

Dubrulle's version of the implicit tridiagonal QR algorithm continues to be important today. We had it in EISPACK. Look at the source for IMTQL1.F. You will see code that is very similar to Austin's note. Except there is a call to a function `PYTHAG` to compute `E(I+1)`.

In LAPACK the *imtql1* and *imtql2* functions are combined into one subroutine named DSTEQR.F Here again you will see code that is very similar to Austin's note. Now `PYTHAG` is named `DLATRG`.

I haven't run any timing experiments myself recently, but I have been told that the method of choice for the real symmetric matrix eigenvalue problem is Cuppen's divide and conquer algorithm, as perfected and implemented in LAPACK's DSTEDC.F. A large problem is recursively broken into smaller ones. When the matrices become small enough, and eigenvectors are desired, DSTEQR is called. This is Dubrulle's version of implicit QR. Then the recursion is unwound to produce the final result.

I could use an appropriate graphic for this post. Here is a picture of the tridiagonal QR algorithm in action, generated by `eigsvdgui` from *Numerical Computing with MATLAB*.

Get
the MATLAB code

Published with MATLAB® R2015a

The 1960 programming language ALGOL 60 influenced the design of many subsequent languages, including C and a miniature language called PL/0. I based the design of the first MATLAB on PL/0.... read more >>

]]>The 1960 programming language ALGOL 60 influenced the design of many subsequent languages, including C and a miniature language called PL/0. I based the design of the first MATLAB on PL/0.

*Algol* is a bright, eclipsing three-star system in the constellation *Perseus*. It provides the name of a family of "algorithmic" programming languages developed very early in the computer era by an international group of, primarily, mathematicians.

I do not have room here to describe ALGOL. Let me just say that ALGOL has profoundly affected my professional life.

In the first published report in 1958, the programming language was originally called IAL, for International Algebraic Language. There were eight authors of the report, including John Backus, an American who headed the IBM team that developed the original Fortran; Friedrich Bauer, a German who was the subject of my previous blog post; and Heinz Rutishauser, a Swiss numerical analyst who invented the LR algorithm.

We had two of versions of ALGOL 58 at Stanford when I was a grad student there in the early '60s. Burroughs had their BALGOL on the B220. We used that for the beginning programming course for about a year. Then two students, Larry Breed and Roger Moore, implemented BALGOL on the new campus IBM mainframe, the 7094. We called it SUBALGOL and used it heavily, including in the programming courses, for a couple of years.

One of my all-time favorite acronyms is JOVIAL, for Jules Own Version of International Algebraic Language. It was developed beginning in 1959 by a group headed by Jules Schwartz at SDC, System Development Corporation. JOVIAL and it successors were in use for almost fifty years by US military projects.

ALGOL, short for Algorithmic Language, was presented in a 1960 report with thirteen authors and edited by Danish computer scientist Peter Naur. The report summarized almost everything that was known about the design of programming languages at the time, which was barely ten years after the development of the first stored program computers. Wikipedia has this to say about ALGOL's influence.

It was the standard method for algorithm description used by the ACM in textbooks and academic sources for more than thirty years. In the sense that most modern languages are "algol-like", it was arguably the most successful of the four high-level programming languages with which it was roughly contemporary: Fortran, Lisp, and COBOL.

The ACM began publishing collections of scientific algorithms in 1960. The collected algorithms appeared in *Communications of the ACM* (*CACM*) until 1975 when *Transactions on Mathematic Software* (*TOMS*) began publication and took over the collection. Initially, all the contributions were required to be in ALGOL 60. Fortran was permitted after it was standardized in 1968. The last ALGOL submission was in 1977. In more recent years, even some MATLAB programs have been accepted.

In the early years of the ACM algorithms section, Stanford played an important role in the editorial process. George Forsythe, who was the founding chairman of Stanford's Computer Science Department and my Ph. D. advisor, was also President of ACM at one time and editor of the algorithms section at another.

My grad school buddy, Bill McKeeman, contributed a couple of early algorithms, in ALGOL. One is algorithm 135, a linear equation solver, "Crout with equilibration and iteration". The Crout algorithm computes the LU factorization from the inner product description of the matrix elements. Another is algorithm 145, "Adaptive numerical integration by Simpson's rule". This one uses recursion even though our BALGOL compiler at Stanford at the time could not handle recursion and Bill could not actually run his program. When we did acquire a Burroughs B5000 with a full ALGOL compiler, Bill verified that his recursive program would work within hours of the machine's installation. (Many years later, one of McKeeman's crowning achievements is the current version of the MATLAB `why` function, which also uses recursion.)

My first ever professional publication is "Remark on Certification of Matrix Inversion Procedures" in the Algorithms section of the July 1963 *CACM*. It is in response to a contribution by Peter Naur, the editor of the ALGOL 60 report. Naur, who was not a numerical analyst, had made the common mistake of comparing the computed inverse of the floating point approximation to the Hilbert matrix to the known inverse of the exact Hilbert matrix. The difficulty is that the effect of the initial floating point approximation is larger than the effect of the rounding errors in the inversion.

For an example of an ALGOL procedure, let's take a look at SOLVE, a program from *Computer Solution of Linear Algebraic Systems*, the 1967 textbook by Forsythe and myself. We have programs in this book in ALGOL, Fortran, and PL/1. This program is preceded by DECOMPOSE, a longer program that computes the *LU* decomposition. Reading these programs today, I'm surprised to see that I used a global array to store the pivot indices. These programs are in the publication format which uses bold face characters for the key words.

In the late 1960s the journal *Numerische Mathematik* ran a series of research papers about numerical linear algebra that includes ALGOL procedures. Altogether there are 29 papers by a combined total of 19 authors. Sixteen of the papers are by J. H. Wilkinson and coauthors. In 1970 all of the papers, with a few modifications, were collected in a volume edited by Wilkinson and C. Reinsch titled "Handbook for Automatic Computation, Volume 2: Linear Algebra". (I'm not sure if Volume 1 of this intended larger project was ever written or published.)

ALGOL was always intended to be a language for research and publication, not for actual execution. IBM dominated the computer market in the '60s and '70s, certainly in the US, and IBM did not support ALGOL. Moreover, in its standard form, ALGOL 60 had no provision for input and output. So it was always expected that the algorithms published by the ACM and *Numerische Mathematik* would have to be modified or translated to some other language before they could actually be used.

The NATS Project, the National Activity to Test Software, was a project at Argonne National Laboratory whose focus was the testing, certification, distribution, and support of mathematical software. Most of the investigators were Argonne staff. I visited Argonne in the summer. The Algol procedures in the *Numerische Mathematik* linear algebra series were organized in two parts. The second part dealt with matrix eigenvalue problems. We translated the procedures in part two to Fortran, carefully preserving their structure and numerical properties. We then tested and timed them at about twenty university and national laboratory sites, and made the collection publically available. "Certification" consisted of a carefully worded offer to respond to any reports of poor or incorrect performance.

The result was EISPACK, for Matrix Eigensystem Routines. We published the Users' Guide in three parts in the 1970s. *Numerische Mathematik* and the Wilkinson/Reinsch *Handbook* had been published by the German publisher Springer-Verlag and we regarded EISPACK as a translation project, so published our guides with them as well.

EISPACK was followed at Argonne by another project, LINPACK, to produce a package for solving linear equations and related problems. Together, the two packages would provide a fairly complete set of routines for computations involving dense matrices. LINPACK was not a NATS project and did not have immediate ALGOL antecedents.

Niklaus Wirth is a Swiss computer scientist who received his PhD from Berkeley in 1963, was on the Stanford faculty at the time the Computer Science department was created in 1965, and returned to the ETH in Switzerland in 1968.

While he was at Stanford, Wirth and Tony Hoare proposed some modifications and simplifications of ALGOL 60 that made it easier to compile and practical for actual execution. The IFIP working group on algorithmic languages felt that the proposal was not a sufficient advance over ALGOL 60, but Wirth implemented it nevertheless. The result became known as ALGOL W, "W" for Wirth, and was used in many programming courses for a number of years.

Wirth is probably best known for his language PASCAL, another descendant of ALGOL. Introduced in 1970, PASCAL is one of the first languages to provide structured programming and data structures. PASCAL proved to be very popular on personal computers in the 1980s.

In 1976 Wirth published a textbook, *Algorithms + Data Structures = Programs*. The last chapter is "Language Structures and Compilers". Wirth gives a formal definition of a miniature programming language that he calls PL/0. He describes how to parse PL/0 and generate code for a hypothetical machine.

PL/0 has only one data type, integer. Wirth defines PL/0 with the seven syntax diagrams for program, block, statement, condition, expression, term, and factor. Here are three of those diagrams.

These definitions are recursive. An expression consists of terms separated by additive operators. A term consists of factors separated by multiplicative operators. And a factor consists of identifiers, numbers, or expressions in parentheses.

The syntax diagram for statement introduces **begin**, **end**, **if**, **then**, **while**, **do**, and **call**. The diagram for block introduces **procedure**. And the diagram for condition involves relational operators. So, despite the simplicity of the definitions, there is a rich structure.

At the time, I was at professor at the University of New Mexico, teaching courses in linear algebra and numerical analysis. I wanted my students to have access to LINPACK and EISPACK without writing Fortran programs. Almost as a hobby, I used Wirth's chapter to develop a simple Matrix Laboratory. Consequently, the syntax of the original MATLAB is that of PL/0. I will discuss this in more detail in my next blog post.

Get
the MATLAB code

Published with MATLAB® R2015a

Fritz Bauer, eminent German computer scientist and last surviving member of the organizing committee of the 1964 Gatlingburg Conference on Numerical Algebra, passed away on March 26 at the age of 90.... read more >>

]]>Fritz Bauer, eminent German computer scientist and last surviving member of the organizing committee of the 1964 Gatlingburg Conference on Numerical Algebra, passed away on March 26 at the age of 90.

Execute these commands in any MATLAB since version 4.0, 1992.

load gatlin image(X) colormap(map) axis off

This is the organizing committee of the 1964 Conference on Numerical Algebra in Gatlinburg, Tennessee. They are J. H. Wilkinson, Wallace Givens, George Forsythe, Alston Householder, Peter Henrici, and Friedrich Bauer. Bauer, on the far right, was the youngest member of the committee.

I obtained the 8-by-10, black-and-white, glossy photo when I attended the conference as a grad student. It was my first big professional meeting and the first of what has become an important series of meetings for me and MATLAB. I kept the photo in my files until 1992 when we could handle images in MATLAB. I scanned in the photo and it became one of our first example images.

All six of the men in the photo have made contributions that have had a lasting impact on numerical analysis, computer science, and, ultimately, MATLAB. And all six have influenced my life personally.

Credit: <http://upload.wikimedia.org/wikipedia/commons/c/c7/FriedrichLudwigBauer.jpg>

Friedrich Bauer, 1924-2015.

My friend Walter Gander has written an obituary of Bauer for the Communications of the ACM that I reference below. Bauer was born on June 10, 1924, in Regensburg, Germany. After World War II, he studied mathematics and physics at Ludwig Maximillians Universitat in Munich, where he received his Ph.D. in 1951.

Bauer was involved in the early development of computers in Germany. These included STANISLAUS, a relay based computer, in 1951, and PERM, a stored program electronic computer, in 1952-56.

Bauer held a professorship in mathematics at Mainz from 1958 until 1963. He then returned to Munich and the Munich Institute of Technology, where he spent the remainder of his professional career. He advised 39 PhD students before he retired in 1989.

Bauer, together with his colleague Klaus Samelson, invented the stack for use in parsing and evaluating algebraic expressions. Today this is also known as a LIFO, for Last In First Out, data structure.

Bauer was an early advocate for the recognition of *computer science* as an independent discipline in German universities. He also advocated the notion of *software engineering* and, in 1972, suggested a definition.

Establishment and use of sound engineering principles to obtain, economically, software that is reliable and works on real machines efficiently.

This definition of software engineering is now universally quoted.

Bauer was one of the principal authors of the reports on the programming languages International Algebraic Language, IAL, also known as Algol 58, and Algorithmic Language 1960, Algol 60. Wilkinson's research on algorithms for matrix eigenvalues was published in *Numerische Mathematik* in Algol. Equally important, Algol led directly to Niklaus Wirth's pedagogical programming language PL/0, which led to the design of MATLAB. So Bauer had a strong influence on the design of the MATLAB language. I want to tell that story in my next blog post.

The numerical linear algebra community represented by our Gatlinburg photo was very closely knit. Bauer visited Oak Ridge several times. He visited Stanford for a quarter in 1967, stayed in Gene Golub's home, and gave a course where he lectured on the theory of norms. Bauer and Householder were close friends. In fact, Householder married Bauer's wife's older sister, Heidi.

One of my favorite results from the numerical analysis that underlies matrix computation is the *Bauer-Fike Theorem*. It tells how close a computed eigenvalue is to an exact one. You need to be able to estimate the condition of the eigenvectors.

**Bauer-Fike Theorem**. Let $A$ be a diagonalizable matrix and let $V$ be its matrix of eigenvectors. Let $\mu$ and $x$ be an *approximate* eigenvalue and eigenvector with corresponding *residual*

$$ r = A x - \mu x $$

Then there exists $\lambda$, an exact eigenvalue of $A$, such that

$$ | \lambda - \mu | \le \kappa(V) \frac{||r||}{||x||} $$

where

$$ \kappa(V) = ||V|| ||V^{-1}|| $$

is the *condition number* of the eigenvector matrix.

There is a proof in Wikipedia.

I didn't know Fritz Bauer well. I only saw him for a few days every three years at the Gatlinburg/Householder meetings. But the memories are vivid.

In 1996, Householder XIII was in Pontresina, Switzerland. The traditional Wednesday afternoon "break" consisted of an excursion to the Morteratsch Glacier. The adventurous among the world's leading numerical analysts took off on a hike down the glacier, back to the base. I did not want to be left out, although the hiking boots I had hauled to Europe were not in good shape. Halfway down the glacier, a few of us were falling behind. Here comes Fritz and a couple of others. They had taken an longer, more difficult route to inspect a waterfall. He was over seventy years old at the time, but he looked great. He is an avid hiker and was in terrific shape. He was wearing the traditional Alpine lederhosen and first-rate hiking boots.

That was almost twenty years ago, but that's how I'll remember him. Fritz Bauer -- Computer Science Renaissance Man and Bavarian Alpinist.

Walter Gander, "The Life of a Computer Pioneer", Blog@CACM, April 13, 2015, <http://cacm.acm.org/blogs/blog-cacm/185577-the-life-of-a-computer-pioneer/fulltext>

Wikipedia, "Bauer-Fike Theorem", <http://en.wikipedia.org/wiki/Bauer%E2%80%93Fike_theorem>

Get
the MATLAB code

Published with MATLAB® R2015a

This is the third in a multi-part series on the MATLAB random number generators. MATLAB has used variants of George Marsaglia's ziggurat algorithm to generate normally distributed random numbers for almost twenty years.... read more >>

]]>This is the third in a multi-part series on the MATLAB random number generators. MATLAB has used variants of George Marsaglia's ziggurat algorithm to generate normally distributed random numbers for almost twenty years.

It's important to have a memorable name for an algorithm. Ziggurats are ancient Mesopotamian terraced temple mounds that, mathematically, are two-dimensional step functions. A one-dimensional ziggurat underlies Marsaglia's algorithm. I'm not sure if we would still be using the algorithm today if Marsaglia had called it the "step function algorithm".

The probability density function, or *pdf*, of the normal distribution is the bell-shaped curve

$$ f(x) = c e^{-x^2/2} $$

where $c = 1/(2\pi)^{1/2}$ is a normalizing constant that we can ignore. If we were to generate random points $(x,y)$, uniformly distributed in the plane, and reject any of them that do not fall under this curve, the remaining $x$'s form our desired normal distribution. The ziggurat algorithm covers the area under the pdf by a slightly larger area with $n$ sections. The following figure has $n = 6$; the actual code we use today has $n = 256$. The choice of $n$ affects the speed, but not the accuracy of the code.

z = zigplot(6)

z = 0 0.8288 1.1713 1.4696 1.7819 2.1761

The top $n-1$ sections of the ziggurat are rectangles. The bottom section is a rectangle together with an infinite tail under the graph of $f(x)$. The right-hand edges of the rectangles are at the points $z_k$ shown with circles in the figure. With $f(z_1) = 1$ and $f(z_{n+1}) = 0$, the height of the $k$ th section is $f(z_k) - f(z_{k+1})$. The key idea is to choose the $z_k$'s so that all $n$ sections, including the unbounded one on the bottom, have the same area.

There are other algorithms that approximate the area under the pdf with rectangles. The distinguishing features of Marsaglia's algorithm are the facts that the rectangles are horizontal and have equal areas.

For a specified number, $n$, of sections, it is possible to solve a transcendental equation to find $z_n$, the point where the infinite tail meets the last rectangular section. In our picture with $n = 6$, it turns out that $z_n = 2.18$. In an actual code with $n = 256$, $z_n = 3.6542$. Once $z_n$ is known, it is easy to compute the common area of the sections and the other right-hand endpoints $z_k$. It is also possible to compute $\sigma_k = z_{k-1}/z_k$, which is the fraction of each section that lies underneath the section above it. Let's call these fractional sections the *core* of the ziggurat. The right-hand edge of the core is the dotted line in our figure. The remaining portions of the rectangles, to the right of the dotted lines in the area covered by the graph of $f(x)$, are the *tips*. The computation of the $z_k$ 's and $\sigma_k$ 's is done only once and the values included in the header of the source code.

With this initialization, normally distributed random numbers can be computed very quickly. The key portion of the code computes a single random integer, $j$, between $1$ and $n$ and a single uniformly distributed random number, $u$, between $-1$ and $1$. A check is then made to see if $u$ falls in the core of the $j$ th section. If it does, then we know that $u z_j$ is the $x$-coordinate of a point under the pdf and this value can be returned as one sample from the normal distribution. The code looks something like this.

j = randi(256); u = 2*rand-1; if abs(u) < sigma(j) r = u*z(j); else r = tip_computation end

Most of the $\sigma_j$'s are greater than 0.98, and the test is true over 98.5% of the time. One normal random number can usually be computed from one random integer, one random uniform, an if-test, and a multiplication. The point determined by $j$ and $u$ will fall in a tip region less than 1.5% of the time. This happens if $j = 1$ because the top section has no core, if $j$ is between $2$ and $n-1$ and the random point is in one of the tips, or if $j = n$ and the point is in the infinite tail. In these cases, the more costly tip computation is required.

It is important to realize that, even though the ziggurat step function only approximates the probability density function, the resulting distribution is exactly normal. Decreasing $n$ decreases the amount of storage required for the tables and increases the fraction of time that extra computation is required, but does not affect the accuracy. Even with $n = 6$, we would have to do the more costly corrections over 25% of the time, instead of less than 1.5%, but we would still get an exact normal distribution.

Details of the ziggurat implementation have varied over the years. The original 1984 paper by Marsaglia and Tsang included a fairly elaborate transformation algorithm for handling the tips. We used this algorithm for several releases of MATLAB and I reproduced its behavior in the program `randntx` in *Numerical Computing with MATLAB*. That method and that code are now obsolete.

The 2000 paper by Marsaglia and Tsang available at the `jstatsoft` link given below has a simpler rejection algorithm for use in the tips. We have been using this in more recent releases of MATLAB, including current ones.

Marsaglia and Tsang were anxious to make their normal generator as fast as their uniform generator. But they were a little too anxious. Their original code made one call to a 32-bit uniform generator. They used the high order 7 bits for the vertical index $j$ into the ziggurat and then reused all 32 bits to get the horizontal abscissa $u$. Later investigators, including Jurgen Doornik, found this correlation led to failures of certain statistical tests.

We now make two calls to the 32-bit Mersenne Twister generator (during sequential computation). We take 8 bits to get $j$ and then 52 of the remaining 56 bits to get a double precision $u$.

How does this affect the timing? Allocate a long vector.

clear m = 25e6 x = zeros(m,1);

m = 25000000

Generate a random uniform vector and a random normal vector. Then compare the two execution times.

tic, x = rand(m,1); tu = toc tic, x = randn(m,1); tn = toc ratio = tu/tn

tu = 0.3416 tn = 0.4520 ratio = 0.7558

So, random uniform execution time is about three-quarters of the random normal execution time.

Thanks to Peter Perkins for his help on the entire series about the MATLAB random number generators.

This post on the ziggurat is adapted from section 9.3 of *Numerical Computing with MATLAB*.

George Marsaglia and W. W. Tsang, "A fast, easily implemented method for sampling from decreasing or symmetric unimodal density functions." *SIAM Journal on Scientific and Statistical Computing* 5, 349-359, 1984. <http://epubs.siam.org/doi/pdf/10.1137/0905026>

George Marsaglia and W. W. Tsang, "The ziggurat method for generating random variables." *Journal of Statistical Software* 5, 1-7, 2000. <http://www.jstatsoft.org/v05/i08>

Jurgen A. Doornik, "An Improved Ziggurat Method to Generate Normal Random Samples." PDF, Nutfield College, Oxford, 2005. <http://www.doornik.com/research/ziggurat.pdf>

Cleve Moler, *Numerical Computing with MATLAB*, SIAM 2004, <http://dx.doi.org/10.1137/1.9780898717952> Also available at: <http://www.mathworks.com/moler/random.pdf>

Get
the MATLAB code

Published with MATLAB® R2014a

This is the second of a multi-part series about the MATLAB random number generators. If you ask for `help rng`, you will get lots of information, including the fact that there are three modern generators.... read more >>

This is the second of a multi-part series about the MATLAB random number generators. If you ask for `help rng`, you will get lots of information, including the fact that there are three modern generators.

'twister' Mersenne Twister 'combRecursive' Combined Multiple Recursive 'multFibonacci' Multiplicative Lagged Fibonacci

My previous post was about `twister`. Today's post is about the other two, `'combRecursive'` and `'multFibonacci'`, which are designed for parallel computing.

I frequently use the card game Blackjack to demonstrate parallel computing. At the same time I can demonstrate the random number generators. I regard Blackjack as a financial instrument, not unlike the stock of a publicly traded company. Simulating the size of an investment as a function of time is a typical application of the Monte Carlo technique.

Begin by opening a pool of workers.

parpool;

Starting parallel pool (parpool) using the 'local' profile ... connected to 2 workers.

Four players each play twenty thousand hands of Blackjack.

p = 4; n = 20000;

Collect the results in this array.

B = zeros(n,p);

A parallel `for` loop executes a Blackjack program that knows nothing about parallel computing.

parfor k = 1:p B(:,k) = cumsum(blackjacksim(n)); end

Plot the results.

plot(B) title(['Blackjack , $ ', int2str(sum(B(end,:)))]) xlabel('Number of hands') ylabel('Stake, ($)') axis([0 n -2500 2500])

The card dealing calls the random number generator. It is essential that the different parallel workers have different, independent streams of random numbers. The default MATLAB generator `twister` does not offer this feature. Simply starting `twister` with different seeds on different workers does not provide statistically independent streams. So we turn to the other generators. To see which one is in use here, run a small `spmd`, "single program, multiple data" block.

spmd r = rng s = r.State' format short x = rand(1,7) end

Lab 1: r = Type: 'combRecursive' Seed: 0 State: [12x1 uint32] s = Columns 1 through 6 1720035765 2052922678 1637499698 3048064580 1173461082 2391850890 Columns 7 through 12 1862757735 2368998908 1385613640 1660833332 146924518 3104031825 x = 0.8789 0.6969 0.0409 0.4609 0.7528 0.2871 0.5241 Lab 2: r = Type: 'combRecursive' Seed: 0 State: [12x1 uint32] s = Columns 1 through 6 323405913 3817048408 3712601073 1070773748 1552739185 3267875480 Columns 7 through 12 1594297407 2533167732 3377045245 3413340742 2651847732 1248925296 x = 0.1072 0.3194 0.1048 0.6623 0.0878 0.3692 0.8035

We see that we have two workers, that they are both using the `combRecursive` generator, that they have the same `seed`, but different `states`, so they are generating different random numbers.

Also known as `mrg32k3a`. A 32-bit combined multiple recursive generator (CMRG), due to Pierre L'Ecuyer, at the Universite de Montreal, and his colleagues, described in the papers referenced below. This generator is similar to the CMRG implemented in the RngStreams package. It has a period of $2^{127}$, and supports up to $2^{63}$ parallel streams, via sequence splitting, and $2^{51}$ substreams each of length $2^{76}$. Here is a link to the C source code. combmrg2.c

The state of the backbone generator is a 2-by-3 array S that evolves at each step according to the linear recurrence expressed succinctly in MATLAB by

m1 = 2^32 - 209; m2 = 2^32 - 22853; x1 = mod(1403580*S(1,2)-810728*S(1,3),m1); x2 = mod(527612*S(2,1)-1370589*S(2,3),m2); S = [x1 S(1,1) S(1,2)) x2 S(2,1) S(2,2)];

A single precision random real `u` is then produced by

z = mod(x1-x2,m1); if z > 0, u = z/(m1+1); else, u = m1/(m1+1); end

The important feature of this generator is that it is possible to create different initial states for each worker in a parallel pool so that the resulting streams of random numbers are statistically independent.

Also known as `mlfg6331_64`. A 64-bit multiplicative lagged Fibonacci generator (MLFG), developed by Michael Mascagni and Ashok Srinivasan at Florida State University. This generator, which has lags $l=63$ and $k=31$, is similar to the MLFG implemented in the SPRNG package. It has a period of approximately $2^{124}$. It supports up to $2^{61}$ parallel streams, via parameterization, and $2^{51}$ substreams each of length $2^{72}$.

The state of this generator is a length 63 vector of 64-bit integers S. The recurrence relation is

$$ x_n = x_{n-k} \times x_{n-l} (mod 2^{64}) $$

Each random double precision value is created using one 64-bit integer from the generator; the possible values are all multiples of $2^{-53}$ strictly within the interval (0,1).

Again, the important feature of this generator is that it is possible to create different initial states for each worker in a parallel pool so that the resulting streams of random numbers are statistically independent.

Which one should you use? Most of the time, stick with the default and you'll be OK. You will get `'twister'` in serial computations and `'combRecursive'` on the workers in a parallel pool. You can use

rng('shuffle')

at the beginning of a session if you want different sequences of random numbers in different sessions. Otherwise, don't worry about setting the generator or the seed.

If you want to experiment, you can use `rng` to try different generators and different starting seeds on your computation. If you find a problem where it makes a significant difference, please let us know.

Thanks again to Peter Perkins.

Pierre L'Ecuyer, R. Simard, E. J. Chen, and W. D. Kelton. "An Objected-Oriented Random-Number Package with Many Long Streams and Substreams." Operations Research, 50(6):1073-1075. 2002. <http://www.iro.umontreal.ca/~lecuyer/myftp/papers/streams00.pdf>

Pierre L'Ecuyer. "Good Parameters and Implementations for Combined Multiple Recursive Random Number Generators." Operations Research 47(1):159-164. 1999. <http://dx.doi.org/10.1287/opre.47.1.159>

Michael Mascagni and Ashok Srinivasan. "Parameterizing Parallel Multiplicative Lagged-Fibonacci Generators." Parallel Computing, 30: 899-916. 2004. <http://www.cs.fsu.edu/~asriniva/papers/mlfg.ps>

Michael Mascagni and Ashok Srinivasan. "SPRNG: A Scalable Library for Pseudorandom Number Generation." ACM Transactions on Mathematical Software, Vol 26 436-461. 2000. <http://www.cs.fsu.edu/~asriniva/papers/sprngacm.ps>

Get
the MATLAB code

Published with MATLAB® R2014b

This is the first of a multi-part series about the MATLAB random number generators.... read more >>

]]>This is the first of a multi-part series about the MATLAB random number generators.

If you issue the following commands at any point in any recent version of MATLAB, you will always get this plot.

```
rng default
hist(randn(10000,1),100)
```

The `rng` command controls the random number generator that is used by the `rand`, `randn`, and `randi` functions. When called with the `default` parameter, `rng` resets the generator to the condition that it has when a fresh MATLAB is started. To see what that condition is, just call `rng` by itself.

rng

ans = Type: 'twister' Seed: 0 State: [625x1 uint32]

We see that the default generator is `'twister'`, the default seed is 0, and that the state is a length 625 vector of unsigned 32-bit integers.

If you ask for `help rng`, you will get lots of information, including the fact that there are three modern generators.

Generator Description ------------------------------------------------------ 'twister' Mersenne Twister 'combRecursive' Combined Multiple Recursive 'multFibonacci' Multiplicative Lagged Fibonacci

The remainder of today's post is about `'twister'`. I will cover the others in future posts.

Mersenne Twister is, by far, today's most popular pseudorandom number generator. It is used by every widely distributed mathematical software package. It has been available as an option in MATLAB since it was invented and has been the default for almost a decade.

Mersenne Twister was developed by professors Makoto Matsumoto and Takuji Nishimura of Hiroshima University almost twenty years ago. Here is their home page. The C source code is available here.

Mersenne primes are primes of the form *2^p - 1* where *p* itself is prime. They are named after a French friar who studied them in the early 17th century. We learn from Wikipedia that the largest known prime number is the Mersenne prime with *p* equal to 57,885,161. The Mersenne Twister has *p* equal to 19937. This is tiny as far as Mersenne primes go, but huge as far as random number generators are concerned.

Matsumoto explains how the name "Mersenne Twister" originated in the Mersenne Twister Home Page.

MT was firstly named "Primitive Twisted Generalized Feedback Shift Register Sequence" by a historical reason. Makoto: Prof. Knuth said in his letter "the name is mouthful." Takuji: ........ a few days later Makoto: Hi, Takkun, How about "Mersenne Twister?" Since it uses Mersenne primes, and it shows that it has its ancestor Twisted GFSR. Takuji: Well. Makoto: It sounds like a jet coaster, so it sounds quite fast, easy to remember and easy to pronounce. Moreover, although it is a secret, it hides in its name the initials of the inventors. Takuji: ....... Makoto: Come on, let's go with MT! Takuji: ....well, affirmative. Later, we got a letter from Prof. Knuth saying "it sounds a nice name." :-)

The integer portion of the Mersenne twister algorithm does not involve any arithmetic in the sense of addition, subtraction, multiplication or division. All the operations are shifts, and's, or's, and xor's.

All the elements of the state, except the last, are unsigned 32 bit random integers that form a cache which is carefully generated upon startup. This generation is triggered by a `seed`, a single integer that initiates the whole process.

The last element of the state is a pointer into the cache. Each request for a random integer causes an element to be withdrawn from the cache and the pointer incremented. The element is "tempered" with additional logical operations to improve the randomness. When the pointer reaches the end of the cache, the cache is refilled with another 623 elements.

The algorithm is analyzed by investigating the group theoretic properties of the permutation and tempering operations. The parameters have been chosen so that the period is the Mersenne prime 2^19937-1. This period is much longer than any other random number generator proposed before or since and is one of the reasons for MT's popularity.

By design, the results generated satisfy an equidistribution property in a 623-dimensional cube.

Here is the function in the Mersenne Twister source code that converts a pair of random uint32s into a random double. You can see that it takes the top 27 bits of one int and the top 26 bits of the other, sticks them together, and multiplies by what MATLAB would call `eps/2`. This is the only place in the code where floating point arithmetic is involved.

double genrand_res53(void) { unsigned long a=genrand_int32()>>5, b=genrand_int32()>>6; return(a*67108864.0+b)*(1.0/9007199254740992.0); }

The result would be zero in the highly unlikely event that both `a` and `b` are zero. If that happens, the MATLAB interface rejects the result and calls this function again. So the smallest double results when `a` equals zero and `b` equals 1. The largest double results when both `a` and `b` are all 1's. Consequently, the output from `rand` is in the closed interval

$$ 2^{-53} \leq x \leq 1-2^{-53} $$

For more about random number seeds, streams, and state, see Peter Perkins, guest blogger in Loren's blog. And, of course, see the documentation.

Thanks to Peter Perkins for the work he has done on our random number suite over the years, and for enlightening me.

M. Matsumoto and T. Nishimura. "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudorandom Number Generator." ACM Transactions on Modeling and Computer Simulation, 8(1):3-30. 1998. Available online at: <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/ARTICLES/mt.pdf>.

Get
the MATLAB code

Published with MATLAB® R2014b