The Quantum Matrix Processor, being announced today, is the world's first viable quantum array processor. Basic matrix operations are done instantaneously, with infinite procession. The programming environment is classic MATLAB.

Contents

Quantum Research Partners

I am pleased to be involved in today's announcement of the QMP, the Quantum Matrix Processor. This new device has been developed by Quantum Research Partners, a startup here in Santa Fe that was formed two years ago by a team from Los Alamos National Lab. I have been a consultant from the beginning of their project because they have based their design on classic MATLAB.

QMP

The QMP is similar conceptually to a GPU, or Graphics Processing Unit. It is a coprocessor that interfaces with the motherboard via an expansion slot such as PCI Express or Accelerated Graphics Port. Matrices are transferred between the memories and registers of the CPU, or central processor, and the coprocessor. Once the data is available in the QMP, quantum computing can take place.

Qureal

Traditional quantum computing involves "qubits", where a quantum circuit represents a single bit and quantum arithmetic operations are similar to digital arithmetic operations. The primary innovative aspect of the QMP is the "qureal" where a single quantum circuit represents a real number. The five basic arithmetic operations of addition, subtraction, multiplication, division, and square root, can be carried out instantaneously, according to the rules of real arithmetic.

Quantities are stored in qureal registers and qureal memory as real numbers. There is no representation error. Arithmetic operations are carried out exactly. There is no roundoff error. Arithmetic operations are done instantaneously. If it could be measured, the execution time of an addition or multiplication operation would be found to be zero.

Uncertainty Principle

Heisenberg's uncertainty principle implies that it is not possible to access an infinitely precise real number in finite time. It would take an infinite amount of time, for example, to access the binary or decimal representation of a real result.

Programming Environment

The programming environment of the QMP is based on MATLAB. The beta release of the environment is an implementation of what is now called "classic MATLAB", the simple interactive matrix calculator that I developed before we started MathWorks. The available features include basic MATLAB matrix management and command processing. Each of the following functions has been implemented as a single QMP instruction, so it is executed essentially instantaneously on a matrix of qureals.

   ABS   ATAN  BASE  CHAR  CHOL  CHOP  COND  CONJ
   COS   DET   DIAG  DIAR  DISP  EIG   EPS   EXEC
   EXP   EYE   FLOP  HESS  HILB  IMAG  INV   KRON
   LINE  LOAD  LOG   LU    MAGI  NORM  ONES  ORTH
   PINV  PLOT  POLY  PRIN  PROD  QR    RAND  RANK
   RAT   RCON  REAL  ROOT  ROUN  RREF  SAVE  SCHU
   SIN   SIZE  SQRT  SUM   SVD   TRIL  TRIU  USER
   CLEA  ELSE  END   EXIT  FOR   HELP  IF    LONG
   RETU  SEMI  SHOR  WHAT  WHIL  WHO   WHY

LINPACK Benchmark

This processor can solve a system of simultaneous linear equations instantaneously. Does that qualify the QMP as the faster computer in the world? My colleague Jack Dongarra and the other folks running the top 500 list will have two concerns when we submit our benchmark results. First, we can't access the computed solution in zero time. And, second my program is written in MATLAB, not Fortran.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Easter Sunday is March 31 this year. Why? How is the date of Easter determined?

Contents

Easter Day

Easter Day is one of the most important events in the Christian calendar. It is also one of the most mathematically elusive. In fact, regularization of the observance of Easter was one of the primary motivations for calendar reform centuries ago.

Easter is linked to the Jewish Passover. The informal rule is that Easter Day is the first Sunday after the first full moon after the vernal equinox. But the ecclesiastical full moon and equinox involved in this rule are not always the same as the corresponding astronomical events, which, after all, depend upon the location of the observer on the earth. The date varies between March 22 and April 25.

There is an easter program in Experiments with MATLAB. Let's check this year.

datestr(easter(2013))

ans =
31-Mar-2013

How about next year?

datestr(easter(2014))

ans =
20-Apr-2014

Don Knuth

The EXM program is based on the algorithm presented in the first volume of the classic series by Donald Knuth, The Art of Computer Programming. Knuth has used it in several publications to illustrate different programming languages. The task has often been the topic of an exercise in computer programming courses.

Knuth says that the algorithm is due to the Neapolitan astronomer Aloysius Lilius and the German Jesuit mathematician Christopher Clavious in the late 16th century and that it is used by most Western churches to determine the date of Easter Sunday for any year after 1582.

Metonic cycle

The earth's orbit around the sun and the moon's orbit around the earth are not in sync. It takes the earth about 365.2425 days to orbit the sun. This is known as a tropical year. The moon's orbit around the earth is complicated, but an average orbit takes about 29.53 days. This is known as a synodic month. The fraction

year = 365.2425;
month = 29.53;
format rat
ratio = year/month

ratio =
    6444/521   

is not the ratio of small integers. However, in the 5th century BC, an astronomer from Athens named Meton observed that the ratio is very close to 235/19.

format short
ratio
meton = 235/19

ratio =
   12.3685
meton =
   12.3684

In other words, 19 tropical years is close to 235 synodic months. This Metonic cycle was the basis for the Greek calendar and is the key to the algorithm for determining Easter.

MATLAB program

Here is the complete MATLAB program.

% addpath ../../exm
type easter

function dn = easter(y)
% EASTER  Date of Easter.
% EASTER(y) is the datenum of Easter in year y.
% Ref: Donald Knuth, The Art of Computer Programming,
%      Fundamental Algorithms, pp. 155-156.

% Golden number in 19-year Metonic cycle.
g = mod(y,19) + 1;

% Century number.
c = floor(y/100) + 1;

% Corrections for leap years and moon's orbit.
x = floor(3*c/4) - 12;
z = floor((8*c+5)/25) - 5;

% Epact.
e = mod(11*g+20+z-x,30);
if (e==25 && g>11 || e==24), e = e + 1; end

