Matrices and recurrences

Flights run between 21 airports: from Almaty to Astana 72 different callsigns fly, from Almaty to Balkhash one. How many different chains of exactly k flights lead from one airport to another, if each callsign counts as a separate way? For k = 1 the answer is read off the table of flights. For k = 2 it has to be gathered over all intermediate airports. And for k = 10 there are more chains than fit in a 64-bit number. Tap an airport and move k.

All these numbers are a row of one matrix: the table of flights C raised to the power k. From Almaty to Astana in two flights there are 1192 chains, and back to Almaty 7358: you can fly out and return through any city Almaty has flights to. The power is computed by fast exponentiation: C¹⁰ = C²·C⁸, and together with the squarings that is four matrix products instead of nine. Switch to “where it takes you”: a traveller each time leaves on a random flight. After the first flight they are certainly not in Almaty, after the second they are back there with probability 51%: almost all flights link cities with the two main airports, and the traveller swings between them. But the swinging dies down, and after a dozen flights it no longer matters where they started.

The idea

Let C_ij be the number of ways to make one step from i to j. Then the number of ways to arrive in two steps is the sum over the intermediate point t of the products C_it·C_tj, and that is exactly an entry of the matrix product C·C. In k steps, an entry of C^k. The step matrix may be a table of flights, a table of an automaton’s transitions or the rule by which the next term of a sequence is obtained from the previous ones; the power in O(log k) products is given by fast exponentiation, which does not care whether it multiplies numbers or matrices.

Adjacency matrices and matrix products are the chapters on graphs and on linear maps of “Mathematics, the Queen of the Sciences”, and powers through eigenvectors, Binet’s formula for the Fibonacci numbers and a Markov chain of rental scooters are the chapter on eigenvectors. Here: how to compute with it, which problems hide a matrix power, what it costs and when it is better not to raise it.

A recurrence is a matrix power too

In how many ways can a 3 × n strip be laid with 2 × 1 dominoes? Go along the strip from left to right, column by column. All you need to remember about the part already laid is which cells of the next column are already taken by lying dominoes sticking out of the previous one. That is a mask of three bits, eight states. The transfer matrix T says in how many ways a column can be filled when entered in state s and left in state t. A tiling of the strip is a walk of n steps through the states from a flat edge to a flat edge, and by the claim above there are (Tⁿ)₀₀ such walks.

2 x 10: 89 tilings, the matrix is 4 x 4
3 x 10: 571 tilings, the matrix is 8 x 8
4 x 10: 18061 tilings, the matrix is 16 x 16
8 x 8 : 12988816 tilings, the matrix is 256 x 256
4 x 10^18: 328515267 mod 1e9+7, 82 matrix products instead of 10^18 steps

A strip of width 2 is tiled in a Fibonacci number of ways: 89 for length 10, the very recurrence F_n = F_n−1 + F_n−2 that is usually written as a 2 × 2 matrix. Here nobody wrote it out: the transfer matrix built itself by trying the fillings of one column. A chessboard 8 × 8 is tiled in 12,988,816 ways, the number found by formula in 1961 by Pieter Kasteleyn and, independently, by Harold Temperley with Michael Fisher; the program gets it by multiplying 256 × 256 matrices. For length 10¹⁸, steps one at a time would take thousands of years even on a fast computer, while exponentiation makes 82 matrix products modulo a prime.

Any linear recurrence of order d, a_n = c₁a_n−1 + … + c_da_n−d, works the same way: the state is a vector of the d last terms, a step is a multiplication by a d × d matrix with the coefficients in the first row and ones below that shift the vector. The cost is O(d³ log n). For recurrences of order in the hundreds and thousands there is a faster way, through remainders of polynomial division, O(d² log n), and with fast polynomial multiplication less still; and if the recurrence is unknown but the first terms are known, it is recovered by the Berlekamp–Massey algorithm from the article on random number generators.

Words that are absent: an automaton instead of a graph

The restriction enzyme EcoRI cuts DNA where it meets the six letters GAATTC. How many such places should a sequence of 48,502 letters have if the letters are random? The mean is easy: each of the 48,497 positions starts a site with probability about 1/4096. But to know the probability that there are at most five sites, you have to count how many sequences avoid them, and these are walks in the graph of the Knuth–Morris–Pratt automaton. Its state is how many letters of the site have already matched; each letter moves the automaton to a new state, and the transitions are set by the same prefix function as in substring search. Add a counter of sites found to the state, weight the transitions by the probabilities of letters, and a row of the matrix to the power 48,502 is the distribution of the number of sites.

The genome of phage lambda, a virus that infects E. coli, contains GAATTC only five times, while a random sequence of the same length with the same letter frequencies has it 11.9 times on average; five or fewer comes up for it with probability about 2%. With the BamHI site, GGATCC, the picture is the same: five against 11.7. For human mitochondria the difference is within chance, and the small phage φX174 has no sites, but only one and a half would be expected. The shortage of restriction enzyme recognition sites in phage genomes was noticed long ago: in 1986 Paul Sharp wrote about it, explaining it by selection, since in bacteria such enzymes cut foreign DNA. But the rule is not universal, and the model “every letter is independent” is crude: in real genomes neighbouring letters are linked. The program computes the same for three genomes and two enzymes.

