Accurately computing a 2×2 determinant

Posted on 31 May 2020 by John

The most obvious way to compute the determinant of a 2×2 matrix can be numerically inaccurate. The biggest problem with computing ad − bc is that if ad and bc are approximately equal, the subtraction could lose a lot of precision. William Kahan developed an algorithm for addressing this problem.

Fused multiply-add (FMA)

Kahan’s algorithm depends on a fused multiply-add function. This function computes xy + z using one rounding operation at the end, where the direct approach would use two.

In more detail, the fused multiply-add behaves as if it takes its the floating point arguments x, y, and z and lifts them to the Platonic realm of real numbers, calculates xy + z exactly, and then brings the result back to the world of floating point numbers. The true value of xy + z may not have an exact representation in the floating point world, and so in general it will be necessary to round the result.

The direct way of computing xy + z would behave as if it first lifts x and y to the Platonic realm, computes xy exactly, then pushes the result back down to the floating point world, rounding the result. Then it takes this rounded version of xy back up to the Platonic realm, possibly to a different value than before, and pushes z up, adds the two numbers, then pushes the sum back down to the world of floating point.

You can get more accuracy if you avoid round trips back and forth to the Platonic realm. If possible, you’d like to push everything up to the Platonic realm once, do as much work as you can up there, then come back down once. That’s what fused multiply-add does.

Kahan’s determinant algorithm

Here is Kahan’s algorithm for computing ad − bc for floating point numbers.

w = RN(bc)
e = RN(w – bc)
f = RN(ad – w)
x = RN(f + e)

The function RN computes its argument exactly, then rounds the result to the nearest floating point value, rounding to even in case of a tie. (This is the default rounding mode for compilers like gcc, so the C implementation of the algorithm will directly follow the more abstract specification above.)

If there were no rounding we would have

w = bc
e = w − bc = 0
f = ad − w = ad − bc
x = f + e = ad − bc + 0 = ad − bc

But of course there is rounding, and so Kahan’s steps that seem redundant are actually necessary and clever. In floating point arithmetic, e is not zero, but it does exactly equal w – bc. This is important to why the algorithm works.

Here is a C implementation of Kahan’s algorithm and a demonstration of how it differs from the direct alternative.

    #include <math.h>
    #include <stdio.h>
    
    double naive_det(double a, double b, double c, double d) {
        return a*d - b*c;
    }
    
    double kahan_det(double a, double b, double c, double d) {
        double w = b*c;
        double e = fma(-b, c, w);
        double f = fma(a, d, -w);
        return (f+e);
    }
    
    int main() {
        double a = M_PI;
        double b = M_E;
        double c = 355.0 / 113.0;
        double d = 23225.0 / 8544.0;
    
        printf("Naive: %15.15g\n", naive_det(a, b, c, d));
        printf("Kahan: %15.15g\n", kahan_det(a, b, c, d));
    }

The values of c and d were chosen to be close to M_PI and M_E so that the naive computation of the determinant would less accurate. See these post on rational approximations to π and rational approximations to e.

The code above prints

    Naive: -7.03944087021569e-07
    Kahan: -7.03944088015194e-07

How do we know that Kahan’s result is accurate? Careful analysis of Kahan’s algorithm in [1] gives bounds on its accuracy. In particular, the absolute error is bounded by 1.5 ulps, units of least precision.

Ratio of area to perimeter

Posted on 30 May 2020 by John

Given a curve of a fixed length, how do you maximize the area inside? This is known as the isoperimetric problem.

The answer is to use a circle. The solution was known long before it was possible to prove; proving that the circle is optimal is surprisingly difficult. I won’t give a proof here, but I’ll give an illustration.

Consider a regular polygons inscribed in a circle. What happens to the ratio of area to perimeter as the number of sides increases? You might suspect that the ratio increases with the number of sides, because the polygon is becoming more like a circle. This turns out to be correct, and it’s not that hard to be precise about what the ratio is as a function of the number of sides.

For a regular polygon inscribed in a circle of radius r,

$\text{Area } = \frac{n}{2} r^2 \sin\left(\frac{2\pi}{n} \right )$

