I could carry less

31 August 2010

The fabled carefree residents of the Carryless Islands in the remote South Pacific have very few possessions, which is just as well, since their notion of arithmetic is ill-suited to accurate record-keeping. When they add or multiply numbers, they follow similar rules to ours, except that there are no carries into other digit positions. Addition and multiplication of single-digit numbers are performed by a process that we would call “reduction mod 10.” Any carry digits are simply ignored. So 9 + 4 = 3, 5 + 5 = 0, 9 × 4 = 6, 5 × 4 = 0, and so on.

With this fable, David Applegate, Marc LeBrun and N. J. A. Sloane introduce a new scheme of arithmetic in a paper newly posted on the arXiv.

And if you think the mathematics sounds trivial, try explaining the structure of this sequence:

21, 23, 25, 27, 29, 41, 43, 45, 47, 49, 51, 52, 53, 54, 56, 57, 58, 59, 61, 63,

which lists the first 20 “carryless primes.”

In the zone

24 August 2010

Mt. Shasta, from a nearby hilltop owned by Coca Cola

Before leaving on a trip to the West Coast, I copied my return flight information onto my Google calendar: SFO to BOS, 12:50 p.m. to 9:30 p.m. Now that I’m nearing the end of my visit to the Mythical State of Jefferson (see local landmark above), I’ve just checked the departure details by calling up the calendar on my cell phone. It tells me the flight departs at 9:50 a.m. and arrives at 6:30 p.m.

It could be worse, of course. As an eastbound traveler all that I risk is wasting three hours at the airport. If I were westbound, I might well miss my flight.

The developers of Google calendar would doubtless argue that automatic time-zone conversion is a feature, not a bug. If the event in question had been a conference call scheduled for 12:50 p.m. Eastern Daylight Time, then 9:50 a.m. Pacific Daylight Time would indeed be the moment to dial in. Or, if I had a series of pills to be taken at fixed intervals, it might be helpful to have the program remind me at the correct times as I wander across continents. But in the case of today’s flight home, the result is just plain wrong.

It’s not only a Google calendar problem. A few weeks ago a friend was entering a schedule of talks into a Drupal web page for a conference that begins November 7, 2010. All the times were mysteriously shifted by an hour. A 1:00 p.m. talk on the input form became a 2:00 p.m. talk on the displayed web page. The key to solving this mystery is knowing that November 7 is the date daylight saving time ends in (most of) the U.S.

I am certainly not the first to encounter such problems. Peter Neumann’s RISKS Digest reports hundreds of computational mishaps involving time zones or daylight saving time, going back over the past 25 years. The issue is known to Google. But what is the right fix?

Some other calendar software offers the option of specifying a time zone for an event. To handle the airline case correctly, the program needs to allow for different zones for the start and the end times. And the Drupal problem suggests we may also need some means of indicating whether or not to adjust for daylight saving time. It gets very messy. Note that an event scheduled for 1:30 a.m. on November 7, 2010, will happen twice. An event at 2:30 a.m. on March 13, 2011, will never take place.

For years my own makeshift solution to these complexities was to live in a single time zone, no matter where I was. I carried a laptop, but I never changed its time-zone setting, even when I was away from home for weeks or months. When conversion was needed, it happened in my head. Even now the Google calendar display on my laptop gives the correct departure time from SFO, because the laptop doesn’t know that it ever left Boston. But in our new world of location-aware devices, pretending to stay home is no longer an option.

I begin to wonder if the whole railroad-age concept of time zones hasn’t outlived its usefulness. But I haven’t time just now to consider the alternatives. Right now it’s time to leave for the airport. I think.

The ormat game

16 August 2010

Here’s the deal. I’m going to give you a square grid, with some of the cells colored and others possibly left blank. We’ll call this a template. Perhaps the grid will be one of these 3×3 templates:

colored 3x3 ormat grids

You have a supply of transparent plastic overlays that match the grid in size and shape and that also bear patterns of black dots:

dot patterns for the six 3x3 permutation matrices

Note that each of these patterns has exactly three dots, with one dot in each row and each column. The six overlays shown are the only 3×3 grids that have this property.

Your task is to assemble a subset of the overlays and lay them on the template in such a way that dots cover all the colored squares but none of the blank squares. You are welcome to superimpose multiple dots on any colored square, but overall you want to use as few overlays as possible. To make things interesting, I’ll suggest a wager. I’ll pay you $3 for a correct covering of a 3×3 template, but you have to pay me $1 for each overlay you use. Is this a good bet?

Before going further, I should mention that not every conceivable template can be covered under these rules. To take an obvious example, no 3×3 template with fewer than three colored squares can possibly be covered by any combination of the six overlays. But I promise to submit only templates that can be covered by some combination of the given dot patterns; if I err about this, I forfeit the bet.

How does the game play out? If I give you the template marked “1″ above, you can easily win; just choose permutations a and b, which together cover all the colored squares and no others. You pay $2 and get $3. Template 2, with all nine squares colored, looks like it might be the toughest challenge. Clearly, it cannot be covered with fewer than three overlays, since we need a total of nine dots; and it turns out that exactly three overlays are required. Indeed, there are two ways of covering the template with three overlays: a + d + e and b + c + f. Thus this template is a breakeven proposition: You earn $3 and pay $3.

Now we come to template 3, which has eight colored squares and one blank. Surely if you can cover the full nine squares with just three overlays, then you should also be able to cover eight squares—no? I invite you to try it. In fact the only covering that works requires four overlays: b + d + e + f. Thus you shouldn’t take my bet, since I can always give you a template with just one blank, and you’ll have a net loss of $1.

Some background. I’ll return to the gaming table momentarily, but first let me explain what this is all about and where it came from. A few weeks ago, I was writing about “ranges of rankings,” which led me into the topic of permutation matrices. To recapitulate:

  • A permutation matrix is a square matrix with a single 1 in each column and each row, and all the rest of the elements 0.
  • An ormat is a superposition of permutation matrices, formed by applying the Boolean OR function to corresponding elements of the permutation matrices. For example: matrix-or-sum.png
  • Not all square matrices with (0,1) entries can be formed by OR-ing permutation matrices, but there’s an efficient algorithm for deciding whether or not a given matrix is an ormat. (I thank some helpful commenters for enlightening me on this point.)
  • Given an ormat, the total number of distinct permutation paths that can be threaded through the 1 entries of the matrix is equal to the permanent of the matrix. Calculating the permanent is known to be a hard computational problem.

In a comment, Barry Cipra posed the following query:

The permanent tells us the maximum number of different permutations that can be OR-summed to produce a given ormat, but what is the corresponding minimum number? Also, in how many different ways can the minimum be achieved?

The connection between ormats and my little game is probably apparent by now. The template of colored and blank squares is an ormat; the dotted overlays represent permutation matrices; to maximize your payoff in the game (or to minimize your loss), you need to answer Barry’s first question, finding the minimum number of permutations that can be combined to yield the given ormat.

For 3×3 matrices, we can solve this problem by exhaustive search, calculating the OR-sums of all possible combinations of the six 3×3 permutation matrices taken 1, 2, 3, …, 6 at a time. I did this with pencil and paper on a recent airplane trip. Here is a summary of the results:

number of ormats generated by various combinations of permutation matrices

Some of the numbers on this card are easy to explain. The six ormats with just three 1 entries are the permutation matrices themselves. There are six of them because there are 3! = 6 permutations of three things. There are no ormats with four 1 entries for a reason that bears thinking about: There can be no permutations that differ from one another in just one element. When you superimpose any of the six overlays shown above, you can wind up with three, five or six dots, but never four.

At the other end of the scale, it’s no surprise that there’s exactly one ormat with nine 1 entries, and that it takes three permutations to produce it. And then there are the nine ormats with eight 1 entries, which each require four permutations to be OR-ed. These are the single-blank patterns like template 3 above.

Based on these results, I began speculating about what I would see in a tabulation of all 4×4 ormats.

guesses about stats for 4x4 ormats