lambda phage         48502  EcoRI GAATTC: found 5, expected 11.86, P(none) = 6.96e-06, P(<= 5) = 0.0221
lambda phage         48502  BamHI GGATCC: found 5, expected 11.68, P(none) = 8.29e-06, P(<= 5) = 0.0246

Where it takes you: a Markov chain

Divide each row of the table of flights by its sum: P_ij is the probability that a traveller in i flies to j if they take any flight from the airport with equal chances. This is a Markov chain: where to fly next depends only on where they are now. The distribution after k flights is a row of P^k, computed by the same power or, if one starting point is needed, by k multiplications of a vector by the matrix. For our network these rows have a limit: Almaty 32.2%, Astana 22.8%, regardless of where the start is. And it almost coincides with the share of flights departing from the airport: 31.8% and 22.1%. This is no coincidence.

The real network is almost symmetric, sometimes with one or two more flights in one direction, so the limit is almost equal to the shares. More important is how fast it is approached. The distance to the limit (half the sum of the absolute differences of probabilities) is 0.33 after the first flight, 0.19 after the second, 0.11 after the third, each time about 0.57 times smaller, while the probability of Almaty jumps: 0, 51%, 22%, 38%, 29%. The factor 0.568 (the program below prints it more precisely) is the absolute value of the second-largest eigenvalue of P, and the fact that this value is itself negative produces the swinging: the network resembles a star where flights go from a big city to small ones, and from small ones back to big ones. On eigenvalues and the power method, see the chapter on eigenvectors and the article on PageRank, where the same limit of a Markov chain ranks pages by importance.

Code: chains and probabilities

 k  to Astana                       back to Almaty
 1  72                            0
 2  1192                          7358
 3  599774                        181406
 4  23398952                      61553540
10  38516548485362965772          49395526650223704252
24  a number of 48 digits      48 digits
from k = 10 some counts no longer fit in a signed 64-bit integer

The number of chains grows as a power of the largest eigenvalue of C, almost a hundredfold with each flight, and with ten flights it no longer fits in 64 bits. Python computes it with big integers by itself, Go needs math/big for it, and contest problems ask for the answer modulo a number: then every addition and multiplication is taken modulo it, and the numbers stay small. The product needs the most care: two remainders up to 10⁹ give up to 10¹⁸, which still fits in a signed 64-bit integer, but a sum of several such products no longer does, so the remainder is taken after each term.

The second half of the program is the random traveller: the probability vector is multiplied by P flight after flight, and the limit is obtained by two thousand such multiplications.

12  0.323   0.227   0.00066
each flight multiplies the distance by about 0.568
the limit: Almaty 0.322, Astana 0.228; their shares of all callsigns: 0.318, 0.221

What it costs

The product of two m × m matrices is m³ multiplications, and the power k takes at most 2 log₂k products: O(m³ log k). But often what is needed is not the whole matrix but one row: chains from one airport, tilings from a flat edge. Then you can go in steps, a vector times a matrix, m² multiplications per step, k·m² in all. Steps pay off while k is below roughly 2m·log₂k. A measurement on the network of 21 airports, modulo 10⁹ + 7:

Steps grow in direct proportion to k, exponentiation as its logarithm. In Python the lines cross between k = 64 and 256, in Go around 256. With k = 10¹⁸ there is no choice: only the power. But if the matrix is large and k modest, as in graph problems with thousands of vertices, exponentiation loses in both time and memory: m³ for m = 10⁴ is 10¹² multiplications per product.

One more cost is the size of the numbers. Without a modulus the entries of C^k grow by a constant number of digits per step and have to be multiplied with big-integer arithmetic; with a modulus each number takes one machine word, but the answer is only a remainder.

Where it works

• Contest problems. “Compute modulo 10⁹ + 7 for n up to 10¹⁸” almost always means a linear recurrence, walks in a small graph or words without forbidden pieces, that is, a matrix power. Both problems below are from this family.
• Markov chains. Models of queues, weather, shopper behaviour, probabilistic models of text; their distribution after k steps is a row of P^k, and the limit is sought by the power method, as in PageRank.
• Combinatorics and physics. The transfer matrix method counts tilings, colourings of strips and partitions, and in statistical physics sums over all states of chains of atoms.
• Bioinformatics. How many times a motif should occur in a random sequence and how surprising it is that it occurred less often: the same automaton matrix raised to the length of the sequence.

History

In 1913 Andrey Markov split the first 20,000 letters of Eugene Onegin into vowels and consonants and counted how often a vowel follows a vowel and a consonant a consonant. It turned out that letters depend on the previous ones, and yet the share of vowels holds steadily at about 43%: the law of large numbers works for dependent events too, if the dependence is arranged as a chain. So appeared the first example, worked out on real data, of the chain that now bears his name. The transfer matrix method came from physics: it computes sums over the states of chains of atoms, and in 1961 two ways to count the domino tilings of a rectangle, Kasteleyn’s and Temperley and Fisher’s, made the chessboard problem a classic. The trick “raise the step matrix to a power” itself is attributed to nobody: it is fast exponentiation applied to the associative multiplication of matrices.

Problems

In the first problem you have to count chains of flights in a real network over a huge number of flights. In the second, how many DNA sequences of a given length never contain the site an enzyme cuts, with an automaton and a power of its matrix.