Math

The imbalance theorem

Posted on by John

The imbalance conjecture is now a theorem. James Alexander Schreib and Yousof Yavari posted a proof last week.

What does the conjecture theorem say? Start with a graph G with no edge between two nodes of the same degree. Then for every edge, calculate the absolute value of the difference of the degree of each end. The imbalance theorem says there exists another graph H whose vertices have degrees corresponding to the differences of degrees in G.

For example, let G be the graph below.

The edges from the top red vertex A to each of the blue vertices around it all have degree difference 5 because A has degree 6 and the vertices a0 to a4 have degree 1. The edge between the two red vertices, A and B, has degree difference 2. The remaining vertices have degree difference 3.

So the multiset of degree differences is

{5, 5, 5, 5, 5, 2, 3, 3, 3}

The imbalance theorem says there exists a graph H whose nodes have these degrees. Here is an example of such an H.

Note that in H, the 5 red nodes have degree 5, the single green node has degree 2, and the three blue nodes have degree 3.

More graph posts

Mean distance to the sun

Posted on by John

Suppose you have a planet in an elliptical orbit around a star. The math is identical for any light object orbiting a heavy object, such as a moon or satellite orbiting a planet, but we’ll call the heavy object a star and the light object a planet.

The center of the star is not quite the center of the orbit. The planet moves along an ellipse with the star at one focus of that ellipse.

Let a be the semi-major axis of planet’s orbit, the maximum distance from the center of the ellipse to a point on the ellipse. Then the distance of a focus to the center of the ellipse is ae where e is the eccentricity of the ellipse. This defines eccentricity. The center of earth’s orbit is between three and four solar radii away from the center of the sun [1].

The planet is farthest from the star when it is along the major axis of the ellipse on the opposite side as the star. The distance is then a + ae, the distance to the center plus the distance from the center to the star. On the opposite side of its orbit, the planet is closest to the star. There the distance is a − ae. In summary the maximum distance to the star is

a(1 + e)

and the minimum distance is

a(1 − e).

If you had to guess the average distance between the planet and its star, a would be a good guess since it’s the average of the maximum and minimum distance. And that’s a good approximation, provided e is small. The mean distance over time is

a(1 + ½e²).

See derivation. The average distance is greater than a because the planet moves faster when nearest the star and slower when further from the star.

The relative error in approximating the mean distance by a is then ½e². When e is small, ½e² is very small. For the earth’s orbit, e = 0.01671, and so the approximation is off by around 0.014%.

The eccentricity of Pluto’s orbit is 0.2488, and so in that case the approximation is off by about 3.1%. The eccentricity of a Molniya orbit, used by some Russian satellites, is 0.74 [2]. For such satellites the error in approximating the mean distance to earth as the semimajor axis is around 27%.

Related posts

[1] For earth’s orbit, e = 0.01671, a = 1.496×1011 m, and the sun’s radius is r = 6.957×108 m. And so ea / r = 3.59.

[2] An object in such a highly elliptical orbit will spend a long time at the far side of its orbit, i.e. over Russia. Sort of a poor man’s geostationary orbit.

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.

 

Probability of correcting errors

Posted on by John

Error correcting codes are most simply described in terms of the errors they can certainly correct. For example, the Hadamard code used for the Mariner 9 probe to Mars encoded each 6-bit pixel to a 32-bit codeword in such a way that the original pixel could be recovered if no more than 7 bits were corrupted in transit.

What is the probability that a pixel could be repaired if corrupted? That depends on your probability model. We will assume that the probability of each bit being flipped is p and that errors are independent.

(Are errors independent, i.e. if a bit flips, is the next bit more or less likely to flip? That would depend on context.)

It’s straight-forward to calculate the probability that 7 or fewer or fewer bits out of 32 flip; this is the cumulative distribution of a binomial random variable. The following Python code will return the probability of k or fewer successes out of n trials, each with probability of success p:

    from scipy.stats import binom
    print(binom.cdf(k, n, p))

