A report about a possible bug in `format bank` and a visit to a local hardware store made me realize that doing decimal arithmetic with binary floating point numbers is like tightening a European bolt with an American socket wrench.... read more >>

A report about a possible bug in `format bank` and a visit to a local hardware store made me realize that doing decimal arithmetic with binary floating point numbers is like tightening a European bolt with an American socket wrench.

It's mid-April and so those of us who file United States income taxes have that chore to do. Years ago, I used MATLAB to help with my taxes. I had a program named `form1040.m` that had one statement for each line on the tax form. I just had to enter my income and deductions. Then MATLAB would do all the arithmetic.

If we're really meticulous about our financial records, we keep track of things to the nearest penny. So line 28 of `form1040` might have been something like

interest = 48.35

interest = 48.3500

I didn't like those trailing zeros in the output. So I introduced

```
format bank
```

into MATLAB. Now

interest

interest = 48.35

is printed with just two digits after the decimal point.

`format bank` has turned out to be useful more broadly and is still in today's MATLAB.

We recently had a user ask about the rounding rule employed by `format bank`. What if a value fails halfway between two possible outputs? Which is chosen and why?

Here's the example that prompted the user's query. Start with

```
format short
```

so we can see four decimal places.

x = (5:10:45)'/1000 y = 1+x z = 2+x

x = 0.0050 0.0150 0.0250 0.0350 0.0450 y = 1.0050 1.0150 1.0250 1.0350 1.0450 z = 2.0050 2.0150 2.0250 2.0350 2.0450

These values appear to fall halfway between pairs of two-digit decimal fractions. Let's see how the ties are broken.

```
format bank
x
y
z
```

x = 0.01 0.01 0.03 0.04 0.04 y = 1.00 1.01 1.02 1.03 1.04 z = 2.00 2.02 2.02 2.04 2.04

Look at the last digits. What mysterious hand is at work here? Three of the `x`'s, `x(1)`, `x(3)`, and `x(4)`, have been rounded up. None of the `y`'s have been rounded up. Two of the `z`'s, `z(2)` and `z(4)`, have been rounded up.

Email circulated internally at MathWorks for a few days after this was reported suggesting various explanations. Is it a bug in `format bank`? In the I/O library? Does it depend upon which compiler was used to build MATLAB? Do we see the same behavior on the PC and the Mac? Has it always been this way? These are the usual questions that we ask ourselves when we see curious behavior.

Do you know what's happening?

Well, none of the suspects I just mentioned is the culprit. The fact is none of the numbers fall exactly on a midpoint. Think binary, not decimal. A value like 0.005 expressed as a decimal fraction cannot be represented exactly as a binary floating point number. Decimal fractions fall in between the binary fractions, and rarely fall precisely half way.

To understand what is going on, the Symbolic Toolbox function

sym(x,'e')

is your friend. The `'e'` flag is provided for this purpose. The `help` entry says

'e' stands for 'estimate error'. The 'r' form is supplemented by a term involving the variable 'eps' which estimates the difference between the theoretical rational expression and its actual floating point value. For example, sym(3*pi/4,'e') is 3*pi/4-103*eps/249.

To see how this works for the situation encountered here.

symx = sym(x,'e') symy = sym(y,'e') symz = sym(z,'e')

symx = eps/2133 + 1/200 3/200 - eps/400 eps/160 + 1/40 (3*eps)/200 + 7/200 9/200 - (3*eps)/400 symy = 201/200 - (12*eps)/25 203/200 - (11*eps)/25 41/40 - (2*eps)/5 207/200 - (9*eps)/25 209/200 - (8*eps)/25 symz = 401/200 - (12*eps)/25 (14*eps)/25 + 403/200 81/40 - (2*eps)/5 (16*eps)/25 + 407/200 409/200 - (8*eps)/25

The output is not as clear as I would like to see it, but I can pick off the sign of the error terms and find

x + - + + - y - - - - - z - + - + -

Or, I can compute the signs with

```
format short
sigx = sign(double(symx - x))
sigy = sign(double(symy - y))
sigz = sign(double(symz - z))
```

sigx = 1 -1 1 1 -1 sigy = -1 -1 -1 -1 -1 sigz = -1 1 -1 1 -1

The sign of the error term tells us whether the floating point numbers stored in `x`, `y` and `z` are larger or smaller than the anticipated decimal fractions. After this initial input, there is essentially no more roundoff error. `format bank` will round up or down from the decimal value accordingly. Again, it is just doing its job on the input it is given.

Photo credit: http://toolguyd.com/

A high-end socket wrench set includes both metric (left) and fractional inch (right) sizes. Again, think decimal and binary. Metric sizes of nuts and bolts and wrenches are specified in decimal fractions, while the denominators in fractional inch sizes are powers of two.

Conversion charts between the metric and fractional inch standards abound on the internet. Here is a link to one of them: Wrench Conversion Chart.

But we can easily compute our own conversion chart. And in the process compute `delta`, the relative error made when the closest binary wrench is used on a decimal bolt.

make_chart

Inch Metric delta 1/64 0.016 1/32 0.031 3/64 0.047 1mm 0.039 -0.191 1/16 0.063 5/64 0.078 2mm 0.079 0.008 3/32 0.094 7/64 0.109 1/8 0.125 3mm 0.118 -0.058 9/64 0.141 5/32 0.156 4mm 0.157 0.008 11/64 0.172 3/16 0.188 13/64 0.203 5mm 0.197 -0.032 7/32 0.219 15/64 0.234 6mm 0.236 0.008 1/4 0.250 9/32 0.281 7mm 0.276 -0.021 5/16 0.313 8mm 0.315 0.008 11/32 0.344 9mm 0.354 0.030 3/8 0.375 13/32 0.406 10mm 0.394 -0.032 7/16 0.438 15/32 0.469 12mm 0.472 0.008 1/2 0.500 9/16 0.563 14mm 0.551 -0.021 5/8 0.625 16mm 0.630 0.008 11/16 0.688 18mm 0.709 0.030 3/4 0.750 13/16 0.813 20mm 0.787 -0.032 7/8 0.875 22mm 0.866 -0.010 15/16 0.938 24mm 0.945 0.008 1 1.000

Let's plot those relative errors. Except for the small sizes, where this set doesn't have enough wrenches, the relative error is only a few percent. But that's still enough to produce a damaging fit on a tight nut.

bar(k,d) axis([0 25 -.05 .05]) xlabel('millimeters') ylabel('delta')

You might notice that my conversion chart, like all such charts, and like the wrenches themselves, has a little bit of floating point character. The spacing of the entries is not uniform. The spacing between the binary values is 1/64, then 1/32, then 1/16. The spacing of the metric values is 1mm at the top of the chart and 2mm later on.

Suppose we want to tighten a 10mm nut and all we have are binary wrenches. The diameter of the nut in inches is

meter = 39.370079; d = 10*meter/1000

d = 0.3937

Consulting our chart, we see that a 13/32 wrench is the best fit, but it's a little too large. 10mm lies between these two binary values.

[floor(32*d) ceil(32*d)] b1 = 12/32; b2 = 13/32; [b1 d b2]'

ans = 12 13 ans = 0.3750 0.3937 0.4063

The fraction of the interval is

f = (d - b1)/(b2 - b1)

f = 0.5984

10mm is about 60% of the way from 12/32 inches to 13/32 inches.

Now let's turn to floating point numbers. What happens when execute this MATLAB statement?

x = 1/10

x = 0.1000

One is divided by ten and the closest floating point number is stored in `x`. The same value is produced by the statement

x = 0.1

x = 0.1000

The resulting `x` lies in the interval between 1/16 and 1/8. The floating point numbers in this interval are uniformly spaced with a separation of

e = eps(1/16)

e = 1.3878e-17

This is `2^-56`. Let

e = sym(e)

e = 1/72057594037927936

This value `e` plays the role that 1/64 plays for my wrenches.

The result of `x = 1/10` lies between these two binary fractions.

b1 = floor(1/(10*e))*e b2 = ceil(1/(10*e))*e

b1 = 7205759403792793/72057594037927936 b2 = 3602879701896397/36028797018963968

The tiny interval of length `e` is

c = [b1 1/10 b2]'

c = 7205759403792793/72057594037927936 1/10 3602879701896397/36028797018963968

In decimal

vpa(c)

ans = 0.099999999999999991673327315311326 0.1 0.10000000000000000555111512312578

Where in this interval does 1/10 lie?

f = double((1/10 - b1)/e)

f = 0.6000

