Experiments with the significance of autocorrelation

Given a sequence of numbers {x_j} of length {L} one may want to look for evidence of its periodic behavior. One way to do this is by computing autocorrelation, the correlation of the sequence with a shift of itself. Here is one reasonable way to do so: for lag values {\ell=1,\dots, \lfloor L/2 \rfloor} compute the correlation coefficient of {(x_1,\dots, x_{L-\ell}} with {(x_{\ell+1},\dots, x_L)}. That the lag does not exceed {L/2} ensures the entire sequence participates in the computation, so we are not making a conclusion about its periodicity after comparing a handful of terms at the beginning and the end. In other words, we are not going to detect periodicity if the period is more than half of the observed time period.

Having obtained the correlation coefficients, pick one with the largest absolute value; call it R. How large does R have to be in order for us to conclude the correlation is not a fluke? The answer depends on the distribution of our data, but an experiment can be used to get some idea of likelihood of large R.

I picked {x_j} independently from the standard normal distribution, and computed {r} as above. After 5 million trials with a sequence of length 100, the distribution of R was as follows:

ac100-5M
Extremal correlation coefficient in a sequence of length 100

Based on this experiment, the probability of obtaining |R| greater than 0.5 is less than 0.0016. So, 0.5 is pretty solid evidence. The probability of {|R| > 0.6} is two orders of magnitude less, etc. Also, |R| is unlikely to be very close to zero unless the data is structured in some strange way. Some kind of correlation ought to be present in the white noise.

Aside: it’s not easy to construct perfectly non-autocorrelated sequences for the above test. For length 5 an example is 1,2,3,2,3. Indeed, (1,2,3,2) is uncorrelated with (2,3,2,3) and (1,2,3) is uncorrelated with (3,2,3). For length 6 and more I can’t construct these without filling them with a bunch of zeros.

Repeating the experiment with sequences of length 1000 shows a tighter distribution of R: now |R| is unlikely to be above 0.2. So, if a universal threshold is to be used here, we need to adjust R based on sequence length.

ac1000-5M
Extremal correlation coefficient in a sequence of length 1000

I did not look hard for statistical studies of this subject, resorting to an experiment. Experimentally obtained p-values are pretty consistent for the criterion {L^{0.45}|R| > 4}. The number of trials was not very large (10000) so there is some fluctuation, but the pattern is clear.
 

Length, L P(L0.45|R| > 4)
100 0.002
300 0.0028
500 0.0022
700 0.0028
900 0.0034
1100 0.0036
1300 0.0039
1500 0.003
1700 0.003
1900 0.0042
2100 0.003
2300 0.0036
2500 0.0042
2700 0.0032
2900 0.0043
3100 0.0042
3300 0.0025
3500 0.0031
3700 0.0027
3900 0.0042

Naturally, all this depends on the assumption of independent normal variables.

And this is the approach I took to computing r in Python:

import numpy as np
n = 1000
x = np.random.normal(size=(n,))
acorr = np.correlate(x, x, mode='same')
acorr = acorr[n//2+1:]/(x.var()*np.arange(n-1, n//2, -1))
r = acorr[np.abs(acorr).argmax()]

Relating integers by differences of reciprocals

Let’s say that two positive integers m, n are compatible if the difference of their reciprocal is the reciprocal of an integer: that is, mn/(m-n) is an integer. For example, 2 is compatible with 1, 3, 4, and 6 but not with 5. Compatibility is a symmetric relation, which we’ll denote {m\sim n} even though it’s not an equivalence. Here is a chart of this relation, with red dots indicating compatibility.

Capture
Compatibility of integers up to 200

Extremes

A few natural questions arise, for any given {n}:

  1. What is the greatest number compatible with {n}?
  2. What is the smallest number compatible with {n}?
  3. How many numbers are compatible with {n}?

Before answering them, let’s observe that {m\sim n} if and only if {|m-n|} is a product of two divisors of {n}. Indeed, for {mn/(m-n)} to be an integer, we must be able to write {|n-m|} as the product of a divisor of {m} and a divisor of {n}. But a common divisor of {m} and {m-n} is also a divisor of {n}.

Of course, “a product of two divisors of {n}” is the same as “a divisor of {n^2}“. But it’s sometimes easier to think in terms of the former.

Question 1 is now easy to answer: the greatest number compatible with {n} is {n+n^2 = n(n+1)}.

But there is no such an easy answer to Questions 2 and 3, because of the possibility of overshooting into negative when subtracting a divisor of {n^2} from {n}. The answer to Question 2 is {n-d} where {d} is the greatest divisor of {n^2} than is less than {n}. This is the OEIS sequence A063428, pictured below.

min_jump10000
The smallest compatible number for n

The lines are numbers with few divisors: for a prime p, the smallest compatible number is p-1, while for 2p it is p, etc.

The answer to Question 3 is: the number of distinct divisors of {n^2}, plus the number of such divisors that are less than {n}. This is the OEIS sequence A146564.

Chains

Since consecutive integers are compatible, every number {n} is a part of a compatible chain {1\sim \cdots\sim n}. How to build a short chain like this?

Strategy A: starting with {n}, take the smallest integer compatible with the previous one. This is sure to reach 1. But in general, this greedy algorithm is not optimal. For n=22 it yields 22, 11, 10, 5, 4, 2, 1 but there is a shorter chain: 22, 18, 6, 2, 1.

Strategy B: starting with {1}, write the greatest integer compatible with the previous one (which is {k(k+1)} by the above). This initially results in the OEIS sequence A007018: 1, 2, 6, 42, 1806, … which is notable, among other things, for giving a constructive proof of the infinitude of primes, even with a (very weak) lower density bound. Eventually we have to stop adding {k^2} to {k} every time; so instead add the greatest divisor of {k^2} such that the sum does not exceed {n}. For 22, this yields the optimal chain 1, 2, 6, 18, 22 stated above. But for 20 strategy B yields 1, 2, 6, 18, 20 while the shortest chain is 1, 2, 4, 20. Being greedy is not optimal, in either up or down direction.

Strategy C is not optimal either, but it is explicit and provides a simple upper bound on the length of a shortest chain. It uses the expansion of n in the factorial number system which is the sum of factorials k! with coefficients less than k. For example {67 = 2\cdot 4! + 3\cdot 3! + 0\cdot 2! + 1\cdot 1!}, so its factorial representation is 2301.

If n is written as abcd in the factorial system, then the following is a compatible chain leading to n (possibly with repetitions), as is easy to check:

1, 10, 100, 1000, 1000, a000, ab00, abc0, abcd

In the example with 67, this chain is

1, 10, 100, 1000, 2000, 2300, 2300, 2301

in factorial notation, which converts to decimals as

1, 2, 6, 24, 48, 66, 66, 67

The repeated 66 (due to 0 in 2301) should be removed.

Thus, for an even number s, there is a chain of length at most s leading to any integer that is less than (s/2+1)!.

As a consequence, the smallest possible length of chain leading to n is {O(\log n/\log \log n)}. Is this the best possible O-estimate?

All of the strategies described above produce monotone chains, but the shortest chain is generally not monotone. For example, 17 can be reached by the non-monotone chain 1, 2, 6, 18, 17 of length 5 but any monotone chain will have length at least 6.

The smallest-chain-length sequence is 1, 2, 3, 3, 4, 3, 4, 4, 4, 4, 5, 4, 5, 5, 4, 5, 5, 4, 5, 4, 5, 5, 5, 4, 5, … which is not in OEIS.  Here is its plot:

stepde
Smallest chain length for integers from 1 to 1000

The sequence 1, 2, 3, 5, 11, 29, 67, 283, 2467,… lists the smallest numbers that require a chain of given length — for example, 11 is the first number that requires a chain of length 5, etc. Not in OEIS; the rate of growth is unclear.

Top 10 xkcd comics according to Stack Overflow

Sorted according to the number of Stack Overflow posts (as of now) in which the comic is linked. The posts themselves can be seen by clicking the “posts” link.

#10: Standards (22 posts)

#9: Compiling (29 posts)

#8: Python (33 posts)

#7: Regular Expressions (36 posts)

#6: ISO 8601 (50 posts)

#5: Password Strength (60 posts)

#4: Wisdom of the Ancients (63 posts)

#3: goto (64 posts)

#2: Random Number (74 posts)

#1: Exploits of a Mom (680 posts)

Need for speed vs bounded position and acceleration

You are driving a car with maximal acceleration (and deceleration) A on a road that’s been blocked off in both directions (or, if you prefer, on the landing strip of an aircraft carrier). Let L be the length of the road available to you.
What is the maximal speed you can reach?

Besides A and L, the answer also depends on your mood: do you want to live, or are you willing to go out in a blaze of glory? In the latter case the answer is obvious: position the car at one end of the interval, and put the pedal to the metal. The car will cover the distance L within the time {\sqrt{2L/A}}, reaching the speed {v=\sqrt{2AL}} at the end. In the former scenario one has to switch to the brake pedal midway through the distance, so the maximal speed will be attained at half the length, {\sqrt{AL}}.

Rephrased in mathematical terms: if {f} is a twice differentiable function and {M_k = \sup|f^{(k)}|} for {k=0,1,2}, then {M_1^2 \le 4M_0M_2} if {f} is defined on a half-infinite interval, and {M_1^2 \le 2M_0M_2} if the domain of {f} is the entire line. To connect the notation, just put {L=2M_0} and {A=M_1} in the previous paragraph… and I guess some proof other than “this is obvious” is called for, but it’s not hard to find one: this is problem 5.15 in Rudin’s Principles of Mathematical Analysis.

Perhaps more interesting is to study the problem in higher dimensions: one could be driving in a parking lot of some shape, etc. Let’s normalize the maximal acceleration as 1, keeping in mind it’s a vector. Given a set E, let S(E) be the square of maximal speed attainable by a unit-acceleration vehicle which stays in E indefinitely. Also let U(E) be the square of maximal speed one can attain while crashing out of bounds after the record is set. Squaring makes these quantities scale linearly with the size of the set. Both are monotone with respect to set inclusion. And we know what they are for an interval of length L: namely, {S = L} and {U=2L}, so that gives some lower bounds for sets that contain a line interval.

When E is a circle of radius 1, the best we can do is to drive along it with constant speed 1; then the centripetal acceleration is also 1. Any higher speed will exceed the allowable acceleration in the normal direction, never mind the tangential one. So, for a circle both S and U are equal to its radius.

On the other hand, if E is a disk of radius R, then driving along its diameter is better: it gives {S\ge 2R} and {U\ge 4R}.

Some questions:

  1. If E is a convex set of diameter D, is it true that {S(E) = D} and {U(E) = 2D}?
  2. Is it true that {U\le 2S} in general?
  3. How to express S and U for a smooth closed curve in terms of its curvature? They are not necessarily equal (like they are for a circle): consider thin ellipses converging to a line segment, for which S and U approach the corresponding values for that segment.

The answer to Question 1 is yes. Consider the orthogonal projection of E, and of a trajectory it contains, onto some line L. This does not increase the diameter or the acceleration; thus, the one-dimensional result implies that the projection of velocity vector onto L does not exceed {\sqrt{D}} (or {\sqrt{2D}} for the crashing-out version). Since L was arbitrary, it follows that {S(E) \le D} and {U(E) \le 2D}. These upper bounds hold for general sets, not only convex ones. But when E is convex, we get matching lower bounds by considering the longest segment contained in E.

I don’t have answers to questions 2 and 3.

Lightness, hyperspace, and lower oscillation bounds

When does a map {f\colon X\to Y} admit a lower “anti-continuity” bound like {d_Y(f(a),f(b))\ge \lambda(d_X(a,b))} for some function {\lambda\colon (0,\infty)\to (0, \infty)} and for all {a\ne b}? The answer is easy: {f} must be injective and its inverse must be uniformly continuous. End of story.

But recalling what happened with diameters of connected sets last time, let’s focus on the inequality {\textrm{diam}\, f(E)\ge \lambda (\textrm{diam}\, E)} for connected subsets {E\subset X}. If such {\lambda} exists, the map f has the LOB property, for “lower oscillation bound” (oscillation being the diameter of image). The LOB property does not require {f} to be injective. On the real line, {f(x)=|x|} satisfies it with {\lambda(\delta)=\delta/2}: since it simply folds the line, the worst that can happen to the diameter of an interval is to be halved. Similarly, {f(x)=x^2} admits a lower oscillation bound {\lambda(\delta) = (\delta/2)^2}. This one decays faster than linear at 0, indicating some amount of squeezing going on. One may check that every polynomial has the LOB property as well.

On the other hand, the exponential function {f(x)=e^x} does not have the LOB property, since {\textrm{diam}\, f([x,x+1])} tends to {0} as {x\to-\infty}. No surprise there; we know from the relation of continuity and uniform continuity that things like that happen on a non-compact domain.

Also, a function that is constant on some nontrivial connected set will obviously fail LOB. In topology, a mapping is called light if the preimage of every point is totally disconnected, which is exactly the same as not being constant on any nontrivial connected set. So, lightness is necessary for LOB, but not sufficient as {e^x} shows.

Theorem 1: Every continuous light map {f\colon X\to Y} with compact domain {X} admits a lower oscillation bound.

Proof. Suppose not. Then there exists {\epsilon>0} and a sequence of connected subsets {E_n\subset X} such that {\textrm{diam}\, E_n\ge \epsilon} and {\textrm{diam}\, f(E_n)\to 0}. We can assume {E_n} compact, otherwise replace it with its closure {\overline{E_n}} which we can because {f(\overline{E_n})\subset \overline{f(E_n)}}.

The space of nonempty compact subsets of {X} is called the hyperspace of {X}; when equipped with the Hausdorff metric, it becomes a compact metric space itself. Pass to a convergent subsequence, still denoted {\{E_n\}}. Its limit {E} has diameter at least {\epsilon}, because diameter is a continuous function on the hyperspace. Finally, using the uniform continuity of {f} we get {\textrm{diam}\, f(E) = \lim \textrm{diam}\, f(E_n) = 0}, contradicting the lightness of {f}. {\quad \Box}

Capture
A homeomorphism lacking LOB

Here is another example to demonstrate the importance of compactness (not just boundedness) and continuity: on the domain {X = \{(x,y)\colon 0 < x < 1, 0 < y < 1\}} define {f(x,y)=(x,xy)}. This is a homeomorphism, the inverse being {(u,v)\mapsto (u, v/u)}. Yet it fails LOB because the image of line segment {\{x\}\times (0,1)} has diameter {x}, which can be arbitrarily close to 0. So, the lack of compactness hurts. Extending {f} to the closed square in a discontinuous way, say by letting it be the identity map on the boundary, we see that continuity is also needed, although it’s slightly non-intuitive that one needs continuity (essentially an upper oscillation bound) to estimate oscillation from below.

All that said, on a bounded interval of real line we need neither compactness nor continuity.

Theorem 2: If {I\subset \mathbb R} is a bounded interval, then every light map {f\colon I\to Y} admits a lower oscillation bound.

Proof. Following the proof of Theorem 1, consider a sequence of intervals {(a_n, b_n)} such that {b_n-a_n\ge \epsilon} and {\textrm{diam}\, f((a_n,b_n))\to 0}. There is no loss of generality in considering open intervals, since it can only make the diameter of the image smaller. Also WLOG, suppose {a_n\to a} and {b_n\to b}; this uses the boundedness of {I}. Consider a nontrivial closed interval {[c,d]\subset (a,b)}. For all sufficiently large {n} we have {[c,d]\subset (a_n,b_n)}, which implies {\textrm{diam}\, f([c,d])\le \textrm{diam}\, f((a_n,b_n))\to 0}. Thus {f} is constant on {[c,d]}, a contradiction. {\quad \Box}

The property that distinguishes real line here is that nontrivial connected sets have nonempty interior. The same works on the circle and various tree-like spaces, but fails for spaces that don’t look one-dimensional.

Continuity and diameters of connected sets

The definition of uniform continuity (if it’s done right) can be phrased as: {f\colon X\to Y} is uniformly continuous if there exists a function {\omega\colon (0,\infty)\to (0,\infty)}, with {\omega(0+)=0}, such that {\textrm{diam}\, f(E)\le \omega (\textrm{diam}\, E)} for every set {E\subset X}. Indeed, when {E} is a two-point set {\{a,b\}} this is the same as {|f(a)-f(b)|\le \omega(|a-b|)}, the modulus of continuity. Allowing general sets {E} does not change anything, since the diameter is determined by two-point subsets.

Does it make a difference if we ask for {\textrm{diam}\, f(E)\le \omega (\textrm{diam}\, E)} only for connected sets {E}? For functions defined on the real line, or on an interval of the line, there is no difference: we can just consider the intervals {[a,b]} and obtain

{|f(a)-f(b)|\le \textrm{diam}\, f([a,b]) \le \omega(|a-b|)}

as before.

However, the situation does change for maps defined on a non-convex domain. Consider the principal branch of square root, {f(z)=\sqrt{z}}, defined on the slit plane {G=\mathbb C\setminus (-\infty, 0]}.

sqrt
Conformal map of a slit domain is not uniformly continuous

This function is continuous on {G} but not uniformly continuous, since {f(-1 \pm i y) \to \pm i } as {y\to 0+}. Yet, it satisfies {\textrm{diam}\, f(E)\le \omega(\textrm{diam}\, E)} for connected subsets {E\subset G}, where one can take {\omega(\delta)=2\sqrt{\delta}}. I won’t do the estimates; let’s just note that although the points {-1 \pm i y} are close to each other, any connected subset of {G} containing both of them has diameter greater than 1.

Capture3
These points are far apart with respect to the inner diameter metric

In a way, this is still uniform continuity, just with respect to a different metric. Given a metric space {(X,d)}, one can define inner diameter metric {\rho} on {X} by letting {\rho(a,b)} be the infimum of diameters of connected sets that contain both {a} and {b}. This is indeed a metric if the space {X} is reasonable enough (i.e., any two points are contained in some bounded connected set). On a convex subset of {\mathbb R^n}, the inner diameter metric coincides with the Euclidean metric {d_2}.

One might think that the equality {\rho=d_e} should imply that the domain is convex, but this is not so. Indeed, consider the union of three quadrants on a plane, say {A = \{(x,y) \colon x > 0\text{ or }y > 0\}}. Any two points of {A} can be connected by going up from whichever is lower, and then moving horizontally. The diameter of a right triangle is equal to its hypotenuse, which is the Euclidean distance between the points we started with.

Capture2
A non-convex domain where inner diameter metric is the same as Euclidean

Inner diameter metric comes up (often implicitly) in complex analysis. By the Riemann mapping theorem, every simply connected domain {G\subset \mathbb C}, other than {\mathbb C} itself, admits a conformal map {f\colon G\to \mathbb D} onto the unit disk {D}. This map need not be uniformly continuous in the Euclidean metric (the slit plane is one example), but it is uniformly continuous with respect to the inner diameter metric on {G}.

Furthermore, by normalizing the situation in a natural way (say, {G \supset \mathbb D} and {f(0)=0}), one can obtain a uniform modulus of continuity for all conformal maps {f} onto the unit disk, whatever the domain is. This uniform modulus of continuity can be taken of the form {\omega(\delta) = C\sqrt{\delta}} for some universal constant {C}. Informally speaking, this means that a slit domain is the worst that can happen to the continuity of a conformal map. This fact isn’t often mentioned in complex analysis books. A proof can be found in the book Conformally Invariant Processes in the Plane by Gregory Lawler, Proposition 3.85. A more elementary proof, with a rougher estimate for the modulus of continuity, is on page 15 of lecture notes by Mario Bonk.

Green envelopes

Green functions of the Laplacian are bounded in one dimensional domains, a reflection of the fact that square integrability of the derivative implies boundedness in one dimension but not in two or more. So, it’s possible to say something about the envelope {\sup_y g(x,y)} obtained by taking the supremum over source location y.

In order to satisfy the distributional equation {-\frac{d^2}{dx^2}g(x,y) = \delta_{x=y}}, Green functions have a corner at y, where the derivative jumps down by 1. Combining this property with the boundary conditions identifies the function. Let’s use the interval (0,1) as the domain.

The Dirichlet boundary condition {g(0,y) = 0 = g(1,y)} leads to the formula {g(x,y)=\min(x,y)-xy}:

green10

The upper envelope is the parabola {x(1-x)}. This fact becomes more evident if we use more functions.

green200
Dirichlet envelope

Mixed Dirichlet-Neumann conditions {g(0,y) = 0 = g_x'(1,y)} yield an even simpler formula {g(x,y) = \min(x,y)}. The envelope is boring in this case:

greendn50
Dirichlet-Neumann envelope

The pure Neumann condition {g_x'(0,y) = 0 = g_x'(1,y)} requires an adjustment to the definition of Green function, since it is incompatible with {-\frac{d^2}{dx^2}g(x,y) = \delta_{x=y}}. Instead, let’s ask for {-\frac{d^2}{dx^2}g(x,y) = \delta_{x=y} - c} where the constant c is the reciprocal of the length of the interval. This ensures that the integral of the second derivative is zero. Also, it is more convenient to use the interval (-1,1) here. The Neumann Green function for this interval is {g(x,y) = (x^2+y^2)/4 - |x-y|/2}, up to an additive constant.

greenn50large
Neumann envelope

So the envelope is also determined up to a constant. On the picture above it is {x^2/2}.

Finally, let’s consider the periodic condition {g(0,y) = g(1,y)}, {g_x'(0,y) = g_x'(1,y)}. As with the Neumann condition, the definition is adjusted to allow the second derivative integral be zero. Using the interval (-1, 1) we get {g(x,y) = (1+x\wedge y)(1-x\vee y)/2 + (x^2+y^2)/4} using the notation {x\wedge y = \min(x,y)} and {x\vee y = \max(x,y)} for brevity. Note that the first term, {(1+x\wedge y)(1-x\vee y)/2}, is the Dirichlet Green function for this interval. The additional quadratic term brings the integral of the second derivative to zero, which ensures the periodic condition for the first derivative. And the upper envelope is, up to a constant…

greenp50
Periodic envelope

a constant. Not exciting but makes perfect sense, considering that we are dealing with the Green function of a circle, where moving the source amounts to shifting the function.