There would have to be 4! = 24 patterns with four 1 entries, and just one pattern with all 1s, generated by OR-ing four permutations. And there should be 16 ormats that require five permutations, namely the 16 matrices with a single 0 element. This last prediction seemed a little less self-evident than the others.

Pocket change and Cheerios. My thoughts about the single-zero (or single-blank) case went something like this. To cover 15 squares with sets of four dots each, we need at least four sets, or else we simply won’t have enough dots. So a useful starting point is one of the optimal arrangements that cover all 16 squares without gaps or overlaps. By this time I had grown tired of drawing zillions of dots, and so I started working with sets of coins.

initial configuration of four permutations of coins

In this arrangement each coin denomination forms a permutation, with no two pennies, nickels, dimes or quarters in the same row or the same column. We have successfully covered all the colored squares, but unfortunately we’ve also covered the blank at the lower right. Thus this pattern of coins is not an acceptable solution, but maybe we can fix it up somehow?

adjusted configuration after one coin is moved

Moving the penny from the blank square to another square in the same column solves one problem but creates another: Now the arrangement of pennies is no longer a permutation. There are two pennies in the third row.

coins after second adjustment to restore permutation

So now we have to shift another penny to restore the one-per-row-and-column property. Inevitably, this leaves a colored square uncovered. The only way we can cover that exposed square is to introduce a fifth permutation. Since I had run out of coin denominations, I chose a popular brand of breakfast toroids. Voila:

coins and cheerios -- the five-permutation solution

There’s nothing special about the particular moves I chose in this sequence. If you try some alternatives, you should be able to persuade yourself that moving the penny that covers the blank to any other square in the fourth column (or in the fourth row) would lead to essentially the same situation. Likewise the game would come out the same if the single blank square were placed anywhere else in the grid. And you could also start with a different set of initial permutations (provided they cover all the squares).

This coin-shuffling exercise demonstrates that we can cover any 4×4 template that has a single blank by combining no more than five permutations, but how do we know that five are actually needed? Maybe there’s some totally different arrangement that would do the job with just four permutations? Well, think about what such an arrangement would look like. It would have to differ at exactly one position from some other layout of permutations that covers the full 16-square grid. But no two permutations can differ at one and only one place. Thus the reason there can be no four-permutation cover of 15 squares is essentially the same as the reason no 4×4 ormat pattern can cover just five squares.

This argument generalizes to k×k matrices: For any integer k, there must be at least k ormat patterns that cannot be covered with fewer than k+1 permutations. But then comes the bigger speculative leap: Perhaps k+1 is an upper bound. Perhaps part of the answer to Barry’s question is that no k×k ormat pattern requires more than k+1 permutations. At one point I even had a “proof” of this conjecture. Then I wrote a program to check it, doing much the same thing I did with the dots on the airplane.

Out of bounds. My program found the expected 16 ormat patterns that require five permutations—and it found many more as well. In all it identified 2,032 4×4 ormats that can’t be composed from fewer than five permutations. And then came a bigger surprise: The program also found 480 patterns that require six permutations. So much for my proposed upper bound.

One of those problematic 480 ormats takes this form:

((1 1 1 1) (1 1 1 1) (1 1 1 1) (1 1 0 0))

Looking over this pattern, I thought I understood where my earlier reasoning had gone awry. This matrix is just like the single-zero pattern, but with two zeros! (I do mean for that statement to make sense. Bear with me.) Suppose we start again with a set of four permutations that completely cover the grid, including the two blanks.

starting configuration of 16 coins on 4x4 template with two blanks

Then we can uncover each blank just as we did in the coin-shuffling procedure above, although we have to be careful the two sets of movements don’t interfere with each other. (Not much point in removing the penny from a blank square, then putting the nickel there.) Here is a strategy for clearing both blank squares while maintaining the one-per-column-and-row permutation property:

four coins and two blanks: first solution

Inevitably, when we uncover the two blank squares, we also remove coins from two colored squares, which now have to be filled in again. The key point is that no single permutation can repair that damage, because the two open colored squares are in the same row. To cover both of those squares we need two additional permutations.

Other ways of reshuffling the coins avoid putting the two open squares in the same row or column, but they still foil all attempts to complete the covering with just five permutations. Try adding four Cheerios to the diagram below. If you cover both of the open blue squares, then either you also cover one of the blank squares or you wind up with two Cheerios in the same row.

four coins and two blanks: second solution

So now it’s clear we need as many as six permutations to cover a 4×4 ormat. Does that suggest that the general upper bound might be k+2 rather than k+1? Or perhaps the appropriate formula is 2k–2? In support of this latter possibility I offer these two ormats, which require 8 and 10 permutations respectively:

2k-2ormats.png

Another wager. Having fooled myself several times about the upper bound on minimal ormat coverings, I feel I should build in a little margin for error before I invite you to make a further wager. We already have direct evidence that covering a k×k ormat can take as many as 2k–2 permutations. So I’ll be generous and offer a full $2k for a proper covering, while charging $1 per permutation. If k=3 or k=4, you can definitely make money on this deal. But is it a good bet for larger k? (Hint: I’d be willing to play the game on these terms for real money.)

•     •     •

Update 2010-08-19: No takers for my bet, eh? Too bad; I had already spent my winnings.

Barry Cipra, who raised the question about minimal ormat covers in the first place, sends this illuminating letter:

I’m going to tiptoe a short ways out on a long long limb and conjecture (really just guess) that the “worst case” behavior, in terms of the minimum number of permutations it takes to produce a given ormat, occurs for ormats of the following form, shown here for k=7:

((1 1 1 1 1 1 1) (1 1 1 1 1 1 1) (0 1 1 1 1 1 1) (0 0 1 1 1 1 1) (0 0 0 1 1 1 1) (0 0 0 0 1 1 1) (0 0 0 0 0 1 1))

So as not to abuse existing matrix terminology, I’ll call any (square) matrix of this type—i.e., whose entries below the main subdiagonal are all 0—”uppity triangular.” I can (and will!) show that this uppity triangular ormat for k=7 requires (at least) 16 permutations—and the number appears to grow concavely upwards from that, so I, for one, will definitely not take you up on your $2k wager.

The trick, I realized, is to view each ormat as the “shadow” of what I’ll call an “addmat.” If you let P1, P2, …, Pr be k×k permutation matrices, their addmat is simply the ordinary result of addition: S = P1 + P2 + … + Pr, whose entries are positive integers wherever one or more of the constituent permutations has a 1 and otherwise 0. The associated ormat is obtained by changing each of these entries to a 1, while leaving the 0’s alone. In this sense, the ormat’s 1’s are the “shadows” of the addmat’s positive entries.

What’s crucial is that addmats have a lovely little property not shared with their shadows: the row and columns sums of the entries of an addmat all equal the number of permutations that produce them, r.

Come now, let us reason together…. The uppity triangular ormat example above (for k=7) must come from an addmat of the form