and

$\text{Perimeter} = 2n r \sin\left(\frac{\pi}{n} \right )$

For a regular n-gon inscribed in a unit circle, we have

$\begin{align*} \frac{\text{Area}}{\text{Perimeter}} &= \frac{(n/2) \sin(2\pi/n)}{2n \sin(\pi/n)} \\ &= \frac{n \sin(\pi/n)\cos(\pi/n)}{2n \sin(\pi/n)} \\ &= \cos(\pi/n)/2 \end{align*}$

We used the double-angle identity for sine in the second line above.

As n increases, the ratio increases toward 1/2, the ratio of the area of a unit circle to its circumference. In a little more detail, the difference between the ratios for a circle and for a regular n-gon goes to zero like O(1/n²), based on a Taylor expansion for cosine.

Here’s a plot of the ratios as a function of the number of sides.

Curse of dimensionality and integration

Posted on 30 May 2020 by John

The curse of dimensionality refers to problems whose difficulty increases exponentially with dimension. For example, suppose you want to estimate the integral of a function of one variable by evaluating it at 10 points. If you take the analogous approach to integrating a function of two variables, you need a grid of 100 points. For a function of three variables, you need a cube of 1000 points, and so forth.

You cannot estimate high-dimensional integrals this way. For a function of 100 variables, using a lattice with just two points in each direction would require 2¹⁰⁰ points.

There are much more efficient ways to approximate integrals than simply adding up values at grid points, assuming your integrand is smooth. But when applying any of these methods to multi-dimensional integrals, the computational effort increases exponentially with dimension.

The methods that first come to mind don’t scale well with dimension, but that doesn’t necessarily mean there aren’t any methods that do scale well.

Are there numerical integration methods whose computational requirement scale slower than exponentially with dimension? In general, no. You can beat the curse of dimension for special integrals. And if you’re content with probabilistic bounds rather than deterministic bounds, you can get around the curse by using Monte Carlo integration, sort of [1].

If you want worst-case error bounds, you cannot escape the curse of dimensionality [2]. You can require that the functions you’re integrating have so many derivatives and that these derivatives are bounded by some constant, but the amount of computation necessary to guarantee that the error is below a specified threshold increases exponentially with dimension.

Nuclear Tornadoes

Posted on 29 May 2020 by John

Tornado with lightning

I’ve recently started reading One Giant Leap, a book about the Apollo missions. In a section on nuclear testing, the book describes what was believed to be a connection between weapon tests and tornadoes.

In 1953 and 1954 there were all-time record outbreaks of tornadoes across the US, accompanied by fierce thunderstorms. The belief that the tornadoes were caused by the atomic testing in Nevada was so common among the public that Gallup conducted a poll on the topic, which showed that 29% of Americans thought that A-bombs were causing tornadoes many states away; 20% couldn’t say whether they were or not. Members of Congress required the military to provide data on any connection.

It wasn’t unreasonable at the time to speculate that there could be a connection between nuclear testing and tornadoes. No one knew much about the effects of nuclear explosions or about the causes of tornadoes in 1954. But here’s the part that was preventable:

In fact it turned out there weren’t really more tornadoes in those years; it was just that tornado reporting and tracking had gotten much more effective and thorough.

Once the public got the impression that tornadoes were becoming more frequent, confirmation bias would firm up that impression with each tornado citing.

I imagine someone knew this in 1954. Someone must have looked at the data and said “Hey, wait a minute.” This person was probably ridiculed by the few who even listened.

Update: I found a newspaper article cited in One Giant Leap: “A-Bomb Link to Tornadoes is Discounted,” Washington Post. November 18, 1954. The article quoted D. Lee Harris of the United States Weather Bureau saying “Improvement in tornado reporting is largely responsible for the great increase in total number of the storms reported in recent years.”

The Gallup poll was prior to the article cited. It would be interesting if there had been a follow-up poll to see how much impact Harris’ statement had on public opinion.

Fundamental Theorem of Arithmetic

Posted on 28 May 2020 by John