So 1/10 is about 60% of the way from `b1` to `b2` and so is closer to `b2` than to `b1`. The two statements

x = 1/10 x = double(b2)

x = 0.1000 x = 0.1000

store exactly the same value in `x`.

The 13/32 wrench is the closest tool in the binary collection to the 10mm nut.

Get
the MATLAB code

Published with MATLAB® R2017a

I've long known that my Erdös Number is 3. This means that the length of the path on the graph of academic coauthorship between me and mathematician Paul Erdös is 3. Somewhat to my surprise, I recently discovered that I can also trace a chain of coauthorship to Donald J. Trump. My Trump number is 5.... read more >>

]]>I've long known that my Erdös Number is 3. This means that the length of the path on the graph of academic coauthorship between me and mathematician Paul Erdös is 3. Somewhat to my surprise, I recently discovered that I can also trace a chain of coauthorship to Donald J. Trump. My Trump number is 5.

Paul Erdös. Photo: http://www.iitvidya.com/paul-erdos.

The *collaborative distance* between two authors is the length of the path of coauthorship of scientific papers, books, and articles connecting the two. If A and B are coauthors, then the collaborative distance between them is 1. Furthermore, if B and C are also coauthors, then the collaborative distance between A and C is 2. And so on. If there is no chain of coauthorship, then the collaborative distance is infinite.

Paul Erdös (1911-1996) was the world's most prolific modern mathematician. He wrote 1,523 papers with 511 distinct coauthors. These 511 people have Erdös number equal to 1. And these authors have, in turn, written papers with over 11,000 other people. That means that over 11,000 people have Erdös number of 2. My estimate is that a few hundred thousand people have Erdös number of 3. I'm one of them.

There are two length three paths of coauthorship between me and Erdös. Both go through my thesis advisor, George Forsythe. I wrote a blog post about Forsythe a few years ago. Forsythe has an Erdös number of 2 in two different ways because he wrote papers with his thesis advisor, William Feller, and with Ernst Straus, both of whom had worked directly with Erdös.

An Erdös Number calculator is available at MathSciNet Collaboration Distance. More than you every wanted to know about Erdös numbers is available at the Oakland University Erdös Number Project.

Steve Johnson is a buddy of mine who also has Erdös number of 3. His path goes through Jeff Ullman and Ron Graham to Paul Erdös. But that's not why I bring him up today. In the 1970's Steve was part of the group at Bell Labs that developed Unix. He wrote the Unix tool Yacc (Yet Another Compiler Compiler), as well as the original C compiler, PCC, (Portable C Compiler). He is coauthor, along with three other Unix guys, of a paper about C in the 1978 issue of the Bell System Technical Journal that was devoted entirely to Unix.

Steve and I wrote a paper about compiling MATLAB, although the references to that paper on the Internet have my named spelled incorrectly. It is through this paper that I was surprised to find that my collaborative distance from Donald J. Trump is only 5.

Finite Trump numbers are possible because Trump's famous book, "The Art of the Deal", was actually ghost-written by a free-lance writer named Tony Schwartz. And Schwartz has coauthored many articles and books with other people. Some of these articles might not exactly be classified as academic papers, but what the heck.

The coauthorship path between me and Trump is through articles by John Winslow Morgan, a professor in the Harvard Business School who specializes in the history of technology. He has written articles with Schwarz and with Dennis Ritchie, one of the originators of Unix.

So, the path of length 5 is Moler - Johnson - Ritchie - Morgan - Schwartz - Trump.

**Erdös Number**

Forsythe, G. E. and Straus, E. G. *On best conditioned matrices.* Proc. Amer. Math. Soc. 6, (1955). 340–345.

Erdös, P., Lovász, L., Simmons, A., and Straus, E. G. *Dissection graphs of planar point sets.* A survey of combinatorial theory. (Proc. Internat. Sympos., Colorado State Univ., Fort Collins, Colo., 1971), 139–149. North-Holland, Amsterdam, 1973.

Feller, William, and George E. Forsythe. *New matrix transformations for obtaining characteristic vectors*. Quarterly of Applied Mathematics 8.4 (1951), 325-331.

Erdös, Paul, William Feller, and Harry Pollard. *A property of power series with positive coefficients.* Bull. Amer. Math. Soc 55.2 (1949): 201-204.

Forsythe, G. E. and C. B. Moler, Computer Solution of Linear Algebraic Systems, (Series in Automatic Computation) XI + 148, Prentice Hall, Englewood Cliffs, N.J. 1967.

Forsythe, George E., Malcolm, Michael A. and Moler, Cleve B., Computer Methods for Mathematical Computations, (Series in Automatic Computation) XI + 259, Prentice Hall, Englewood Cliffs, N.J. 1977.

**Trump Number**

Johnson, S. C. and C. Mohler (Moler), *Compiling MATLAB*, Proceedings of the USENIX Symposium on Very High Level Languages (VHLL), 119-27, Santa Fe, New Mexico, October 1994. USENIX Association.

Ritchie, D. M., Johnson, S. C., Lesk, M. E. and Kernighan, B. W., *UNIX Time-Sharing System: The C Programming Language.* Bell System Technical Journal, 57: 1991–2019, 1978.

Ritchie, D. M. and Morgan, J. W, *The Origins of UNIX*, Harvard Business Review, 48: 28-35, 1985.

Morgan, John Winslow and Schwartz, Tony, *Does UNIX Have A Future?*, MIT Technology Review, 53: 1-8, 1992.

Trump, Donald J. and Schwartz, Tony, The Art of the Deal, Ballantine Books, (paperback), 384 pp., 2004.

Get
the MATLAB code

Published with MATLAB® R2017a

A binary tree is an elegant way to represent and process Morse code. The new MATLAB graph object provides an elegant way to manipulate binary trees. A new app, `morse_tree`, based on this approach, is now available in version 2.40 of Cleve's Laboratory.... read more >>

A binary tree is an elegant way to represent and process Morse code. The new MATLAB graph object provides an elegant way to manipulate binary trees. A new app, `morse_tree`, based on this approach, is now available in version 2.40 of Cleve's Laboratory.

I have a chapter on Morse Code and binary trees in *Experiments with MATLAB*. Some of this blog post is taken from that chapter. But today's implentation is much more, as I said, elegant. In *EXM* I represented a binary tree as a cell array of cell arrays. Now I am using clever indexing into a linear character string, together with the MATLAB graph object. The resulting code is much shorter and more readable.

Here's all we have to do to get a picture of a binary tree. First we create the adjacency matrix and view its spy plot.

n = 31; j = 2:n; k = floor(j/2); A = sparse(k,j,1,n,n); spy(A)

Then we create and plot the directed graph. The MATLAB `digraph` object, and its supporting `plot` method, are doing all of the work.

G = digraph(A); Gp = plot(G); set(gca,'xtick',[],'ytick',[])

This is the complete binary tree with five levels, and consequently $2^5-1$ nodes. (In defiance of nature, computer scientists put the *root* of a tree at the top.)

Follow the nodes in numerical order, from top to bottom and left to right. You are doing a *breadth first traversal* of the structure. The *successors* of the node with index *j* have indices *2j* and *2j+1*. This makes old-fashioned linear indexing possible.

I want to save the coordinates of the nodes for use in the layout in my next plot.

x = Gp.XData; y = Gp.YData;

Morse code is no longer important commercially, but it still has some avid fans among hobbyists. Morse code was invented over 150 years ago, not by Samuel F. B. Morse, but by his colleague, Alfred Vail. It has been in widespread use ever since. The code consists of short *dots*, '.', and longer *dashes*, '-', separated by short and long spaces. You may be familiar with the international distress signal, '... --- ...', the code for "SOS", abbreviating "Save Our Ships" or perhaps "Save Our Souls". But did you notice that some modern cell phones signal '... -- ...', the code for "SMS", indicating activity of the "Short Message Service".

Until 2003, a license to operate an amateur radio in the US required minimal proficiency in Morse code. When I was in junior high school, I learned Morse code to get my ham license, and I've never forgotten it.

According to Wikipedia, in 2004, the International Telecommunication Union formally added a code for the ubiquitous email character, @, to the international Morse code standard. This was the first addition since World War I.

I could provide a table showing that '.-' is the code for A, '-...' the code for B, and so on. But I'm not going to do that, and my MATLAB program does not have such a table. As far as I am concerned, Morse code is *defined* by a MATLAB statement containing the 26 upper case letters of the English language, an asterisk, and four blanks.

```
morse = '*ETIANMSURWDKGOHVF L PJBXCYZQ ';
```

