## Boundary value problems: not all explicit solutions are useful

Consider this linear differential equation: ${y''(t) + 4y'(t) + 2t y(t) = 7}$ with boundary conditions ${y(0)=1}$ and ${y(1)=0}$. Nothing looks particularly scary here. Just one nonconstant coefficient, and it’s a simple one. Entering this problem into Wolfram Alpha produces the following explicit solution:

I am not sure how anyone could use this formula for any purpose.

Let us see what simple linear algebra can do here. The differential equation can be discretized by placing, for example, ${4}$ equally spaces interior grid points on the interval: ${t_k = k/5}$, ${k=1, \dots, 4}$. The yet-unknown values of ${y}$ at these points are denoted ${y_1,\dots, y_4}$. Standard finite-difference formulas provide approximate values of ${y'}$ and ${y''}$:

${\displaystyle y'(t) \approx \frac{y(t+h) - y(t-h)}{2h}}$

${\displaystyle y''(t) \approx \frac{y(t+h) - 2y(t) + y(t-h)}{h^2}}$

where ${h}$ is the step size, ${1/5}$ in our case. Stick all this into the equation: we get 4 linear equations, one for each interior point. Namely, at ${t_1 = 1/5}$ it is

${\displaystyle \frac{y_2 - 2y_1 + 1}{(1/5)^2} + 4 \frac{y_2 - 1}{2/5} + 2\cdot \frac15 y_1 = 7 }$

(notice how the condition ${y(0)=1}$ is used above), at ${t_2 = 2/5}$ it is

${\displaystyle \frac{y_3 - 2y_2 + y_1}{(1/5)^2} + 4 \frac{y_3 - y_1}{2/5} + 2 \cdot \frac25 y_2 = 7 }$

and so on. Clean up this system and put it in matrix form:

${\displaystyle \begin{pmatrix} -49.6 & 35 & 0 & 0 \\ 15 & -49.2 & 35 & 0 \\ 0 & 15 & -48.8 & 35 \\ 0 & 0 & 15 & -48.4 \end{pmatrix} \vec y = \begin{pmatrix} -8 \\ 7 \\ 7 \\ 7 \end{pmatrix} }$

This isn’t too hard to solve even with pencil and paper. The solution is

${\displaystyle y = \begin{pmatrix} -0.2859 \\ -0.6337\\ -0.5683\\ -0.3207 \end{pmatrix}}$

It can be visualized by plotting 4 points ${(t_k, y_k)}$:

Not particularly impressive is it? And why are all these negative y-values in a problem with boundary condition ${y(0)=1}$? They do not really look like they want to approach ${1}$ at the left end of the interval. But let us go ahead and plot them together with the boundary conditions, using linear interpolation in between:

Or better, use cubic spline interpolation, which only adds another step of linear algebra (see Connecting dots naturally) to our computations.

This begins to look believable. For comparison, I used a heavier tool: BVP solver from SciPy. Its output is the red curve below.

Those four points we got from a 4-by-4 system, solvable by hand, pretty much tell the whole story. At any rate, they tell a better story than the explicit solution does.

## The distribution of pow(2, n, n)

For a positive integer ${n}$ let ${f(n)}$ be the remainder of the division of ${2^n}$ by ${n}$. This is conveniently computed in Python as pow(2, n, n). For example, the first 20 values are returned by

[pow(2, n, n) for n in range(1, 21)]

and they are:

[0, 0, 2, 0, 2, 4, 2, 0, 8, 4, 2, 4, 2, 4, 8, 0, 2, 10, 2, 16]

Obviously, ${f(n)=0}$ if and only if ${n}$ is a power of ${2}$. It also looks like ${f}$ is always even, with powers of 2 dominating the list of its values. But ${f}$ does take on odd values, although this does not happen often: only 6 times in the first 100 integers.

[(n, pow(2, n, n)) for n in range(1, 101) if pow(2, n, n) % 2]

returns

[(25, 7), (45, 17), (55, 43), (91, 37), (95, 13), (99, 17)]

It is conjectured that the range of ${f}$ consists of all nonnegative integers except 1. Here is a proof that ${f(n)\ne 1}$, provided by Max Alexeyev on the OEIS page for A036236

Suppose ${2^n \equiv 1 \bmod n}$ for some ${n>1}$. Let ${p}$ be the smallest prime divisor of ${n}$. Then ${2^n\equiv 1 \bmod p}$. This means that the order of ${2}$ in the multiplicative group ${(\mathbb Z/ p \mathbb Z)^*}$ is a divisor of ${n}$. But this group has ${p-1}$ elements, and the only divisor of ${n}$ that is smaller than ${p}$ is ${1}$. Thus, the order of ${2}$ is ${1}$, which means ${2\equiv 1 \bmod p}$, which is absurd.

The fact that the smallest ${n}$ with ${f(n) = 3}$ is 4700063497 deserves to be mentioned here. But I would like to consider the most frequent values of ${f}$ instead. Experimentally, they are found like this:

from collections import Counter
freq = Counter(pow(2, n, n) for n in range(1, 100000000 + 1))
print(freq.most_common(20))

which prints 20 pairs (value, frequency):

(2, 5763518),
(4, 3004047),
(8, 2054167),
(16, 1569725),
(64, 1076293),
(256, 824438),
(512, 737450),
(1024, 668970),
(32, 638063),
(4096, 569705),
(128, 466015),
(65536, 435306),
(262144, 389305),
(1048576, 351723),
(2097152, 326926),
(16384, 246533),
(16777216, 243841),
(32768, 232037),
(2048, 153537),
(8192, 131614)

These are all powers of 2… but not exactly in the order one might expect. I repeated the experiment for the first ${10^8k}$ integers with ${k=1, 2, 3, 4, 5, 6}$. The frequency order remained stable at the top. These are the most common 11 values, with their frequency (in %) listed for the aforementioned six ranges.

2, [5.8, 5.5, 5.4, 5.3, 5.3, 5.2]
4, [3.0, 2.9, 2.8, 2.8, 2.7, 2.7]
8, [2.1, 2.0, 1.9, 1.9, 1.9, 1.8]
16, [1.6, 1.5, 1.5, 1.4, 1.4, 1.4]
64, [1.1, 1.0, 1.0, 1.0, 1.0, 1.0]
256, [0.8, 0.8, 0.8, 0.8, 0.7, 0.7]
512, [0.7, 0.7, 0.7, 0.7, 0.7, 0.7]
1024, [0.7, 0.6, 0.6, 0.6, 0.6, 0.6]
32, [0.6, 0.6, 0.6, 0.6, 0.6, 0.6]
4096, [0.6, 0.5, 0.5, 0.5, 0.5, 0.5]
128, [0.5, 0.4, 0.4, 0.4, 0.4, 0.4]

It seems that 32 and 128 are consistently less common than one might expect.

The most common value, 2, is contributed by primes and pseudoprimes to base 2. The value of 4 appears when ${n = 2p}$ with ${p}$ prime (but not only then). Still, the primes having zero density makes it difficult to explain the frequency pattern based on primality. A more convincing reason could be: when ${n=2k}$ is even, we are computing ${4^k \bmod 2k}$ and that is likely to be a power of ${4}$. This boosts the frequency of the even (generally, composite) powers of 2 in the sequence.

## Numerical integration visualized

scipy.integrate.quad is a popular method of numerical integration with Python. Let’s see how it chooses the points at which to evaluate the function being integrated. Begin with a simple example, the exponential function.

The blue dots indicate the evaluation points, their y-coordinate being the order of evaluation. So, it begins with x=0, continues with something close to -1, etc. The function, drawn in red, is not to its true vertical scale.

The placement of dots gives away the method of integration: it is the Gauss-Kronrod quadrature with 10 Gauss nodes and 21 Kronrod nodes, abbreviated (G10, K21).  The Gauss nodes are included in the Kronrod nodes, and of course the function is not evaluated there again. The evaluation process is slightly out of order in that 0 is a Kronrod node but comes first, followed by 10 Gauss nodes, followed by 10 remaining Kronrod nodes. The process ends there, as the comparison of G10 and K21 results shows the necessary precision was reached.

The square root on [0, 1] is not such a nice function, so (G10, K21) does not reach the required precision at once. The interval is bisected again and again until it does.

Surely the cube root is even worse. Or is it?

The nodes are symmetric about the midpoint of the interval of integration. So for any odd function on an interval symmetric about 0 we get G10 = 0 and K21 = 0, and the process stops at once. To see the effect of the cube root singularity, one has to use a non-symmetric interval such as [-1, 2].

That’s still fewer subdivisions than for the square root: cancellation between the left and right neighborhoods of 0 still helps. Let’s look at smooth functions next.

Rapid oscillations forces subdivisions here. There are also other reasons for subdivision, such as the loss of analyticity:

Although exp(-1/x) is infinitely differentiable on this interval, the fact that it is not analytic at 0 makes it a more difficult integrand that exp(x). Finally, an analytic function with no oscillation which still needs a bunch of subintervals:

This is the standard example used to illustrate the Runge phenomenon. Although we are not interpolating here, numerical integration is also influenced by the function having a small radius of convergence of its Taylor series. After all, G10 can be thought of as degree-9 interpolation of the given function at the Gauss nodes (the zeros of a Legendre polynomial),  with the formula returning the integral of the interpolating polynomial.

The code used to plot these things:

import numpy as np
import matplotlib.pyplot as plt
global eval_points
f = lambda x: np.exp(x)
def f_int(x):
eval_points.append(x)
return f(x)

eval_points = []
a, b = -1, 1
n = len(eval_points)
t = np.linspace(a, b, 1000)
y = f(t)
yy = n*(y - y.min())/(y.max() - y.min())
plt.plot(t, yy, 'r')
plt.plot(eval_points, np.arange(0, n), '.')
plt.show()


## Critical points of a cubic spline

The choice of piecewise polynomials of degree 3 for interpolation is justifiably popular: even-degree splines are algebraically awkward to construct, degree 1 is simply piecewise linear interpolation (not smooth), and degree 5, while feasible, entails juggling too many coefficients. Besides, a cubic polynomial minimizes the amount of wiggling (the integral of second derivative squared) for given values and slopes at the endpoints of an interval. (Recall Connecting dots naturally.)

But the derivative of a cubic spline is a quadratic spline. And one needs the derivative to find the critical points. This results in an awkward example in SciPy documentation, annotated with “(NB: sproot only works for order 3 splines, so we fit an order 4 spline)”.

Although not implemented in SciPy, the task of computing the roots of a quadratic spline is a simple one. Obtaining the roots from the internal representation of a quadratic spline in SciPy (as a linear combination of B-splines) would take some work and reading. But a quadratic polynomial is determined by three values, so sampling it at three points, such as two consecutive knots and their average, is enough.

Suppose we know the values of a quadratic polynomial q at -1, 0, 1, and wish to find if it has roots between -1 and 1. Let’s normalize so that q(0)=1, and let x = q(-1), y = q(1). If either x or y is negative, there is definitely a root on the interval. If they are positive, there is still a chance: we need the parabola to be concave up, have a minimum within [-1, 1], and for the minimum to be negative. All of this is easily determined once we note that the coefficients of the polynomial are a = (x+y)/2 – 1, b = (y-x)/2, and c = 1.

The inequality ${(x-y)^2 \ge 8(x+y-2)}$ ensures the suitable sign of the discriminant. It describes a parabola with vertex (1, 1) and focus (2, 2), contained in the first quadrant and tangent to the axes at (4, 0) and (0, 4). Within the orange region there are no real roots.

The line x+y=2, tangent to the parabola at its vertex, separates convex and concave parabolas. While concavity in conjunction with x, y being positive definitely precludes having roots in [-1, 1], slight convexity is not much better: it results in real roots outside of the interval. Here is the complete picture: green means there is a root in [-1, 1], orange means no real roots, red covers the rest.

## Back to splines

Since the derivative of a spline is implemented in SciPy (B-splines have a nice formula for derivatives), all we need is a root-finding routine for quadratic splines. Here it is, based on the above observations but using built-in NumPy polynomial solver np.roots to avoid dealing with various special cases for the coefficients.

def quadratic_spline_roots(spl):
roots = []
knots = spl.get_knots()
for a, b in zip(knots[:-1], knots[1:]):
u, v, w = spl(a), spl((a+b)/2), spl(b)
t = np.roots([u+w-2*v, w-u, 2*v])
t = t[np.isreal(t) & (np.abs(t) <= 1)]
roots.extend(t*(b-a)/2 + (b+a)/2)
return np.array(roots)

A demonstration, which plots the spline (blue), its critical points (red), and original data points (black) as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import InterpolatedUnivariateSpline

x = np.arange(7)
y = np.array([3, 1, 1, 2, 2, 4, 3])
f = InterpolatedUnivariateSpline(x, y, k=3)

t = np.linspace(x[0], x[-1], 500)
plt.plot(t, f(t))
plt.plot(x, y, 'kd')
plt.plot(crit_pts, f(crit_pts), 'ro')
plt.show()

## The sum of pairwise distances and the square of CDF

Suppose we have ${n}$ real numbers ${x_0,\dots, x_{n-1}}$ and want to find the sum of all distances ${|x_j-x_k|}$ over ${j < k}$. Why? Maybe because over five years ago, the gradient flow of this quantity was used for "clustering by collision" (part 1, part 2, part 3).

If I have a Python console open, the problem appears to be solved with one line:

>>> 0.5 * np.abs(np.subtract.outer(x, x)).sum()

where the outer difference of x with x creates a matrix of all differences ${x_i-x_j}$, then absolute values are taken, and then they are all added up. Double-counted, hence the factor of 0.5.

But trying this with, say, one million numbers is not likely to work. If each number takes 8 bytes of memory (64 bits, double precision), then the array x is still pretty small (under 8 MB) but a million-by-million matrix will require over 7 terabytes, and I won’t have that kind of RAM anytime soon.

In principle, one could run a loop adding these values, or store the matrix on a hard drive. Both are going to take forever.

There is a much better way, though. First, sort the numbers in nondecreasing order; this does not require much time or memory (compared to quadratic memory cost of forming a matrix). Then consider the partial sums ${s_k = x_0+\dots+x_k}$; the cost of computing them is linear in time and memory. For each fixed ${k}$, the sum of distances to ${x_j}$ with ${j is simply ${kx_k - s_{k-1}}$, or, equivalently, ${(k+1)x_k - s_k}$. So, all we have to do is add these up. Still one line of code (after sorting), but a much faster one:

>>> x.sort()
>>> (np.arange(1, n+1)*x - np.cumsum(x)).sum()

For example, x could be a sample from some continuous distribution. Assuming the distribution has a mean (i.e., is not too heavy tailed), the sum of all pairwise distances grows quadratically with n, and its average approaches a finite limit. For the uniform distribution on [0, 1] the computation shows this limit is 1/3. For the standard normal distribution it is 1.128… which is not as recognizable a number.

As ${n\to \infty}$, the average distance of a sample taken from a distribution converges to the expected value of |X-Y| where X, Y are two independent variables with that distribution. Let’s express this in terms of the probability density function ${p}$ and the cumulative distribution function ${\Phi}$. By symmetry, we can integrate over ${x> y}$ and double the result:

${\displaystyle \frac12 E|X-Y| = \int_{-\infty}^\infty p(x)\,dx \int_{-\infty}^x (x-y) p(y)\,dy}$

Integrate by parts in the second integral: ${p(y) = \Phi'(y)}$, and the boundary terms are zero.

${\displaystyle \frac12 E|X-Y| = \int_{-\infty}^\infty p(x)\,dx \int_{-\infty}^x \Phi(y)\,dy}$

Integrate by parts in the other integral, throwing the derivative onto the indefinite integral and thus eliminating it. There is a boundary term this time.

${\displaystyle \frac12 E|X-Y| = \Phi(\infty) \int_{-\infty}^\infty \Phi(y)\,dy - \int_{-\infty}^\infty \Phi(x)^2\,dx}$

Since ${\Phi(\infty) = 1}$, this simplifies nicely:

${\displaystyle \frac12 E|X-Y| = \int_{-\infty}^\infty \Phi(x) (1-\Phi(x))\,dx}$

This is a lot neater than I expected: ${E|X-Y|}$ is simply the integral of ${2\Phi(1-\Phi)}$. I don’t often see CDF squared, like here. Some examples: for the uniform distribution on [0,1] we get

${\displaystyle E|X-Y| = \int_0^1 2x(1-x)\,dx = \frac13}$

and for the standard normal, with ${\Phi(x) = (1+\mathrm{erf}\,(x/\sqrt{2}))/2}$, it is

${\displaystyle \int_{-\infty}^\infty \frac12 \left(1-\mathrm{erf}\,(x/\sqrt{2}) ^2 \right)\,dx = \frac{2}{\sqrt{\pi}}\approx 1.12838\ldots }$

The trick with sorting and cumulative sums can also be used to find, for every point ${x_k}$, the sum (or average) of distances to all other points. To do this, we don’t sum over ${k}$ but must also add ${|x_j-x_k|}$ for ${j>k}$. The latter sum is simply ${S - s_k - (n-k-1)x_k}$ where ${S}$ is the total sum. So, all we need is

>>> (2*np.arange(1,n+1)-n)*x - 2*np.cumsum(x) + x.sum()

Unfortunately, the analogous problems for vector-valued sequences are not as easy. If the Manhattan metric is used, we can do the computations for each coordinate separately, and add the results. For the Euclidean metric…

## Discrete Cosine Transforms: Trapezoidal vs Midpoint

The existence of multiple versions of Discrete Cosine Transform (DCT) can be confusing. Wikipedia explains that its 8 types are determined by how one reflects across the boundaries. E.g., one can reflect 1, 2, 3 across the left boundary as 3, 2, 1, 2, 3 or as 3, 2, 1, 1, 2, 3, and there are such choices for the other boundary too (also, the other reflection can be odd or even). Makes sense enough.

But there is another aspect to the two most used forms, DCT-I and DCT-II (types I and II): they can be expressed in terms of the Trapezoidal and Midpoint rules for integration. Here is how.

The cosines ${\cos kx}$, ${k=0,1,\dots}$ are orthogonal on the interval ${[0,\pi]}$ with respect to the Lebesgue measure. The basis of discrete transforms is that these cosines are also orthogonal with respect to some discrete measures, until certain frequency is reached. Indeed, if cosines are orthogonal with respect to measure ${\mu}$ whose support consists of ${n}$ points, then we can efficiently use them to represent any function defined on the support of ${\mu}$, and those are naturally identified with sequences of length n.

How to find such measures? It helps that some simple rules of numerical integration are exact for trigonometric polynomials up to some degree.

For example, the trapezoidal rule with n sample points exactly integrates the functions ${\cos kx}$ for ${k=0,\dots, 2n-3}$. This could be checked by converting to exponential form and summing geometric progressions, but here is a visual explanation with n=4, where the interval ${[0,\pi]}$ is represented as upper semi-circle. The radius of each red circle indicates the weight placed at that point; the endpoints get 1/2 of the weight of other sample points. To integrate ${\cos x}$ correctly, we must have the x-coordinate of the center of mass equal to zero, which is obviously the case.

Replacing ${\cos x}$ by ${\cos kx}$ means multiplying the polar angle of each sample point by ${k}$. This is what we get:

In all cases the x-coordinate of the center of mass is zero. With k=6 this breaks down, as all the weight gets placed in one point. And this is how it goes in general with integration of sines and cosines: equally spaced points work perfectly until they don’t work at all, which happens when the step size is equal to the period of the function. When ${k=2n-2}$, the period ${2\pi/k}$ is equal to ${\pi/(n-1)}$, the spacing of points in the trapezoidal rule.

The orthogonality of cosines has to do with the formula ${\cos kx \cos j x = \frac12(\cos (k-j)x) + \frac12(\cos (k+j)x)}$. Let ${\tau_n}$ be the measure expressing the trapezoidal rule on ${[0,\pi]}$ with ${n}$ sample points; so it’s the sum of point masses at ${0, \pi/(n-1), \dots , \pi}$. Then ${\{\cos k x \colon k=0,\dots, n-1\}}$ are orthogonal with respect to ${\tau_n}$ because any product ${\cos k x\cos j x }$ with ${k >j}$ taken from this range will have ${k-j, k+j \le 2n-3}$. Consequently, we can compute the coefficients of any function ${f}$ in the cosine basis as

${\displaystyle c_k = \int f(x)\cos kx\,d\tau_n(x) \bigg/ \int \cos^2 kx\,d\tau_n(x)}$

The above is what DCT-I (discrete cosine transform of type 1) does, up to normalization.

The DCT-II transform uses the Midpoint rule instead of the Trapezoidal rule. Let ${\mu_n}$ be the measure expressing the Midpoint rule on ${[0,\pi]}$ with ${n}$ sample points; it gives equal mass to the points ${(2j-1)\pi/(2n)}$ for ${j=1,\dots, n}$. These are spaced at ${\pi/n}$ and therefore the midpoint rule is exact for ${\cos kx}$ with ${k=0,\dots, 2n-1}$ which is better than what the trapezoidal rule does. Perhaps more significantly, by identifying the given data points with function values at the midpoints of subintervals we stay away from the endpoints ${0,\pi}$ where the cosines are somewhat restricted by having to have zero slope.

Let’s compare DCT-I and DCT-II on the same data set, ${y=(0, \sqrt{1}, \sqrt{2}, \dots, \sqrt{8})}$. There are 9 numbers here. Following DCT-I we place them at the sample points of the trapezoidal rule, and expand into cosines using the inner product with respect to ${\tau_n}$. Here is the plot of the resulting trigonometric polynomial: of course it interpolates the data.

But DCT-II does it better, despite having exactly the same cosine functions. The only change is that we use ${\mu_n}$ and so place the ${y}$-values along its support.

Less oscillation means the high-degree coefficients are smaller, and therefore easier to discard in order to compress information. For example, drop the last two coefficients in each expansion, keeping 6 numbers instead of 8. DCT-II clearly wins in accuracy then.

Okay, so the Midpoint rule is better, no surprise. After all, it’s in general about twice as accurate as the Trapezoidal rule. What about Simpson’s rule, would it lead to some super-efficient form of DCT? That is, why don’t we let ${\sigma_n}$ be the discrete measure that expresses Simpson’s rule and use the inner product ${\int fg\,d\sigma_n}$ for cosine expansion? Alas, Simpson’s rule on ${n}$ points is exact only for ${\cos kx}$ with ${k=0,\dots, n-2}$, which is substantially worse than either Trapezoidal or Midpoint rules. As a result, we don’t get enough orthogonal cosines with respect to ${\sigma_n}$ to have an orthogonal basis. Simpson’s rule has an advantage when dealing with algebraic polynomials, not with trigonometric ones.

Finally, the Python code used for the graphics; I did not use SciPy’s DCT method (which is of course more efficient) to keep the relation to numerical integration explicit in the code. The method trapz implements the trapezoidal rule, and the midpoint rule is just the summation of sampled values. In both cases there is no need to worry about factor dx, since it cancels out when we divide one numerical integral by the other.

import numpy as np
import matplotlib.pyplot as plt
#   Setup
y = np.sqrt(np.arange(9))
c = np.zeros_like(y)
n = y.size
#   DCT-I, trapezoidal
x = np.arange(n)*np.pi/(n-1)
for k in range(n):
c[k] = np.trapz(y*np.cos(k*x))/np.trapz(np.cos(k*x)**2)
t = np.linspace(0, np.pi, 500)
yy = np.sum(c*np.cos(np.arange(9)*t.reshape(-1, 1)), axis=1)
plt.plot(x, y, 'ro')
plt.plot(t, yy)
plt.show()
#   DCT-II, midpoint
x = np.arange(n)*np.pi/n + np.pi/(2*n)
for k in range(n):
c[k] = np.sum(y*np.cos(k*x))/np.sum(np.cos(k*x)**2)
t = np.linspace(0, np.pi, 500)
yy = np.sum(c*np.cos(np.arange(9)*t.reshape(-1, 1)), axis=1)
plt.plot(x, y, 'ro')
plt.plot(t, yy)
plt.show()

## 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:

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.

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()]


## Recursive randomness of reals: summing a random decreasing sequence

In a comment to Recursive randomness of integers Rahul pointed out a continuous version of the same problem: pick ${x_0}$ uniformly in ${[0,1]}$, then ${x_1}$ uniformly in ${[0,x_0]}$, then ${x_2}$ uniformly in ${[0, x_1]}$, etc. What is the distribution of the sum ${S=\sum_{n=0}^\infty x_n}$?

The continuous version turns out to be easier to analyse. To begin with, it’s equivalent to picking uniformly distributed, independent ${y_n\in [0,1]}$ and letting ${x_n = y_0y_1\cdots y_n}$. Then the sum is

${\displaystyle y_0+y_0y_1 + y_0y_1y_2 + \cdots }$

which can be written as

${\displaystyle y_0(1+y_1(1+y_2(1+\cdots )))}$

So, ${S}$ is a stationary point of the random process ${X\mapsto (X+1)U}$ where ${U}$ is uniformly distributed in ${[0,1]}$. Simply put, ${S}$ and ${(S+1)U}$ have the same distribution. This yields the value of ${E[S]}$ in a much simpler way than in the previous post:

${E[S]=E[(S+1)U] = (E[S] + 1) E[U] = (E[S] + 1)/2}$

hence ${E[S]=1}$.

We also get an equation for the cumulative distribution function ${C(t) = P[S\le t]}$. Indeed,

${\displaystyle P[S\le t] = P[(S+1)U \le t] = P[S \le t/U-1]}$

The latter probability is ${\int_0^1 P[S\le t/u-1]\,du = \int_0^1 C(t/u-1)\,du}$. Conclusion: ${C(t) = \int_0^1 C(t/u-1)\,du}$. Differentiate to get an equation for the probability density function ${p(t)}$, namely ${p(t) = \int_0^1 p(t/u-1)\,du/u}$. It’s convenient to change the variable of integration to ${v = t/u-1}$, which leads to

${\displaystyle p(t) = \int_{t-1}^\infty p(v)\,\frac{dv}{v+1}}$

Another differentiation turns the integral equation into a delay differential equation,

${\displaystyle p'(t) = - \frac{p(t-1)}{t}}$

Looks pretty simple, doesn’t it? Since the density is zero for negative arguments, it is constant on ${[0,1]}$. This constant, which I’ll denote ${\gamma}$, is ${\int_{0}^\infty p(v)\,\frac{dv}{v+1}}$, or simply ${E[1/(S+1)]}$. I couldn’t get an analytic formula for ${\gamma}$. My attempt was ${E[1/(S+1)] = \sum_{n=0}^\infty (-1)^n M_n}$ where ${M_n=E[S^n]}$ are the moments of ${S}$. The moments can be computed recursively using ${E[S^n] = E[(S+1)^n]E[U^n]}$, which yields

${\displaystyle M_n=\frac{1}{n} \sum_{k=0}^{n-1} \binom{n}{k}M_k}$

The first few moments, starting with ${M_0}$, are 1, 1, 3/2, 17/6, 19/3, 81/5, 8351/80… Unfortunately the series ${\sum_{n=0}^\infty (-1)^n M_n}$ diverges, so this approach seems doomed. Numerically ${\gamma \approx 0.5614}$ which is not far from the Euler-Mascheroni constant, hence the choice of notation.

On the interval (1,2) we have ${p'(t) = -\gamma/t}$, hence

${p(t) = \gamma(1-\log t)}$ for ${1 \le t \le 2}$.

The DDE gets harder to integrate after that… on the interval ${[2,3]}$ the solution already involves the dilogarithm (Spence’s function):

${\displaystyle p(t) = \gamma(1+\pi^2/12 - \log t + \log(t-1)\log t + \mathrm{Spence}\,(t))}$

following SciPy’s convention for Spence. This is as far as I went… but here is an experimental confirmation of the formulas obtained so far (link to full size image).

To generate a sample from distribution S, I begin with a bunch of zeros and repeat “add 1, multiply by U[0,1]” many times. That’s it.

import numpy as np
import matplotlib.pyplot as plt
trials = 10000000
terms = 10000
x = np.zeros(shape=(trials,))
for _ in range(terms):
np.multiply(x+1, np.random.uniform(size=(trials,)), out=x)
_ = plt.hist(x, bins=1000, normed=True)
plt.show()

I still want to know the exact value of ${\gamma}$… after all, it’s also the probability that the sum of our random decreasing sequence is less than 1.

### Update

The constant I called “${\gamma}$” is in fact ${\exp(-\gamma)}$ where ${\gamma}$ is indeed Euler’s constant… This is what I learned from the Inverse Symbolic Calculator after solving the DDE (with initial value 1) numerically, and calculating the integral of the solution. From there, it did not take long to find that

Oh well. At least I practiced solving delay differential equations in Python. There is no built-in method in SciPy for that, and although there are some modules for DDE out there, I decided to roll my own. The logic is straightforward: solve the ODE on an interval of length 1, then build an interpolating spline out of the numeric solution and use it as the right hand side in the ODE, repeat. I used Romberg’s method for integrating the solution; the integration is done separately on each interval [k, k+1] because of the lack of smoothness at the integers.

import numpy as np
from scipy.integrate import odeint, romb
from scipy.interpolate import interp1d
numpoints = 2**12 + 1
solution = [lambda x: 1]
integrals = [1]
for k in range(1, 15):
y0 = solution[k-1](k)
t = np.linspace(k, k+1, numpoints)
rhs = lambda y, x: -solution[k-1](np.clip(x-1, k-1, k))/x
y = odeint(rhs, y0, t, atol=1e-15, rtol=1e-13).squeeze()
solution.append(interp1d(t, y, kind='cubic', assume_sorted=True))
integrals.append(romb(y, dx=1/(numpoints-1)))
total_integral = sum(integrals)
print("{:.15f}".format(1/total_integral))


As a byproduct, the program found the probabilities of the random sum being in each integer interval:

• 56.15% in [0,1]
• 34.46% in [1,2]
• 8.19% in [2,3]
• 1.1% in [3,4]
• 0.1% in [4,5]
• less than 0.01% chance of being greater than 5

## The Kolakoski-Cantor set

A 0-1 sequence can be interpreted as a point in the interval [0,1]. But this makes the long-term behavior of the sequence practically invisible due to limited resolution of our screens (and eyes). To make it visible, we can also plot the points obtained by shifting the binary sequence to the left (Bernoulli shift, which also goes by many other names). The resulting orbit  is often dense in the interval, which doesn’t really help us visualize any patterns. But sometimes we get an interesting complex structure.

The vertical axis here is the time parameter, the number of dyadic shifts. The 0-1 sequence being visualized is the Kolakoski sequence in its binary form, with 0 and 1 instead of 1 and 2. By definition, the n-th run of equal digits in this sequence has length ${x_n+1}$. In particular, 000 and 111 never occur, which contributes to the blank spots near 0 and 1.

Although the sequence is not periodic, the set is quite stable in time; it does not make a visible difference whether one plots the first 10,000 shifts, or 10,000,000. The apparent symmetry about 1/2 is related to the open problem of whether the Kolakoski sequence is mirror invariant, meaning that together with any finite word (such as 0010) it also contains its complement (that would be 1101).

There are infinitely many forbidden words apart from 000 and 111 (and the words containing those). For example, 01010 cannot occur because it has 3 consecutive runs of length 1, which implies having 000 elsewhere in the sequence. For the same reason, 001100 is forbidden. This goes on forever: 00100100 is forbidden because it implies having 10101, etc.

The number of distinct words of length n in the Kolakoski sequence is bounded by a power of n (see F. M. Dekking, What is the long range order in the Kolakoski sequence?). Hence, the set pictured above is covered by ${O(n^p)}$ intervals of length ${2^{-n}}$, which implies it (and even its closure) is zero-dimensional in any fractal sense (has Minkowski dimension 0).

The set KC apparently does not have any isolated points; this is also an open problem, of recurrence (whether every word that appears in the sequence has to appear infinitely many times). Assuming this is so, the closure of the orbit is a totally disconnected compact set without isolated points, i.e., a Cantor set. It is not self-similar (not surprising, given it’s zero-dimensional), but its relation to the Bernoulli shift implies a structure resembling self-similarity:

Applying the transformations ${x\mapsto x/2}$ and ${x\mapsto (1+x)/2}$ yields two disjoint smaller copies that cover the original set, but with some spare parts left. The leftover bits exist because not every word in the sequence can be preceded by both 0 and 1.

Applying the transformations ${x\mapsto 2x}$ and ${x\mapsto 2x-1}$ yields two larger copies that cover the original set. There are no extra parts within the interval [0,1] but there is an overlap between the two copies.

The number ${c = \inf KC\approx 0.146778684766479}$ appears several times in the structure of the set: for instance, the central gap is ${((1-c)/2, (1+c)/2)}$, the second-largest gap on the left has the left endpoint ${(1-c)/4}$, etc. The Inverse Symbolic Calculator has not found anything about this number. Its binary expansion begins with 0.001 001 011 001 001 101 001 001 101 100… which one can recognize as the smallest binary number that can be written without doing anything three times in a row. (Can’t have 000; also can’t have 001 three times in a row; and 001 010 is not allowed because it contains 01010, three runs of length 1. Hence, the number begins with 001 001 011.) This number is obviously irrational, but other than that…

In conclusion, the Python code used to plot KC.

import numpy as np
import matplotlib.pyplot as plt
n = 1000000
a = np.zeros(n, dtype=int)
j = 0
same = False
for i in range(1, n):
if same:
a[i] = a[i-1]
same = False
else:
a[i] = 1 - a[i-1]
j += 1
same = bool(a[j])
v = np.array([1/2**k for k in range(60, 0, -1)])
b = np.convolve(a, v, mode='valid')
plt.plot(b, np.arange(np.size(b)), '.', ms=2)
plt.show()

## Pisot constant beyond 0.843

In a 1946 paper Charles Pisot proved a theorem involving a curious constant ${\gamma_0= 0.843\dots}$. It can be defined as follows:

${\gamma_0= \sup\{r \colon \exists }$ monic polynomial ${p}$ such that ${|p(e^z)| \le 1}$ whenever ${|z|\le r \}}$

Equivalently, ${\gamma_0}$ is determined by the requirement that the set ${\{e^z\colon |z|\le \gamma_0\}}$ have logarithmic capacity 1; this won’t be used here. The theorem is stated below, although this post is really about the constant.

Theorem: If an entire function takes integer values at nonnegative integers and is ${O(e^{\gamma |z|})}$ for some ${\gamma < \gamma_0}$, then it is a finite linear combination of terms of the form ${z^n \alpha^z}$, where each ${\alpha }$ is an algebraic integer.

The value of ${\gamma_0}$ is best possible; thus, in some sense Pisot’s theorem completed a line of investigation that began with a 1915 theorem by Pólya which had ${\log 2}$ in place of ${\gamma_0}$, and where the conclusion was that ${f}$ is a polynomial. (Informally speaking, Pólya proved that ${2^z}$ is the “smallest” entire-function that is integer-valued on nonnegative integers.)

Although the constant ${\gamma_0}$ was mentioned in later literature (here, here, and here), no further digits of it have been stated anywhere, as far as I know. So, let it be known that the decimal expansion of ${\gamma_0}$ begins with 0.84383.

A lower bound on ${\gamma_0}$ can be obtained by constructing a monic polynomial that is bounded by 1 on the set ${E(r) = \{e^z \colon |z|\le r \}}$. Here is E(0.843):

It looks pretty round, except for that flat part on the left. In fact, E(0.82) is covered by a disk of unit radius centered at 1.3, which means that the choice ${p(z) = z-1.3}$ shows ${\gamma_0 > 0.82}$.

How to get an upper bound on ${\gamma_0}$? Turns out, it suffices to exhibit a monic polynomial ${q}$ that has all zeros in ${E(r)}$ and satisfies ${|q|>1}$ on the boundary of ${E(r)}$. The existence of such ${q}$ shows ${\gamma_0 < r}$. Indeed, suppose that ${p}$ is monic and ${|p|\le 1}$ on ${E(r)}$. Consider the function ${\displaystyle u(z) = \frac{\log|p(z)|}{\deg p} - \frac{\log|q(z)|}{\deg q}}$. By construction ${u<0}$ on the boundary of ${E(r)}$. Also, ${u}$ is subharmonic in its complement, including ${\infty}$, where the singularities of both logarithms cancel out, leaving ${u(\infty)=0}$. This contradicts the maximum principle for subharmonic functions, according to which ${u(\infty)}$ cannot exceed the maximum of ${u}$ on the boundary.

The choice of ${q(z) = z-1.42}$ works for ${r=0.89}$.

So we have ${\gamma_0}$ boxed between 0.82 and 0.89; how to get more precise bounds? I don’t know how Pisot achieved the precision of 0.843… it’s possible that he strategically picked some linear and quadratic factors, raised them to variable integer powers and optimized the latter. Today it is too tempting to throw some optimization routine on the problem and let it run for a while.

But what to optimize? The straightforward approach is to minimize the maximum of ${|p(e^z)|}$ on the circle ${|z|=r}$, approximated by sampling the function at a sufficiently fine uniform grid ${\{z_k\}}$ and picking the maximal value. This works… unspectacularly. One problem is that the objective function is non-differentiable. Another is that taking maximum throws out a lot of information: we are not using the values at other sample points to better direct the search. After running optimization for days, trying different optimization methods, tolerance options, degrees of the polynomial, and starting values, I was not happy with the results…

Turns out, the optimization is much more effective if one minimizes the variance of the set ${\{|p(\exp(z_k))|^2\}}$. Now we are minimizing a polynomial function of ${p(\exp(z_k)}$, which pushes them toward having the same absolute value — the behavior that we want the polynomial to have. It took from seconds to minutes to produce the polynomials shown below, using BFGS method as implemented in SciPy.

As the arguments for optimization function I took the real and imaginary parts of the zeros of the polynomial. The symmetry about the real axis was enforced automatically: the polynomial was the product of quadratic terms ${(z-x_k-iy_k) (z-x_k+iy_k)}$. This eliminated the potentially useful option of having real zeros of odd order, but I did not feel like special-casing those.

### Three digits

Real part: 0.916, 1.186, 1.54, 1.783
Imaginary part: 0.399, 0.572, 0.502, 0.199

Here and below, only the zeros with positive imaginary part are listed (in the left-to-right order), the others being their conjugates.

Real part: 0.878, 1.0673, 1.3626, 1.6514, 1.8277
Imaginary part: 0.3661, 0.5602, 0.6005, 0.4584, 0.171

### Four digits

Real part: 0.8398, 0.9358, 1.1231, 1.357, 1.5899, 1.776, 1.8788
Imaginary part: 0.3135, 0.4999 ,0.6163, 0.637, 0.553, 0.3751, 0.1326

Real part: 0.8397, 0.9358, 1.1231, 1.3571, 1.5901, 1.7762, 1.879
Imaginary part: 0.3136, 0.5, 0.6164, 0.6372, 0.5531, 0.3751, 0.1326

No, I didn’t post the same picture twice. The polynomials are just that similar. But as the list of zeros shows, there are tiny differences…

### Five digits

Real part: 0.81527, 0.8553, 0.96028, 1.1082, 1.28274, 1.46689, 1.63723, 1.76302, 1.82066, 1.86273
Imaginary part: 0.2686, 0.42952, 0.556, 0.63835, 0.66857, 0.63906, 0.54572, 0.39701, 0.23637, 0.08842

Real part: 0.81798, 0.85803, 0.95788, 1.09239, 1.25897, 1.44255, 1.61962, 1.76883, 1.86547, 1.89069
Imaginary part: 0.26631, 0.4234, 0.54324, 0.62676, 0.66903, 0.65366, 0.57719, 0.44358, 0.26486, 0.07896

Again, nearly the same polynomial works for upper and lower bounds. The fact that the absolute value of each of these polynomials is below 1 (for lower bounds) or greater than 1 (for upper bounds) can be ascertained by sampling them and using an upper estimate on the derivative; there is enough margin to trust computations with double precision.

Finally, the Python script I used. The function “obj” is getting minimized while function “values” returns the actual values of interest: the minimum and maximum of polynomial. The degree of polynomial is 2n, and the radius under consideration is r. The sample points are collected in array s. To begin with, the roots are chosen randomly. After minimization runs (inevitably, ending in a local minimum of which there are myriads), the new starting point is obtained by randomly perturbing the local minimum found. (The perturbation is smaller if minimization was particularly successful.)

import numpy as np
from scipy.optimize import minimize

def obj(r):
rc = np.concatenate((r[:n]+1j*r[n:], r[:n]-1j*r[n:])).reshape(-1,1)
p = np.prod(np.abs(s-rc)**2, axis=0)
return np.var(p)

def values(r):
rc = np.concatenate((r[:n]+1j*r[n:], r[:n]-1j*r[n:])).reshape(-1,1)
p = np.prod(np.abs(s-rc), axis=0)
return [np.min(p), np.max(p)]

r = 0.84384
n = 10
record = 2
s = np.exp(r * np.exp(1j*np.arange(0, np.pi, 0.01)))
xr = np.random.uniform(0.8, 1.8, size=(n,))
xi = np.random.uniform(0, 0.7, size=(n,))
x0 = np.concatenate((xr, xi))

while True:
res = minimize(obj, x0, method = 'BFGS')
if res['fun'] < record:
record = res['fun']
print(repr(res['x']))
print(values(res['x']))
x0 = res['x'] + np.random.uniform(-0.001, 0.001, size=x0.shape)
else:
x0 = res['x'] + np.random.uniform(-0.05, 0.05, size=x0.shape)