For example, if there is a 10% chance that each bit will flip, there’s a 98.8% chance that 7 or fewer bits out of 32 will flip.

However this only gives a lower bound on the probability of correcting an error. If eight bits flip in transit, we cannot tell with certainty which codeword was sent, but there will be a couple possibilities that stand out. We’ll have to guess, but we’ve narrowed down the possibilities. With even more flipped bits, there’s always a chance of recovering the original data. Still, the lower bound captures most of the probability of recovery.

Now suppose you’re given a desired error recovery rate and have to determine what value of p it can sustain. For example, someone might say they want a 98.8% chance of recovering a pixel correctly, and you could come back and say p must be less than or equal to 0.1. This would be a conservative answer because as discussed above, p = 0.1 gives a pixel recovery probability of something more than 98.8, though it’s messy to calculate how much more.

You could solve for p by trial and error, or you could use some more sophisticated math to compute p directly. Given a probability F, you can solve for p such that the probability of up to k successes out of n trials using the inverse of the regularized incomplete beta function.

from scipy.special import betaincinv
p = 1 - betaincinv(n - k, k + 1, F)

Calculating F given n, k, and p could be a homework exercise in an introductory probability course. Solving for p given F, n, and k either requires some numerical programming or special functions and so would be a more challenging problem.

Related posts

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.

Hadamard Codes and Sphere Packing

Posted on by John

Yesterday Levent Alpöge announced that he and his colleagues had discovered a new Hadamard matrix using Claude AI. That motivated a post I wrote this morning on how to construct Hadamard matrices. I mentioned in that post that these matrices arise in applications.

This evening I gave an example, describing how NASA used a Hadamard matrix of order 32 to transmit photos from the Mariner 9 spacecraft in 1971. This post will give another application: sphere packing.

Conway and Sloane [1] give a correspondence between binary codes and sphere packings that they call Construction A. Given an (n, M, d) binary code C, center a sphere on a point x if and only if x is a congruent (mod 2) to codeword in C.

Here (n, M, d) means an error correcting code that encodes M bits of data as strings of n bits, with a minimum Hamming distance between code words of d, i.e. all codewords differ in at least d bits.

The previous post described how to create a (32, 6, 16) code by stacking a Hadamard matrix H of order 32 on top of −H and turning −1’s into 0’s. The analogous construction for a (8, 4, 4) Hadamard code gives E8, the densest packing in ℝ8.

We start with the Hadamard matrix

H_8 = \begin{pmatrix} 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 & 1 & -1 & -1 & 1 & -1 & 1 \\ 1 & 1 & -1 & -1 & -1 & -1 & 1 & 1 \\ 1 & -1 & -1 & 1 & -1 & 1 & 1 & -1 \end{pmatrix}

and obtain the matrix