% Full moon.
n = 44 - e;
if n < 21, n = n + 30; end

% Find a Sunday.
d = floor(5*y/4) - x - 10;

% Easter is a Sunday in March or April.
d = n + 7 - mod(d+n,7);
dn = datenum(y,3,d);

References

[1] Donald E. Knuth, The Art of Computer Programming, Volume 1: Fundamental Algorithms (3rd edition), pp. 155-156, Addison-Wesley, 1997, ISBN 0-201-89683-4.

[2] Wikipedia, Primary article on Easter. <http://en.wikipedia.org/wiki/Easter>

[3] Wikipedia, Computus, details on calculation of Easter. <http://en.wikipedia.org/wiki/Computus>

[4] Wikipedia, Metonic cycle. <http://en.wikipedia.org/wiki/Metonic_cycle>

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Wilkinson's polynomials are a family of polynmials with deceptively sensitive roots.

Contents

$P_n$

The Wilkinson polynomial of degree $n$ has as roots the integers from $1$ to $n$.

$$ P_n(x) = \prod_{k=1}^{n} (x-k) $$

$P_{20}$

The usual value of $n$ is 20.

format compact
n = 20
syms x
P20 = prod(x-(1:n))

n =
    20
P20 =
(x - 1)*(x - 2)*(x - 3)*(x - 4)*(x - 5)*(x - 6)*(x - 7)*(x - 8)*(x - 9)*(x - 10)*(x - 11)*(x - 12)*(x - 13)*(x - 14)*(x - 15)*(x - 16)*(x - 17)*(x - 18)*(x - 19)*(x - 20)

Expressed in this form, the polynomial is well behaved.

Coefficients

It gets interesting when $P_n(x)$ is expanded in the monomial basis, $\{x^j\}$.

P = expand(P20)

P =
x^20 - 210*x^19 + 20615*x^18 - 1256850*x^17 + 53327946*x^16 - 1672280820*x^15 + 40171771630*x^14 - 756111184500*x^13 + 11310276995381*x^12 - 135585182899530*x^11 + 1307535010540395*x^10 - 10142299865511450*x^9 + 63030812099294896*x^8 - 311333643161390640*x^7 + 1206647803780373360*x^6 - 3599979517947607200*x^5 + 8037811822645051776*x^4 - 12870931245150988800*x^3 + 13803759753640704000*x^2 - 8752948036761600000*x + 2432902008176640000

The coefficients are huge. The constant term is $20!$, and that's not even the largest coefficient.

