Combinatorics

Balanced Incomplete Block Designs

Posted on by John

Suppose you have n different fertilizers you want to test. The simplest experiment would be to take a big plot of land, divide it into n sections, and apply each fertilizer to a section. That’s fine if you have one big enough plot of land, and the land itself doesn’t matter.

If the fertilizer supplies a nutrient that is already abundant in a particular soil, you won’t see as much benefit from that fertilizer on that soil. Plots of land will vary in other ways besides soil composition—some may receive more sun than others, some may have better drainage, etc.—that may interact with the effect of the fertilizer. So you’d like to test your fertilizers on several different plots of land. Call the number of plots b.

Suppose each of your plots is large enough to test q < n fertilizers, i.e. the plots are not large enough to test every fertilizer on every plot. You would like to design your experiment so that each fertilizer is tested the same number of times, and so that every pair of fertilizers is tested somewhere on the same plot of land for a direct comparison.

The resulting experimental design is called a balanced incomplete block design or BIBD. The blocks are the plots of land. The design is incomplete because you cannot test everything on every plot. The design is balanced because every pair of fertilizers is tested together the same number of times. BIBDs are important in the design of experiments.

A BIBD is a special case of the combinatorial designs discussed in the previous post with parameters (n, q, 2, λ), i.e. the parameter r is set to 2 because you want every pair to appear in some number of blocks.

A design with parameters (n, q, r, λ) is sometimes called an r-(n, q, λ) design [1]. Steiner systems correspond to λ = 1 and BIBDs correspond to r = 2.

While the number of blocks b is of interest, it is not an independent parameter. For combinatorial designs in general the number of blocks is determined by the other parameters:

b = \lambda \frac{\binom{n}{r}}{\binom{q}{r}}

Similarly, you may see the number of blocks each treatment appears in given as a parameter, but this count is also determined by the other parameters, namely

\lambda \frac{\binom{n-1}{r-1}}{\binom{q-1}{r-1}}

Examples

Here is a design to test n = 6 treatments on b = 10 plots, testing q = 3 treatments on each plot, with each pair of treatments (i.e. r = 2) appearing in λ = 2 plots, i.e. a (6, 3, 2) design.

    0000011122
    1123423433
    2345554545

The way to read this is that each treatment is a number 0 through 5, and each column indicates the treatments tested in a particular block (plot).

Here is another example, an (8, 4, 3) design: n = 8 treatments, blocks of size q = 4, with each pair of treatments being tested in λ = 3 blocks.

    00000001111222
    11123342334334
    22554666455455
    34675777677766

[1] Combinatorial designs are sometimes called t-designs, which makes sense if you denote the size of the set you’re looking for in each block by t but not if you denote it by r. You’ll see both conventions. And while we’re mentioning notational variations, the number of treatments is sometimes denoted v for variety, going back to agricultural applications, testing crop varieties.

Combinatorial designs

Posted on by John

The last several posts have been building up to this one, moving from tangible to more abstract.

The posts on finite geometry gave examples of Steiner triples and more generally Steiner systems, introduced in the previous post. Combinatorial designs take this one step further.

The previous post gave the following definition:

A Steiner triple system, written STS(n), is a collection of three-element subsets, called blocks, of a set X with n elements, such that every two-element subset of X belongs to exactly one block.

The post pointed out that the smallest projective plane is STS(7) and the smallest affine plane is STS(9).

The post then defined a general Steiner system:

A Steiner system S(t, k, n) is a collection of k-element subsets, called blocks, of a set X with n elements, such that every t-element subset of X belongs to exactly one block.

Moving from Steiner triples to Steiner systems we generalized the size of the blocks from 3 to k, and generalized the size of sets contained in a unique block from 2 to t.

Combinatorial designs take this one step further by generalizing the number of blocks that a set of t elements falls in from 1 to λ.

Unfortunately when we go from Steiner systems to combinatorial designs the notation changes and the order of the parameters changes.

A design with parameters (n, q, r, λ) is a collection of q-element subsets of an n-element set X such that every r-element subset of X is contained in exactly λ blocks.

So a Steiner system S(t, k, n) corresponds to a (n, k, t, 1) design.

Existence of designs