M = \begin{pmatrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 0 & 1 & 0 & 1 & 0 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 1 & 0 & 0 & 1 \\ 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 0 & 0 & 0 & 1 & 1 \\ 1 & 0 & 0 & 1 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 & 0 & 0 & 1 & 1 \\ 0 & 1 & 1 & 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & 1 & 1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \\ 0 & 1 & 1 & 0 & 1 & 0 & 0 & 1 \end{pmatrix}

whose centers form the sphere packing.

This doesn’t look like the E8 sphere packing as it is usually presented, but it’s isomorphic.

[1] J. H. Conway and N. J. A. Sloane. Sphere Packings, Lattices and Groups. Springer. 1999.

 

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

Inverse differential equations

Posted on by John

In science and engineering classes, you might describe a system using Newton’s laws and end up with a differential equation. You then solve the differential equation, analytically or numerically, to see how the solutions behave.

You might also do the opposite, especially in a mathematics class: look at what differential equation a set of functions satisfy in order to understand those functions.

Bessel functions came out of solving differential equations from astronomy. But then they turned out to be useful, not just in other areas of science, but in pure mathematics as well. So there are Bessel function users who came to the functions first and haven’t seen the differential equation they came from.

You can learn a lot about Bessel functions, and other special functions, by looking at their defining differential equation even if you’re not directly interested in the differential equation or the physical problem that motivated it.

Bessel functions satisfy

x² y″ + x y′ + (x² − n²) y = 0.

You can tell a lot about Bessel functions just by inspecting this equation without solving it. If we divide by x² and write the equation in the form

y″ + (p(x)/x) y′ + (q(x)/x²) y = 0

then p(x) = 1 and q(x) = x² − n². The indicial equation

r(r − 1) + p(0) r + q(0) = 0

reduces to

r² = n²

and so r = ±n. That alone tells us there are two solutions, one analytic at zero and one singular at zero. These are Jn and Yn respectively. It also tells us the behavior of these functions as x goes to zero and as x goes to infinity. To find out more, look up “method of Frobinius.”

How not to calculate cosine

Posted on by John

Calculus professors with no experience in numerical computing will tell students that computers calculate trig functions with power series. They don’t. I worked on the implementation of trig functions in hardware, and I can assure you we didn’t just use power series.

Power series are an excellent way to calculate functions near the center of the series, such as computing sine for small angles. But the further you get from the center, the less useful power series are.

Let’s suppose you want to calculate cos(200) using the power series for cosine. The nth term of that series is

(−1)n x2n / (2n)!

This is an alternating series, and so the error in truncating the series after n terms is bounded by the size of the n+1 term, if you’ve gone far enough out in the series that the terms are monotonically decreasing in absolute value.

To calculate cos(200) to machine precision, i.e. with an error of less than 2−52, we’d need to sum the series up to n where

| 2002n+2 / (2n + 2)! | < 2−52

Actually, that will ensure that the absolute error is small enough, but not that the relative error is small enough; if the value of cos(200) is small, we’d need more terms. Let’s ignore that and assume we’re only concerned with absolute error.

Turns out we’d need 287 terms. That’s a lot of terms. But you might say “That’s fine. I’m not in a hurry, and it’s just more work for the computer, not for me.” OK, so let’s try.

from math import *

s = 0
for n in range(288):
    s += (-1)**n * 200**(2*n) / factorial(2*n)
print(s)

This prints -3.6840358571084123e+67. You may suspect the answer is incorrect since values of cosine are on the order of 1, not on the order of 1067. Something went spectacularly bad. On closer inspection, it’s remarkable the code didn’t crash.

If you changed 200 to 200.0 above, the code would crash. Calculating 200.0**(2*n) overflows when n = 67. But when we calculate 200**(2*n), the result is an integer. And we’re dividing by factorial(2*n), which is also an integer. Both of these integers become too large to fit in a float, but their ratio has a maximum value of around 1080, smaller than the maximum float, which is on the order of 10308.

When we don’t overflow, we have a different problem: catastrophic cancellation. You can’t calculate a number between −1 and 1 as an alternating sum of numbers as large as 1080. You’d need more than 80 + 16 = 96 decimal places of precision to compute the sum accurately, and floating point only gives you between 15 and 16 decimal places of precision.

So how would you calculate cos(200)? The first step would be to use some sort of range reduction on 200. You could reduce 200 mod 2π to get a smaller number to work with.

>>> from math import cos, pi
>>> x = 200 % (2*pi)
>>> x
5.221255477432827

Using a power series to compute the cosine of 5.221255477432827 is feasible, but not optimal. There’s also another problem: the naive range reduction above loses some precision.

>>> cos(x)
0.48718767500701254
>>> cos(x) - cos(200)
6.661338147750939e-15

The error is small, but it’s still an order of magnitude larger than machine precision. You can’t simply reduce n mod 2π with ordinary float division because the integer part of n / 2π pushes some digits of precision off the right end. I intend to write about how range reduction works in future posts.

Update: See this post for a simple range reduction algorithm that is fine for values of x such as 200, but not adequate for much larger values.

Calculating log(1000!)

Posted on by John

The previous post pointed out that the following code such as the following unexpectedly works.

>>> from math import log, factorial
>>> log(factorial(1000))
5912.128178488163

If you don’t find this unexpected, note that if you replace math.log with numpy.log the code will fail [1]. Functions like natural logarithm operate on real numbers. Real numbers are represented as floating point numbers in programming languages, and 1000! factorial is too large to represent as a standard floating point number. (More on that here.)

In this post I’d like to look at how you might calculate log(1000!) with less capable software, and even without software.

One approach would be to sum the logarithms of the numbers 1 through 1000. This will give essentially the same result as above, with a little difference in the last couple decimal places due to rounding error.

If you have a way to calculate 1000! but not a way to cast it to a floating point number, you could do this manually.

>>> s = str(factorial(1000))
>>> s[:16]
'4023872600770937'
>>> len(s)
2568

This tells us 1000! = 4.023872600770937 × 102567. Therefore

log(1000!) = log(4.023872600770937) + 2567 log(10)

which only requires working with numbers of modest size.

Calculating by hand

Now suppose it’s 1964. You don’t have a computer, or even a calculator, but you do have a copy of the recently published Handbook of Mathematical Functions by Abramowitz and Stegun (A&S). You turn to Table 6.6 “Factorials for large arguments.” This has values of factorial for 100, 200, 300, …, 1000, so you can simply look up your answer to 20 decimal places.

That was too easy; I didn’t expect that to be there when I started writing this post. If you wanted to compute log(950!), for example, you’d have to work harder. You could find A&S equation 6.1.41 (Stirling’s series) which says

\begin{align*} \ln \Gamma(z) &\sim (z - \tfrac{1}{2}) \ln z - z + \tfrac{1}{2} \ln 2\pi + \frac{1}{12z} - \frac{1}{360z^3} \\ &+ \frac{1}{1260z^5} - \frac{1}{680z^7} + \cdots \end{align*}

So how would you use this formula to calculate log(1000!)? Since n! = Γ(n + 1), you set z = 1001.

You’d need to decide how many terms you need to use. Assuming the error is on the order of the first term you leave out, you’d reason that you could probably stop with the 1/12z term because the next term is between 10−11 and 10−12.

You find Table 4.2 has natural logarithms, but not for 1001. You can look up log(1.001), however, and at the bottom of the same page is log(10) to 16 decimal places, and you can find log(10) to 24 decimal places in Table 1.1. So you calculate

log(1001) = log(1.001 × 10³) = log(1.001) + 3 log(10).

You can find log(2) and log(π) in Table 1.1, and average them to find ½ log(2π).

Here’s Python code to simulate the hand calculations.

log2     = 0.6931_47180_55994_53094_172321 # Table 1.1
log10    = 2.3025_85092_99404_56840_179915 # Table 1.1
logpi    = 1.1447_29885_84940_01741_43427  # Table 1.1
log1_001 = 0.00099_95003_330835            # Table 4.2

z = 1001
logz = log1_001 + 3*log10
s = (z - 0.5)*logz - z + (log2 + logpi)/2 + 1/(12*z)

print(s)

This result differs from the one at the top of the post only in the last decimal place.

Related posts

Doing calculations with tables is not as simple as “just look it up.” It takes a bit of skill.

[1] The code will also fail if you replace math.log with math.cos. Both logarithm and cosine return moderate sized real numbers when given enormous inputs like 1000!, so representing the output as a float is not the problem. But logarithms of huge numbers can be computed with ordinary precision functions, as above. But computing the cosine of a huge number requires extended precision.

Update: The next post expands on why computing the cosine of a large number is more difficult than computing the log.