The asterisk at the start serves as the root of the tree. The four blanks will correspond to four missing nodes in the bottom level of the tree. So the binary tree is not going to be complete. The automatic layout in the `graph/plot` method does not draw a perfect tree. This is why I saved the node coordinates from the plot of the full tree. Let's remove the rows and columns of these potentially unlabeled nodes from the adjacency matrix, and the coordinates.

```
m = find(morse == ' ');
A(:,m) = [];
A(m,:) = [];
x(m) = [];
y(m) = [];
```

Convert the character array into a cell array of individual chars.

```
nodes = num2cell(morse);
nodes(m) = '';
```

Create the directed graph with these node labels.

G = digraph(A,nodes);

Plot the graph, using the layout saved from the full tree.

plot(G,'xdata',x,'ydata',y, ... 'showarrows','off'); set(gca,'xtick',[],'ytick',[])

If you follow the links downward and emit a *dot* when you go left and a *dash* when you go right, you have Morse code. For example start at the root, go left once, then right once. You will reach A and will have generated '.-', the code for A.

Morse code is already present in the defining character vector `morse`, even before we create and plot a graph. The character with index `5` is `A`. The characters with indices `2*5` and `2*5+1` are `R` and `W`.

j = 5; morse(j) morse(2*j:2*j+1)

ans = 'A' ans = 'RW'

Consequently, the Morse code for R and W is '.-.' and '.--'. This works everywhere in the first half of `morse`. The characters in the second half don't (yet) have successors.

Morse code can be extended by replacing the blanks in the basic tree with characters that are not in the 26 letter English alphabet and adding two more levels to the tree. Several different extensions have been in use at different times and different regions of the world. The Wikipedia page for Morse code includes a binary tree with some commonly used extensions. The fifth level includes the 10 decimal digits, several non-Latin chacters, and a few punctuation marks. The sixth level includes several punctuation marks.

Here is a screen shot from the new app, `morse_tree`, that is available in version 2.40 of Cleve's Laboratory. Clicking on the toggle labeled "extend" shows two more levels of the tree and some extensions, although not many as Wikipedia's tree. Clicking on a node in our tree highlights that node while the sound of the corresponding code is played. This screen shot shows A highlighted.

Text entered in the box below the tree will be encoded and the sound played.

The speed of the generated code is set by the control currently labeled "5 wpm", for five words per minute.

morse_tree_extended

Here's a snap quiz: What is Morse code for the @ sign? Download `morse_tree` to check your answer.

Get
the MATLAB code

Published with MATLAB® R2017a

The Lake Arrowhead Coauthor Graph came out of the Householder XII conference in 1993 at the UCLA conference center in the mountains north of San Bernardino. John Gilbert now remembers it as one of the first computational social network analyses he had ever seen. Today I revisit it using the new MATLAB graph object.... read more >>

]]>The Lake Arrowhead Coauthor Graph came out of the Householder XII conference in 1993 at the UCLA conference center in the mountains north of San Bernardino. John Gilbert now remembers it as one of the first computational social network analyses he had ever seen. Today I revisit it using the new MATLAB graph object.

I wrote a blog post about the Lake Arrowhead Coauthor Graph a couple of years ago.

At the Householder XII conference in 1993 Nick Trefethen posted a flip chart with Gene Golub's name in the center. He invited everyone present to add their name to the chart and draw lines connecting their name with the names of all their coauthors. The diagram grew denser throughout the week. At the end of the conference it was a graph with 104 vertices (or people) and 211 edges.

Nick framed the original chart and has it on the wall of his office at Oxford. Another Nick, Nick Hale, took this photo.

John Gilbert, Rob Schreiber and I entered the names and coauthor connections into MATLAB, creating an adjacency matrix `A`. For a long time the data was available as a `.mat` file from John's old web site at Xerox PARC. But John left PARC years ago, so I am now making the data available in human-readable form at the end of this blog post. There are two files.

A = arrowhead_adjacency; names = arrowhead_names;

The ordering of the authors in the original data is determined by how we happened to read them off the flip chart. Golub is first, and has the most connections to other authors. A symmetric reverse Cuthill-McGee ordering of the rest of the nodes puts coauthors near each other and reduces the band width of the adjacency matrix.

r = [1 symrcm(A(2:end,2:end))+1]; A = A(r,r); names = names(r); spy(A)

Let's make a MATLAB graph and plot it with the `circle` layout. By default, the `graph/plot` method omits the node names if there are more than 100 of them. This removes clutter and is usually a good idea. The resulting unannotated circle layout is a pretty picture, but not very informative.

The RCM ordering places most of the graph edges, except those connecting to Golub, near the circumference of the circle.

G = graph(A,names); plot(G,'layout','circle') axis square

Move Golub to the center.

p = plot(G,'Layout','circle'); idGolub = findnode(G,'Golub'); p.XData(idGolub) = 0; p.YData(idGolub) = 0; axis square

Placing the rest of the names around the outside of the circle requires some custom processing.

axis off t = text(p.XData,p.YData,names); theta = linspace(0,360,numel(t)+1); for k = 1:numel(t) t(k).Rotation = theta(k); t(k).String = [blanks(4) t(k).String]; t(k).FontSize = 8; end

This plot can now be compared with the plot we originally produced at Householder XII.

In this two-dimensional circle plot it's nearly impossible to follow the links from one author to his or her coauthors. 'Moler' is at about 5 o'clock on the circle. What is the degree of my node? Who are my coauthors?

A three-dimensional layout and judicious use of the mouse in conjunction with `cameratoolbar` makes it possible to examine this graph in detail.

p = plot(G,'layout','force3','nodelabelmode','auto'); % Setting the nodelabelmode says include all of the names. axis off vis3d set(gca,'clipping','off') zoom(2) gif_frame('arrowhead_name.gif') gif_frame(2)

Let's highlight me and my coauthors.

me = findnode(G,'Moler') friends = neighbors(G,me)' color = get(gca,'colororder'); highlight(p,[me friends],'nodecolor',color(2,:),'markersize',6) highlight(p,me,friends,'edgecolor',color(2,:),'linewidth',2) gif_frame(2)

me = 88 friends = 61 62 63 68 75 77 89 98

Make my node the center of this little universe.

set(p,'xdata',p.XData - p.XData(me), ... 'ydata',p.YData - p.YData(me), ... 'zdata',p.ZData - p.ZData(me)) gif_frame(2)

Now zoom and rotate. This is much better done in person with `cameratoolbar`. The best we can do in this blog is an animated gif.

% Zoom for k = 1:4 zoom(sqrt(2)) gif_frame end % Rotate [a,e] = view; for k = 1:4 view(a,e-12*k) gif_frame end gif_frame(5)

Here is the animated .gif centered on me and my coauthors.

And here we zoom in on Paul VanDooren by replacing

me = findnode(G,'Moler')

with

me = findnode(G,'VanDooren')

