2018-02-05.ipynb

Kernel: SageMath 8.1

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

February 5, 2018: Combinatorics (part 1 of 3)

Administrivia:

• We are continuing to monitor the course waitlist. If you are on it and still want to join the course, please continue to keep up with lectures and homework, and watch your email for further instructions.

• Here (again) is the precise link for the documentary screening later today.

• There was an error in the HW 2 solution set (now fixed) which may have affected the grading of problem 3a. If you believe you were affected, contact course staff.

Added in class:

• HW 4 is now posted.

Postscript: algorithms for determinants

How does a computer algebra system compute the determinant of an $n \times n$ matrix $A$? Let's consider some options.

• Recall that there is an explicit formula $$\det(A) = \sum_{\sigma \in S_n} (-1)^{s(\sigma)} A_{1\sigma(1)} A_{2\sigma(2)} \cdots A_{n \sigma(n)}$$ where $\sigma$ runs over all permutations of $\{1,\dots,n\}$ and $s(\sigma)$ is the sign of $\sigma$ (if you factor $\sigma$ as a sequence of transpositions, the sign is +1 or -1 depending on whether this sequence has even or odd length). The sum involves $n!$ terms, which means that already for $n=100$ this is unusable.

• There is a recursive algorithm given by expansion in minors. For example, if I expand in minors along the first row, I get $$\det(A) = \sum_{i=1}^n (-1)^i A_{1 i}\det(B_i)$$ where $B_i$ is the submatrix obtained by removing row 1 and column $i$. However, once you unwind the recursion, this also involves $O(n!)$ computations.

• The most reasonable generic method is Gaussian elimination, which requires $O(n^3)$ arithmetic operations.

• For very large matrices (say n > 50000), there are some more advanced methods that do a bit better, e.g., ones which incorporate Strassen's algorithm for matrix multiplication.

However, as you learned on the last problem set, Sage does not use Gaussian elimination by itself to do this kind of computation.

A = random_matrix(ZZ, 1000, x=-10, y=10);  # random 1000x1000 matrix over ZZ with entries in [-10..10]
%time A.det()   # should take around 10 seconds!

CPU times: user 11.4 s, sys: 103 ms, total: 11.5 s
Wall time: 14 s

-2849968122179466491223698417283711974654232562105638734135332928564615805911430153031967249575087725892310490051042419471717894992961842001133415054034911330822286669353386639811952127516767300995155636837047585686502163841436856805724169807621473993572911899933436702420366693356236453089979256580760853492429812813760531854473545059513417880393740455219897267424285863433965173236885194738566888569932655683273797930395898602289852551154772348233601169483056081345162609676576815777486604892784764781640606200833689263971844928574269297522632738182318760922390322455043624625953680956044526741416404441521711047221342112239739534603859288652897834288678325753080689615125253049263901393370046748145364643195483355613302737670336301176329863041771721960306372831552120826495631014979934545248269998741818477982505334609095865241562192529191583467579422971576053646816534675389032764172210902445701532474499220639052741868842030190513252742426969767065635995077175382925105237545362901650409883374958318141287486467606057572100726039023391395750488529842890708784251016971030391039486318334724577779692184388629878954890479105626929326783236764512408094476057444651336753518657899579298264910854867489367784231375510532464745046673463016300700436964001218760290685428526512059714444207769629782232983409781010667672044894201263196138112856900917569444809455865861831433777522630685185857871752888069849407734076755359774619645093924520445508817088202528511890262931852472284871432296782646505657304548807234304097977512764282538877423604241687510154601108761594268029935202732162070593736948424000674645198671134288975698820810101314121312477742667845922714206611973369898660294603878343726897738220383846990786815286278967782668014939556926541006505156788821087456211888502368849579215846030160256233409907984503398777118633819194620392733459614095899678993551738510667829919455440137138908901786211378757087315635233741232560329331903409923637693029107704672610059596015776176597670334760483634478515270855631966537980421397088023500941455872