Given a set of parameters (n, q, r, λ) there may or may not exist a design with these parameters. It’s easy to show that a few consistency conditions are necessary. It was long conjectured that these conditions were sufficient, except possibly for a finite number of values of n. Peter Keevash proved in 2024 that the conjecture is correct.

Theorem [1]: Fix q > r ≥ 1 and λ ≥ 1. There is an N, depending on q, r, and λ, such that for every n > N with

\left. \binom{q-i}{r-i} \;\middle|\; \lambda\binom{n-i}{r-i} \right.
for all 0 ≤ i ≤ r, a design with parameters (n, q, r, λ) exists.

[1] Peter Keevash. The existence of designs. https://arxiv.org/abs/1401.3665v4

Steiner Systems

Posted on by John

Steiner triples

Several examples from the recent posts on finite geometry are special cases of Steiner triples.

A Steiner triple system, written STS(n), is a collection of three-element subsets, called blocks, of a set X with n elements, such that every two-element subset of X belongs to exactly one block.

STS(7) is the smallest finite projective plane, the Fano plane, and STS(9) is the smallest finite affine plane. In both cases you interpret the elements of S as points and the blocks as lines.

Several questions come up naturally. First, for which values of n does an STS(n) exist? A necessary and sufficient condition for existence is that n = 1 or 3 mod 6.

For a given n, can there be more than one STS(n), i.e. are there Steiner triple systems of the same size that are not isomorphic? Yes indeed. The number of non-isomorphic instances grows rapidly with n.

When does (some instance of) STS(n) correspond to a projective or affine space? When n = 2k − 1 for k ≥ 3. For n ≥ 13 there are multiple isomorphism classes of STS(n), so for k > 3 there will be an instance of STS(2k − 1) that corresponds to a projective space and other instances that do not. The instance that does correspond to a projective plane will be isomorphic to a finite field analog of the real projective plane.

Steiner systems

Steiner triple systems are a special case of Steiner systems.

A Steiner system S(t, k, n) is a collection of k-element subsets, called blocks, of a set X with n elements, such that every t-element subset of X belongs to exactly one block. Steiner triples are S(2, 3, n) systems.

Steiner systems often correspond to finite geometries with the blocks interpreted as lines. As with the special case of triples, there can be non-isomorphic instances of Steiner systems, and so when I say below that S(t, k, n) corresponds to a certain geometry, I mean that at least one instance corresponds to that geometry.

Projective

The Steiner system

S(2, q + 1, (qd + 1 − 1)/(q − 1))

with q a prime power greater than 2 is a d-dimensional projective space.

Affine

The space

S(2, q, qd)

is an affine space. Again q > 2 is a prime power.

Hyperbolic

The finite hyperbolic plane given as an example in an earlier post corresponds to S(2, 2, 6). Many other Steiner systems are hyperbolic, but I’m not aware of a tidy classification similar to the ones above.

Steiner systems outside geometry

Steiner systems were named after Swiss mathematician Jakob Steiner because he published a highly cited article asking questions about what we now call Steiner triples. His interest in these systems was combinatorial. Others developed the connections to geometry later.

Steiner systems arise in several areas of mathematics outside of geometry, including coding theory, design of experiments, and graph theory.

Harmonic triangle

Posted on by John

Pascal’s triangle is the most famous triangular arrangement of numbers. I’ve written about other numerical triangles, such as Bell’s triangle and the trinomial triangle. There’s another well-known triangle that somehow I never heard of until recently: the harmonic triangle, first explored by Leibniz around 1700.

To create the harmonic triangle, start with the harmonic series

1, 1/2, 1/3, 1/4, 1/5, …

This is the first column of the triangle. Next take the differences of the terms:

(1 − 1/2), (1/2 − 1/3), (1/3 − 1/4), … = 1/2, 1/6, 1/12, …

This, shifted down one space, is the second column of the triangle. Taking differences again gives the third column, and so forth. Column c of the triangle consists of the cth differences of the harmonic series.

    1
    1/2  1/2
    1/3  1/6   1/3
    1/4  1/12  1/12  1/4
    1/5  1/20  1/30  1/20  1/5
    …

Each entry is the difference of the entries to its west and northwest. If we let L(r, c) be the element in the rth row and cth column, then

L(r,c) = L(r-1, c-1) - L(r, c - 1)