function A = arrowhead_adjacency k = [ 1 1 1 1 1 3 3 1 1 1 1 2 3 3 38 1 2 ... 32 42 3 32 42 43 44 44 1 46 1 1 37 13 49 1 32 1 ... 44 53 1 53 44 45 54 1 1 20 41 44 1 39 52 60 15 58 ... 61 3 63 41 43 1 15 35 51 61 41 43 59 65 53 3 58 63 ... 42 1 41 44 1 46 47 1 10 15 32 35 62 66 2 19 39 58 ... 61 66 69 37 53 1 9 72 42 8 21 32 42 77 1 7 15 46 ... 47 49 6 42 66 1 3 25 41 43 59 65 67 80 33 52 59 27 ... 28 46 74 31 51 46 1 22 23 24 25 35 41 52 81 3 21 55 ... 60 44 58 62 4 17 12 19 23 72 77 89 1 23 86 1 23 86 ... 91 4 72 90 19 90 48 1 18 25 77 82 86 17 52 66 86 89 ... 16 32 34 58 66 80 81 1 34 74 77 14 82 77 79 49 102 26 ... 32 39 44 48 51 52 56 63 69 74 78 81 83 90]; j = [ 2 3 5 10 11 29 30 32 34 35 36 36 38 40 40 41 42 ... 42 43 44 44 44 44 45 46 47 47 48 49 49 50 50 51 51 53 ... 54 54 55 55 56 56 57 58 59 59 59 59 60 60 60 61 62 62 ... 62 63 64 65 65 66 66 66 66 66 67 67 67 67 68 69 69 69 ... 70 71 71 71 72 72 72 73 73 73 73 73 73 73 74 74 74 74 ... 74 74 74 75 75 76 76 76 77 78 78 78 78 78 79 79 79 79 ... 79 79 80 80 80 81 81 81 81 81 81 81 81 81 82 82 82 83 ... 83 83 83 84 84 85 86 86 86 86 86 86 86 86 86 87 87 87 ... 87 88 88 88 89 89 90 90 90 90 90 90 91 91 91 92 92 92 ... 92 93 93 93 94 94 95 96 96 96 96 96 96 97 97 97 97 97 ... 98 98 98 98 98 98 98 99 99 99 99 100 100 101 102 103 103 104 ... 104 104 104 104 104 104 104 104 104 104 104 104 104 104]; U = sparse(k,j,1,104,104); A = logical(U + U'); end

function names = arrowhead_names names = {'Golub', 'Wilkinson', 'TChan', 'He', 'Varah', 'Kenney', ... 'Ashby', 'LeBorne', 'Modersitzki', 'Overton', 'Ernst', 'Borges', ... 'Kincaid', 'Crevelli', 'Boley', 'Anjos', 'Byers', 'Benzi', 'Kaufman', ... 'Gu', 'Fierro', 'Nagy', 'Harrod', 'Pan', 'Funderlic', 'Edelman', ... 'Cullum', 'Strakos', 'Saied', 'Ong', 'Wold', 'VanLoan', ... 'Chandrasekaran', 'Saunders', 'Bojanczyk', 'Dubrulle', 'Marek', ... 'Kuo', 'Bai', 'Tong', 'George', 'Moler', 'Gilbert', 'Schreiber', ... 'Pothen', 'NTrefethen', 'Nachtigal', 'Kahan', 'Varga', 'Young', ... 'Kagstrom', 'Barlow', 'Widlund', 'Bjorstad', 'OLeary', 'NHigham', ... 'Boman', 'Bjorck', 'Eisenstat', 'Zha', 'VanHuffel', 'Park', 'Arioli', ... 'MuntheKaas', 'Ng', 'VanDooren', 'Liu', 'Smith', 'Duff', 'Henrici', ... 'Tang', 'Reichel', 'Luk', 'Hammarling', 'Szyld', 'Fischer', ... 'Stewart', 'Bunch', 'Gutknecht', 'Laub', 'Heath', 'Ipsen', ... 'Greenbaum', 'Ruhe', 'ATrefethen', 'Plemmons', 'Hansen', 'Elden', ... 'BunseGerstner', 'Gragg', 'Berry', 'Sameh', 'Ammar', 'Warner', ... 'Davis', 'Meyer', 'Nichols', 'Paige', 'Gill', 'Jessup', 'Mathias', ... 'Hochbruck', 'Starke', 'Demmel'}; end

Special thanks to Christine Tobler and Cosmin Ionita at MathWorks for help with today's post.

Get
the MATLAB code

Published with MATLAB® R2017a

The adjacency matrix of a hypercube demonstrates the new MATLAB graph object.... read more >>

]]>The adjacency matrix of a hypercube demonstrates the new MATLAB graph object.

My previous blog post was about the Patience Chinese Rings Puzzle. The puzzle's state lives in a six-dimensional space. There are 64 vertices, numbered 0 from to 63. Each vertex is connected to six other vertices, the ones whose numbers, expressed in binary, differ from its own by one bit. That connectivity arrangement -- those vertices and the connections between them -- form a *graph* known as a *hypercube*.

Recent versions of MATLAB include two new objects, `graph` and `digraph`. Here is an abbreviated `help` entry for `graph`.

>> help graph graph Undirected Graph G = graph(A) uses the square symmetric matrix A as an adjacency matrix and constructs a weighted graph with edges corresponding to the nonzero entries of A. If A is logical then no weights are added.

G = graph(A,NAMES) additionally uses the cell array of strings NAMES as the names of the nodes in G. NAMES must have as many elements as size(A,1).

And here is a function that generates the adjacency matrix for a hypercube.

```
type hypercube
```

function A = hypercube(d) % hypercube(d) is a logical adjacency matrix for a d-dimensional hypercube. pow2 = 2.^(0:d-1); f2b = @(j) floor(rem(j./pow2,2)); % flint to binary b2f = @(b) b*pow2'; % binary to flint n = 2^d; A = zeros(n,n,'logical'); for j = 0:n-1 % Convert column index to binary. b = f2b(j); % Flip bits to get corresponding row indices. for i = 1:d b(i) = ~b(i); k = b2f(b); b(i) = ~b(i); A(k+1,j+1) = 1; end end end

I first encountered hypercubes thirty years ago. For a couple of years in the late 1980s I was involved with one of the world's first commercial parallel computers, the IPSC, the *Intel Personal Super Computer*. Here is a link to my blog post in 2013.

The IPSC was constructed from single board computers, similar to the mother boards in the PCs of the day. We had three models -- *d5*, *d6* and *d7* -- with 32, 64 or 128 processors. Since it was impractical to connect every processor directly to every other processor, the connection network formed a hypercube. For a machine with $n = 2^d$ processors, that's $n \log_2{n}$ instead of $n^2$ connections. For the d6 that's 384 connections instead of 4096.