The reason is that the use of Gaussian elimination creates some numbers in the middle of the computation which are about twice as long as this. To avoid this, Sage uses a bit of elementary number theory, in the form of the Chinese remainder theorem. (I'll demonstrate this for determinants, but similar considerations apply to other complex operations in linear algebra: reduced row echelon forms, matrix inverses, etc.)

n = 5
set_random_seed(42)
A = random_matrix(QQ, n, n, num_bound=10, den_bound=1)
show(A)

Let's step through Gaussian elimination in this example.

for i in range(1,n):
    A.add_multiple_of_row(i, 0, -A[i,0]/A[0,0])
show(A)

for i in range(2,n):
    A.add_multiple_of_row(i, 1, -A[i,1]/A[1,1])
show(A)

Skipping ahead...

for j in range(2,n):
    for i in range(j+1, n):
        A.add_multiple_of_row(i, j, -A[i,j]/A[j,j])
show(A)

det(A)

37033

One way to get around this would be to be more careful about division, to get rid of the denominators. However, this requires some care and is not quite the best approach.

Alternate strategy:

1. Choose some random prime $p$.

2. Reduce our matrix modulo a prime number $p$ to get $A\pmod{p}$.

3. Perform row reduction to compute $\det(A) \pmod{p}$.

4. (*) Lift the result back to the rational numbers to get the right answer.

This is potentially better, since until the last step we only deal with integers from 0 to $p-1$.

n = 5
set_random_seed(42)
A = random_matrix(QQ, n, n, num_bound=10, den_bound=1)
show(A)

next_prime(500)

503

previous_prime(500)

499

B = A.change_ring(GF(389))
show(B)
B.parent()

Full MatrixSpace of 5 by 5 dense matrices over Ring of integers modulo 389

for i in range(1,n):
    B.add_multiple_of_row(i, 0, -B[i,0]/B[0,0])
show(B)

for i in range(2,n):
    B.add_multiple_of_row(i, 1, -B[i,1]/B[1,1])
show(B)

for j in range(2, n):
    for i in range(j+1, n):
        B.add_multiple_of_row(i, j, -B[i,j]/B[j,j])
show(B)

det(B)

78

Notice that the numbers are not exploding in size. Imagine doing all this with a 200x200 matrix (or worse)!

Unfortunately, the prime $p=389$ was not big enough: this computation tells me that $\det(A) \equiv 78 \pmod{389}$ but this is not enough to deduce that $\det(A) = 37033$. To fix this, I need to compute a bound $N$ on $|\det(A)|$ (e.g., using the explicit formula in terms of permutations) and then take $p > 2|N|$. (Why does this suffice?)

This creates its own difficulties, because if $p$ is really large, then we again end up working with large integers and this creates computational overhead. (CPUs can typically only handle 64 bits at a time, so large integers have to be broken into smaller chunks in order to deal with them. Think about the effect on multiplication, for example.)

Enter the Chinese remainder theorem. Let's instead choose primes $p_1,\dots,p_k$ such that $p_1 \cdots p_k > 2 |N|$; we can then compute $\det(A) \pmod{p_i}$ for each $i$ and use this knowledge to reconstruct $\det(A)$.

show(5^(2.5)*10^5) # Hadamard bound

l = prime_range(25)
show(l)
show(prod(l))

m = [det(A.change_ring(GF(p))) for p in l]
show(m)

crt?

crt(m, l) ## This won't work

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-59-8883128bdf31> in <module>()
----> 1 crt(m, l) ## This won't work

/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/arith/misc.pyc in crt(a, b, m, n)
   2756     """
   2757     if isinstance(a, list):
-> 2758         return CRT_list(a, b)
   2759     if isinstance(a, integer_types):
   2760         a = Integer(a) # otherwise we get an error at (b-a).quo_rem(g)
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/arith/misc.pyc in CRT_list(v, moduli)
   2841     from sage.arith.functions import lcm
   2842     for i in range(1, len(v)):
-> 2843         x = CRT(x,v[i],m,moduli[i])
   2844         m = lcm(m,moduli[i])
   2845     return x%m
/ext/sage/sage-8.1/local/lib/python2.7/site-packages/sage/arith/misc.pyc in crt(a, b, m, n)
   2760         a = Integer(a) # otherwise we get an error at (b-a).quo_rem(g)
   2761     g, alpha, beta = XGCD(m, n)
-> 2762     q, r = (b - a).quo_rem(g)
   2763     if r != 0:
   2764         raise ValueError("No solution to crt problem since gcd(%s,%s) does not divide %s-%s" % (m, n, a, b))
/ext/sage/sage-8.1/src/sage/structure/element.pyx in sage.structure.element.Element.__sub__ (build/cythonized/sage/structure/element.c:11444)()
   1361             return (<Element>left)._sub_(right)
   1362         if BOTH_ARE_ELEMENT(cl):
-> 1363             return coercion_model.bin_op(left, right, sub)
   1364 
   1365         try:
/ext/sage/sage-8.1/src/sage/structure/coerce.pyx in sage.structure.coerce.CoercionModel_cache_maps.bin_op (build/cythonized/sage/structure/coerce.c:10736)()
   1137         # We should really include the underlying error.
   1138         # This causes so much headache.
-> 1139         raise bin_op_exception(op, x, y)
   1140 
   1141     cpdef canonical_coercion(self, x, y):
TypeError: unsupported operand parent(s) for -: 'Ring of integers modulo 3' and 'Ring of integers modulo 2'

m[-1].parent()

Ring of integers modulo 23

lift(m[-1])

3

crt(map(lift, m), l) ## Need to lift mod-p reductions back to integers first.

37033

P.S. If for some reason you really wanted to work modulo some very large prime $p$, you would be better off lifting to integers, then doing this CRT approach. Sage doesn't yet do this, but maybe someday...

Combinatorics

Sage has very well-developed infrastructure for combinatorics, of which we will only scratch the surface here. For PhD-level research in the subject, it is probably the best software available, open-source or otherwise. (And remember that some basic functionality is available in Python via the itertools module, as described last time.)

Reference: the Sage combinatorics tutorial http://doc.sagemath.org/html/en/reference/combinat/sage/combinat/tutorial.html

To begin with, recall that the number of ways to make an unordered choice of $k$ distinguishable objects from among $n$ given objects is the binomial coefficient $$\binom{n}{k} = \frac{n!}{k! (n-k)!}$$ where as usual $n! = 1 \times 2 \times \cdots \times n$ denotes the factorial function.

## Count unordered pairs chosen from a set of 6 objects.
binomial(6, 2)

15

multinomial(2,3,5)

2520

You can use the binomial coefficients to count subsets without generating them. However, if you do actually want to generate the list of $k$-element subsets of a given set, you can do that easily!

l = ["ERC", "Marshall", "Revelle", "Muir", "Warren", "Sixth"]
S = Subsets(l, 2)
show(S)

print(S.cardinality())

15

S[0]

{'Marshall', 'ERC'}

S[3:5]

[{'Warren', 'ERC'}, {'Sixth', 'ERC'}]

S[1]

{'Revelle', 'ERC'}

list(S)

[{'Marshall', 'ERC'},
 {'Revelle', 'ERC'},
 {'Muir', 'ERC'},
 {'Warren', 'ERC'},
 {'Sixth', 'ERC'},
 {'Marshall', 'Revelle'},
 {'Marshall', 'Muir'},
 {'Marshall', 'Warren'},
 {'Sixth', 'Marshall'},
 {'Revelle', 'Muir'},
 {'Revelle', 'Warren'},
 {'Sixth', 'Revelle'},
 {'Muir', 'Warren'},
 {'Sixth', 'Muir'},
 {'Sixth', 'Warren'}]

The Subsets function is an example of a structured set construction. Another important example is the Cartesian product; let's see an example based on playing cards.

suits = ["Clubs", "Diamonds", "Hearts", "Spades"]
ranks = ["Ace", "2", "3", "4", "5", "6", "7", "8", "9", "10", "Jack", "Queen", "King"]
deck = cartesian_product([suits, ranks])

type(deck)

<class 'sage.sets.cartesian_product.CartesianProduct_with_category'>

deck[1]

('Clubs', '2')

Another example is the set of permutations of a list.

print(list(Permutations([1,2,3])))

[[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]

print(list(Permutations([1,2,1,3])))

[[1, 1, 2, 3], [1, 1, 3, 2], [1, 2, 1, 3], [1, 2, 3, 1], [1, 3, 1, 2], [1, 3, 2, 1], [2, 1, 1, 3], [2, 1, 3, 1], [2, 3, 1, 1], [3, 1, 1, 2], [3, 1, 2, 1], [3, 2, 1, 1]]

Exercise for right now: "Deal" a random hand of 13 cards from the deck. You can do this either by "shuffling" the whole deck, or by directly selecting a suitable subset.

t = Permutations(deck).random_element()

t[:13]

[('Clubs', '8'),
 ('Hearts', '6'),
 ('Hearts', 'King'),
 ('Diamonds', '6'),
 ('Spades', '3'),
 ('Clubs', 'Ace'),
 ('Diamonds', '7'),
 ('Spades', 'Jack'),
 ('Hearts', '8'),
 ('Spades', '6'),
 ('Hearts', 'Ace'),
 ('Hearts', '4'),
 ('Hearts', '2')]

S = Subsets(deck, 13)
t = S.random_element()
show(t)

Another example is the set of partitions of a given nonnegative integer. Let's illustrate with an example:

print(list(Partitions(5)))

[[5], [4, 1], [3, 2], [3, 1, 1], [2, 2, 1], [2, 1, 1, 1], [1, 1, 1, 1, 1]]

print(list(Partitions(7)))

[[7], [6, 1], [5, 2], [5, 1, 1], [4, 3], [4, 2, 1], [4, 1, 1, 1], [3, 3, 1], [3, 2, 2], [3, 2, 1, 1], [3, 1, 1, 1, 1], [2, 2, 2, 1], [2, 2, 1, 1, 1], [2, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1]]

There is a method for counting partitions without enumerating them. This is due to Rademacher, and builds on the phenomenal work of Hardy and Ramanujan where they found an asymptotic series.

%time n = Partitions(10^12).cardinality()

print(N(log(n, 10))) ## This number has over 1 million digits! Let's not print it.

CPU times: user 139 µs, sys: 3 µs, total: 142 µs
Wall time: 150 µs
1.11399578738969e6

At a somewhat smaller scale, we can sample random partitions:

Partitions(10^4).random_element()

[424, 379, 287, 245, 240, 228, 228, 189, 181, 173, 172, 171, 152, 151, 150, 147, 144, 137, 129, 126, 126, 111, 111, 105, 101, 101, 98, 94, 92, 89, 85, 80, 78, 77, 75, 75, 75, 75, 70, 70, 70, 68, 65, 65, 64, 64, 63, 63, 58, 58, 57, 57, 57, 52, 52, 51, 51, 51, 51, 51, 51, 51, 49, 49, 48, 48, 48, 47, 46, 46, 46, 46, 46, 44, 44, 41, 40, 40, 39, 38, 36, 36, 33, 32, 32, 32, 30, 30, 30, 30, 30, 30, 30, 30, 26, 26, 26, 26, 26, 26, 26, 24, 24, 22, 22, 22, 22, 21, 20, 19, 19, 19, 19, 19, 19, 19, 17, 17, 17, 17, 17, 17, 17, 17, 17, 17, 16, 15, 15, 15, 15, 14, 14, 14, 12, 12, 12, 11, 11, 11, 11, 11, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 8, 8, 8, 8, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 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, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

One can also enumerate partitions where the parts have fixed size, as in my examples of last time, making change:

WeightedIntegerVectors(37, [1,5,10,25]).list()

[[2, 0, 1, 1],
 [2, 2, 0, 1],
 [7, 1, 0, 1],
 [12, 0, 0, 1],
 [2, 1, 3, 0],
 [7, 0, 3, 0],
 [2, 3, 2, 0],
 [7, 2, 2, 0],
 [12, 1, 2, 0],
 [17, 0, 2, 0],
 [2, 5, 1, 0],
 [7, 4, 1, 0],
 [12, 3, 1, 0],
 [17, 2, 1, 0],
 [22, 1, 1, 0],
 [27, 0, 1, 0],
 [2, 7, 0, 0],
 [7, 6, 0, 0],
 [12, 5, 0, 0],
 [17, 4, 0, 0],
 [22, 3, 0, 0],
 [27, 2, 0, 0],
 [32, 1, 0, 0],
 [37, 0, 0, 0]]

... or solving the (original) Chicken McNuggets problem: https://en.wikipedia.org/wiki/Coin_problem. (I'm punting on my American football example because both Super Bowl teams have proved that failing to convert the point-after-touchdown is not such a rare occurrence!)

WeightedIntegerVectors(92, [6,9,20]).list()

[[2, 0, 4], [0, 8, 1], [3, 6, 1], [6, 4, 1], [9, 2, 1], [12, 0, 1]]

Many combinatorics problems lead naturally to integer sequences. These can be looked up in the On-line Encyclopedia of Integer Sequences (OEIS); but better yet, this database is integrated into Sage! (Whenever you encounter an unknown integer sequence, always consult OEIS. If it's not there yet, maybe it should be added.)

oeis([1,2,3,5,8,13], max_results=4) # requires internet access

0: A000045: Fibonacci numbers: F(n) = F(n-1) + F(n-2) with F(0) = 0 and F(1) = 1.
1: A027926: Triangular array T read by rows: T(n,0) = T(n,2n) = 1 for n >= 0; T(n,1) = 1 for n >= 1; T(n,k) = T(n-1,k-2) + T(n-1,k-1) for k = 2..2n-1, n >= 2.
2: A001129: Iccanobif numbers: reverse digits of two previous terms and add.
3: A005347: First differences of A005579.

Maybe a less familiar example: let's count set partitions of {1, ..., n}.

list(SetPartitions([1..4]))

[{{1, 2, 3, 4}},
 {{1}, {2, 3, 4}},
 {{1, 3, 4}, {2}},
 {{1, 2, 4}, {3}},
 {{1, 2, 3}, {4}},
 {{1, 2}, {3, 4}},
 {{1, 3}, {2, 4}},
 {{1, 4}, {2, 3}},
 {{1}, {2}, {3, 4}},
 {{1}, {2, 4}, {3}},
 {{1}, {2, 3}, {4}},
 {{1, 4}, {2}, {3}},
 {{1, 3}, {2}, {4}},
 {{1, 2}, {3}, {4}},
 {{1}, {2}, {3}, {4}}]

l = [SetPartitions([1..n]).cardinality() for n in [1..10]]
print l

[1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975]

oeis(l, max_results=4)

0: A000110: Bell or exponential numbers: number of ways to partition a set of n labeled elements.
1: A164864: Number of ways of placing n labeled balls into 10 indistinguishable boxes; word structures of length n using a 10-ary alphabet.
2: A276726: Number of set partitions of [n] such that for each block b the smallest integer interval containing b has at most ten elements.
3: A287281: Number of set partitions of [n] such that for each block all absolute differences between consecutive elements are <= nine.

Aha, let's read about that first one some more.

t = oeis('A000110')

#What methods?
t.

t.name()

'Bell or exponential numbers: number of ways to partition a set of n labeled elements.'

t.comments()[:5]

0: The leading diagonal of its difference table is the sequence shifted, see Bernstein and Sloane (1995). - _N. J. A. Sloane_, Jul 04 2015
1: Also the number of equivalence relations that can be defined on a set of n elements. - Federico Arboleda (federico.arboleda(AT)gmail.com), Mar 09 2005
2: a(n) = number of nonisomorphic colorings of a map consisting of a row of n+1 adjacent regions. Adjacent regions cannot have the same color. - _David W. Wilson_, Feb 22 2005
3: If an integer is squarefree and has n distinct prime factors then a(n) is the number of ways of writing it as a product of its divisors. - _Amarnath Murthy_, Apr 23 2001
4: Consider rooted trees of height at most 2. Letting each tree 'grow' into the next generation of n means we produce a new tree for every node which is either the root or at height 1, which gives the Bell numbers. - _Jon Perry_, Jul 23 2003