C = flipud(coeffs(P)')

C =
                     1
                  -210
                 20615
              -1256850
              53327946
           -1672280820
           40171771630
         -756111184500
        11310276995381
      -135585182899530
      1307535010540395
    -10142299865511450
     63030812099294896
   -311333643161390640
   1206647803780373360
  -3599979517947607200
   8037811822645051776
 -12870931245150988800
  13803759753640704000
  -8752948036761600000
   2432902008176640000

Symbolic form

Let's have the Symbolic Toolbox generate LaTeX and cut and paste the result into the source for the blog.

L = latex(P);
% edit(L)

$$ x^{20} - 210\, x^{19} + 20615\, x^{18} \\ - 1256850\, x^{17} + 53327946\, x^{16} - 1672280820\, x^{15} \\ + 40171771630\, x^{14} - 756111184500\, x^{13} + 11310276995381\, x^{12} \\ - 135585182899530\, x^{11} + 1307535010540395\, x^{10} - 10142299865511450\, x^9 \\ + 63030812099294896\, x^8 - 311333643161390640\, x^7 + 1206647803780373360\, x^6 \\ - 3599979517947607200\, x^5 + 8037811822645051776\, x^4 - 12870931245150988800\, x^3 \\ + 13803759753640704000\, x^2 - 8752948036761600000\, x + 2432902008176640000 $$

With the monomial basis the Wilkinson polynomial does not look so innocent. The spacing between the roots is tiny relative to these coefficients. But the toolbox can easily find the roots $1$ to $20$ with no error.

Z = sort(solve(P))'

Z =
[ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]

Double precision

Convert the symbolic form to double precision floating point.

format long e
p = sym2poly(P)';
c = flipud(coeffs(poly2sym(p))')

c =
                     1
                  -210
                 20615
              -1256850
              53327946
           -1672280820
           40171771630
         -756111184500
        11310276995381
      -135585182899530
      1307535010540395
    -10142299865511450
     63030812099294896
   -311333643161390656
   1206647803780373248
  -3599979517947607040
   8037811822645051392
 -12870931245150988288
  13803759753640704000
  -8752948036761600000
   2432902008176640000

Five of the coefficients cannot be represented in double precision format. They have been perturbed.

delta = C-c

delta =
    0
    0
    0
    0
    0
    0
    0
    0
    0
    0
    0
    0
    0
   16
  112
 -160
  384
 -512
    0
    0
    0

Double precision roots

The roots of the polynomial with the double precision coefficients are no longer 1:20. The largest perturbations occur to the roots at 16 and 17.

z = sort(roots(p));
fmt = '%25.16f\n';
fprintf(fmt,z)

       0.9999999999997972
       2.0000000000475513
       2.9999999973842115
       4.0000000600349486
       4.9999992200080809
       6.0000064898423089
       6.9999645183068031
       8.0001214537649652
       8.9997977879058215
       9.9997884273944386
      11.0021078422257670
      11.9944660326764140
      13.0076639325972590
      13.9958171375909190
      14.9955431043450050
      16.0114122070902880
      16.9888815075391190
      18.0060272186900080
      18.9981637365860050
      20.0002393259700360

Wilkinson's perturbation

Wilkinson made a different perturbation, a deliberate roundoff error on his machine to the second coefficient, the $-210$. He changed this coefficient by $2^{-23}$ and discovered that several of the roots were driven into the complex plane.

I am not sure $^\dagger$ about the sign of Wilkinson's perturbation, so let's do both. Here is a movie, an animated GIF, of the root locus in the complex plane produced by perturbations like his. It shows the trajectories of the roots from $9$ to $20$ of

$$ P_{20}(x) - \alpha x^{19} $$

as we vary $\alpha$ over the range

$$ \alpha = \pm 2^{-k}, k = 23, ..., 36 $$ The roots $1$ to $8$ stay real for perturbations in this range. Wilkinson's result is at the end of either the red or the black trajectories.

Sensitivity

With a little calculus, we can get an analytic expression for the sensitivity of the roots with Wilkinson's perturbation. Regard each root $x$ as function of $\alpha$ and differentiate the following equation with respect to $\alpha$.

$$ P_{20}(x) - \alpha x^{19} = 0 $$

At $\alpha = 0$, we have

$$ \frac{dx}{d\alpha} = \frac{x^{19}}{P_{20}'(x)} $$

This is easy to evaluate.

pprime = sym2poly(diff(P20));
xdot = zeros(n,1);
for k = 1:n
   xdot(k) = k^19/polyval(pprime,k);
end
format short e
xdot

xdot =
  -8.2206e-18
   8.1890e-11
  -1.6338e-06
   2.1896e-03
  -6.0774e-01
   5.8248e+01
  -2.5424e+03
   5.9698e+04
  -8.3916e+05
   7.5994e+06
  -4.6378e+07
   1.9894e+08
  -6.0434e+08
   1.3336e+09
  -2.1150e+09
   2.4094e+09
  -1.9035e+09
   9.9571e+08
  -3.0901e+08
   4.3100e+07

The sensitivities vary over 27 orders of magnitude, with the largest values again at $16$ and $17$. Here is a log plot.

This is a different perturbation than the one from symbolic to double, but the qualitative effect is the same.

Note added March 10, 2013.

$^\dagger$ I am back in the office where I have access to Wilkinson's The Algebraic Eigenvalue Problem. This polynomial is discussed, among other places, on pages 417 and 418. The perturbation he makes is negative, so the coefficient of $x^{19}$ becomes $-210 - 2^{-23}$ and the resulting roots are at the ends of our red trajectories.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

I have already made several posts about gatlin, the image distributed with MATLAB of the organizing committee for the Gatlinburg III conference. From the MATLAB point of view, the first man on the left in the photo, J. H. Wilkinson, is by far the most important. From the early days of computers in the 1950s until his death in 1986, Wilkinson was the world's authority on matrix computation. His research on eigenvalue algorithms and their implementation in Algol led directly to EISPACK, the mathematical foundation for the first MATLAB.

Contents

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

Early Life

James Hardy Wilkinson was born in 1919 in Strood, England, a small town east of London, on the way towards Canterbury and Dover. The family owned a small dairy.

When Jim was only eleven years old he received a scholarship to Sir Joseph Williamson's Mathematical School in nearby Rochester. When he was barely 17, two years earlier than average, he entered Cambridge on a Trinity Major Scholarship in Mathematics. At Cambridge, his tutors included the famous mathematicians G. H. Hardy, J. E. Littlewood and A. S. Besicovitch. Besicovitch characterized Jim's performance on the Tripos exams, which are the culminatation of the undergraduate mathematics program at Cambridge, as "easily at the top of the First Class".

World War II

In the early years of World War II Wilkinson joined a team of mathematicians in an office at Cambridge. They computed trajectories of artillery shells using mechanical desk calculators. Later he transferred to a laboratory at Fort Halstead to work on the thermodynamics of explosions.

At Fort Halsead, he met Heather Ware, a mathematician from King's College. They soon married and lived happily together for over 40 years.

At the end of the war, Jim transferred to the National Physical Laboratory, where he would remain for the rest of his career. NPL is located just west of London, not far from Heathrow airport, and is like a British version of the American National Bureau of Standards.

Leslie Fox

One of Jim's colleagues at NPL, who became his lifelong friend, was Leslie Fox. Leslie was a year older than Jim, had been educated at Oxford, and returned there as a professor after eleven years at NPL. After Jim's death, Leslie wrote a 40 page biographical memoir for the British Royal Society. I am using that memoir as a reference for this blog.

I bring up Fox at this point because I was so impressed by his remark in the memoir that in their early days at NPL the two of them, together with E. T. Goodwin, solved an 18-by-18 set of simultaneous linear equations using desk computing equipment. It took them two weeks! It was one of the reasons they decided to build an automatic computer.

Pilot Ace

The famous mathematician Alan Turing was at NPL at the time and proposed the construction of a elaborate machine known as the Automatic Computing Engine, or ACE. Wilkinson became an assistant to Turing. After some disagreements with the Laboratory management, Turing departed for the University of Manchester and Jim was put in charge of the project. The team built only a portion of Turing's design, the Pilot ACE. Jim actually did some of the soldering. The Pilot ACE proved to be a useful machine on its own and was used for several years. A commercial version was produced by the English Electric Company under the name DEUCE.

Here is Jim at the console of the Pilot Ace, probably at its public debut in November 1950.

Ann Arbor

Beginning in the late 1950's, the University of Michigan offered a summer short course in computing. At first it was just one course, but it soon grew to be several courses in different aspects of computing -- programming, compilers, systems, numerical methods.

Bob Bartels was both a professor of mathematics and director of the University of Michigan computing center. Bob invited eminent numerical analysts to give the lectures in the numerical methods short course. Starting around 1958, Alston Householder and Jim Wilkinson were regular lecturers in the course and came to Ann Arbor for two weeks every summer.

In 1965, Bartels decided that running the computing center had become a full time job and that he had to give up his regular teaching duties. So the university recruited someone to replace him and I became an Assistant Professor of Mathematics at the University of Michigan in the fall of 1966. I also took over organizing the summer numerical methods short course. I was pleased that Alston and Jim continued to participate. Jim's notes for the Michigan summer course became the basis for his now classic 1965 book, The Algebraic Eigenvalue Problem.

Wilkinson and Householder knew all the restaurants in Ann Arbor, including the Pretzel Bell and the Old German, where the short course lecturers would usually order beer by the pitcher. Olga Taussky-Todd, who accompanied her husband John Todd to these dinners, had a fairly reserved manner. One day she planned to take a visiting friend to the Pretzel Bell for lunch and she asked Jim, "Is it possible to order beer in more modest quantities?" Jim replied, "Yes, but it has never been one of my ambitions."

Stanford

I took my first numerical analysis course in grad school at Stanford from George Forsythe in 1962-63. There were only a handful of students in the class. Forsythe was editor of the Prentice-Hall Series in Automatic Computation and was preparing Wilkinson's book Rounding Errors in Algebraic Processes for publication. He gave us the galley proofs to both read as course notes and check for typographical errors. Forsythe was a fastidious editor, so I suspect the class did not find any typos that he had not already caught.

Sometime after I left Stanford in 1965 Wilkinson was appointed to a part time professorship there. He visited perhaps one quarter a year. I'm not sure if he ever taught a regular course. He was a terrific lecturer for a seminar or a short course and wonderful in one-on-one discussions, but he was never anxious to be a full time academic.

Handbook

Much of Wilkinson's research was published in a series of papers in the journal Numerische Mathematik in the late 1960's. The papers described algorithms for matrix computation and included programs in Algol 60, the International Algorithmic Language. A number of other authors contributed papers and programs to the series, which was known as the Handbook Series.

All of the papers in the Handbook Series were eventually collected in an actual handbook, edited by Wilkinson and Christian Reinsch, and published by Springer-Verlag in 1971. Even today this Wilkinson and Reinsch handbook remains an important resource for work in matrix computation. The Algol codes provide a readable description of the techniques that lie at the heart of the matrix library in MATLAB.

Argonne

Argonne National Laboratory is a located in a forest preserve 25 miles southwest of Chicago. A herd of white deer roam freely on the grounds. The laboratory carries out a program of primarily unclassified research for the U.S. Department of Energy on basic science, energy, transportation, environmental issues, computing and national security.

Very few people could actually run the programs in the Handbook because very few computer centers had Algol compilers. Around 1970, the Mathematics and Computer Science Division, MCS, at Argonne began a project to translate a portion of the Handbook into Fortran, and to study the subsequent testing, distribution and support of the software. The project involved the subset of the Handbook devoted to eigenvalue computation because most computer center libraries already had satisfactory linear equation solvers, but few had the new eigenvalue methods that Wilkinson and colleagues had developed.

The project was called NATS, National Activity to Test Software. It was supported by the National Science Foundation and what was then called the Atomic Energy Commission. The stated objective of the project was research into the software development, testing and distribution process. The software itself, which became known as EISPACK, for Eigensysem Package, was just a byproduct. Simply writing good software wasn't regarded as supportable research. It wasn't in 1970, and still isn't in 2013.

Jim was an important participant in the NATS/EISPACK project. He visited Argonne at least once a year, usually after the Michigan summer course. He didn't write any Fortran code, but he was an active consultant in the process of rewriting what was mostly his Algol. And, the testing process in those days was more complicated because there were so many machines with different arithmetic behavior.

If you were involved in mathematical software back then, Argonne was the place to visit. I moved to the University of New Mexico in 1972 and came to Argonne for most of every summer. (When I walked in, Jim Boyle would say, "Moler is here, it must be June.") There were always lots of international visitors. In addition to Wilkinson, there were many other Brits. MCS had a close connection with NAG, the Numerical Algorithms Group in the UK. At the annual picnic one summer, we got the NAG chaps to play softball. That resulted in a challenge to try cricket and they brought the necessary equipment with them the next summer. Jim was a cricket connoisseur. He gave us a lecture on cricket strategy at the same blackboard he is using in the photo above. A day or two later at the next picnic, he continued the lecture. Jack Dongarra found this faded color photo in one of his shoe boxes.

The EISPACK project was followed by the LINPACK project, also centered at Argonne. These two projects motivated the first MATLAB. Wilkinson wasn't directly involved in MATLAB, but he certainly was the primary source of the algorithms that formed its original mathematical core.

References

[1] Leslie Fox, James Hardy Wilkinson, Biographical Memoirs of Fellows of the Royal Society 33 (1987): 669-708. <http://rsbm.royalsocietypublishing.org/content/33/670>

[2] Nick Higham's Gallery of Wilkinson Photos, <http://www.maths.manchester.ac.uk/~higham/photos/wilkinson/index.htm>.

[3] History of Householder Symposia, <http://sites.uclouvain.be/HHXIX/SymposiumHistory>.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

The inverse Hilbert matrix, invhilb, has recently made surprise appearances in Cody, the programming game on MATLAB Central, and one of Ned's posts in the MATLAB Spoken Here blog. Inverse Hilbert matrices had nearly been forgotten in MATLAB. Their comeback is due to the sign pattern of their entries. But I want to take you back to their original role demonstrating ill conditioning in numerical calculation.

Contents

   T = invhilb(8);
   imagesc(sign(T))
   axis image, axis off
   colormap(copper)

Hilbert Matrix

The Hilbert matrix is the first matrix I ever knew by name. I met it in my first numerical analysis course, when I was a junior at Caltech in 1959. The matrix comes from least squares approximation by monomials, $x^{j}$, on the unit interval $[0,1]$. The elements of $H_n$ are inner products

$$ h_{i,j} = \int_0^1 x^{i-1} x^{j-1} dx = \frac{1}{i+j-1} $$

Here is the code to generate the 6-by-6 Hilbert matrix in a single array operation. Similar code is used in the function hilb that is distributed with MATLAB. I am using single precision so we can see roundoff error in matrices that print nicely in the width of one blog page.

   format compact
   format long
   n = 6
   J = 1:n;
   J = J(ones(n,1),:);
   I = J';
   E = single(ones(n,n));
   H = E./(I+J-1)

n =
     6
H =
   1.0000000   0.5000000   0.3333333   0.2500000   0.2000000   0.1666667
   0.5000000   0.3333333   0.2500000   0.2000000   0.1666667   0.1428571
   0.3333333   0.2500000   0.2000000   0.1666667   0.1428571   0.1250000
   0.2500000   0.2000000   0.1666667   0.1428571   0.1250000   0.1111111
   0.2000000   0.1666667   0.1428571   0.1250000   0.1111111   0.1000000
   0.1666667   0.1428571   0.1250000   0.1111111   0.1000000   0.0909091

Badly Conditioned

The monomials $x^{j}$ are nearly linearly dependent on $[0,1]$. So the Hilbert matrices are nearly singular. The condition of $H_6$ is comparable with single precision roundoff error and the condition of $H_{12}$ is comparable with double precision roundoff error.

   format short e
   [1/cond(hilb(6)) eps(single(1))]
   [1/cond(hilb(12)) eps(double(1))]

ans =
  6.6885e-08  1.1921e-07
ans =
   5.7268e-17   2.2204e-16

The $l_2$ condition number, $\kappa$, of $H_n$ grows exponentially.

$$ \kappa(H_n) \approx 0.01133e^{3.49n} $$

Inverse Hilbert Matrix

A Hilbert matrix qualifies as a Cauchy matrix, which is a matrix whose entries are of the form

$$ a_{i,j} = \frac{1}{x_i - y_j} $$

A classic Knuth homework problem or the Wikipedia entry on Cauchy matrices (see References) shows how it is possible to express the elements of the inverse of a Cauchy matrix in terms of products involving the $x_i$'s and $y_j$'s. For a Hilbert matrix, those products become binomial coefficients.

Here is a program that generates the inverse Hilbert matrix using doubly nested for loops and many scalar evaluations of binomial coefficients.

type invh

function T = invh(n)
   for i = 1:n
      for j = 1:n
         T(i,j) = (-1)^(i+j)*(i+j-1)*nchoosek(n+i-1,n-j)* ...
            nchoosek(n+j-1,n-i)*nchoosek(i+j-2,i-1)^2;
      end
   end
end

An Old Program

One of the first programs that I ever wrote generated the inverse Hilbert matrix from these binomial coefficients. It was written in absolute numeric machine language for the Burroughs 205 Datatron at Caltech around 1959. I needed to avoid time consuming floating point arithmetic, so I used recursion relationships among the coefficients. Here is the core of the MATLAB version of that old program. For a long time this has been distributed with MATLAB as invhilb. This is the function that the Cody participants discovered.

type invhilb_code

   T = zeros(n,n);
   p = n;
   for i = 1:n
       r = p*p;
       T(i,i) = r/(2*i-1);
       for j = i+1:n
           r = -((n-j+1)*r*(n+j-1))/(j-1)^2;
           T(i,j) = r/(i+j-1);
           T(j,i) = r/(i+j-1);
       end
       p = ((n-i)*p*(n+i))/(i^2);
   end

We can use either code to generate the exact inverse of $H_6$.

T = invhilb(6)

T =
          36        -630        3360       -7560        7560       -2772
        -630       14700      -88200      211680     -220500       83160
        3360      -88200      564480    -1411200     1512000     -582120
       -7560      211680    -1411200     3628800    -3969000     1552320
        7560     -220500     1512000    -3969000     4410000    -1746360
       -2772       83160     -582120     1552320    -1746360      698544

Since the elements of a Hilbert matrix are rational fractions, the elements of its inverse must also be rational fractions. But it turns out that the denominators of all the fractions are equal to one, so, as you can see, the inverse has integer entries. The integers are large, hence the large condition number.

Roundoff Error

Let's invert T and see how close we get to H.

   format long
   inv(single(T))

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND =  3.380340e-08. 
ans =
   1.0102043   0.5085138   0.3406332   0.2563883   0.2056789   0.1717779
   0.5085138   0.3404373   0.2560914   0.2053309   0.1714057   0.1471225
   0.3406332   0.2560914   0.2052232   0.1712378   0.1469208   0.1286576
   0.2563883   0.2053309   0.1712378   0.1468577   0.1285564   0.1143121
   0.2056789   0.1714057   0.1469208   0.1285564   0.1142728   0.1028457
   0.1717779   0.1471225   0.1286576   0.1143121   0.1028457   0.0934704

We can just recognize one or two digits of H. Since the condition number of T, is comparable with single precision roundoff level, we might have expected to lose all the digits. Roundoff error observed in actual matrix computation is usually not as bad as estimates based on condition numbers predict.

My Project

Over 50 years ago, after my first numerical analysis course, where I was introduced to matrix computation by Prof. John Todd, I did an individual undergrad research project under his direction. Part of the project involved studying the roundoff error in matrix inversion. First, I wrote a matrix inversion program. Then, in order to assess roundoff error, I generated Hilbert matrices, inverted them with my program, and compared the computed inverses with the exact inverses generated by the program that I've just described. If I'd had MATLAB at the time, the project would have gone something like this.

   n = 6;
   H = single(hilb(n));
   X = inv(H);
   T = single(invhilb(n));
   relerr = norm(X-T,inf)/norm(T)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND =  3.570602e-08. 
relerr =
   0.0470157

(My old matrix inversion program would not have estimated condition and issued a warning message.)

The Datatron floating point arithmetic had eight decimal digits, so I used values of n up to 7 and reported the observed relerr's as the result of roundoff error from matrix inversion by Gaussian elimination for this particular badly conditioned matrix. Case closed. Project grade: A.

Project Reexamined

I went on to grad school at Stanford and met George Forsythe and later Jim Wilkinson. As I came to appreciate the work that Wilkinson was doing, I realized that in my project at Caltech I had made a subtle but important mistake. The same mistake is still being made today by others. I had not actually inverted the Hilbert matrix. I had inverted the floating point approximation to the Hilbert matrix. It turns out that the initial step -- converting the fractions that are the elements of the Hilbert matrix to floating point numbers -- has a larger effect on the computed inverse than the inversion process itself. Even if my inversion program could somehow compute the exact inverse of the matrix it is given, it could not match the result generated by my theoretical inverse Hilbert matrix program.

I should have inverted the inverse, because I would not have to make any roundoff errors to generate the input. Since the elements of the starting matrix are integers, they can be represented exactly in floating point, at least if n is not too large. I should have interchanged the roles of T and H in my project.

   n = 6;
   T = single(invhilb(n));
   X = inv(T);
   H = single(hilb(n));
   relerr = norm(X-H,inf)/norm(H)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND =  3.380340e-08. 
relerr =
   0.0266826

This confirms what Wilkinson had been saying, and proving. The roundoff error caused by inverting a matrix using Gaussian elimination with partial pivoting is comparable with -- in this case actually smaller than -- the error caused by rounding the data in the first place.

Checkerboard

You can see the +/- sign pattern of the elements in the inverse Hilbert matrix by examining T or the (-1)^(i+j) multiplying the binomial coefficients in invh. That's what caught the attention of the folks playing Cody. I have no idea how they came across it.

References

[1] Donald E. Knuth, The Art of Computer Programming, vol. 1, Problems 1.2.3-41 and 1.2.3-45, pp. 36-37, Addison-Wesley, 1968.

[2] Wikipedia entry on Cauchy matrix, <http://en.wikipedia.org/wiki/Cauchy_matrix>.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

I investigated the reduced penultimate remainder algorithm in an undergraduate research project under professor John Todd at Caltech in 1961. I remember it today for two reasons. First, I learned what penultimate means. And second, it is the most obscure, impractical algorithm that I know. I suspect none of my readers have ever heard of it.

Contents

Penultimate

Penultimate means "next to last". November is the penultimate month of the year. In the fifty-plus years since I learned its meaning, I've never had a chance to use penultimate in everyday conversation, or in my writing. I think of it as British English. Is it?

Reduced Penultimate Remainder

Start with two polynomials, both with leading coefficient equal to one. When you divide one polynomial into the other using long division, you usually stop when the degree of the remainder is one less than the degree of the divisor. But instead, stop one step earlier, with the penultimate remainder. Then scale the leading coefficient of the remainder to be equal to one. This is the reduced penultimate remainder. If the original divisor had actually been an exact factor of the dividend, then this reduced remainder would be equal to the divisor. Iterating the process with exact divisors exhibits fixed points.

An Example

Our example involves dividing $q(x) = x^2 - 10x + 5$ into $p(x) = x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120$. Here is a scan of one step of the long division carried out by hand. (I don't know enough LaTeX to be able to do the typesetting.)

The reduced penultimate remainder is $r(x) = x^2 + 1.24x - 1.20$. Since $r(x)$ is not close to $q(x)$, we can conclude that $q(x)$ is not close to a factor of $p(x)$ and consequently that the roots of $q(x)$ are not close to the roots of $p(x)$. I will return to this example later.

Reduced Penultimate Remainder Algorithm

If the divisor happens to be a factor of the dividend, then the reduced penultimate remainder process is a fixed point iteration. In 1943, Shi-nge Lin suggested that starting with an initial approximation that is nearly a factor and repeating the RPR process, the iteration might converge to a factor and hence could be used to find roots of polynomials. I do not know who Shi-nge Lin was and have not actually seen his original paper.

A. C. Aitken

A. C. Aitken(1895-1967) wrote papers about Lin's idea in 1952 and 1955. Aitken was born in New Zealand, educated in Scotland, and became a professor of mathematics at the University of Edinburgh. He made many contributions to many areas of mathematics. In numerical analysis he is best known for the delta-squared sequence acceleration technique. His prodigious memory made some aspects of his life unpleasant. He served with the New Zealand forces in World War I and his memories of the battle at Gallipoli haunted him throughout his life.

John Todd was the professor at Caltech who introduced me to numerical analysis, matrices, and computing. I believe that Todd knew Aitken personally and he was certainly familiar with his work. So he suggested I investigate the method and I read Aitken's papers. I wrote programs in machine language for the Burroughs 205 Datatron. Unfortunately, about all I can remember about the project is the name of the algorithm.

Householder's Comments

Alston Householder(1904-1993) is one of my heroes. I knew him quite well. He's the fourth guy from the left in our photo of the organizing committee of the 1964 Gatlinburg Conference on Numerical Algebra

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

I just learned that Householder wrote a paper about the Reduced Penultimate Remainder algorithm in 1970, ten years after I did my undergraduate research project. I didn't know about that paper until I came across it while writing this blog.

Householder makes some very candid comments about the earlier papers by Lin and Aitken. He says "Lin's development is complicated and unsuggestive". He says "Aitken's derivation of the conditions for convergence is very sketchy and leaves much to the reader." Householder's paper was published in the SIAM Journal of Applied Math in 1970. I don't think we could make comments like these in SIAM papers today.

Aitken had introduced a device called the catalytic multiplier and Householder tries to explain it. I do not understand his paper at all. Near the end of the paper he says "In notation elsewhere used by the author ..." That is not very helpful. I'm afraid that is this particular paper Alston is no clearer than the authors he is criticizing. I can say this in a blog.

MATLAB Code

type rpr

function q = rpr(p,q)
% RPR Reduced Penultimate Remainder.
% r = rpr(p,q), for polynomials p and q.
% Assume p(1) = 1, q(1) = 1.
for c = 1:15
   % Polynomial long division.
   % Stop when degree of r = degree of q.
   n = length(p);
   m = length(q);
   r = p(1:m);
   for k = m+1:n
      r = r - r(1)*q;
      r = [r(2:end) p(k)];
   end
   q = r/r(1);
   disp(q)
end

Example, continued

The example I used in hand calculation above is the polynomial whose roots are the integers from one to five.

p = charpoly(diag(1:5))

p =
     1   -15    85  -225   274  -120

The starting guess is a deliberately bad approximation to a quadratic factor.

q = [1 -10 5]

q =
     1   -10     5

Let's run 15 iterations

format long
format compact
rpr(p,q)

   1.000000000000000   1.240000000000000  -1.200000000000000
   1.000000000000000  -1.067114979620489   0.318854992571954
   1.000000000000000  -1.723550954295960   0.821587108699062
   1.000000000000000  -2.062229096576079   1.106543069892922
   1.000000000000000  -2.272884387512085   1.294527998591067
   1.000000000000000  -2.417143826065064   1.428233837285976
   1.000000000000000  -2.522010435111559   1.527880806063684
   1.000000000000000  -2.601434375494343   1.604615444108659
   1.000000000000000  -2.663426234246922   1.665180563853953
   1.000000000000000  -2.712940049546721   1.713920778923779
   1.000000000000000  -2.753213571017107   1.753767759646395
   1.000000000000000  -2.786455813893868   1.786771702201661
   1.000000000000000  -2.814227004785662   1.814408345663866
   1.000000000000000  -2.837661124583663   1.837765841878627
   1.000000000000000  -2.857602508367706   1.857663278238952
ans =
   1.000000000000000  -2.857602508367706   1.857663278238952

It's converging, slowly, to the factor $x^2 - 3x + 2$ with roots $1$ and $2$. It takes about 300 iterations to reach this factor to double precision accuracy.

Another Example

One of my favorite polynomials is $p(x) = x^3 - 2x - 5$. It's only real root is near $2.09$. Let's try to find it.

p = [1 0 -2 -5]
q = [1 -2]
rpr(p,q)

p =
     1     0    -2    -5
q =
     1    -2
   1.000000000000000  -2.500000000000000
   1.000000000000000  -1.176470588235294
   1.000000000000000   8.117977528089890
   1.000000000000000  -0.078245352175702
   1.000000000000000   2.507676417722349
   1.000000000000000  -1.165924862052266
   1.000000000000000   7.804948516596162
   1.000000000000000  -0.084864830447043
   1.000000000000000   2.509035084827172
   1.000000000000000  -1.164074683720088
   1.000000000000000   7.752777799980693
   1.000000000000000  -0.086050279678108
   1.000000000000000   2.509290208665588
   1.000000000000000  -1.163727809437372
   1.000000000000000   7.743083431952472
ans =
   1.000000000000000   7.743083431952472

No luck, it's not going to converge. The real root repels any attempt to find it with this iteration.

Conclusions

Finding polynomial roots is not a very important task and this algorithm is not very good at it anyway. Its convergence properties are too unreliable. So for me, this is just a quirky memory. If anybody else has heard of the reduced penultimate remainder algorithm, please contribute a comment.

References

[1] A. C. AITKEN, Studies in practical mathematics. VII. On the theory of factorizing polynomials by iterated division, Proc. Roy. Soc. Edinburgh Sec..A, 63 (1952), pp. 326-335.

[2] A. C. AITKEN, Note on the acceleration of Lin's process of iterated penultimate remainder, Quart. J. Mech. Appl. Math., 8 (1955), pp. 251-258.

[3] SHI-NGE LIN, A method for finding roots of algebraic equations, J. Math. Phys. Mass. Inst. Tech., 22 (1943), pp. 60-77.

[4] ALSTON S. HOUSEHOLDER, On the penultimate remainder algorithm and the catalytic multiplier, SIAM J. Appl. Math., 19 (1970), pp. 668-671.

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b

Tuesday, January 8, 2013, would have been George Forsythe's 96th birthday. He passed away in 1972 at the age of 55. A pioneer in the establishment of computer science as an intellectual discipline, he was my Ph.D. thesis advisor, colleague, and friend.

Contents

Early Life

George Elmer Forsythe was born in State College, Pennsylvania, in 1917 and moved to Ann Arbor, Michigan, when he was a small boy. He went to Swarthmore as an undergraduate, where he majored in mathematics, and then to Brown for a Ph.D., also in mathematics, which he received in 1941.

Forsythe served in the Air Force in World War II. Working there as a mathematician and meteorologist, he developed some of the first methods of numerical weather prediction. He is the coauthor of an early textbook on numerical meteorology.

Institute for Numerical Analysis

After the war, the National Bureau of Standards, with funding from the Office of Naval Research, sponsored the development of two electronic computers. The one on the west coast, at UCLA, was known as the SWAC, for Standards Western Automatic Computer.

The INA, the Institute for Numerical Analysis, was established at UCLA in 1947 to support the SWAC. Forsythe was one of its principal members.

Our photo, taken at UCLA sometime in the early 1950s, shows Forsythe in the center. Olga Taussky-Todd is the only woman in the photo and her husband, John Todd, is in the back row peering out between George and Olga.

The other institute members in the front row are Mark Kac, who would later ask if one could hear the shape of a drum; J. Barkley Rosser, who created the amazing rosser test matrix in MATLAB; Wolfgang Wasow, who would later coauthor a seminal book with Forsythe about the numerical solution of PDEs; and Magnes Hestenes, coinventor of the conjugate gradient method.

The Todds and Caltech

The INA was dissolved in 1953 as collateral damage from political scandals at the Bureau of Standards involving McCarthyism and congressional hearings about battery additives. The Todds moved a few miles across Los Angeles to Pasadena and Caltech and Forsythe moved a few hundred miles north to Palo Alto and Stanford.

I enrolled in Caltech in 1957 and a couple of years later was taking courses from John Todd. He introduced me to matrices, numerical analysis and computing and set me in the direction that brought me to where I am today. Olga would occasionally guest lecture in John's classes on her specialty, matrix theory.

One quarter in my junior year I got a B- in advanced calculus and an A+ in numerical analysis. I never looked back. When it was time to go to grad course, Todd recommended Stanford and his friend George Forsythe. I didn't apply anywhere else.

Stanford Computer Science

In 1961 Forsythe was a math professor at Stanford. Computing was just a master's degree specialization within the math department. I am not sure if it was even yet called "Computer Science". But Forsythe and others could see the need for a full fledged academic program, department, and discipline. They wanted to see Stanford appoint faculty and offer courses in new fields such as artificial intelligence and programming languages that certainly were not research mathematics.

So in January, 1965, the Stanford Computer Science Department was born. Forsythe was its founding department chairman and remained head until his death seven years later.

Today Stanford Computer Science is famous as the incubator for Google, Yahoo, and many other not-so-famous Silicon Valley commercial successes. But at their 50th anniversary celebration Stanford's president John Hennessy, himself a CS professor, reminded his audience that the department also deserves to be known for its some of its world renowned academics who did not start companies, including John McCarthy and Donald Knuth. These professors were recruited by George Forsythe.

Computer Science Discipline

Forsythe not only started Computer Science at Stanford, he was one of its advocates as an academic and intellectual discipline nationally and internationally. He served as president of the ACM early in its history. He was an editor of several journals. He was asked to advise other universities. He wrote influential articles and letters.

Shortly after Forsythe's death, Donald Knuth wrote, "It is generally agreed that he [Forsythe], more than any other man, is responsible for the rapid development of computer science in the world's colleges and universities."

Gatlinburg and Householder Meetings

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

This photo is one of the very first images distributed with MATLAB. The photo, taken in 1964, shows the organizing committee of the Gatlinburg conferences on numerical algebra. All of these men were to have some influence on the development of MATLAB. From left to right are:

Householder was head of computing at Oak Ridge National Lab and a professor at the University of Tennessee in Knoxville. He organized the first two conferences on numerical algebra in the early 1960s. The meetings were held in Gatlinburg, Tennessee, which was then a lovely, small resort town in the Smoky Mountains.

These were intensive, week-long, research workshops with attendance by invitation only from this committee. When it came to organizing the third Gatlinburg conference in 1964, Forsythe wanted to invite one of his graduate students -- me. But the other committee members objected, students had not been invited to the previous meetings. I understand that Forsythe said if Moler could not come, neither would he. So I was invited. I was the only student there, among 40 or 50 of the most important people in my field. It was one of the most significant events of my professional career. And I was only 25 years old.

An Oak Ridge photographer took 8-by-10 black-and-white glossy photos. I kept the one of the organizing committee. Twenty-five years later we scanned it in to become a MATLAB, and Internet, classic.

After Householder retired, the committee evolved, and the locations changed. The meetings are now called the Householder Conferences and are held every four years. The locations alternate between North American and Europe. Attendance is still by invitation only, but anybody can apply. The next conference will be Householder XIX in June 2014, in Belgium, hosted by Universite catholique de Louvain and the Katholieke Universiteit Leuven.

Prentice-Hall Series

Forsythe served as advising editor of the prestigious Prentice-Hall Series in Automatic Computation. He solicited manuscripts from many influential authors. With more than 75 titles published in the early years, this series helped define the Computer Science discipline. The distinctive red and white dust jackets can still be seen on many library shelves.

Our Books

Two of the books in the Prentice-Hall series were by Forsythe and Moler, and by Forsythe, Malcolm and Moler. Before I was known for MATLAB, I could meet people at conferences and say, "You know those books by Forsythe and somebody, well I'm that somebody."

The Forsythe and Moler book, "Computer Solution of Linear Algebraic Systems", evolved from notes for the second quarter of a two-quarter basic numerical analysis course that Forsythe taught in the early 1960's and that I taught in 1965. The first quarter used a text by Peter Henrici that covered all of the material in our syllabus except matrix computation. The notes that George had written for the second quarter eventually became the book.

I wrote programs DECOMP and SOLVE that implement column-oriented Gaussian elimination in three languages, Algol 60, Fortran, and PL/1. When I was working on the book, the year was 1966 and Algol, the International Algorithmic Language, was the generally accepted standard for publishing mathematical material. Fortran was already over 10 years old and we didn't expect it to last much longer. PL/1 had been developed recently by IBM and we thought it might represent the future. You can see how accurate those predictions turned out to be.

The Forsythe, Malcolm, and Moler book, "Computer Methods for Mathematical Computations", had a longer incubation period. After I left grad school at Stanford, Forsythe started a more elementary course in numerical methods intended for engineers and other students outside of math and computer science. The course centered around a collection of well written mathematical software. Mike Malcolm was his teaching assistant.

Forsythe fell ill before he and Malcolm could refine their course notes into a book and they asked me to help complete the project. The reviewers of our first draft said that it read like it was written by three different authors who never read each other's chapters. They were right. It took Mike and I five years to smooth out the exposition, by which time we had rewritten much of the software. And, oh yes, it was 1977 and the clear choice for programming language was Fortran.

Forsythe was a great believer that published computer programs should be a means of scientific discourse between humans, as well as a means of controlling machines. He would read and edit code as carefully as he would prose, which was very carefully. 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.

The Forsythe Tree

One of Forsythe's most important legacies is his students and his students' students. The Forsythe Tree has George at the root and links that are Ph.D. theses. There are 17 first-level nodes, one of which is me. The version on the Web is 10 years old and woefully out of date. Even this version has a depth of five or six and over 200 nodes. There are dozens of university professors and several deans, provosts and presidents. I suspect that a current version might be almost twice as deep and have several hundred nodes. It reaches far into academia and industry throughout the world. We should bring the tree up to date.

Postscript

George never got to see MATLAB. He passed away before its birth. He would have loved it.

Reference

<http://www.stanford.edu/dept/ICME/docs/history/forsythe_knuth.pdf.> Knuth, Donald E. (1972). "George Forsythe and the Development of Computer Science". Communications of the ACM 15 (8).

Get the MATLAB code

Published with MATLAB® R2012b

Published with MATLAB® R2012b