Every numerator is 1, though it’s not immediately obvious why this should be the case. And every number is a value of the beta function, the function that gives the normalizing constant for the beta distribution [1].

L(r,c) = \frac{1}{r\binom{r-1}{c-1}} = \beta(c, r-c+1)

Related posts

[1] Ayoub B. Ayoub. The Harmonic Triangle and the Beta Function. Mathematics Magazine, Vol. 60, No. 4 (Oct., 1987), pp. 223-225

Proportion of 1s in a Hadamard matrix

Posted on by John

The first post in the recent series of posts on Hadamard matrices describes a way of constructing new Hadamard matrices from two other Hadamard matrices by taking their Kronecker product.

Starting with a Hadamard matrix H0 and a Hadamard matrix G, you can construct a sequence of Hadamard matrices by

Hn+1 = G ⊗ Hn

for  positive integers n. This is known as the generalized Sylvester method.

Let pn be the proportion of 1s in Hn and let q be the proportion of 1s in G. Then you can show that the recurrence holds

pn+1 = q pn + (1 − q)(1 − pn).

You can solve the recurrence to show that

limn → ∞ pn = ½

and so as the iterations proceed, the ratio of number of 1s to the number of −1s approaches 1.

This doesn’t say anything Hadamard matrices in general, but it does apply to all Hadamard matrices created by repeatedly applying the generalized Sylvester method.

If you set G and H equal to the matrix

\begin{bmatrix} 1 & 1\\ 1 & -1 \end{bmatrix}

then p0 = q = ¾. Then for n = 1, 2, 3, …, 8 the values of pn are

0.625
0.5625
0.53125
0.515625
0.5078125
0.50390625
0.501953125
0.5009765625.

 

Compressing a Hadamard matrix

Posted on by John

Hadamard matrices are in the news following the recent announcement of a newly discovered Hadamard matrix. I’ve written three posts on Hadamard matrices recently, one as a sort of introduction and two on applications: the error correcting code used in the Mariner 9 probe and constructing sphere packings.

A Hadamard matrix is an orthogonal matrix with all entries equal to ±1. Jacques Hadamard conjectured that there exist Hadamard matrices of order 4n for all positive integers n. It’s necessary that the order be divisible by 4, and Hadamard conjectured that this is sufficient [1].

How could you compactly represent a Hadamard matrix? Since the entries are all either 1 or − 1 each entry could be represented by a single bit, and n² bits could store an n × n Hadamard matrix. But we can do better.

Methodical matrices

If the matrix can be produced by an algorithm, you only need to store the name of the algorithm and the argument to the algorithm. So, for a 1024 × 1024 matrix applied by iterating Sylvester’s algorithm could be stored by saying “Apply Sylvester’s algorithm 10 times” rather than storing a megabyte of data.

Paley’s method can create a Hadamard matrix corresponding to every prime power. So you could determine a Paley type matrix by storing the prime and the exponent.

Next in complexity would be hybrid algorithms, such as start with the Paley method applied to 376 and then apply Sylvester’s method 3 times.

There are more methods of creating Hadamard matrices than Sylvester’s method and Paley’s method, though they’re harder to describe and parameterize.

Sporadic matrices

If a Hadamard matrix cannot be constructed using an algorithm, you can still store the matrix in fewer than n² bits. Since the rows are orthogonal, the last row of the matrix is determined by all the previous rows, up to sign. So you could store a Hadamard matrix using n(n − 1) + 1 bits.

Some Hadamard matrices are symmetric or skew. A symmetric matrix is determined by its diagonal and the elements above the diagonal. So a symmetric Hadamard matrix could be represented by n(n + 1)/2 bits.

A skew Hadamard matrix isn’t quite skew-symmetric. A matrix M is skew symmetric if

MT = −M.

This implies the diagonal elements are 0, and Hadamard matrices cannot contain 0s. A Hadamard matrix H is called skew if

H + HT = 2I.

This implies the diagonal elements are all 1s and the elements below the diagonal have the opposite sign of the elements above the diagonal. Since the elements on the diagonal are determined, a skew Hadamard matrix can be sorted using n(n − 1)/2 bits.

Incidentally, there is a conjecture that there exist skew Hadamard matrices of order 4n for all positive n.

 

[1] There are Hadamard matrices of order 1 and 2, but larger orders must be divisible by 4.