It’s hard to understand anything from just one example. One of the reason for studying other planets is that it helps us understand Earth. It can even be helpful to have more examples when the examples are purely speculative, such as xenobiology, or even known to be false, i.e. counterfactuals, though here be dragons.

The fundamental theorem of arithmetic seems trivial until you see examples of similar contexts where it isn’t true. The theorem says that integers have a unique factorization into primes, up to the order of the factors. For example, 12 = 2² × 3. You could re-order the right hand side as 3 × 2², but you can’t change the list of prime factors and their exponents.

I was unimpressed when I first heard of fundamental theorem of arithmetic. It seemed obvious, and hardly worth distinguishing as the fundamental theorem of arithmetic. (If not this, what would the fundamental theorem of arithmetic be? Any candidates? I didn’t have one, and still don’t.)

In order to appreciate the unique factorization property of the integers, mathematicians created algebraic structures analogous to the integers, i.e. rings, in which unique factorization may or may not hold. In the language of ring theory, the integers are a “unique factorization domain” but not all rings are.

An easy way to create new rings is to extend the integers by adding new elements, analogous to the way the complex numbers are created by adding i to the reals.

If you extend the integers by appending √−5, i.e. using the set of all numbers of the form a + b√−5 where a and b are integers, you get a ring that is not a unique factorization domain. In this ring, 6 can be factored as 2 × 3 but also as (1 + √−5)(1 − √−5).

However, you can extend the integers by other elements and you may get a unique factorization domain. For example, the Eisenstein [1] integers append (−1 + i√ 3)/2 to the ordinary integers, and the result is a unique factorization domain.

The Fundamental Theorem of Algebra

Posted on 27 May 2020 by John

This post will take a familiar theorem in a few less familiar directions.

The Fundamental Theorem of Algebra (FTA) says that an nth degree polynomial over the complex numbers has n roots. The theorem is commonly presented in high school algebra, but it’s not proved in high school and it’s not proved using algebra! A math major might see a proof midway through college.

You’re most likely to see a proof of the Fundamental Theorem of Algebra in a course in complex analysis. That is because the FTA depends on analytical properties of the complex numbers, not just their algebraic properties. It is an existence theorem that depends on the topological completeness of the complex numbers, and so it cannot be proved from the algebraic properties of the complex numbers alone.

(The dividing lines between areas of math, such as between algebra and analysis, are not always objective or even useful. And for some odd reason, some people get touchy about what belongs where. But the categorization of math subjects implicit in the creation of course catalogs puts the proof of the FTA on the syllabus of complex analysis courses, sometimes topology courses, but not algebra courses.)

Gauss offered a proof of the FTA in 1799, but his proof was flawed because he implicitly assumed the Jordan curve theorem, a seemingly obvious theorem that wasn’t rigorously proven until 1911. In keeping with our the theme above, the hole was due to overlooking a topological subtlety, not due to an algebraic mistake.

The proof of the FTA using complex analysis is a quick corollary of Liouville’s theorem: If a function is analytic everywhere in the complex plane and bounded, then it must be constant. Assume a non-constant polynomial p(z) is nowhere zero. Then 1/p(z) is analytic everywhere in the complex plane. And it must be bounded because p(z) → ∞ as |z| → ∞. We have a contradiction, so p(z) must be zero for some z.

The usual statement of the FTA says that an nth degree polynomial has n roots, but the proof above only says it must have one root. But you can bootstrap your way to the full result from there. The analyst tells the algebraist “I’ve done the hard part for you. You can take it from there.”

In some sense the rest of the theorem properly belongs to algebra, because the bookkeeping of how to count zeros is more of an algebraic matter. In order for every nth degree polynomial to have exactly n roots, some roots must be counted more than once. For example, (z − 42)² has a double root at 42. From a strictly analytical view point, there is one point where (z − 42)² is zero, namely at z = 42, but it is useful, even to analysts, to count this as a double root.