{{

where a, b, and all the *’s are positive integers. In particular, each * is at least 1. Since all row and column sums must be equal, the sum a+b must equal the sum of b and all 6 *’s above it. Hence a is at least 6. Likewise a+b must equal the sum of a and all 6 *’s above it, so b is also at least 6. Hence a+b is at least 12, which means the OR-sum that produced the given ormat involves at least 12 permutations.

This clearly generalizes to arbitrary k, which is more than “direct evidence” that covering a k×k ormat can take as many as 2k–2 permutations, it’s rigorous proof! But we can immediately do better, at least on a case-by-case basis. If we try to get by with just 12 permutations for this uppity triangular ormat, we quickly run into trouble. We obviously must have a = b = 6, and it follows that all the *’s above them are 1’s (to make those column sums 12). That is, we have the addmat

{{

where I now wish to call your attention to the entry labeled “@”. To make its row-sum equal 12, we need @ = 10. But that means its column sum (with the 5 *’s above it) is at least 15, which cannot be! So we are forced to try larger values of a and/or b—which is to say, we need more permutation matrices to produce this addmat.

It turns out you can’t satisfy the row and column sum condition until you get to a = b = 8. I won’t take you through all the steps, but just give you a taste with the penultimate possibility, a = 7, b = 8. The best you can hope for in this case is

{{

Note that I put as much of the “weight” in the last two columns as close to the 7 and 8 as possible, so that I could use the smallest possible value (10) as the entry with 5 positive entries above it. This makes the last three rows, and the right three columns all have the same sum, 15, but now we see a problem in the 12’s column: Its sum is at least 16. So once again, we’re screwed. It’s only with the next attempt that we avoid contradiction:

{{8, 3, 1, 1, 1, 1, 1}, {8, 3, 1, 1, 1, 1, 1}, {0, 10, 2, 1, 1, 1,<br />
  1}, {0, 0, 12, 1, 1, 1, 1}, {0, 0, 0, 12, 2, 1, 1}, {0, 0, 0, 0, 10, 3, 3}, {0, 0, 0, 0, 0, 8, 8}}

This matrix finally has all its row and column sums equal. Please note, this may or may not be an actual addmat of a set of 16 permutation matrices—I suspect it probably is, but I haven’t bothered to check. All we know is that it satisfies a necessary condition of being an addmat, namely that its row and column sums are all equal. (It’d be nice if that were also a sufficient condition, but something tells me it isn’t.)

This example, which can clearly be played out for larger values of k, suggests that not only are you safe with a $2k wager, but with a $(2k+2) wager and higher I’ve played around with this a bit, and persuaded myself that the number of permutation matrices will go to 2k + a lot—for k=10, if I did things correctly, you need 24 permutations (or possibly more, if the uppity triangular matrix the analysis leads to is not an actual addmat). I am entirely convinced that some additional careful thought can streamline the analysis into a nice, slick proof. I’m just not sure I haven’t already made a mistake, and built an elaborate house of cards….

Does any of this jibe with what you’ve already found to be the case?

It does indeed jibe.

First of all, to answer a small question Barry left open, here is a set of 16 permutations that will successfully cover his 7×7 “uppity triangular” matrix:

{1,2,3,4,5,6,7} {2,1,4,3,6,7,5} {1,3,2,5,4,7,6} {1,3,4,2,6,5,7}
{2,3,1,5,6,4,7} {1,2,3,5,6,7,4} {1,2,4,5,3,6,7} {1,2,4,5,6,3,7}
{1,2,4,5,6,7,3} {1,3,4,5,2,6,7} {1,3,4,5,6,2,7} {1,3,4,5,6,7,2}
{2,3,4,1,5,6,7} {2,3,4,5,1,6,7} {2,3,4,5,6,1,7} {2,3,4,5,6,7,1}

This was found with a simple greedy search.

My own attempts to find an upper bound have focused not on uppity triangular matrices but on matrices I’ve been calling “flags,” like this 7×7 case:

{{0, 0, 0, 1, 1, 1, 1}, {0, 0, 0, 1, 1, 1, 1}, {0, 0, 0, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1, 1}, {1, 1, 1, 1, 1, 1, 1}}

This matrix also requires 16 permutations for a proper covering. To see why, try threading permutations through the columns of the matrix, starting at the left edge and in each column choosing a 1 element (never a 0) from a different row. Because of the block of zeros at the upper left, the first three elements of every permutation must lie in rows 4 through 7. Thus each permutation “uses up” three of the last four rows in the first three columns, and the rest of the permutation can revisit this range of rows only once. It follows that each permutation can touch only one element in the 4×4 block of 1s in the lower right corner of the matrix, and at least 16 permutations are needed to cover all the 1s in the matrix. Showing that 16 are sufficient is not hard.

This kind of analysis works for any odd k, and thus we know that such matrices can require as many as

{\biggl\lceil\frac{k}{2}\biggr\rceil}^2

permutations. (For even k the situation is a little less symmetrical, and I haven’t worked out the exact details.)

These results give us a lower bound on the upper bound on the number of permutations that may be needed to cover a k×k ormat. But we haven’t proved it’s the true upper bound. Are there other ormats that require even more permutations? My guess is no, but keep in mind that almost all my conjectures along these lines have turned out to be wrong.

The state of the spamosphere

8 August 2010

BP has finally stoppered the Macondo well with a plug of mud and cement, but the gusher of spam continues to pollute inboxes everywhere. Maybe we need a relief well?

spam statistics Aug 2008 through July 2010

It’s been six months since my last spam update. The good news, I suppose, is that last summer’s huge spurt of spam has subsided. But I’m still getting 2,000 inanities a month. That’s six or seven times the rate I was seeing when I first started monitoring my spam intake in 2003. For the two years covered by the graph above, the cumulative total is 111,945 unwanted emails received.

Needless to say, most of those messages are utterly unremarkable. But as I reviewed the most recent batch of dreck, one series of emails caught my eye. In the past 24 hours I’ve received five copies of a spam with the subject line: “Solve this if you could…!!!!”:

For all Math’s champs, Accounting Experts & Number Numbssss… (including future champs too)

Find the 6th Number

1, 2, 6, 42, 1806, ____?

6th number is the password of the attachment.

The attachment mentioned in the last line is a zip archive that unpacks to reveal an .exe file, which I would not run even if I could. But I confess that I did stop to solve the little puzzle sequence. It’s not difficult, although it is a little harder than any of the series I use in the spambot filter on the bit-player comment form. Maybe I should add it to my repertory.

Four questions about fuzzy rankings

24 July 2010

The National Research Council is getting ready to release a new assessment of graduate-education programs in the U.S. The previous study, published in 1995, gave each Ph.D.-granting department a numerical score between 0 and 5, then listed all the programs in each discipline in rank order. For example, here’s the top-10 list for doctoral programs in mathematics (as presented by H. J. Newton of Texas A&M University):

 rank    school                 score   
    1    Princeton               4.94
    2    Cal Berkeley            4.94
    3    MIT                     4.92
    4    Harvard                 4.90
    5    Chicago                 4.69
    6    Stanford                4.68
    7    Yale                    4.55
    8    NYU                     4.49
    9    Michigan                4.23
   10    Columbia                4.23

Note that the scores of the first two schools are identical (to two decimal places), and the first four scores differ by less than 1 percent. Given the uncertainties in the data, it seems reasonable to suppose that the ranking could have turned out differently. If the whole survey had been repeated, the first few schools might have appeared in a different order. Doctoral candidates in mathematics are presumably sophisticated enough to understand this point. Nevertheless, the spot at top of the list still carries undeniable prestige, even when you know that the distinction could be merely an artifact of statistical noise.

The committee appointed by the NRC to conduct the new graduate-school study wants to avoid this “spurious precision problem.” They’ve adopted some jazzy statistical methods—mainly a technique called resampling—to model the uncertainty in the data, and they’ve also decreed that the results will be presented differently. There will be no sorted master list showing overall ranks in descending order. Instead the programs in each discipline will be listed alphabetically, and each program will be given a range of possible ranks. For example, a program might be estimated to rank between fifth place and ninth place. Let’s call such a range of ranks a rank-interval, and denote it {5, 6, 7, 8, 9} or {5–9}.

For a hypothetical set of 10 institutions, A through J, here’s what a set of rank-intervals might look like.

bar graph showing ranges of rankings for schools A through J.png

Acknowledging the uncertainty in your findings is commendable. But let’s be realistic. If you actually want to make use of these results—for example, if you’re a student choosing a grad-school program—the first thing you’re going to do is sort those bars into some sort of rank order, trying to figure out which school is best and how they all stack up against one another. In other words, you’re going to undo all the elaborate efforts the NRC committee has put into obscuring that information.

Below is one possible ordering of the bars. I have sorted first on the top of the rank-intervals, then, if two columns have the same top rank, I’ve sorted on the bottom rank. Other sorting rules give similar but not identical results. For example, sorting on the midpoints of the intervals would interchange columns B and F.

bar graphs showing rank-ranges sorted into one canonical order.png

Question 1. Does sorting a set of rank-intervals by one of these simple rules yield a consistent and meaningful total ordering of the data? To put it another way, can you trust this attempt to reconstruct a ranking?

I hasten to add that this is not really a practical question about finding the best grad school. If you’re facing such a choice in real life, the NRC rank-intervals are not the only available source of information. But, for the sake of the mathematical puzzle, let’s pretend that all we know about schools A through J is embodied in those ranges of rankings.

It turns out that rank-intervals have some fairly peculiar behavior. Ranges of ratings are not a problem. If the NRC merely gave each school a fuzzy rating on the 0-to-5 scale, no one would have much trouble interpreting the results. But when you turn fuzzy ratings into fuzzy rankings, there are hidden constraints. For example, not all sets of rank-intervals are well-formed.

two impossible sets of rank-ranges

The set at left is impossible because there’s no one in last place. (We can’t all be above average.) The example at right is also nonsensical because D has no ranking at all. For a set of rank-intervals to be valid, there has to be at least one entry in each row and each column.

That’s a necessary condition, but not a sufficient one, as the two graphs below illustrate.

two more impossible rank-intervals

Do you see the problem with the example at left? Column B has a rank-interval of {1–2}, but in fact B can never rank first because A has no alternative to being first. The case at right is conceptually similar but a little subtler: If B is ranked third, then either first place or second place will have to remain vacant.

The underlying issue here is the presence of constraints or linkages within a set of rankings. Suppose you have calculated ratings and rankings of several schools, and then some new information turns up about one school. You can change the rating of that school without any need to adjust other ratings, but not so the ranking. If a school goes from third place to fourth place, the old fourth-place school has to move to some other rung of the ladder, and somebody has to fill the vacancy in third place. These interdependencies are obvious in a non-fuzzy ranking, but they also exist in the fuzzy case. You can’t just assign arbitrary rank-intervals to the items in a set and assume they’ll all fit together. This observation leads to a second question:

Question 2. What are the admissible sets of rank-intervals? How do we characterize them?

I have a partial answer to this question. It goes like this. Any ranking of k things must be a permutation of the integers from 1 through k. A permutation can be embodied in a permutation matrix—a square k × k matrix in which every row has a single 1, every column has a single 1, and all the other entries are 0. For example, here are the six possible 3 × 3 permutation matrices:

3x3-permutation-matrices.png

They correspond to the rankings (1, 2, 3), (1, 3, 2), (2, 1, 3), (3, 1, 2), (2, 3, 1) and (3, 2, 1).

Since a permutation matrix represents a specific (non-fuzzy) ranking, we can build up a set of rank-intervals by taking the OR-sum of two or more permutation matrices. What do I mean by an OR-sum? It’s just the element-by-element sum of the matrices using the boolean OR operator, ∨, instead of ordinary addition. OR has the following addition table:

                      0 ∨ 0 = 0
                      0 ∨ 1 = 1
                      1 ∨ 0 = 1
                      1 ∨ 1 = 1

For the first two 3 × 3 matrices shown above the arithmetic sum is:

matrix-addition.png

whereas the OR-sum looks like this:

matrix-or-sum.png

Every valid set of rank-intervals must correspond to an OR-sum of permutation matrices, simply because a set of rank-intervals is in fact a collection of permutations. The converse also holds: Any OR-sum of permutation matrices yields an admissible set of rank-intervals. Thus the OR-sums of permutation matrices—let’s call them ormats for brevity—are in one-to-one correspondence with the admissible sets of rank-intervals. (There’s just one catch when applying this idea to the NRC study. The columns of an ormat may well have “gaps,” as in the column pattern (0 1 1 0 0 1 1), which corresponds to the rank-interval {2–3, 6–7}. Will the NRC allow such discontinuous ranges in their grad-school assessments? Perhaps the issue will never come up in practice. In any case, I’m ignoring it here.)

Arithmetic sums of permutation matrices form an open-ended, infinite series; in contrast, there are only finitely many distinguishable OR-sums. The reason is easy to see: Ormats have k2 entries, each of which can take on only two possible values, and so there can’t be more than \(2^{k^{2}}\) distinct matrices. Because of the various constraints on the arrangement of the entries, the actual number of ormats is smaller. For example, at k = 3 the \(2^{k^{2}}\) upper bound allows for 512 ormats, but there are only 49:

the-49-3-by-3-or-sums.png

Thus we come to the next question.

Question 3. For each k ≥ 1, how many distinct ormats can we build by OR-ing subsets of k × k permutation matrices? Is there a closed-form expression for this number?

I have answers only for puny values of k.

   k       upper bound        # of ormats  
   1                 1                  1
   2                16                  3
   3               512                 49
   4            65,536              7,443
   5        33,554,432          6,092,721
   6    68,719,476,736                  ?

The tallies of ormats were calculated by direct enumeration, which is not a promising approach for larger k. (I note—to spare folks the bother of looking—that the sequence 1, 3, 49, 7443, 6092721 does not yet appear in the OEIS.)

To extend this series, we might try to exploit the internal structure and symmetries of the ormats. By sorting the columns and rows of the matrices, we can reduce the 49 3×3 ormats to just six equivalence classes, with the following exemplars:

exemplars of six ormat equivalence classes

Enumerating just these reduced sets of matrices should make it possible to reach larger values of k, but I have not pursued this idea. (Furthermore, the two-dimensional sorting of matrices looks to be a curiously challenging task in itself.)

By the way, I think the number of ormats will approach the \(2^{k^{2}}\) upper bound asymptotically as k increases. Many of the features that disqualify a matrix from ormathood—such as all-zero rows or columns—become rarer when k is large. I have tested this conjecture by generating random (0,1) matrices and then counting how many of them turn out to be ormats.

fraction-of-ormats.png

For k = 1 through 5 the results are in close agreement with the actual counts of ormats, and up to k = 10 the trend is clearly upward. But continuing this inquiry to larger values of k will depend on a positive answer to the next question.

Question 4. Given a square matrix with (0,1) entries, is there an efficient algorithm for deciding whether or not it is an OR-sum of permutation matrices, and thus an admissible set of rank-intervals?

The question asks for a recognition predicate—a procedure that will return true if a matrix is an ormat and otherwise false. If efficiency doesn’t matter, there’s no question such an algorithm exists. At worst, we can generate all the k × k ormats and see if a given matrix is among them. But that’s like saying we can factor integers by producing a complete multiplication table. It just won’t do in practice. Isn’t there a quick and easy shortcut, some distinctive property of ormats that will let us recognize them at a glance?

If we could replace the OR-sum with the ordinary arithmetic sum, the answer would be yes. Permutation matrices have the handy property that all rows and columns sum to 1. An arithmetic sum of r permutation matrices has rows and columns that all sum to r. (It is a semi-magic square.) The converse is also true (though harder to prove): If a matrix of nonnegative integers has rows and columns that all sum to r, it is a sum of r permutation matrices. This fact yields a simple test: Sum the rows and the columns and check for equality.

Unfortunately, the trick won’t work for ormats, because the boolean OR operation throws away even more information than summing does. Because 0 ∨ 1 = 1 ∨ 0 = 1 ∨ 1, infinitely many sets of operands map into the same result, and there’s no obvious way to recover the operands or even to determine how many permutation matrices entered into the OR-sum.

Maybe there’s some other clever trick for recognizing ormats, but I haven’t found it. Let me make the question more concrete. Below are three (0,1) square matrices. Two of them are ormats but the third is not. Can you tell the difference?

three-puzzle-matrices.png

If it’s so hard to recognize an ormat, how did I count the ormats among a bunch of randomly generated (o,1) matrices? By hard work: I reconstructed the set of permutations allowed by each matrix. Visualize a permutation as a path threading its way through the matrix from left to right, connecting only non-zero elements and touching each column and each row just once. When you have drawn all possible permutation paths, check to see if every non-zero element is included in at least one path; if so, then the matrix is an ormat. Note that this is not an efficient recognition procedure. In the worst case (namely, an all-ones matrix), there are k! permutations, so this method has exponential running time. But k! is better than \(2^{k^2}\); and, besides, for sparse matrices the number of permutations is much smaller than k!. The 10 × 10 matrix presented as an example at the start of this post gives rise to 580 permutations, a manageable number. Here’s what they look like, plotted as a spider web of red paths across the bar chart.

ranges-with-paths

Every nonzero site is visited by at least one permutation path, so this set of rank-intervals is indeed valid.

This process of lacing permutations through a matrix finally brings me back to Question 1, about how to make sense of the NRC’s fuzzy ranking scheme. Let’s take a small example:

probability-example-1.png

Examining the graph above shows that A must rank either first or second—but which is more likely? In the absence of more-detailed information, it seems reasonable to assume the two cases are equally likely; we assign them each a probability of 1/2. Similarly, B has the rank-interval {1–3}, and so we might suppose that each of these three cases has probability 1/3. Continuing in the same way, we assign probabilities to every element of the matrix.

probability-example-2.png

But wait! This can’t be right; our probabilities have sprung a leak. Any proper set of probabilities has to sum to 1. Our procedure assures that each column obeys this rule, but there is no such guarantee for the rows. In row 1, we’re missing one-sixth of our probability, and in row 2 we have an excess of 1/2; row 4 comes up short by 1/3.

Is there any self-consistent assignment of probabilities for the elements of this matrix? Sure. As a matter of fact, there are infinitely many such assignments, including this one:

probability-example-3.png

I’ll return in a moment to the question of how I plucked those particular numbers out of the air, but note first what they imply about the ranking of items A through D. For item A, with the rank-interval {1–2}, the odds are two-to-one that it ranks first rather than second. B has the behavior we expected from the outset, with probability uniformly distributed over the three cases. But if you pick either C or D, each with the rank-interval {2–4}, your chance of getting second place is only 1/6, and half the time you’ll be in last place.

Where do these numbers come from? Instead of starting with the assumption that probability is uniformly distributed over each rank-interval, assume that each possible permutation of the ranks is equiprobable. For this matrix there are six allowed permutations: (1, 2, 3, 4), (1, 2, 4, 3), (1, 3, 2, 4), (1, 3, 4, 2), (2, 1, 3, 4) and (2, 1, 4, 3). Observe that four of the six ordering put A first, and only two permutations place A second. We can also tally up such “occupation numbers” for all the other matrix elements:

probability-example-4.png

Dividing these numbers by the total number of permutations, 6, yields the probabilities given above.

We can do the same computation for the 10 × 10 example matrix, which turns out to allow 580 permutations:

ranges-with-path-weights.png

If you care to check, you’ll find that each column and each row sums to 580; dividing all the entries by this number yields a probability matrix with columns and rows that sum to 1 (also known as a doubly stochastic matrix).

This process of tabulating permutation paths recovers some of the information we would have gotten from the arithmetic sum of the permutation matrices—information that was lost in the OR-ing operation. But we get back only some of the information because we have to assume that each permutation included in the OR-sum appears only once. (This is just another way of saying that the allowed permutations are equiprobable.) There’s no particularly good reason to make this assumption, but at least it leads to a feasible probability matrix.

Is there any way of calculating the entries in the doubly stochastic matrix without explicitly tracing out all the permutation paths? I’m sure there is. I think the construction of the matrix can be approached as an integer-programming problem, and perhaps through other kinds of optimization technology. What seems less likely is that there’s some simple and efficient shortcut algorithm. But I could be wrong about that; there’s a lot of mathematics connected with this subject that I don’t understand well enough to write about (e.g., the Birkhoff polytope). I hope others will fill in the gaps.

Getting back to the assessment of grad schools—have we finally found the right way to understand those rank-intervals that the NRC promises to publish any day now? My sense is that a semi-magic square (or, equivalently, a doubly stochastic matrix) will give a less-misleading impression than a simple eyeball sorting on the spans or midpoints of the rank-intervals. But what a lot of bother to get to that point! How many prospective grad students are going to repeat this analysis?

Acknowledgment: Thanks to Geoff Davis of PhDs.org for introducing me to this story. PhDs.org will have the new ratings as soon as the NRC releases them, and may even find a way to make them intelligible! Disclaimer: I’ve done paid work for the PhDs.org web site (but this is not a paid endorsement).

Update 2010-07-27: If you’ve gotten this far, please read the comments as well. A number of commenters have provided important insights and context, which have helped me understand what’s going on in the matrices I’ve been calling ormats. But I’m still a bit murky about the best way to recognize and count them. I’m not sure that publishing my still-murky thoughts is terribly helpful, but maybe someone else will read what follows and give us a dazzling, gemlike synthesis.

For the ormat-recognition problem (Question 4 above), three basic approaches have been mentioned: enumerating the permutation paths through the matrix, examining matrix minors, and looking for perfect matchings in a bipartite graph defined by the matrix. It seems to me that all of these methods are doing the same thing.

Start with Barry Cipra’s method of minors. The basic operation is to choose a nonzero matrix element, then delete the row and the column in which that element occurs. You then apply the same operation to the remaining, smaller matrix.

In tracing permutation paths, we’re looking for sequences of nonzero elements, drawing one element from each column and each row. A way of organizing this search is to choose a nonzero element and then, after recording its location, delete the corresponding column and row, so that no other elements can be chosen from that column or row.

In the method based on Hall’s theorem, as explained by John R., we view the ormat as the adjacency matrix of a bipartite graph, where every nonzero element designates an edge connecting a row vertex to a column vertex. To find a matching, we delete an edge, along with the two vertices it connects (and also all the other edges incident on those vertices). Then we recurse on the smaller remaining graph. (See further update below.) If you translate this operation on the graph back into the language of matrices, deleting an edge and its endpoints amounts to deleting a row and a column of the adjacency matrix.

I am not asserting that these three algorithms are all identical, but they all rely on the same underlying operation. To say more, we would need to consider the control structure of the algorithms—how the basic operations are organized, how the recursion works, all the details of the bookkeeping. I don’t trust myself to make those comparisons without trying to implement the three methods, which I have not yet done. However, at this point I just don’t see how any method can guarantee correct results without something resembling backtracking (or else exhaustive search through an exponential space). After all, we’re not looking for just one matching in the graph, or one decomposition into matrix minors, or one permutation path; we have to examine them all.

Here’s a further hand-wavy argument for the essential difficulty of the task. For a (0,1) matrix, the number of permutation paths that avoid all zero entries is equal to the permanent of the matrix. Computing the permanent of such a matrix is known to be #P-complete.

Update 2010-07-31: With lots of help from my friends, I think I finally get it. Although there could be as many as k! permutation paths in a k × k matrix, you don’t need to examine all of the paths to decide whether or not the matrix is an ormat. It’s enough to establish that one such path passes through each nonzero element. This is what the algorithm based on Hall’s theorem does. As Frans points out in a comment below, I misunderstood the essential nature of that algorithm (in spite of having it explained to me several times). There is no recursive deconstruction into progressively smaller matrix minors; instead, we just loop over all the nonzero elements of the matrix, find the minor associated with each such element, then check for a perfect matching in the minor. (Still more refinements are possible—but already we have a polynomial algorithm.)

With this efficient recognizer predicate, it’s easy to measure the proportion of ormats in random matrices at larger values of k:

fraction-of-ormats-k25.png

As expected, the fraction of ormats approaches 1 beyond about k = 20.

So much for identifying ormats. I am still unable to extend the series of exact counts beyond k = 5. The tabulations for random (0,1) matrices suggest that for k = 6 there should be about 20 billion ormats, and counting that high is just too painful. I need to work out the symmetries of the problem.

As far as I can tell, assigning exact probabilities to the nonzero matrix elements requires a full enumeration of all the permutation paths, and thus a calculation equivalent to the permanent. There may be a useful approximation.

Barry Cipra asks a really good question: The permanent tells us the maximum number of permutations that could possibly be included in a given ormat, but what is the minimum number? A naive upper bound is the number of 1s in the matrix, but I don’t see an easy path to an exact count. But enough for now.

The thrill of the chase

11 July 2010

How I love to go out hunting on a bright Sunday morning—though it’s not my style to shoot furry/feathery/finny animals. My game is to get up early and stalk a wily factoid.

A posting from Mat Roberts, whose blog I’ve recently discovered, sent me out this morning to chase down a passage in How Long Is a Piece of String, a book by Rob Eastaway and Jeremy Wyndham:

passage from Eastaway-Wyndham, page 160

The concept here seemed familiar, but the term “Lincoln Index” was new to me. Lincoln who? What index?

Google offered some useful clues. (Also a generous helping of false scents—books about Honest Abe that happen to have an index.) Without even clicking on a link I had the general context:

The Lincoln Index provides a way to measure population sizes of individual animal species. It is based on a capture/mark/ recapture method…

So we’re talking ecology and population biology. The original idea was not to catch the same typo twice but to catch the same furry/feathery/finny creature twice. Interesting. However, the first couple of web pages that Google sent me to (here and here) told me nothing about Lincoln. And, oddly, I found no Wikipedia entry for “Lincoln Index.” If it’s not in Wikipedia, does it exist?

With a little more poking around, I stumbled upon another clue that seemed promising: a mention of “the Lincoln-Pearson equation for estimating population size.” I was still in the dark about Lincoln, but Pearson is quite a familiar figure. Surely that’s Karl Pearson, the pioneering statistician, who did much of his work in the biological sciences and might very well have come up with a scheme for estimating population sizes.

Back at Google, though, searching for “Lincoln-Pearson” turned up nothing pertinent other than the page I’d come from (though I did learn that Karl Pearson “read in chambers in Lincoln’s Inn” during his early years studying law).

More beating the bushes. Eventually I realized I had wandered into a blind alley. Somebody needs to hire a pair of proofreaders: The formula is not “Lincoln-Pearson” but “Lincoln-Petersen.” Try those names at Google and you’ll get an abundance of useful pointers. (You’ll also learn that Abraham Lincoln died in Petersen’s Boarding House, across the street from Ford’s Theater. Google is not just a search engine but also a coincidence engine.)

The particular web page where I finally got the correct names (notes for a course at North Carolina State University) explains that capture-mark-recapture methods

are used extensively to estimate populations of fish, game animals, and many non-game animals. The approach was first used by Petersen (1896) to study European plaice in the Baltic Sea and later proposed by Lincoln (1930) to estimate numbers of ducks. Petersen’s and Lincoln’s method is often referred to as the Lincoln-Petersen Index, even though it is not an index but a method to estimate actual population sizes. (Should it not be the Petersen-Lincoln Estimate?)

I decided to pursue Petersen first—and immediately ran into a few further bibliographic brambles. Some citations spell the name “Petersen” and others “Peterson.” Some give the initials “C. G. T.” and others “C. G. J.” or “C. J. G.” The date might be 1895 or 1896 or 1897. Here’s what I believe to be a correct citation:

Petersen, C. G. J. 1896. The yearly immigration of young plaice into the Limfjord from the German Sea. Report of the Danish Biological Station to the Home Department 6:1–48.

Wikipedia identifies our elusive author as Carl Georg Johannes Petersen (1860-1928). He was a founder of the Danish Biological Station, which was not in fact a station but a mobile laboratory—a decommissioned naval vessel that was moved around from year to year. In 1895, Petersen took the station to the Limfjord, a chain of bays, lakes and channels cutting across the Jutland peninsula in northern Denmark. There he studied the plaice fishery. (Back to Wikipedia: “The European plaice is a right-eyed flounder belonging to the Pleuronectidae family.” But let’s not get started on right-eyed and left-eyed flatfish, or we’ll never get to the end of this.)

Petersen’s report is available online, scanned from a copy belonging to the library of the Marine Biological Laboratory and Woods Hole Oceanographic Institution, and hosted by the Biodiversity Heritage Library of the Internet Archive. A second surprise: The report is written in English. But on reading through it I find only vague and murky connections between the work Petersen reports and the mark-recapture method of estimating populations. There’s nothing resembling the E1E2/S formula.

Petersen does describe a series of capture/mark/recapture experiments. A few hundred plaice were caught and marked by attaching numbered buttons, then put back in the water. Fishermen who recaught the labeled fish in later months were asked to report them. But the purpose of this study was not to estimate the total population; instead, Petersen used before-and-after measurements of the marked fish to estimate their growth rate.

In a much larger experiment, some 82,580 plaice (somebody must have counted them!) were transplanted into the fjord, and 10,900 of the fish were marked by having a hole punched in their dorsal fin. The number of marked fish was recorded as the plaice were caught during the coming year. It’s not clear whether the aim of this project was to estimate the total population, but in any case it didn’t work. The fraction of marked fish in the transplanted batch was about 1/7, but the marked fraction in the subsequent catches was 1/5. Petersen remarks, “This result is very strange,” and I have to agree.

When Petersen did try to estimate the plaice population, he didn’t rely on a recapture scheme. He went out with seine nets designed to dredge up every bottom fish in a measured plot, then extrapolated from the density of fish per unit area.

The whole report is fascinating fishy stuff, but it leaves me wondering just how Petersen came to be given credit for the resampling idea. As far as I can tell, it’s not to be found in this paper.

Having chased down Petersen, I turned back to Mr. Lincoln. Without much trouble I was able to identify the work in question:

Lincoln, F. C. 1930. Calculating waterfowl abundance on the basis of banding returns. United States Department of Agriculture Circular 118:1–4.

portrait of Frederick C. Lincoln in his office, with stuffed duck.The author was Frederick C. Lincoln, who was bird-bander-in-chief in the U.S. for some 25 years. The agency he founded has since migrated from the Department of Agriculture to the U.S. Geological Survey and become the Bird Banding Laboratory.

Google returns hundreds of works that cite Lincoln’s paper (including some quite far afield from population biology). But tracking down the USDA document itself was not so easy. If the USDA has it online, I wasn’t able to locate it. But a search of WorldCat eventually turned up an archive in the Hathi Trust Digital Library where you can page through Lincoln’s pamphlet in a copy scanned by Google at the University of Minnesota library.

Lincoln gives only a brief and informal account of the recapture idea, but the basic principle is stated clearly enough:

If in one season 5,000 ducks were banded and yielded 600 first-season returns, or 12 percent, and if during that same season the total number of ducks killed and reported by sportsmen was about 5,000,000, then this number would be equivalent to approximately 12 per cent of the waterfowl population for that year, which would be about 42,000,000.

It’s not hard to translate this formula from the language of duck hunters into the language of proofreaders. The first reader finds 5,000 typos and the second spots 5 million; 600 of these errors are common to both lists, and so the total number of typos is:

\frac{5\,000 \times 5\,000\,000}{600} = 41\,666\,667

So that’s my reward for a morning spent out hunting: 42 million typos.

Does Frederick Lincoln deserve credit for the Lincoln Index? I’d say he has a good claim, except that Pierre Simon de Laplace had the same idea more than a century earlier. In 1802 Laplace applied his method to estimating the (human) population of France. But maybe that’s a story for another Sunday morning.

Epilogue. This is not really a story about typos, or about fish and ducks. It’s about finding things—about the phenomenal ease of chasing facts on the world wide web. Does a marked fish have any hope of escaping recapture there?

 

A twist of fate

11 June 2010

The school of philosophy called Antipodianism briefly flourished on the fringes of the Hellenistic world more than 2,000 years ago. The sect held that every person has an opposite number, a mirror image who inverts all our beliefs, feelings, actions and attitudes. If I smile, my antipodian counterpart frowns; when I wake, she sleeps. If I’m a Mac, she’s a PC. For every liberal Democrat there’s an antipodian Tea Party Republican. In this way the universe is held in balance. It’s an enforced equilibrium, which none of us has the power to upend, hard as we might try.

The earliest of the Antipodians believed that every such matched pair (commonly designated A and ∀) live at diametrically opposite points on the surface of the earth. This arrangement ensures that A and ∀ can never meet—thereby averting a cosmic catastrophe. A later quantum-field version of Antipodianism relaxed the geographic constraint by allowing for the creation and annihilation of A∀ pairs, but that idea never really caught on.

One day an Antipodian master was teaching an exchange student from New Zealand. The child was crafty.

“Is it not true,” she asked, “that A and ∀ always do the opposite thing?”

“Yes, antisymmetry demands it,” the master replied.

“If A walks north, ∀ must walk south?” the child asked.

Again the master assented.

“If A goes east, ∀ must go west?”

“Yes.”

“If A turns to the right, ∀ must turn to the left, no?”

The master agreed, although he sensed trouble coming.

“I’m afraid the universe is out of joint,” said the child. “If A goes north and turns to the right, while ∀ goes south and turns to the left, afterwards they are both walking east. They are doing the same thing.”

Needless to say, this was a moment of crisis in Antipodian doctrine. The Pythagoreans, you may recall, resolved a similar impasse by resorting to violence. When some upstart challenged their precept that “all is number” by showing that no known number can be the square root of 2, the Pythagoreans tossed the troublemaker out of the boat. But in this case the Antipodian master kept his calm.

“Ah my little Kiwi,” he said to the student. “You are clever but not wise. Your own statements refute your claim. Did you not begin by saying that A and ∀ always do the opposite thing? When A walks 10 paces north, ∀ walks 10 paces south. When A turns right, ∀ turns left. But now you would have us believe they both take a step forward, contradicting the most basic law of their nature. What really happens is that A walks forward and ∀ walks backward. Thus A goes eastward and ∀ westward, and all is well with the world.”

Through this brittle sophistry the master extricated himself from the classroom—though he may have had to walk backwards to make good his escape. He never taught again. The Kiwi student went on to a brilliant career studying the weak interactions of neutral K mesons. As for Antipodianism, it vanished without a trace.

Or maybe it left a tiny trace. I’ve never visited the antipodes, but I hear that corkscrews Down Under turn the other way.

 

Disentangling Gaussians

10 June 2010

The printed program for the recent STOC meeting in Cambridge announced the following talk:

Talk by Kalai, Moitra, Valiant and Erdos on printed program, STOC 2010

Unfortunately, the fourth author could not be present (and he is not listed as an author on the paper itself), so the talk was given by Gregory Valiant.

Here is a motivating example. Take a tape measure, and go record the heights of a few thousand adults chosen at random. You’ll come back with a distribution that looks something like this:

combined-distribution.png

How come the curve is so lumpy and lopsided? Isn’t height a variable with an approximately normal distribution? Of course it is. The problem is that we have mixed up two subpopulations—men and women—with different height distributions:

three-distributions.png

Question: Given the lumpy combined distribution, and the knowledge that it represents a mixture of exactly two normal distributions, can we somehow recover the component distributions? Specifically, can we determine the means and standard deviations and the Gaussian curves, as well as their relative weights, or contributions to the total?

The answer is yes. In a series of papers published in the 1950s and 60s, Henry Teicher (then at Purdue, now emeritus at Rutgers) proved a kind of unique-factorization theorem for distributions. He showed that no mixed distribution can be decomposed into two or more different sets of normal distributions, just as no composite whole number can be formed as the product of two or more distinct sets of primes. (Teicher’s theorem generalizes to many distributions other than normal ones. There are also a few caveats; for example, none of the component distributions can be pointlike, with a standard deviation of zero.) Teicher’s work implies that a mixed distribution is identifiable: If you can break it down into a set of primitive distributions, then that decomposition is unique.

If STOC were a mathematics meeting, that might be the end of the story. But STOC is the Symposium on the Theory of Computing, and the question here is not “Does a solution exist?” but rather “Can you solve it in polynomial time?” Teicher’s result offers no such guarantee. It turns out that his proof of identifiability depends on the behavior of the distributions far out in the tails. Measuring data frequencies in these sparsely populated outlying regions requires an exponentially large number of samples. But other approaches to the problem don’t run into this snag; Valiant and his colleagues show that “robust polynomial identifiability” is indeed possible.

A key idea in their proof goes back at least as far as the 1890s, to work done by the indefatigable statisticians W. F. R. Weldon and Karl Pearson. In 1892 Weldon spent a summer on the Bay of Naples, measuring various features on the carapaces of crabs. One such feature yielded a distinctively lumpy distribution much like the height curve shown above. Weldon thought that the asymmetry might signal the incipient splitting of the crab population into two races or species, each of which if taken individually would have a normal distribution. For help with the analysis of his data he turned to Pearson, who was able to identify two component Gaussian curves that sum up to the observed distribution. (The figure comes from Weldon’s 1893 paper; see below for references.)

Figure 3 distribution from Weldon 1893.

Pearson’s method was based on calculating the first six statistical moments of the given distribution. (The nth moment of a distribution is the expectation value of \((x-\bar{x})^n\), where \(x\) is a random value drawn from the distribution and \(\bar{x}\) is the mean value.) The process called for solving a ninth-degree polynomial, which was a heroic feat in 1893.

Valiant et al. show that a procedure like Pearson’s will always succeed in the following sense: If two mixed distributions can be distinguished at all, then they can be distinguished by examining their first six moments, and those moments characterize the component Gaussians. Calculating the moments requires only a polynomial number of samples, and the running time of the overall algorithm is a polynomial function of various parameters such as the required accuracy. Moreover, the process can be extended from one-dimensional distributions to multidimensional Gaussians.

Although the algorithm has polynomial running time, that’s not a guarantee of practicality. (After all, this is a theory conference.) One stage of the process is essentially a brute-force search through a very large (though polynomially bounded) space of parameter values for the six moments.

Sources:

The STOC proceedings with the Kalai-Moitra-Valiant paper are now available online, but only by subscription. A PDF preprint is posted on Valiant’s web site; also an expanded version.

For Henry Teicher’s work see: “Identifiability of Mixtures,” The Annals of Mathematical Statistics 32(1):244–248 (1961) and “On the Mixture of Distributions,”  The Annals of Mathematical Statistics 31(1):55–73 (1960).

For Weldon and Pearson see: W. F. R. Weldon: “On Certain Correlated Variations in Carcinus maenas,“ Proceedings of the Royal Society of London 54:318–329 (1893) and Karl Pearson: “Contributions to the Mathematical Theory of Evolution,” Philosophical Transactions of the Royal Society of London A 185:71- 110 (1894).

My thanks to Virginia Gold and Irene Frawley at ACM, Lance Fortnow, Paul Oka, Jennifer Chayes and Christian Borgs, all of whom helped make it possible for me to attend some of the STOC sessions.

The snarXiv

7 June 2010

The snarXiv is a ran­dom high-energy the­ory paper gen­er­a­tor incor­po­rat­ing all the lat­est trends, entropic rea­son­ing, and excit­ing mod­uli spaces. The arXiv is sim­i­lar, but occa­sion­ally less ran­dom.

Inspiring! Soon bit-player, too, will be generated by a context-free grammar, and no one will know the difference.

 

A hole in the bottom of the ocean

4 June 2010

The explosion and fire that destroyed the drilling rig Deepwater Horizon on the night of April 20 was a run-of-the-mill industrial accident. In saying this, I don’t mean to make light of the disaster, in which 11 workers perished. Nevertheless, it’s important to recognize that events like this one have happened before. They happen all the time. Just two weeks earlier, an explosion in a West Virginia coal mine killed 29 workers. A refinery fire in Anacortes, Washington, killed seven on April 2. In February six workers died in a natural gas explosion at a power plant under construction in Connecticut. In January, leaking phosgene gas killed an operator at a West Virginia chemical plant. Looking back a few years further, a 2005 fire at a BP refinery in Texas City, Texas, killed 15 workers and injured 170 more. There was another fatality at the same refinery in 2008. Texas City was also the site of a fertilizer explosion in 1947 that destroyed much of the port and left almost 600 dead. Then there’s the Bhopal catastrophe, where toxic fumes from an insecticide plant suffocated thousands of residents of nearby neighborhoods.

Still another accident that belongs in this sad catalogue is the destruction of the Piper Alpha platform in the North Sea oil fields. Piper Alpha was not a drilling rig but a production platform; it pumped oil and natural gas from completed wells to a terminal in the Orkney Islands. On the night of July 6, 1988, a huge explosion shattered the main part of the platform, and a subsequent fire consumed the rest; 167 crew members died. The apparent cause was a miscommunication: The night shift tried to start up a pump, not knowing that the day shift had removed a crucial valve for maintenance, leading to a massive leak of flammable gas.

It’s a bitter truth: Industrial accidents are business as usual. Death on the job is a cost that we (as a society) are evidently willing to accept and pay. When a mine collapses or a sugar refinery explodes, we read the story in the newspaper and we watch the film at 11, and then we grit our teeth and move on. The Deepwater Horizon sinking is just another in this long series of mishaps.

And yet it’s also totally different. This accident has consequences that reach beyond the workers and their families. A hole in the bottom of the ocean has been spewing hydrocarbons for six weeks. The rate of loss could well be a million gallons a day. Efforts to control the well are still under way, but if they fail the spill might go on for another three months. All that makes for a nightmare that doesn’t drop out of the news cycle after the first 24 hours.

The Deepwater Horizon oil slick as seen by the NASA Terra satellite, 24 May 2010

In this respect, Deepwater Horizon may turn out to be not another Piper Alpha but another Three Mile Island. The 1979 accident at that nuclear power station had lasting consequences: We haven’t built another nuclear plant in the U.S. in the 30 years since. (There are other reasons for the long lassitude of the nuclear industry, but TMI was a major factor.) Likewise, in the aftermath of Deepwater Horizon, it seems a fair guess that we won’t be drilling another deep offshore well anywhere near the U.S. coastline for years to come, and maybe decades. Maybe never.

A moratorium on offshore drilling, enforced by economic as well as legal and political pressure, looks like a perfectly sensible response to the current situation. But we deserve more. By all means, let’s do whatever necessary to avoid repeating this particular mistake, but at the same time let’s address the broader issue of why things all over the industrial landscape keep blowing up and falling down.

Will we ever learn what happened on the Deepwater Horizon? There are certainly a lot of investigators trying to find out. Three Congressional committees have heard testimony already. A board of inquiry formed by the Coast Guard and the Minerals Management Service is holding its own series of hearings. So is a state panel in Louisiana. On May 11 President Obama asked the National Academy of Engineering to look into the cause of the accident, then on May 22 he appointed a special commission to carry out yet another investigation. Meanwhile the Justice Department is studying possible civil or criminal penalties.

For all that effort, we sure haven’t learned much yet. The well is leaking, but BP isn’t. It’s amazing how they’ve managed to keep the public and the press at a distance; the one art they seem to have mastered is secrecy. Even friendly trade publications (e.g., Oil and Gas Journal, Offshore) have so far failed to penetrate the security cordon. In the end, though, the story will come out. We’re going to learn why that blowout preventer didn’t prevent a blowout. But I’m not so confident we’ll learn how to prevent the next blowout.

Blowout stack beneath a Chesapeake Energy drilling rig, Marlow, Oklahoma, 2004.

The only blowout preventer I’ve ever seen up close and personal was on a drilling rig near Oklahoma City in 2004. At the time, I thought it was quite a brawny-looking piece of gear, with all those rings of torqued bolts clamping the flanges together, and the big red hydraulic rams like pincers clasping the pipe. Apparently, this unit is puny compared with the BOP stack installed on the Deepwater Horizon well (“five stories tall” in news accounts). Still, the principle of operation is the same: If the well “kicks”—meaning that gas and oil begin to push their way toward the surface—the rams close off the well and keep everything sealed tight. Some rams are meant to choke off the annular space around the drill pipe; some are “blind rams” used when no drill pipe is present in the bore; the last resort is a pipe ram or shear ram meant to crush or cut off the drill pipe.

BOP-controls-7980.jpg

It sounds foolproof, but evidently not. Controls can stick, valves can leak, parts can break. I’ve read speculation that the shear ram on the BOP beneath the Deepwater Horizon might have failed because it happened to strike a thicker-walled section of drill pipe where two lengths of tubing are threaded together; I really hope that’s not a serious concern, because it would mean that every BOP on the planet has a failure probability of 3 percent. (Drill pipe comes in 30-foot lengths, and the thickened sections are about a foot long.)

A 1997 article in Offshore lists a fascinating variety of other failure modes for deep-water wells. The article is by Larry H. Flak of Boots and Coots—a company often called in to deal with crises like this one. Flak is writing for an audience of oil-patch insiders, and I don’t follow all of his jargon. For example, there’s this wonderfully opaque passage:

Broached blowouts could happen with casing failure. Recently, an ultra-deepwater operator swabbed in a kick resulting in over 9,000 psi on the subsea BOPs. Fortunately, the casing was sound and set just on top of the sand. This allowed safe kick bullheading.

Even for those of us who have never bullheaded or swabbed in a kick, the message comes through: A lot can go wrong. A well is not just a hole in the ground; it’s a fairly complicated structure with multiple concentric layers of casing and drill pipe, which create several interconnected annular spaces; fluids under pressure can find many pathways to the surface. Flak points out that even the massive steel blades of a hydraulic ram can be torn to pieces by high-velocity streams of abrasive fluids, such as drilling mud. He concludes: “Blowout control options in ultra-deepwater are very limited. Blowout prevention is of paramount importance.”

When the story of the Deepwater Horizon accident is finally told, we’re going to hear at least two interpretations. In one telling, BP and its contractors were incompetent or negligent or criminally greedy; they screwed up. But the basic technology of offshore drilling is sound, and if only the companies had followed established industry practices and common sense, none of this mess would have happened. In the other version, BP and its contractors were incompetent or negligent or criminally greedy; they screwed up. But even if they had followed approved industry practices, a disaster like this was waiting to happen, because the technology of offshore drilling is fatally flawed. I really wish I knew which of these stories to believe. I worry that we lack the institutional means to decide between them.

•     •     •

I’m a great admirer of the National Transportation Safety Board. When an airliner crashes or a couple of commuter trains collide, the NTSB mounts a focused, scientific effort to understand the cause of the accident. Investigating accidents is all they do, and this turns out to be a key to their success; they are insulated from the responsibilities and temptations of regulating the industry, adjudicating legal culpability, assessing penalties or even enforcing their own safety recommendations. Separating the investigative and regulatory functions has worked brilliantly, and the board’s findings command wide respect. But the NTSB won’t be investigating the Deepwater Horizon blowout, because the event doesn’t fall within their statutory purview; it’s not a transportation accident. There’s also a Chemical Safety and Hazard Investigation Board (CSB), which has had a similar role in the chemicals industry since 1998. You might think that petroleum would count as a chemical substance, but as far as I can tell the CSB will also be sitting this one out. Their mission statement says: “The CSB conducts root cause investigations of chemical accidents at fixed industrial facilities.” The Deepwater Horizon wasn’t “fixed.”

So here’s my call for action: We need a single agency, modeled on the NTSB and the CSB, with authority to investigate all accidents involving industry or infrastructure—everything from well and refinery fires to water-main breaks, power blackouts, mine cave-ins, bridge and dam failures, boiler explosions and sewage spills. I believe we’re under a moral imperative to learn all we can from every such accident, in the hope of preventing a recurrence. A permanent, independent agency with investigative authority, in-house expertise and technical resources looks to me like the best way to reach this goal.

After the 2005 refinery fire in Texas City, BP hired James Baker III (the former Secretray of State) to convene a special of inquiry into the company’s safety practices and culture. The Baker panel report observed that BP had an admirable record on personal safety but cited deficiencies in process safety. Personal safety says: Wear your hardhat, your eye protection and your steel-toed boots when you go out in the plant to open a valve. Process safety says: Make sure you open the right valve, not the one that’s going to create a fireball that engulfs the whole plant. The report remarks: “BP mistakenly interpreted improving personal injury rates as an indication of acceptable process safety performance.” On April 20, BP executives flew out to the Deepwater Horizon to celebrate a personal safety milestone on the rig: no loss-of-work accidents in seven years. The executives were still on board when somebody turned the wrong valve.

Note: I have posted a PDF of the oil and gas chapter from my book Infrastructure: A Field Guide to the Industrial Landscape on the book’s web site.

Update 2010-06-21: A superb report on the failure modes of blowout preventers was published this morning in The New York Times. The article is by David Barstow, Laura Dodd, James Glanz,  Stephanie Saul and Ian Urbina. Nothing they say is reassuring.