The d6 model has the same graph as the patience puzzle that I wrote about in my previous post. (Actually, the IPSC's graph is a directed graph, a *digraph*, because the connections are two-way. But undirected graphs have simpler plots.)

A 3-dimensional hypercube is not a conventional cube. The 3-d hypercube is only the vertices and edges of the cube. The cube itself includes its surfaces and interior. If you had a wire-frame cube made from rubber bands and squished it around without destroying the connectivity, you wouldn't have a cube any more, but you would still have a model of a hypercube.

Here is the adjacency matrix of a 3-d hypercube. (It's so small that it's not necessary to make use of the sparsity.)

A = hypercube(3)

A = 8×8 logical array 0 1 1 0 1 0 0 0 1 0 0 1 0 1 0 0 1 0 0 1 0 0 1 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 1 0 0 0 1 0 1 1 0

One view of this matrix is its `spy` plot.

spy(A)

To get other views, let's use the new MATLAB graph object. We'll label the nodes with '0' through '7'.

nodes = num2cellstr(0:7) G = graph(A,nodes)

nodes = 1×8 cell array '0' '1' '2' '3' '4' '5' '6' '7' G = graph with properties: Edges: [12×1 table] Nodes: [8×1 table]

By itself, a graph doesn't have any geometric structure. In order to plot it you specify a *layout*, or coordinates, for the nodes. Coming up with good layouts is an art, as well as a science. Here are the currently available layouts.

'auto' (Default) Automatic choice of layout method based on the structure of the graph. 'circle' Circular layout. 'force' Force-directed layout. Uses attractive and repulsive forces on nodes. 'layered' Layered layout. Places nodes in a set of layers. 'subspace' Subspace embedding layout. Uses projection onto an embedded subspace. 'force3' 3-D force-directed layout. 'subspace3' 3-D Subspace embedding layout.

The default layout of the 3-d hypercube graph is a perspective view of an ordinary cube.

```
plot(G)
axis square
```

Let's move up to four dimensions.

A = hypercube(4); spy(A)

I think the default layout of the 4-d hypercube is attractive.

```
nodes = num2cellstr(0:15)';
G = graph(A,nodes);
plot(G)
axis square
```

But `'circle'` is another attractive layout. Notice that '3' is not connected to '4' because 0011 and 0100 differ by more than a single bit. The blue dot in the center is not a node -- it is just the point where many edges cross.

plot(G,'layout','circle') axis square

Here is yet another view.

plot(G,'layout','layered') axis square

I'm still not satisfied. I want one cube inside of another. I can supply my own coordinates for the vertices.

```
type cubelet
```

function x = cubelet % x = cubelet % 8-by-3 array of vertices of a unit cube. x = [-1 -1 -1 1 -1 -1 -1 1 -1 1 1 -1 -1 -1 1 1 -1 1 -1 1 1 1 1 1]; end

x = cubelet; x = [3*x; x]; plot(G,'xdata',x(:,1),'ydata',x(:,2),'zdata',x(:,3)) axis equal off view(-20,11)

This served as a logo for the IPSC.

The Laplacian of a graph is a matrix that has the degree of the nodes on the diagonal and -1's off the diagonal marking the edges. For the 3-d hypercube the degrees are all 4. So here the Laplacian.

L = full(laplacian(G)); frmt = [' ' repmat('%3d',1,16) '\n']; fprintf('L = \n') fprintf(frmt,L)

L = 4 -1 -1 0 -1 0 0 0 -1 0 0 0 0 0 0 0 -1 4 0 -1 0 -1 0 0 0 -1 0 0 0 0 0 0 -1 0 4 -1 0 0 -1 0 0 0 -1 0 0 0 0 0 0 -1 -1 4 0 0 0 -1 0 0 0 -1 0 0 0 0 -1 0 0 0 4 -1 -1 0 0 0 0 0 -1 0 0 0 0 -1 0 0 -1 4 0 -1 0 0 0 0 0 -1 0 0 0 0 -1 0 -1 0 4 -1 0 0 0 0 0 0 -1 0 0 0 0 -1 0 -1 -1 4 0 0 0 0 0 0 0 -1 -1 0 0 0 0 0 0 0 4 -1 -1 0 -1 0 0 0 0 -1 0 0 0 0 0 0 -1 4 0 -1 0 -1 0 0 0 0 -1 0 0 0 0 0 -1 0 4 -1 0 0 -1 0 0 0 0 -1 0 0 0 0 0 -1 -1 4 0 0 0 -1 0 0 0 0 -1 0 0 0 -1 0 0 0 4 -1 -1 0 0 0 0 0 0 -1 0 0 0 -1 0 0 -1 4 0 -1 0 0 0 0 0 0 -1 0 0 0 -1 0 -1 0 4 -1 0 0 0 0 0 0 0 -1 0 0 0 -1 0 -1 -1 4

The Laplacian could also be obtained from the original adjacency matrix.

L = diag(sum(A)) - A;

The coordinates of possible layouts for the plot of the graph can be obtained by picking three of the eigenvectors of the Laplacian. Here are all of the eigenvalues and three of the eigenvectors. (The eigenvalues are all integers.)

```
[V,E] = eig(L);
evals = round(diag(E)');
fprintf('evals = \n')
fprintf(frmt,evals)
evecs3 = V(:,1:3)
```

evals = 0 2 2 2 2 4 4 4 4 4 4 6 6 6 6 8 evecs3 = -0.2500 0.0216 0.0025 -0.2500 0.1620 0.4067 -0.2500 0.0245 -0.1449 -0.2500 0.1649 0.2593 -0.2500 0.2546 -0.2521 -0.2500 0.3950 0.1521 -0.2500 0.2575 -0.3996 -0.2500 0.3979 0.0046 -0.2500 -0.3979 -0.0046 -0.2500 -0.2575 0.3996 -0.2500 -0.3950 -0.1521 -0.2500 -0.2546 0.2521 -0.2500 -0.1649 -0.2593 -0.2500 -0.0245 0.1449 -0.2500 -0.1620 -0.4067 -0.2500 -0.0216 -0.0025

And here is the resulting plot.

G = graph(A,nodes); plot(G,'xdata',V(:,1),'ydata',V(:,2),'zdata',V(:,3)) axis square off set(gca,'clipping','off') zoom(2) view(-75,45)

The state space of both the puzzle that got me started on this blog post and one of the models of the parallel computer that consumed me thirty years ago is six dimensional.

A = hypercube(6); spy(A)

A plot of a cube inside of a cube inside of a cube inside of a cube is a bit messy. The circle layout is an interesting design, but it is hard to see which nodes are connected to which.

nodes = num2cellstr(0:63)'; G = graph(A,nodes); plot(G,'layout','circle') axis square

Here is a three-dimensional layout.

plot(G,'layout','subspace3') axis off square set(gca,'clipping','off') zoom(2) view(-108,-21)

The Bucky Ball provides an elegant example of a graph. Here's the `help` entry.

```
help bucky
```

BUCKY Connectivity graph of the Buckminster Fuller geodesic dome. B = BUCKY is the 60-by-60 sparse adjacency matrix of the connectivity graph of the geodesic dome, the soccer ball, and the carbon-60 molecule. [B,V] = BUCKY also returns xyz coordinates of the vertices.

When we put `bucky` into MATLAB years ago we went to a whole lot of trouble to compute those vertex coordinates. (To see how much trouble, `type bucky`.) But now, with the three dimensional `force3` layout, the `graph/plot` function does that work. Putting a little spin on it makes it interesting to watch.

```
type bucky_gif
```

B = graph(bucky); plot(B,'layout','force3') axis square off vis3d set(gca,'clipping','off') zoom(2) [a,e] = view; gif_frame('bucky_spin.gif') for d = 10:10:360 view(a+d,e+d) gif_frame end gif_frame('reset')

It's more fun to use the `cameratoolbar` and put the spin on the ball yourself.

Get
the MATLAB code

Published with MATLAB® R2017a

MIT's Professor Daniel Frey recently introduced me to an ancient mechanical puzzle known as "Chinese Rings", "Patience", or "Baguenaudier." I have modified Daniel's simulator to produce a new app. The state space of the puzzle forms a hypercube.... read more >>

]]>MIT's Professor Daniel Frey recently introduced me to an ancient mechanical puzzle known as "Chinese Rings", "Patience", or "Baguenaudier." I have modified Daniel's simulator to produce a new app. The state space of the puzzle forms a hypercube.

There are dozens of YouTube videos showing solutions and speed contests for the tavern puzzle called "Patience". The best video is by Ty Yann, "Baguenaudier Chinese Rings." Here is the link to his puzzle tutorial.

Wikipedia tell us that *Baguenaudier* is French for "time-waster".

Prof. Frey has written a nice MATLAB simulator for his course this spring at MIT, 2.086, "Numerical Computation for Mechanical Engineers." My new app is based on his program.

The puzzle has six steel hooks, often made from nails, attached to a base. A steel ring is attached to each hook. A long movable piece known as the shuttle is initially threaded around the nails and through the rings. The objective is to free the shuttle. It ain't easy.

Image credit: Daniel Frey, MIT

My new interactive app, `patience`, is included with version 2.2 of Cleve's Laboratory in MATLAB Central. I hope you will try to solve the puzzle yourself before you resort to the "no patience" automated solve button. The following animated gif shows several steps at the start of the solution, then skips many intermediate steps before showing the last several.

The position of the six hooks, up or down, provides a six-bit binary number that characterizes the state of the mechanism. Initially, the rings are all engaged, so the state in binary is `000000`, which, of course, is a decimal 0. The objective is to free all the rings, giving a binary `111111`, or decimal 63. The app shows the decimal value at each step in the figure title.

The hooks are numbered from 1 through 6 starting on the left. So, we will use "big endian" binary, with the least significant bit at the left.

This one line in `patience.m`, which refers to the state vector `S` and a hook number `n`, controls the possible moves.

dbtype 52:53 patience

52 % Is the move allowed ? 53 ok = (n==1) || (n==2)&&~S(1) || (n>2)&&(~S(n-1)&&all(S(1:n-2)));

The rule says

- Hook 1 can always be toggled.
- Hook 2 can be toggled if hook 1 is up.
- Any hook numbered 3 through 6 can be toggled if the hook immediately to its left is up and all the others before it are down.

That's all you need to know to operate the puzzle.

It turns out that at any step at most two, and sometimes only one, of the six hooks is a legal move.

For example, if the state is

S = [1 1 0 0 1 0]

there are two possible moves. You can click on hook 1 to give

S = [0 1 0 0 1 0]

Or, you can click on hook 4 to give

S = [1 1 0 1 1 0]

The optimum path from a given starting state, such as 0, to the goal, 63, can be found by a brute force recursive search that tries all possible moves at each step. The code is at the end of this post. Here is the shortest path from 0 to 63.

path = shortest_path(0)

path = Columns 1 through 13 0 2 3 11 10 8 9 13 12 14 15 47 46 Columns 14 through 26 44 45 41 40 42 43 35 34 32 33 37 36 38 Columns 27 through 39 39 55 54 52 53 49 48 50 51 59 58 56 57 Columns 40 through 43 61 60 62 63

Notice that a shortest path never repeats any state because if we find ourselves reaching a state that we've already visited, we could shorten that path by skipping the extra steps.

It turns out that the only starting state that visits all vertices is `31 = [1 1 1 1 1 0]`, which corresponds to starting with the sixth hook up and all the others down.

Each step changes the state by a power of two. You can recover the sequence of mouse clicks that would generate this path by taking the log2 of the changes.

hooks = log2(abs(diff(path))) + 1

hooks = Columns 1 through 13 2 1 4 1 2 1 3 1 2 1 6 1 2 Columns 14 through 26 1 3 1 2 1 4 1 2 1 3 1 2 1 Columns 27 through 39 5 1 2 1 3 1 2 1 4 1 2 1 3 Columns 40 through 42 1 2 1