Constructing Hadamard matrices

Posted on by John

A Hadamard matrix is an orthogonal matrix whose entries are all either 1 or − 1. For example

\begin{bmatrix} 1 & 1\\ 1 & -1 \end{bmatrix}

is a Hadamard matrix of order 2. True to Stigler’s law of eponymy, James Joseph Sylvester investigated Hadamard matrices before Jacques Hadamard. Sylvester saw how to bootstrap the example above into more examples. If H is a Hadamard matrix, then the partitioned matrix

\begin{bmatrix} H & H\\ H & -H \end{bmatrix}

Sylvester’s construction can be generalized as follows. If Hm is a Hadamard matrix of order m and Hn is a Hadamard matrix of order n, the the Kronecker product Hm⊗Hn is a Hadamard matrix of order mn. That is, you can form a new Hadamard matrix by taking the matrix Hm and replacing ±1 with the matrix ±Hn.

Let S be the set of all possible Hadamard matrix orders. By the construction above, this set is closed under multiplication. Since 2 is in S, every power of 2 is in S. Hadamard proved that all n ≥ 4 in S are multiples of 4. That is, the condition 4 | n is necessary. He conjectured that it was also sufficient, though that has not been proven.

So the big question is what is the set S. Is there some multiple of 4 not in S? Until that question is answered, what is the smallest multiple of 4 not known to be in S? Hadamard matrices are useful in applications, so constructing Hadamard matrices of various orders is useful even while Hadamard’s conjecture remains open. For example, see the next post for how NASA used Hadamard matrices to transmit photographic images back from Mars.

Paley’s method

Raymond Paley came up with a way of constructing Hadamard matrices of size q + 1 if q is a prime power congruent to 3 mod 4, and of size 2(q + 1) if q is a prime power congruent to 1 mod 4. Let’s see what we can squeeze out of this.

If p is a prime congruent to 1 mod 4, every power of p is also congruent to 1 mod 4, and so there exist Hadamard matrices of order 2(pk + 1) for every k.

If p is a prime with p = 3 mod 4, then even powers of p are congruent to 1 mod 4 and odd powers of p are congruent to 3 mod 4. So there are Hadamard matrices of order 2(p2k + 1) and of order p2k+1 + 1.

Let’s run a script to see what we can learn from this.

from sympy import primerange

s = set()

for p in primerange(20):
    if p % 4 == 1:
        s.update([2*(p**k + 1) for k in range(1, 10)])
    if p % 4 == 3:
        s.update([2*(p**(2*k) + 1) for k in range(1, 6)])
        s.update([p**(2*k + 1) + 1 for k in range(1, 6)])
print(sorted(s)[:20])

This prints

12, 20, 28, 36, 52, 100, 164, 244, 252, 340]

We can add 16 to the list because it’s a power of 2, and we can add 24 because it’s 2 × 12, etc. But there doesn’t seem to be any way to get 44. There is a way to create a Hadamard matrix of order 44, but it doesn’t follow from anything we’ve seen so far.

New records

I have a book published in 1996 that says Hadamard’s conjecture had been verified for n up to 428. Until yesterday, the smallest multiple of 4 for which nobody had found a corresponding Hadamard matrix was 668. Then Levent Alpöge announced that he and his and collaborators found an example of size 668 and filled in all remaining gaps below 2000.

So now the set S is known to contain {1, 2, 4, 8, 12, 16, …, 2000}. It also contains all orders that can be obtained by Paley’s method and other methods. And it contains all products of its elements. But it is not yet known to contain 2004.

Related posts

Counting rooted trees

Posted on by John

Combinatorial problems can be interesting for their own sake, but they are more interesting when there is a connection to a problem outside combinatorics, and the more unexpected the connection the better.

Counting the number of unlabeled rooted trees [1] with n nodes is a pure mathematics problem. Designing numerical methods for solving differential equations is an applied mathematics problem. And yet the two are closely linked.

Let t(n) be the number of distinct unlabeled rooted trees with n nodes. The diagram below shows that the first few terms of this sequence are 1, 1, 2, and 4.

Recursive calculation

The values of t(n) can be computed recursively using

\begin{align*} g_k &= \sum_{d\mid k} d t_d \\ t_1 &= 1 \\ t_n &= \frac{1}{n-1} \sum_{k=1}^{n-1} g_k t_{n-k} \text{ for } n > 1<br /> \end{align*}<br />

