Approximating rapidly divergent integrals

Posted on by John

A while back I ran across a paper [1] giving a trick for evaluating integrals of the form

I(M) = \int_a^M \exp(f(x)) \, dx

where M is large and f is an increasing function. For large M, the integral is asymptotically

A(M) = \frac{\exp(f(M))}{f'(M)}.

That is, the ratio of A(M) to I(M) goes to 1 as M goes to infinity.

This looks like a strange variation on Laplace’s approximation. And although Laplace’s method is often useful in practice, no applications of the approximation above come to mind. Any ideas? I have a vague feeling I could have used something like this before.

There is one more requirement on the function f. In addition to being an increasing function, it must also satisfy

\lim_{x\to\infty} \frac{f''(x)}{(f'(x))^2} = 0.

In [1] the author gives several examples, including using f(x) = x². If we wanted to approximate

\int_0^100 \exp(x^2)\, dx

the method above gives

exp(10000)/200 = 4.4034 × 104340

whereas the correct value to five significant figures is 4.4036 × 104340.

Even getting an estimate of the order of magnitude for such a large integral could be useful, and the approximation does better than that.

[1] Ira Rosenholtz. Estimating Large Integrals: The Bigger They Are, The Harder They Fall. The College Mathematics Journal, Vol. 32, No. 5 (Nov., 2001), pp. 322-329

Best approximation of a catenary by a parabola

Posted on by John

A parabola and a catenary can look very similar but are not the same. The graph of

y = x²

is a parabola and the graph of

y = cosh(x) = (ex + e–x)/2

is a catenary. You’ve probably seen parabolas in a math class; you’ve seen a catenary if you’ve seen the St. Louis arch.

Depending on the range and scale, parabolas and catenaries can be too similar to distinguish visually, though over a wide range enough range the exponential growth of the catenary becomes apparent.

For example, for x between -1 and 1, it’s possible to scale a parabola to match a catenary so well that the graphs practically overlap. The blue curve is a catenary and the orange curve is a parabola.

The graph above looks orange because the latter essentially overwrites the former. The relative error in approximating the catenary by the parabola is about 0.6%.

But when x ranges over -10 to 10, the best parabola fit is not good at all. The catenary is much flatter in the middle and much steeper in the sides. On this wider scale the hyperbolic cosine function is essentially e|x|.

Here’s an intermediate case, -3 < x < 3, where the parabola fits the catenary pretty well, though one can easily see that the curves are not the same.

Now for some details. How are we defining “best” when we say best fit, and how do we calculate the parameters for this fit?

I’m using a least-squares fit, minimizing the L² norm of the error, over the interval [M, M]. That is, I’m approximating

cosh(x)

with

c + kx²

and finding c and k that minimize the integral

\int_{-M}^M (\cosh(x) - c - kx^2)^2\, dx

The optimal values of c and k vary with M. As M increases, c decreases and k increases.

It works out that the optimal value of c is

-\frac{3 \left(M^2 \sinh (M)+5 \sinh (M)-5 M \cosh (M)\right)}{2 M^3}

and the optimal value of k is

\frac{15 \left(M^2 \sinh (M)+3 \sinh (M)-3 M \cosh (M)\right)}{2 M^5}

Here’s a log-scale plot of the L² norm of the error, the square root of the integral above, for the optimal parameters as a function of M.

More on catenaries

Five places the Sierpiński triangle shows up

Posted on by John

The Sierpiński triangle is a fractal that comes up in unexpected places. I’m not that interested in fractals, and yet I’ve mentioned the Sierpiński triangle many times on this blog just because I run into it while looking at something else.

The first time I wrote about the Sierpiński triangle was when it came up in the context of a simple random process called the chaos game.

Unbiased chaos game results

Next I ran into Sierpiński in the context of cellular automata, specifically Rule 90. A particular initial condition for this rule leads to the image below. With other initial conditions you don’t get such a clean Sierpiński triangle, but you do get similar variations on the theme.

Rule 90 with one initial bit set

Next I ran into Sierpiński in the context of low-level programming. The following lines of C code prints an asterisk when the bit-wise and of two numbers is non-zero.

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++)
            printf("%c", (i&j) ? ' ' : '*');
        printf("\n");
    }

A screenshot of the output shows our familiar triangle.

screen shot that looks like Sierpinski triangle

Then later I wrote a post looking at constructible n-gons, n-sided figures that can be constructed using only a straight edge and a compass. These only exist for special values of n. If you write these special values in binary, and replace the 1’s with a black square and the 0’s with a blank, you get yet another Sierpiński triangle.

Finally, if you look at the odd numbers in Pascal’s triangle, they also form a Sierpiński triangle.

Evolute of an egg

Posted on by John

The set of lines perpendicular to a curve are tangent to a second curve called the evolute. The lines perpendicular to the ellipse below are tangent to the curve inside called an astroid.

If we replace the ellipse with an egg, we get a similar shape, but less symmetric.

The equation for the egg is described here with parameters a = 3, b = 2, and k = 0.1. The ellipse above has the same a and b but k = 0.

I made the lines slightly transparent, setting alpha = 0.4, so the graph would be darker where many lines cross.

Related post: Envelopes of epicycloids

Sample size calculation

Posted on by John

If you’re going to run a test on rabbits, you have to decide how many rabbits you’ll use. This is your sample size. A lot of what statisticians do in practice is calculate sample sizes.

A researcher comes to talk to a statistician. The statistician asks what effect size the researcher wants to detect. Do you think the new thing will be 10% better than the old thing? If so, you’ll need to design an experiment with enough subjects to stand a good chance of detecting a 10% improvement. Roughly speaking, sample size is inversely proportional to the square of effect size. So if you want to detect a 5% improvement, you’ll need 4 times as many subjects as if you want to detect a 10% improvement.

You’re never guaranteed to detect an improvement. The race is not always to the swift, nor the battle to the strong. So it’s not enough to think about what kind of effect size you want to detect, you also have to think about how likely you want to be to detect it.

Here’s what often happens in practice. The researcher makes an arbitrary guess at what effect size she expects to see. Then initial optimism may waver and she decides it would be better to design the experiment to detect a more modest effect size. When asked how high she’d like her chances to be of detecting the effect, she thinks 100% but says 95% since it’s necessary to tolerate some chance of failure.

The statistician comes back and says the researcher will need a gargantuan sample size. The researcher says this is far outside her budget. The statistician asks what the budget is, and what the cost per subject is, and then the real work begins.

The sample size the negotiation will converge on is the budget divided by the cost per sample. The statistician will fiddle with the effect size and probability of detecting it until the inevitable sample size is reached. This sample size, calculated to 10 decimal places and rounded up to the next integer, is solemnly reported with a post hoc justification containing no mention of budgets.

Sample size is always implicitly an economic decision. If you’re willing to make it explicitly an economic decision, you can compute the expected value of an experiment by placing a value on the possible outcomes. You make some assumptions—you always have to make assumptions—and calculate the probability under various scenarios of reaching each conclusion for various sample sizes, and select the sample size that leads to the best expected value.

More on experimental design

[1] There are three ways an A/B test can turn out: A wins, B wins, or there isn’t a clear winner. There’s a tendency to not think enough about the third possibility. Interim analysis often shuts down an experiment not because there’s a clear winner, but because it’s becoming clear there is unlikely to be a winner.

Binomial coefficients mod primes

Posted on by John

Imagine seeing the following calculation:

{95 \choose 57} = {19\cdot 5 \choose 19\cdot 3} = {5 \choose 3} = \frac{5\cdot 4}{2\cdot 1} = 10

The correct result is

{95 \choose 57} = 487343696971437395556698010

and so the first calculation is off by 25 orders of magnitude.

But there’s a variation on the calculation above that is correct! A theorem by Édouard Lucas from 1872 that says for p prime and for any nonnegative integers m and n,

{pm \choose pn} = {m \choose n} \bmod p

So while the initial calculation was grossly wrong as stated, it is perfectly correct mod 19. If you divide 487343696971437395556698010 by 19 you’ll get a remainder of 10.

A stronger versions of Lucas’ theorem [1] says that if p is at least 5, then you can replace mod p with mod p³. This is a stronger conclusion because it says not only is the difference between the left and right side of the congruence divisible by p, it’s also divisible by p² and p³.

In our example, not only is the remainder 10 when 487343696971437395556698010 is divided by 19, the remainder is also 10 when dividing by 19² = 361 and 19³ = 6859.

More on binomial coefficients

[1] V. Brun, J. O. Stubban, J. E. Fjeldstad, L. Tambs, K. E. Aubert, W. Ljunggren, E. Jacobsthal. On the divisibility of the difference between two binomial coefficients, Den 11te Skandinaviske Matematikerkongress, Trondheim, 1949, 42–54.

Surface of revolution with minimum area

Posted on by John

Suppose you’re given two points (x1, y1) and (x2, y2) with y1 and y2 positive. Find the smooth positive curve f(x) that passes through the two points such that the area of the surface formed by rotating the graph of f around the x-axis is minimized.

You can state this as a problem in calculus of variations, which leads to a differential equation, which leads to the solution

f(x) = c cosh((x + d)/c).

In other words, the surface area is minimized when the graph of f is a piece of a catenary [1].

This is interesting because the answer is not something you’re likely to guess, unlike say the isoperimetric problem, where the it’s easy to guess (but hard to prove) that the solution is a circle.

There’s also some interesting fine print to the solution. It’s not quite right to say that the solution is a catenary. To be more precise we should say that if there is a unique catenary that passes through both specified points, then it is the smooth curve with minimal area when rotated about the x-axis. But there are a couple things that could go wrong.

It’s possible that two catenaries pass through the given points, and in that case one of the catenaries leads to minimal surface area. But it’s also possible that there is no catenary passing through the given points.

My first thought would be that you could always find values of c and d so that the function f passes through the points (x1, y1) and (x2, y2), but that’s not true. Often you can, but if the difference in the y‘s is very high relative to the difference in the x‘s it might not be possible.

Suppose the graph of f passes through (0, 1) and (1, y2).

Since the graph passes through the first point, we have

c cosh(d/c) = 1.

Since cosh(x) ≥ 1, we must also have c ≤ 1. And since our curve is positive, we must have c > 0. We can maximize

c cosh((1 + d)/c)

for 0 < c ≤ 1 subject to the constraint

c cosh(d/c) = 1

to find the maximum possible value of y2. If we ask Mathematica

    NMaximize[
        {   c Cosh[(1 + d)/c], 
            {0 < c <= 1}, 
            {c Cosh[d/c] == 1}
        }, 
        {c, d}
    ]

we get

    {6.45659*10^8, {c -> 0.0352609, d -> -0.142316}}

meaning the largest possible value of y2 is 6.45659 × 108, and it occurs when c = 0.0352609, d = -0.142316.

Update: See the comment by Bill Smathers below arguing that the maximum should be unbounded. If the argument is correct, this would imply the code above ran into a numerical limitation.

Related posts

[1] See Calculus of Variations by I. M. Gelfand and S. V. Fomin.

Chemical element frequency in writing

Posted on by John

How do the frequencies of chemical element names in English text compare to the abundance of elements in Earth’s crust? Do we write most frequently about the elements that appear most frequently?

It turns out the answer is “not really.” The rarest elements rarely appear in writing. We don’t have much to say about dysprosium, thulium, or lutetium, for example. But overall there’s only a small correlation between word frequency and chemical frequency. (The rank correlation is substantially higher than ordinary linear correlation.)

We write often about things like oxygen and iron because they’re such a part of the human experience. On the other hand, we care about some things like silver and gold precisely because they are rare.

Here are the most common elements according to text usage.

|------------+--------+-----------+---------+------------|
| element    | word % | word rank | earth % | earth rank |
|------------+--------+-----------+---------+------------|
| lead       |  15.50 |         1 |   0.001 |         36 |
| gold       |  11.64 |         2 |   0.000 |         75 |
| iron       |  11.14 |         3 |   5.612 |          4 |
| silver     |   7.38 |         4 |   0.000 |         68 |
| carbon     |   5.15 |         5 |   0.012 |         17 |
| oxygen     |   5.13 |         6 |  45.956 |          1 |
| copper     |   4.61 |         7 |   0.006 |         26 |
| hydrogen   |   3.51 |         8 |   0.139 |         10 |
| sodium     |   3.38 |         9 |   2.352 |          6 |
| calcium    |   2.84 |        10 |   4.137 |          5 |
| nitrogen   |   2.79 |        11 |   0.002 |         34 |
| mercury    |   2.22 |        12 |   0.000 |         67 |
| tin        |   2.13 |        13 |   0.000 |         51 |
| potassium  |   1.94 |        14 |   2.083 |          8 |
| zinc       |   1.70 |        15 |   0.007 |         24 |
| silicon    |   1.12 |        16 |  28.112 |          2 |
| nickel     |   1.08 |        17 |   0.008 |         23 |
| phosphorus |   1.05 |        18 |   0.104 |         11 |
| magnesium  |   0.98 |        19 |   2.322 |          7 |
| sulfur     |   0.84 |        20 |   0.035 |         16 |
|------------+--------+-----------+---------+------------|

This is based on the Google book corpus summarized here. There’s some ambiguity; I imagine most used of “lead” are the verb and not the element name. Some portion of the uses of “iron” refer to a device for smoothing wrinkles out of clothes.

Word percentage is relative to the set of chemical element names. Earth percentage is relative to the Earth’s crust.

The percentages above have been truncated for presentation’ obviously the abundance of gold, silver, mercury, and tin is not zero, though it is when rounded to three decimal places. The full data for the first 111 elements is available here.

Convex function of diagonals and eigenvalues

Posted on by John

Sam Walters posted an elegant theorem on his Twitter account this morning. The theorem follows the pattern of an equality for linear functions generalizing to an inequality for convex functions. We’ll give a little background, state the theorem, and show an example application.

Let A be a real symmetric n×n matrix, or more generally a complex n×n Hermitian matrix, with entries aij. Note that the diagonal elements aii are real numbers even if some of the other entries are complex. (A Hermitian matrix equals its conjugate transpose, which means the elements on the diagonal equal their own conjugate.)

A general theorem says that A has n eigenvalues. Denote these eigenvalues λ1, λ2, …, λn.

It is well known that the sum of the diagonal elements of A equals the sum of its eigenvalues.

\sum_{i=1}^n a_{ii} = \sum_{i=1}^n \lambda_i

We could trivially generalize this to say that for any linear function φ: R → R,

\sum_{i=1}^n \varphi(a_{ii}) = \sum_{i=1}^n \varphi({\lambda_i})

because we could pull any shifting and scaling constants out of the sum.

The theorem Sam Walters posted says that the equality above extends to an inequality if φ is convex.

\sum_{i=1}^n \varphi(a_{ii}) \leq \sum_{i=1}^n \varphi({\lambda_i})

Here’s an application of this theorem. Assume the eigenvalues of A are all positive and let φ(x) = – log(x). Then φ is convex, and

-\sum_{i=1}^n \log(a_{ii}) \leq -\sum_{i=1}^n \log({\lambda_i})

and so

\prod_{i=1}^n a_{ii} \geq \prod_{i=1}^n \lambda_i = \det(A)

i.e. the product of the diagonals of A is an upper bound on the determinant of A.

This post illustrates two general principles:

  1. Linear equalities often generalize to convex inequalities.
  2. When you hear a new theorem about convex functions, see what it says about exp or -log.

More linear algebra posts

Bit flipping to primes

Posted on by John

Someone asked an interesting question on MathOverflow: given an odd number, can you always flip a bit in its binary representation to make it prime?

It turns out the answer is no, but apparently it is very often the case an odd number is just a bit flip away from being prime. I find that surprising.

Someone pointed out that 2131099 is not a bit flip away from a prime, and that this may be the smallest example [1]. The counterexample 2131099 is itself prime, so you could ask whether an odd number is either a prime or a bit flip away from a prime. Is this always the case? If not, is it often the case?

The MathOverflow question was stated in terms of Hamming distance, counting the number of bits in which two bit sequences differ. It asked whether odd numbers are always Hamming distance 1 away from a prime. My restatement of the question asks whether the Hamming distance is always at most 1, or how often it is no more than 1.

You could ask more generally about the Hamming distance to the nearest prime. Is it bounded, if not by 1, then by another finite number? If so, what is the smallest such bound? What is the probability that its value is 1? Etc.

This ties into a couple of other things I’ve blogged about. A few weeks ago I wrote about new work on the problem of finding the proportion of odd numbers that can be written as the sum of a power of 2 and a prime. That’s a little different problem since bit flipping is taking the XOR (exclusive or) and not always the same as addition. It also leaves out the possibility of flipping a bit beyond the most significant bit of the number, i.e. adding to a number n a power of 2 greater than n.

Another related post is on the Rowhammer attack on public key cryptography. By flipping a bit in the product of two primes, you can produce a number which is much easier to factor.

These two posts suggest a variation on the original problem where we disallow flipping bits higher than the most significant bit of n. So giving a k-bit number n, how often can we flip one of its k bits and produce a prime?

[1] Note that the bit flipped may be higher than the most significant bit of the number, unless ruled out as in the paragraph above. Several people have asked “What about 85?” It is true that flipping any of the seven lowest bits of 85 will not yield a prime. But flipping a zero bit in a more significant position will give a prime. For example, 1024 + 85 is prime. Bur for 2131099 it is not possible to add any larger power of 2 to the number and produce a prime.