It may seem pedantic, or even Orwellian, so say that some roots count more than others. But a great deal of theorems are simpler to state if we agree to count roots according to multiplicity. This idea is taken further in algebraic geometry. For example, Bézout’s theorem says that if X is an mth degree curve and Y is an nth degree curve, then X and Y intersect in mn points. It takes a lot of prep work for this theorem to have such a simple statement. You have to count intersections with multiplicity, and you have to add points at infinity.

Liouville’s theorem’s is useful for more than proving the FTA. It is also useful in studying elliptic functions. (I believe Liouville proved his theorem for this purpose. Someone into math history can check me on that.) If an analytic function is periodic in two distinct directions in the complex plane (more on that here), it must have a singularity because otherwise it would be bounded.

The FTA may seem less profound than it is. You might object that of course every nth degree polynomial has n roots: you’ve cheated by adding in the roots you need. But that’s not true. The complex numbers are created by adding in one root to one particular polynomial, namely z² = −1. The impressive part is that adding one root of one particular quadratic polynomial is enough to force all quadratic polynomials and even all higher degree polynomials also to also have roots.

This does not happen in finite fields. A finite field cannot be algebraically complete. If you extend a finite field so that every quadratic equation has a root, then only quadratic equations have roots. It does not follow that cubic polynomials, for example, have roots. You could extend it further so that all cubic equations have roots, but then 4th degree polynomials will not all have roots. If you stop the extension process after any finite number of steps, there will be some polynomials without roots.

Fundamental theorem of calculus generalized

Posted on 23 May 2020 by John

The first fundamental theorem of calculus says that integration undoes differentiation.

The second fundamental theorem of calculus says that differentiation undoes integration.

This post looks at the fine print of these two theorems, in their basic forms and when generalized to Lebesgue integration.

Second fundamental theorem of calculus

We’ll start with the second fundamental theorem because it’s simpler. In its basic form, it says that if f is a continuous function on an open interval I, and a is a point in I, then the function F defined by

$F(x) = \int_a^x f(t)\, dt$

is an antiderivative for f on the interval I, i.e.

for all x in I. In that sense differentiation undoes integration.

If we remove the requirement that f be continuous, we still have F‘ = f almost everywhere as long as f is absolutely integrable, i.e. the integral of |f| over I is finite. In more detail,

at every Lebesgue point x, i.e. every point x that satisfies

$\lim_{\varepsilon \to 0} \, \frac{1}{2\varepsilon}\int_{x - \varepsilon}^{x + \varepsilon} |f(x) - f(y)| \ dy = 0$

First fundamental theorem of calculus

The first fundamental theorem of calculus says that if the derivative of F is f and f is continuous on an interval [a, b], then

$\int_a^b f(x) \, dx = F(b) - F(a)$

So if F has a continuous derivative, then integration undoes differentiation. What if F is continuous and but differentiable at almost every point rather than at every point? Then the theorem doesn’t necessarily hold. But the theorem does hold if we require F to be absolutely continuous rather than just continuous.

A function is absolutely continuous if it maps sets of measure zero to sets of measure zero. It’s not easy to imagine continuous functions that are not absolutely continuous, but Cantor’s function, a.k.a. the Devil’s staircase, takes the Cantor set, a set of measure zero, to a set of measure one.

The usual definition of absolute continuity, equivalent to the one above, takes the ordinary definition of continuity and chops into n pieces. That is, for every ε > 0 and for every n, there exists a δ > 0 such that for any collection of n intervals of total length less than δ, the sum of the variation in f over all the intervals is less than ε. If n = 1 this is the definition of uniform continuity, so absolute continuity is a more demanding criterion than uniform continuity.

Square wave, triangle wave, and rate of convergence

Posted on 22 May 2020 by John

There are only a few examples of Fourier series that are relatively easy to compute by hand, and so these examples are used repeatedly in introductions to Fourier series. Any introduction is likely to include a square wave or a triangle wave [1].

By square wave we mean the function that is 1 on [0, 1/2] and −1 on [1/2, 1], extended to be periodic.

square wave

Its nth Fourier (sine) coefficient is 4/nπ for odd n and 0 otherwise.

By triangle wave we mean 1 − |2x − 1|, also extended to be periodic.

triangle wave