You can implement this in Python as follows.

from sympy import divisors

def t(n):
    if n <= 1:
        return 1 if n == 1 else 0
    return sum(g(k) * t(n - k) for k in range(1, n)) // (n - 1)

def g(k):
    return sum(d * t(d) for d in divisors(k))

This code is correct, but it will run more efficiently if you cache function values to avoid calculating the same values over and over. You can do this by adding

from functools import lru_cache

and writing @lru_cache(maxsize=None) above both function definitions.

Connection to Runge-Kutta

In an earlier post I showed that designing a 4-stage explicit Runge-Kutta method required solving a system of 8 equations in 10 unknowns, leaving two degrees of freedom in the solutions.

The number of constraints c(s) needed to design an s-stage explicit RK method is equal to the number of rooted trees with up to s nodes:

c(s) = t(1) + t(2) + t(3) + … + t(s)

This is because there is a one-to-one correspondence between constraints on the nth derivative of an RK formula and rooted trees, and an s stage method has to satisfy the constraints of all stages up to s. In the example of the 4th order RK method, we have

c(4) = t(1) + t(2)  + t(3) + t(4) = 1 + 1 + 2 + 4 = 8.

The first few values [2] of t(n) are

1, 1, 2, 4, 9, 20, 48, 115, 286, 719, 1842, 4766, 12486, 32973, …

and so you can see that t(n) grows quickly. In fact, it grows exponentially [3].

However, the number of parameters in an s stage RK method is s(s + 1)/2. The number of equations grows exponentially and the number of variables grows only quadratically, so at some point you have more equations than variables. That’s already the case for s = 5 because you have 17 constraints on 15 variables. The system has a solution because symmetry considerations render some of the equations redundant.

A 10th order RK method requires 17 stages. (See the previous post for why the number of stages exceeds the order when the order is greater than 4.) Designing such a method would require solving over a million equations in 153 variables, and yet it can be done. [4]

Related posts

[1] This is a slightly contradictory term. Unlabeled means that we don’t distinguish the nodes. But we do distinguish one node, namely the root.

[2] See OEIS A000081.

[2] Richard Otter proved in 1948 that the number of unlabeled rooted trees with n nodes is asymptotically C αn / n−3/2 where C = 0.4399… and α = 2.9557…. The cumulative sum is at least this large since Otter’s estimate gives the size of the last term in the sum.

[3] E. Hairer. A Runge-Kutta Method of Order 10. J. Inst. Maths Applics (1978) 21, 47-59

Cryptographic Keys and Decks of Cards

Posted on by John

The previous post looked at the idea of storing a cryptographic key in the order of a deck of cards. A deck of 52 cards can store 225 bits of data because

⌊log2(52!)⌋ = 225.

Here ⌊x⌋ is x rounded down to the nearest integer.

If we want to store bigger keys, we’re going to need a bigger deck of cards.

Bitcoin

A Bitcoin key has 256 bits, which would require a deck of 58 cards. There is a card game called Zwicker that uses a deck of 58 cards, the usual 52 cards plus six jokers. So you could store a Bitcoin key in the permutation of a Zwicker deck.

You could also use a deck of 52 cards, plus 2 jokers, if you also consider orientation. 30 cards are rotationally symmetric, 22 are not, and neither are jokers. So, including two asymmetric jokers, you could add 24 additional bits. Permutations of a 54 card deck can encode 237 bits, and with 24 orientation bits, this is a total of 261 bits.

RSA

RSA key sizes vary, but 2048 and 3072 are common. A 2048-bit key would require a deck of 301 cards. Casinos often use a shoe of 312 cards, combining six decks of 52 cards, to deal Baccarat or Blackjack. However, casinos combine identical decks. If you were to combine six unique decks, you could store a 2048-bit key.

Storing a 3072-bit key would require a deck of 422 cards. You could make a deck of 432 cards by combining 8 distinguishable packs of 54 cards (52 + 2 jokers).

ML-KEM

ML-KEM is a proposed quantum-resistant replacement for RSA. As with RSA, key sizes for ML-KEM vary, the smallest being ML-KEM-512 with a key size of 1632 bytes, which equals 13056 bits. This would require a deck of 1442 cards. You could combine 28 distinct packs of 52 cards, but that’s unwieldy.