Let's plot the decimal value of the state for the shortest path from 0 to 63. Here is the progress of this value over the 42 steps. The even numbered steps involve toggling the first hook, which changes the value of the state by plus or minus one. The sixth hook is moved only once, at step 11, changing the value by +32. The fifth hook contributes its +16 at step 28.

plot_path

Here are the lengths of the shortest paths from each starting position.

L = zeros(64,1); for n = 0:63 L(n+1) = length(shortest_path(n)); end bar(L)

The puzzle's state lives in a six-dimensional space. There are 64 vertices, numbered 0 from to 63. Each vertex is connected to the six vertices whose numbers, expressed in binary, differ from its own by one bit. That connectivity arrangement -- those vertices and the connections between them -- form a *graph* known as a *hypercube*.

Recent versions of MATLAB include methods for a new object, the *graph*. I plan to discuss hypercubes and graphs in my next blog post.

Here is the code for the recursive search that finds the shortest path. The code for the complete `patience` app is now included with Cleve's Laboratory in MATLAB Central.

```
type shortest_path
```

function path = shortest_path(state,m,path) % shortest_path, or shortest_path(state), where state is a decimal state % value between 0 and 63, is the shortest path of allowable puzzle moves % from that state to the objective, 63. % % shortest_path(S,m,path) is the recursive call. if nargin == 0 state = 0; % Default end if nargin <= 1 % Initial call m = 0; % Previous n path = []; end if state == 63 % Terminate recursion path = state; return end if state == 31 m = 0; end % Search for shortest path pow2 = 2.^(0:5); lmin = inf; for n = [1:m-1 m+1:6] Sn = mod(floor(state./pow2),2); % Convert state to binary Sn(n) = ~Sn(n); % Flip one bit ok = (n==1) || (n==2)&&~Sn(1) || ... (n>2)&&(~Sn(n-1)&&all(Sn(1:n-2))); % Allowable move if ok staten = Sn*pow2'; % Recursive call pn = shortest_path(staten,n,[path staten]); if length(pn) < lmin lmin = length(pn); pbest = pn; end end end path = [state pbest]; end % shortest_path

Thanks to Daniel Frey.

Get
the MATLAB code

Published with MATLAB® R2017a

"ULP" stands for "unit in the last place." An *ulps plot* samples a fundamental math function such as $\sin{x}$, or a more esoteric function like a Bessel function. The samples are compared with more accurate values obtained from a higher precision computation. A plot of the accuracy, measured in ulps, reveals valuable information about the underlying algorithms.... read more >>

"ULP" stands for "unit in the last place." An *ulps plot* samples a fundamental math function such as $\sin{x}$, or a more esoteric function like a Bessel function. The samples are compared with more accurate values obtained from a higher precision computation. A plot of the accuracy, measured in ulps, reveals valuable information about the underlying algorithms.

`libm` is the library of elementary math functions that supports the C compiler. `fdlibm` is "freely distributable" source for `libm` developed and put into the public domain over 25 years ago by K. C. Ng and perhaps a few others at Sun Microsystems. I wrote about `fdlibm` in our newsletter in 2002, The Tetragamma Function and Numerical Craftsmanship.

Mathematically `fdlibm` shows immaculate craftsmanship. We still use it today for our elementary transcendental functions. And I suspect all other mathematical software projects do as well. If they don't, they should.

`ulps(x)` is the distance from `x` to the next larger floating point number. It's the same as `eps(x)`.

To assess the accuracy of a computed value

y = f(x)

compare it with the more accurate value obtained from the Symbolic Math Toolbox

```
Y = f(sym(x,'f'))
```

The `'f'` flag says to convert `x` to a `sym` exactly, without trying to guess that it is an inverse power of 10 or the `sqrt` of something. The relative error in `y`, measured in units in the last place, is then

u = (y - Y)/eps(abs(Y))

Since this is *relative* error, it is a stringent measure near the zeros of `f(x)`.

If `y` is the floating point number obtained by correctly rounding `Y` to double precision, then

-0.5 <= u <= 0.5

This is the best that can be hoped for. Compute the exact mathematical value of `f(x)` and make a single rounding error to obtain the final result. Half-ulp accuracy is difficult to obtain algorithmically, and too expensive in execution time. All of the functions in MATLAB that are derived from `fdlibm` have better than one-ulp accuracy.

-1.0 < u < 1.0

Each of the following plots involves 100,000 random arguments `x`, uniformly distributed in an appropriate interval.

We see about 0.8 ulp accuracy from this sample. That's typical.

Argument reduction is the first step in computing $\sin{x}$. An integer multiple $n$ of $\pi/2$ is subtracted from the argument to bring it into the interval

$$ -\frac{\pi}{4} \le x - n \frac{\pi}{2} \le \frac{\pi}{4} $$

Then, depending upon whether $n$ is odd or even, a polynomial approximation of degree 13 to either $\sin$ or $\cos$ gives the nearly correctly rounded result for the reduced argument, and hence for the original argument. The ulps plot shows a slight degradation in accuracy at odd multiples of $\pi/4$, which are the extreme points for the polynomial approximations.

It is important to note that the accuracy is better than one ulp even near the end-points of the sample interval, $0$ and $2\pi$. This is where $\sin{x}$ approaches $0$ and the approximation must follow carefully so that the relative error remains bounded.

Again, roughly 0.8 ulp accuracy.

Similar argument reduction results in similar behavior near odd multiples of $\pi/4$. In between these points, at $\pi/2$ and $3\pi/2$, $\tan{x}$ has a pole and the approximation must follow suit. The algorithm uses a reciprocal and the identity

$$ \tan x = -1/\tan{(x+\frac{\pi}{2})} $$

This comes close to dividing by zero as you approach a pole, but the resulting approximation remains better than one ulp.

Good to within slightly more than 0.8 ulp. The underlying approximation is a piecewise polynomial with breakpoints at a few multiples of 1/16 that are evident in the plot and marked on the axis.

About the same accuracy as the previous plots.

The argument reduction involves the key value

$$ r = \ln{2} \approx 0.6931 $$

and the identity

$$ \exp{x} = 2^n \exp{(x-n r)} $$

The resulting ulps plot shows the extremes of the error at odd multiples of $r/2$.

Now for two functions that are not in `fdlibm`. If you follow this blog, you might have noticed that I am a big fan of the Lambert W function. I blogged about it a couple of years ago, The Lambert W Function. The Wikipedia article is excellent, Lambert W Function. And, you cn just Google "lambert w function" for many more interesting links.

The Lambert W function is not available in MATLAB itself. The Symbolic Math Toolbox has an @double method that accesses the symbolic code for type double arguments. Or, my blog post mentioned above has some simple (but elegant, if I do say so myself) code.

Neither code is one-ulp accurate. The primary branch of the function has a zero at the origin. As we get near that zero, the relative error measured in *ulps* is unbounded. The absolute accuracy is OK, but the relative accuracy is not. In fact, you might see a billion ulps error. That's a *gigaulp*, or *gulp* for short.

As `lambertw(x)` crosses a power of two, the unit in the last place, `eps(lambertw(x))`, jumps by a factor of two. Three of these points are marked by ticks on the x-axis in the ulps plot.

for a = [1/8 1/4 1/2] z = fzero(@(x) lambertw(x)-a,.5) end

z = 0.1416 z = 0.3210 z = 0.8244

Our codes for Bessel functions have fairly good, although not one-ulp, absolute accuracy. But they do not have high relative accuracy near the zeros. Here are the first five zeros, and an ulps plot covering the first two, for $J_0(x)$, the zero-th order Bessel function of the first kind.

for a = (1:5)*pi z = fzero(@(x) besselj(0,x), a) end

z = 2.4048 z = 5.5201 z = 8.6537 z = 11.7915 z = 14.9309

Here is the inverse of the error function. It looks very interesting, but I haven't investigated it yet.

I have recently updated Cleve's Laboratory in the MATLAB Central file exchange. The latest version includes `ulpsapp.m`, which generates the ulps plots in this blog.

Get
the MATLAB code

Published with MATLAB® R2017a

A headline in the *New York Times* at the end of 2016 said "Growth of U.S. Population Is at Slowest Pace Since 1937". This prompted me to revisit an old chestnut about fitting and extrapolating census data. In the process I have added a couple of nonlinear fits, namely the logistic curve and the double exponential Gompertz model.... read more >>