Its nth Fourier (cosine) coefficient is −4/(nπ)² for odd n and 0 otherwise.

This implies the Fourier coefficients for the square wave are O(1/n) and the Fourier coefficients for the triangle wave are O(1/n²). (If you’re unfamiliar with big-O notation, see these notes.)

In general, the smoother a function is, the faster its Fourier coefficients approach zero. For example, we have the following two theorems.

If a function f has k continuous derivatives, then the Fourier coefficients are O(1/n^k).
If the Fourier coefficients of f are O(1/n^k+1+ε) for some ε > 0 then f has k continuous derivatives.

There is a gap between the two theorems so they’re not converses.

If a function is continuous but not differentiable, the Fourier coefficients cannot decay any faster than 1/n², so the Fourier coefficients for the triangle wave decay as fast as they can for a non-differentiable function.

Clipped sine waves

Posted on 20 May 2020 by John

One source of distortion in electronic music is clipping. The highest and lowest portions of a wave form are truncated due to limitations of equipment. As the gain is increased, the sound doesn’t simply get louder but also becomes more distorted as more of the signal is clipped off.

Clipping 0.2

For example, here is what a sine wave looks like when clipped 20%, i.e. cut off to be between -0.8 and 0.8.

Sine clipped at 0.8

A simple sine wave has only one Fourier component, itself. But when we clip the sine wave, we move energy into higher frequency components. We can see that in the Fourier components below.

Fourier coefficients of sine clipped at 0.8

You can show by symmetry that the even-numbered coefficients are exactly zero.

Clipping 0.6

Here are the corresponding plots for 60% clipping, i.e. the absolute value of the signal is cut off to be 0.4. First the signal

Sine clipped at 0.8

and then its Fourier components.

Fourier coefficients of sine clipped at 0.8

Here are the first five sine waves with the amplitudes given by the Fourier coefficients.

Fourier components

And here we see how the of the sines above do a pretty good job of reconstructing the original clipped sine. We’d need an infinite number of Fourier components to exactly reconstruct the original signal, but the first five components do most of the work.

Adding up the first five Fourier components

Continuous range of clipping

Next let’s look at the ratio of the energy in the 3rd component to that of the 1st component as we continuously vary the amount of clipping.

Ratio of energy in 3rd harmonic to fundamental

Now for the 5th harmonic. This one is interesting because it’s not strictly increasing but rather has a little bump before it starts increasing.

Ratio of energy in 5th harmonic to fundamental

Finally, here’s the ratio of the energy in all higher frequencies to the energy in the fundamental.

Ratio of energy in all higher frequences combined to fundamental

Inverse optimization

Posted on 18 May 2020 by John

This morning Patrick Honner posted the image below on Twitter.

The image was created by Robert Bosch by solving a “traveling salesman” problem, finding a nearly optimal route for passing through 12,000 points.

I find this interesting for a couple reasons. For one thing, I remember when the traveling salesman problem was considered intractable. And in theory it still is! But in practice it is not difficult now to find nearly optimal solutions for very large traveling salesman problems.

Another reason I find this interesting is that there is a higher-level problem in the background, the problem of constructing a problem whose solution gives you a desired image.

Robert Bosch is showing us a solution to two problems. The image above is the solution to an optimization problem, and the problem it solves is itself the solution to another problem in the background. He’s doing a sort of inverse optimization, searching for an optimization problem. He goes into some of his techniques in his recent book Opt Art which I wrote about here.

This gets back to an argument I’ve had with statisticians who are resistant to innovative methods. They’ll say “But our method is optimal. You can’t do any better by definition. End of discussion!” But every method is optimal by some criteria. The question is whether the method is optimal by criteria appropriate to the problem at had or whether it is only optimal according to criteria that were mathematically convenient at the time it was formulated.

Robert Bosch has shown beautifully that any image can be viewed as the solution to an optimization problem. Of course the problem has to be crafted for the image. If he were solving a naturally occurring problem, such as planning a tour for an actual sales rep, and the solution had an uncanny resemblance to a piece of art by Michelangelo, that would be astounding.

Month: May 2020