This illustrates one of the difficult trade-offs with post-quantum cryptography: key sizes are much bigger. If you wanted to create a deck of 1442 cards, you’d probably want to make your “cards” something other than standard playing cards. You’d want to use permutations of something else.

Verification

The following Python code verifies the calculations above.

from math import log2, factorial, floor

def capacity(cards):
    return floor(log2(factorial(cards)))

def verify(bits, cards):
    return capacity(cards) >= bits and capacity(cards-1) < bits

print(verify(237, 54))
print(verify(256, 58))
print(verify(2048, 301))
print(verify(3072, 422))
print(verify(1632*8, 1442))

For more on how I came up with the deck sizes, see the next post on computing the inverse factorial.

Hiding data in permutations

Posted on by John

The latest issue of Paged Out! has an article by Stephen Hewitt “An off-line backup of your cryptographic key using playing cards.” The idea is to use a deck of 52 to store a 128-bit cryptographic key. To erase the key, shuffle the deck. Hewitt gives his algorithm for embedding a key, one that can be carried out manually but isn’t maximally efficient.

You could store a 225-bit key as a permutation of 52 cards because

log2(52!) = 225.581.

But then how would you number permutations so you could go from a number to a particular permutation and later decode the permutation to a number? Is this even practical? For a small number n, you could encode a number k < n by enumerating the first k permutations of a set of n items, and you could decode by enumerating permutations until you find the one you have. But this is completely impractical for large n, such as n = 52.

The process of mapping permutation to an integer is called ranking, and the mapping from an integer to a permutation is called unranking. How efficiently can rankings and unrankings be calculated?

Let n be the number of symbols being permuted. Then there are simple algorithms for ranking and unranking with respect to lexicographical order that have complexity O(n²) and more sophisticated algorithms that have complexity O(n log n). There are also O(n) algorithms that do not preserve lexicographical order.

The Permutations class in SymPy has methods unrank_lex and rank to unrank and rank permutations according to lexicographical order.

The notation the Permutations class uses requires a little explanation. For example, suppose we unrank 2026.

>>> from sympy.combinatorics import Permutation
>>> Permutation.unrank_lex(52, 2026)
Permutation(45, 47, 51, 48, 46, 50)

The output is not a full list of 52 numbers in permuted order; it is only a cycle. The notation refers to the permutation that sends 45 to 47, 47 to 51, …, 50 to 45 and leaves everything else fixed.

If we rank the permutation given above, we get 2026 back.

>>> Permutation.rank(Permutation(45, 47, 51, 48, 46, 50))
2026

Note that we didn’t say how many elements (45, 47, 51, 48, 46, 50) is a permutation of. Because of lexicographical order, the rank would be the same whether we viewed this as a permutation of 52 objects or of more objects.

Now let’s do something larger. Let’s generate a 220-bit number and encode it as a permutation.

>>> n = random.getrandbits(225)
>>> a = Permutation.unrank_lex(52, n)
>>> n
40234719030664563684489051530416964877785781669439875437823431388841
>>> a
Permutation(0, 25, 32, 15, 8, 28)(1, 48, 34, 14, 10, 51, 38, 31, 21, 5, 42, 47, 29, 26, 46, 30, 50, 49, 37, 22, 18, 23)(2, 45, 17, 20, 36, 40, 11, 4, 7, 41, 33, 3, 43, 44, 19, 16, 35, 39, 12, 6, 9)
>>> Permutation.rank(a) == n
True

Now just for fun, let’s display the permutation above applied to a standard (French) deck of 52 cards. As explained here, symbols associated with these cards have a range of Unicode values. By printing these values, we can visualize the permuted deck.

Here’s the code that made the image above.

spades = list(range(0x1F0A1, 0x1F0AF))
spades.remove(0x1F0AC) # take out the knight
cards = [s + 16*i for s in spades for i in range(4)]

a = Permutation.unrank_lex(52, n)
p = a(cards)

for i in range(4):
    for j in range(13):
        print(chr(p[13*i + j]), end="")
    print()

The code above is plenty fast, but Permutation has methods rank_nonlex and unrank_nonlex that run in O(n) time, which could be useful for n much larger than 52.