A headline in the *New York Times* at the end of 2016 said "Growth of U.S. Population Is at Slowest Pace Since 1937". This prompted me to revisit an old chestnut about fitting and extrapolating census data. In the process I have added a couple of nonlinear fits, namely the logistic curve and the double exponential Gompertz model.

This experiment is older than MATLAB. It started as an exercise in *Computer Methods for Mathematical Computations*, by Forsythe, Malcolm and Moler, published by Prentice-Hall 40 years ago. We were using Fortran back then. The data set has been updated every ten years since. Today, MATLAB graphics makes it easier to vary the parameters and see the results, but the underlying mathematical principles are unchanged:

- Using polynomials of even modest degree to predict the future by extrapolating data is a risky business.

Recall that the famous computational scientist Yogi Berra said

- "It's tough to make predictions, especially about the future."

The data are from the decennial census of the United States for the years 1900 to 2010. The population is in millions.

1900 75.995 1910 91.972 1920 105.711 1930 123.203 1940 131.669 1950 150.697 1960 179.323 1970 203.212 1980 226.505 1990 249.633 2000 281.422 2010 308.746

The task is to extrapolate beyond 2010. Let's see how an extrapolation of just six years to 2016 matches the Census Bureau announcement. Before you read any further, pause and make your own guess.

Here's is the opening screen of the January 2017 edition of my `censusapp`, which is included in Cleve's Laboratory. The plus and minus buttons change the extrapolation year in the title. If you go beyond 2030, the plot zooms out.

The pull-down menu offers these models. Forty years ago we had only polynomials.

models = {'census data','polynomial','pchip','spline', ... 'exponential','logistic','gompertz'}'

models = 7×1 cell array 'census data' 'polynomial' 'pchip' 'spline' 'exponential' 'logistic' 'gompertz'

The Census Bureau news release that prompted the story in *The Times* said the official population in 2016 was 323.1 million. That was on Census Day, which is April 1 of each year. The Census Bureau also provides a dynamic Population Clock that operates continuously. But let's stick with the 323.1 number.

Polynomials like to wiggle. Constrained to match data in a particular interval, they go crazy outside that interval. Today, there are 12 data points. The app lets you vary the polynomial degree between 0 and 11. Polynomials with degree less than 11 approximate the data in a least squares sense. The polynomial of degree 11 interpolates the data exactly. As the degree is increased, the approximation of the data becomes more accurate, but the behavior beyond 2010 (or before 1900) becomes more violent. Here are degrees 2 and 7, 9, 11, superimposed on one plot.

The quadratic fit is the best behaved. When evaluated at year 2016, it misses the target by six million. Of course, there is no reason to believe that the US population grows like a second degree polynomial in time.

The interpolating polynomial of degree 11 tries to escape even before it gets to 2010, and it goes negative late in 2014.

MATLAB has two piecewise cubic interpolating polynomials. The classic `spline` is smooth because it has two continuous derivatives. Its competitor `pchip` sacrifices a continuous second derivate to preserve shape and avoid overshoots. I blogged about splines and pchips a few years ago.

Neither is intended for extrapolation, but we will do it anyway. Their behavior beyond the interval is determined by their end conditions. The classic `spline` uses the so-called *not-a-knot* condition. It is actually a single cubic in the last two subintervals. That cubic is also used for extrapolation beyond the endpoint. `pchip` uses just the last three data points to create a different shape-preserving cubic for use in the last subinterval and beyond.

Let's zoom in on the two. Both are predicting a decreasing rate of growth beyond 2010, just as the Census Bureau is observing. `pchip` gets lucky and comes within 0.2 million of the announcement for 2016.

As I said, there is no good reason to model population growth by a polynomial, piecewise or not. But because the rate of growth can be expected to be proportional to the size of the population, there is good reason to use an exponential.

$$ p(t) \approx a e^{bt} $$

There have been many proposals for ways to modify this model to avoid its unbounded growth. I have just added two of these to `censusapp`. One is the logistic model.

$$ p(t) \approx \frac{a}{1+b e^{-ct}} $$

And the other is the Gompertz double exponential model, named after Benjamin Gompertz, a 19th century self-educated British mathematician and astronomer.

$$ p(t) \approx a e^{-b e^{-ct}} $$

In both of these models the growth is limited because the approximating term approaches $a$ as $t$ approaches infinity.

In all three of the exponential models, the parameters $a$, $b$, and possibly $c$, appear nonlinearly. In principle, I could use `lsqcurvefit` to search in two or three dimensions for a least squares fit to the census data. But I have an alternative. By taking one or two logarithms, I have a *separable least squares* model where at most one parameter, $a$, appears nonlinearly.

For the exponential model, take one logarithm.

$$ \log{p} \approx \log{a} + bt $$

Fit the logarithm of the data by a straight line and then exponentiate the result. No search is required.

For the logistic model, take one logarithm.

$$ \log{(a/p-1)} \approx \log{b} - ct $$

For any value of $a$, the parameters $\log{b}$ and $c$ appear linearly and can be found without a search. So use a one-dimensional minimizer to search for $a$. I could use `fminbnd`, or its textbook version, `fmintx`, from *Numerical Methods with MATLAB*.

For the Gompertz model, take two logarithms.

$$ \log{\log{a/p}} \approx \log{b} - ct $$

Again, do a one-dimensional search for the minimizing $a$, solving for $\log{b}$ and $c$ at each step.

Here are the three resulting fits, extrapolated over more than 200 years to the year 2250. The pure exponential model reaches 5 billion people by that time, and is growing ever faster. I think that's unreasonable.

The value of $a$ in the Gompertz fit turns out to be 4309.6, so the population will be capped at 4.3 billion. But it has only reached 1.5 billion two hundred years from now. Again unlikely.

The value of $a$ in the logistic fit turns out to be 756.4, so the predicted US population will slightly more than double over the next two hundred years. Despite the Census Bureau's observation that our rate of growth has slowed recently, we are not yet even half the way to our ultimate population limit. I'll let you be the judge of that prediction.

I have recently updated Cleve's Laboratory in the MATLAB Central file exchange. One of the updates changes the name of `censusgui` to `censusapp` and adds the two exponential models. If you do install this new version of the Laboratory, you can answer the following question.

The fit generated by `pchip` defines a cubic for use beyond the year 2000 that predicts the population will reach a maximum in the not too distant future and decrease after that. What is that maximum and when does it happen?

Get
the MATLAB code

Published with MATLAB® R2017a

I don't recall where I found this seasonal fractal. And I can't explain how it works. So please submit a comment if you can shed any light on either of these questions.... read more >>

]]>I don't recall where I found this seasonal fractal. And I can't explain how it works. So please submit a comment if you can shed any light on either of these questions.

I have featured a fractal at this time of the year in two previous years. I see Christmas trees or perhaps a holly decoration.

All these figures are obtained by varying the parameters in a formula that generates complex numbers $z$ from the partial sums.

$$ z = \sum_n{\exp{(\phi n^p+\sigma) \pi i}} $$

A vector of points in the complex plane is produced by taking $n$ to be a vector of consecutive integers and using the MATLAB cumulative summation function `cumsum` to compute the partial sums. There are 8600 points in the figure above and 100,000 points in the figures below.

The default value of the parameter $\phi$ is my old friend the `golden ratio`.

$$ \phi = \frac{1+\sqrt{5}}{2} $$

In previous posts I've taken $\phi$ to be other rational and irrational numbers, but today I am sticking to this value.

The parameter $\sigma$ controls the angular orientation. Taking $\sigma$ near $1/8$ makes the large Christmas tree vertical.

While trying how to understand how this thing works I've varied the power $p$ from its usual value of 2 and taken hundreds of thousands of points. This produces today's pictures. Different values of $p$ produce wildly different results.

For a real variable $x$, the expression $\exp (x \pi i)$ is periodic and lies on the unit circle in the complex plane. So we are plotting the cumulative sum of values taken from around the unit circle. At first glance, this appears to be a complex valued random number generator. But it is a lousy generator because we can see Christmas trees in the output.

```
type greetings_gifs
```

function greetings_gifs % Generate animated season's greeting gifs. % Generate the fractal phi = (1+sqrt(5))/2; s = 1/8; n = 100000; for gifnum = 1:4 switch gifnum case 1, p = 2; case 2, p = 2/3; case 3, p = 1.25; case 4, p = 4; end w = exp((phi*(0:n).^p+s)*pi*1i); z = cumsum(w); % Find local extrema ks = extrema(z); % Make an animated gif plotit(z,p,ks,gifnum) end % gifnum % ------------------------ function ks = extrema(z) n = length(z)-1; m = n/40; ks = []; for j = 0:m:n-m zj = z(j+1:j+m); w = zj - mean(zj); k = find(abs(w) == max(abs(w))) + j; ks = [ks k]; end end % extrema function plotit(z,p,ks,gifnum) % Make an animated gif shg plot(z,'.') axis square ax = axis; gif_frame(['greetings' int2str(gifnum) '.gif']) clf axis(ax) axis square gif_frame for j = 1:length(ks) k = ks(j); plot(z(1:k),'k.','markersize',0.5) axis(ax) axis square hold on plot(z(ks(1:j)),'g.','markersize',18) plot(z(k),'r.','markersize',24) hold off title(sprintf('p = %4.2f',p)) gif_frame end gif_frame(5) gif_frame('reset') end % plotit end % greetings_gifs

*Happy New Year.*

Get
the MATLAB code

Published with MATLAB® R2016a

When I was in high school in the 1950's, I didn't know anything about matrices. But I nearly encountered one when I wrote a paper for my physics class about the color scheme that allowed new color TV broadcasts to be compatible with existing black-and-white TV receivers.... read more >>

]]>When I was in high school in the 1950's, I didn't know anything about matrices. But I nearly encountered one when I wrote a paper for my physics class about the color scheme that allowed new color TV broadcasts to be compatible with existing black-and-white TV receivers.

During my high school years in Utah, I worked part-time and summers as a camera man and projectionist at what was then called KTVT, Channel 4, and the local NBC affiliate in Salt Lake City. The TV cameras at the time employed vacuum tubes, not transistors. They were the size of suit cases, weighed over 100 pounds, and were mounted on dollies.

The first coast-to-coast TV broadcast in color was the Tournament of Roses parade on January 1, 1954. KTVT had a prototype color receiver, perhaps the only one in Utah at the time. We set up bleachers in the studio and invited guests to see the new technology.

The NTSC compatible color television standard, adopted in 1953, allows color information to be encoded in a broadcast in such a way that existing black-and-white receivers will still function without modification. The standard remained in use in the US until broadcasters ceased analog transmissions in 2009.

I wrote a paper for my high school physics class about the NTSC compatible color scheme. It involves a *luminance-chrominance* conversion between the RGB, red-green-blue, signal used in color studio equipment and home receivers, and a YIQ encoding. Y represents the brightness, or *luma*, of the signal, and is used by itself for black-and-white. The letters I and Q stand for *in-phase* and *quadrature*. These two components combine with the brightness to produce color.

Here are the conversion formulas that I had in my high school paper.

$$ Y = 0.2989 R + 0.5870 G + 0.1140 B $$

$$ I = 0.5959 R - 0.2744 G - 0.3216 B $$

$$ Q = 0.2115 R - 0.5229 G + 0.3114 B $$

I now know how to write the conversion formulas in matrix notion. If $x$ is the 3-vector containing RGB and $y$ is the 3-vector of YIQ, then

$$ y = A x $$

where $A$ is the matrix

$$ A = \left[ \begin{array}{rrr} 0.2989 & 0.5870 & 0.1140 \\ 0.5959 & -0.2744 & -0.3216 \\ 0.2115 & -0.5229 & 0.3114 \end{array} \right] $$

In retrospect, this was my very first matrix.

Of course, conversion at the receiving end from YIQ to RGB is done by

$$ x = A^{-1} y $$

Fast forward 25 years so we can use MATLAB.

A = [0.2989 0.5870 0.1140 0.5959 -0.2744 -0.3216 0.2115 -0.5229 0.3114] T = inv(A)

A = 0.2989 0.5870 0.1140 0.5959 -0.2744 -0.3216 0.2115 -0.5229 0.3114 T = 1.0002 0.9560 0.6211 1.0001 -0.2720 -0.6470 1.0000 -1.1060 1.7030

The first column of $A^{-1}$ is all 1's, so each of R, G and B gets a full dose of Y.

Although I had nothing to do with putting them there, these two matrices are central to the functions `ntsc2rgb` and `rgb2ntsc` in the `colorspaces` folder of the Image Processing Toolbox. In the following code `A` is a 3D full color RGB image.

Here is the core of the NTSC to RGB function. The image is reshaped to a tall, skinny 2D array and the transformation applied on the right by the transpose of `T`.

dbtype 36:36 ntsc2rgb dbtype 41:43 ntsc2rgb dbtype 48:48 ntsc2rgb

36 T = [1.0 0.956 0.621; 1.0 -0.272 -0.647; 1.0 -1.106 1.703]; 41 m = size(A,1); 42 n = size(A,2); 43 A = reshape(A(:),m*n,3)*T'; 48 A = reshape(A,m,n,3);

The conversion from RGB to NTSC uses *matrix division* on the right by `T` transpose. This may be the only place in MATLAB where the forward slash matrix operator appears.

dbtype 30:31 rgb2ntsc dbtype 35:35 rgb2ntsc

30 T = [1.0 0.956 0.621; 1.0 -0.272 -0.647; 1.0 -1.106 1.703].'; 31 [so(1),so(2),thirdD] = size(A); 35 A = reshape(reshape(A,so(1)*so(2),thirdD)/T,so(1),so(2),thirdD);

Since analog NTSC broadcasts are history, these conversions are probably not very common any more. But the single formula

$$ Y = 0.2989 R + 0.5870 G + 0.1140 B $$

is the basis for the `rgb2gray` function.

NBC introduced a peacock logo to promote the debut of its "in living color" broadcasts. The peacock has lost a few feathers over the years, but it still provides an excellent color test pattern.

(NBC SVG by FOX 52 [Public domain], via Wikimedia Commons)

Here are the red, green and blue components. The bright portions represent high values. All three have relatively bright backgrounds and bright peacocks in the center; these combine to produce white. The R component has bright feathers on the left to contribute to the yellow, orange and red feathers in the full logo. The G component has bright components on the bottom to contribute to the yellow and green in the full. And B has bright feathers on the upper right to give blue and purple.

Convert to YIQ. The Y component is the grayscale rendition to drive old black-and-white TVs. The I component goes from orange to blue and the Q component goes from purple to green.

Talking about color, here is one of my favorite graphics. Let's start with an animated .gif of `colorcubes(3)`, which contains 27 cubelets. If you run `colorcubes` yourself, you will drive the rotation with your mouse.

You can see white at one vertex and black at the opposite vertex. Red, a bright green, and blue are at three vertices, while their complements -- cyan, magenta, and yellow -- oppose them. Actually, the green vertex is too bright. It is often called "lime". A darker green is along the bottom edge, halfway between lime and black. Here's a little quiz: What color is the hidden cubelet in the center?

Now a look at the construction of `colorcubes(5)`.

Here is the `colorcubes` code.

```
type colorcubes
```

function colorcubes(n,w) % COLORCUBES A cube of cubes in the RGB color space. % COLORCUBES, with no arguments, shows 5^3 = 125 cubes with % colors equally spaced in the RGB color space. % COLORCUBES(n) shows n-by-n-by-n colors. % COLORCUBES(2) shows 8 colors: R, G, B, C, M, Y, W, K (black). % COLORCUBES(n,w) uses cubes of width w. Default is w = 0.85. % Rotate the cube with the mouse or arrow keys. % Copyright 2016 The MathWorks, Inc. if nargin < 1, n = 5; end if nargin < 2, w = 0.85; end initgraphics(n) [x,y,z] = cube(w); m = n-1; for i = m:-1:0 for j = m:-1:0 for k = 0:m r = k/m; g = 1-j/m; b = 1-i/m; surface(i+x,j+y,k+z, ... 'facecolor',[r g b], ... 'facelighting','gouraud'); drawnow end %k end %j end %i % ------------------------ % INITGRAPHCS Inialize the colorcubes axis. % INITGRAPHICS(n) for n-by-n-by-n display. function initgraphics(n) clf reset shg set(gcf,'color','white') axis([0 n 0 n 0 n]); axis off axis vis3d rotate3d on end %initgraphics function [x,y,z] = cube(w) % CUBE Coordinates of the faces of a cube. % [x,y,z] = cube(w); surface(x,y,z) % plots a cube of with w. u = [0 0; 0 0; w w; w w]; v = [0 w; 0 w; 0 w; 0 w]; z = [w w; 0 0; 0 0; w w]; s = [nan nan]; x = [u; s; v]; y = [v; s; u]; z = [z; s; w-z]; end %cube end % colorcubes

Get
the MATLAB code

Published with MATLAB® R2016a