Why fitting a logistic is nearly impossible from early data

Posted on 18 September 2026 by John

Nothing grows exponentially forever. What appears to be an exponential curve often turns out to be some sort of S curve, such as a logistic curve. logistic curve with extrapolations

Suppose you’re collecting data on the left side of the curve. If there’s even a small amount of error in your data, you won’t be able to predict the asymptotic value with any accuracy. But if you have data on both sides of the inflection point, you can make a good prediction of the limiting value.

I’ve written about this before, explaining that the problem is hard, but I didn’t say why it’s hard. Here I’d like to give an idea why it’s hard.

Suppose you want to fit a logistic equation

$y(t) = \frac{L}{1 + \exp(-k(t - t_0))}$

to three distinct values of t and the corresponding values of y. There is a unique solution, but in general you cannot find a solution in closed form. However, if the values of t are evenly spaced

there is a method [1] to solve for the parameters L, k, and t₀. For this post we’re only interested in the limiting value L, and it can be found by

$L = \frac{y_1^2(y_0 + y_2) - 2y_0 y_1 y_2}{y_1^2 - y_0 y_2}$

independent of h.

To find out how small changes in the y‘s change the estimate of L, we take the partial derivatives of L with respect to the y‘s and find

$\frac{\partial L}{\partial y_0} = \frac{\partial L}{\partial y_2} = \frac{y_1^2\, h^2}{\left(y_1^2 - y_0 y_2\right)^2}$

and

$\frac{\partial L}{\partial y_1} = \frac{-2\, y_0 y_2\, h^2}{\left(y_1^2 - y_0 y_2\right)^2}$
All three derivatives have the same expression in the denominator: y₁² − y₀ y₂.

If the function y(t) were an exponential, this expression would be exactly zero [2]. The function y(t) is not exactly exponential, but it is approximately exponential when the t‘s are in the left or right tail of the logistic curve. The further out in either tail the t‘s are, the closer the expression is to zero.

So when all the t‘s come from the same side of the inflection point, y(t) is nearly exponential the partial derivatives are huge and so the fitted value of L is extremely sensitive to changes in the y‘s.

As a concrete example, set L = k = 1 and t₀ = 0. Evaluate y(t) at −2, −1.5, and −1. Then the values of y are

y₀ = 0.11920292
y₁ = 0.18242552
y₂ = 0.26894142

If you forecast L using exactly these three values you’ll get L = 1.

But if you change y₀ to 0.12374097, the forecasted value of L is infinite. Values of y₀ in the interval [0.11920292, 0.12374097] predict values of K in [1, ∞].

[1] Raymond Pearl and Lowell J. Reed. On the Rate of Growth of the Population of the United States Since 1790 and its Mathematical Representation. Proceedings of the National Academy of Sciences of the United States of America, Vol. 6, No. 6 (Jun. 15, 1920), pp. 275-288

[2] exp(x + h)² = exp(x)² exp(h)² = exp(x) exp(x + 2h)

Empirical fractal

Posted on 17 September 2026 by John

There’s a common saying in discussion of fractals that the length of a coastline depends on how small a device you use to measure it. I thought this was a hypothetical, say as applied to the steps in the construction of the Koch snowflake. But the saying has its roots in actually surveying.

Lewis Fry Richardson (1881–1953) noticed that the length of the coast of Scotland depended on the size of segments used to measure it. More specifically, he found that the length followed a power law, i.e. that there’s a linear relation between the log of the coastline length and the log of the ruler length.

Here’s a reproduction of Richardson’s plot, taken from [1].

Mandelbrot built on Richardson’s observation and defined the idea of fractal dimension.

I was under the impression that fractals were invented as mathematical novelties that researchers later found applications for. But as is often the case, the applications came first. Or at least some applications came first.

Ideally there’s always a feedback cycle where applications lead to theory and theory leads to applications. As Donald Knuth put it, “The best theory is inspired by practice. The best practice is inspired by theory.”

[1] Eoghan Bradley and Mark McCartney. Four hundred years of the fractal coastline of Scotland. The Mathematical Gazette, November 2019, Vol. 103, No. 558 (November 2019), pp. 518-521

Converting between cosine similarity and concentration ratio

Posted on 16 September 2026 by John

I’ve written three posts on cosine similarity lately. The first looked at interpreting cosine similarity. The second looked at an approximation related to the first. The third looked at how ranking according to cosine similarity works better than cosine similarity itself.

Normalized word vectors are points on a high dimensional sphere, and geometry in high dimensions is counterintuitive. See the first post in this series for an explanation.

The set of points within a given angular distance of a point on a hypersphere is called a spherical cap. The ratio of the area of this spherical cap to that of the whole sphere is called cap fraction or concentration ratio. Concentration ratio explains why a modest cosine similarity value corresponds to a tiny portion of the area of the sphere and should be interpreted as a close match.

For this post, I wanted to share a plot of concentration ratio as a function of cosine similarity.

This shows that moderate values of cosine similarity correspond to infinitesimal concentration ratios. And yet, as the third post linked at the top showed, word vectors are very unevenly distributed, and even extremely small regions of the sphere can contain multiple word vectors.

I only included cosine similarity values up to 0.8 because the function plotted above takes a nosedive for larger values, even on a logarithmic scale.

Here’s the Python code to make the plot, using the function cap_fraction from here.

s = np.linspace(0, 0.8, 500)
plt.plot(s, cap_fraction(np.acos(s), 200))
plt.yscale("log")
plt.xlabel("cosine similarity")
plt.ylabel("concentration ratio")
plt.show()

Fibonacci product

Posted on 16 September 2026 by John

The product of four consecutive Fibonacci numbers equals the product of two consecutive integers.

For example,

3 × 5 × 8 × 13 = 39 × 40.

I ran across this theorem in a note [1] that says “The product of any four consecutive Fibonacci numbers is twice a triangular number.” Since triangular numbers have the form n(n + 1)/2, twice a triangular number is the product of two consecutive integers.

The note also gives a way to find the numbers on the right hand side. We have

F_n F_n+1 F_n+2 F_n+3 = m(m + 1)

where m equals

F_n+1 F_n+2

if n is odd and

F_n F_n+3

if n is even.

In the example at the top, 3 is the 4th Fibonacci number, so n = 4. Since 4 is even, m is the product of the 4th and 7th Fibonacci numbers, i.e. m = 3 × 13 = 39.

Simple approximation for spherical cap area

Posted on 15 September 2026 by John

The previous post looked at how to interpret cosine similarity, or equivalently angles between word vectors. In a high-dimensional space, randomly chosen vectors are likely nearly perpendicular, and so relatively large angles, such as 50°, indicate very closely related words.

Another way to look at this, as explained in the previous post, is that in high dimensions, a spherical cap of angular radius θ represents a small portion of a sphere, even for moderately large θ.

The proportion of the area inside the spherical cap, given here, involves the “regularized incomplete beta function” and so it’s hard to have an intuition for the value.

For large dimension n, the approximation

n^−1/2 sin^{n − 1}(θ)

gives the proportion of the area inside the cap to within an order of magnitude. It’s easy to see that this function goes to zero quickly as n increases, provided |θ| < π/2.

If you have the cosine similarity c = cos θ rather than θ itself, the approximation becomes

n^−1/2 (1 − c²)^{(n − 1)/2}.

Python script

Let’s try it on the example from the previous post, in which n = 200 and θ = 49°.

import numpy as np
from scipy.special import betainc

# Fraction of S^{n-1} inside a spherical cap of angular radius theta
# theta is measured from the pole
# Assume 0 < theta < pi/2

def cap_fraction(theta, n):
    x = np.sin(theta) ** 2
    return 0.5 * betainc(0.5 * (n - 1), 0.5, x)

def cap_fraction_approx(theta, n):
    return n**(-0.5) * np.sin(theta)**(n-1)

theta = np.deg2rad(49)
print(cap_fraction(theta, 200)) 
print(cap_fraction_approx(theta, 200))

This prints 2.03e-26 and 3.37e-26. The order of magnitude is correct as advertised.

Navier-Stokes in the news

Posted on 8 September 2026 by John

There are rumors that a long-standing math problem, one of the Millennium Prize problems, has been solved.

The problem concerns technical properties of solutions to the Navier-Stokes equations [1], a set of equations that describe the dynamics of fluid flow. Popular accounts of the problem are often oversimplified and misleading.

Some reports will speak of the problem as “solving the Navier-Stokes equations.” The task is not to write down a closed-form solution, which can’t be done, or solve the equations numerically, which has been done for decades. The problem is to prove theoretical properties of solutions which are of little interest in practice.

There has been progress toward settling the Navier-Stokes problem. Terence Tao wrote a post on this yesterday.

What’s also interesting is the intrigue around the possible solution. A post this morning says

If I am reading this correctly, Tristan Buckmaster is alleging OAI has a resolution of Navier-Stokes … which maybe used info from Buckmaster and Levent Alpöge’s private Codex sessions.

Buckmaster asked OpenAI whether they used his private sessions and they have not responded. Update: Statement from OpenAI.

Personal note

This topic connects parts of my career spanning decades. My graduate work was in PDEs and I had some interest in the Navier-Stokes equations. Here are some notes I wrote back in the day.

Now I work more with privacy than with PDEs. The question of whether OpenAI uses private data, contradicting their stated policy, is more relevant to my current work than whether the Navier-Stokes equations have global regular solutions.

[1] I never know whether to say equation or equations. You’ll hear both. You could think of Navier-Stokes as one vector-valued PDE or three scalar-valued equations.

Proof of the rank-trace theorem

Posted on 5 September 2026 by John

The previous post discussed the motivation for and application of the rank-trace theorem. This post will give a proof.

Suppose A is a real symmetric matrix. The rank-trace inequality says

$\operatorname{rank}(A)\ge\frac{(\operatorname{tr} A)^2}{\operatorname{tr}(A^2)}$

where tr is the trace operator, the sum of the elements along the diagonal of the matrix.

Terse proof

Here’s the proof in a nutshell: diagonalize A and use the Cauchy-Schwarz inequality.

Detailed proof

Now let’s unpack that. Any real symmetric matrix A is similar to a matrix D with the eigenvalues of A along the diagonal.

$A = PDP^{-1}$

The trace of a matrix stays the same under a similarity transformation, i.e. multiplying by P on one side and its inverse on the other side. So without loss of generality we may as well assume A is diagonal.

The rank of a matrix equals the number of non-zero eigenvalues, so a vector containing the non-zero eigenvalues of A

$v = [\lambda_1, \lambda_2, \ldots, \lambda_r]$

has length r where r is the rank of A. Define w to be the vector of dimension r consisting of all 1’s.

$w = [1, 1, \ldots, 1]$

Then by the Cauchy-Schwarz inequality we have

$\operatorname{tr}(A)^2 = \langle v, w \rangle^2 \leq \langle v, v \rangle \, \langle w, w \rangle = r \operatorname{tr}(A^2)$

Cyclic trace property

Why should a matrix A and its diagonalization D have the same trace?

The trace of a matrix product AB equals the trace of the product BA. To prove this, write out matrix products and the traces, then note that the two expressions are equal.

$\begin{align*} \operatorname{tr}(AB) &= \sum_i(AB)_{ii}=\sum_i\sum_k A_{ik}B_{ki} \\ \operatorname{tr}(BA) &= \sum_j(BA)_{jj}=\sum_j\sum_k B_{jk}A_{kj} \end{align*}$

Therefore

$\operatorname{tr}(A) = \operatorname{tr}((PD)P^{-1}) = \operatorname{tr}(P^{-1}(PD)) = \operatorname{tr}(D)$

More generally, trace has the cyclic property

$\operatorname{tr}(ABC) = \operatorname{tr}(CAB) = \operatorname{tr}(BCA)$

However, not all permutations preserve the trace. For example, let

$A=\begin{pmatrix}0&1\\0&0\end{pmatrix},\quad B=\begin{pmatrix}0&0\\1&0\end{pmatrix},\quad C=\begin{pmatrix}1&0\\0&0\end{pmatrix}.$

Then

$\operatorname{tr}(ABC) = \operatorname{tr}\begin{pmatrix}1&0\\0&0\end{pmatrix} = 1$

but

$\operatorname{tr}(ACB) = \operatorname{tr}\begin{pmatrix}0&0\\0&0\end{pmatrix} = 0$

Computing a lower bound on matrix rank

Posted on 4 September 2026 by John

Suppose you want to know the rank of an n × n matrix A, the number of linearly independent rows of A, or equivalently the number of linearly independent columns. There are at least three difficulties.

Difficulties in computing rank

First of all, rank is not a continuous function of a matrix. Since rank is an integer, an arbitrarily small change in the matrix could cause a discrete change in the rank [1]. A small error in computing A could produce a matrix with a different rank.

Second, finding the rank takes O(n³) operations, which may or may not be an issue depending on context.

Third, you may not have the matrix A in an explicit form. Maybe you’re able to compute products Av for vectors v but it’s not practical to form the entire matrix A.

Rank-trace inequality

If you don’t need to know the rank of A per se, but only need to know whether it is above a certain size, a lower bound on the rank may enough.

Suppose A is a Hermitian matrix. If A is real, this means A is symmetric. If A is complex, this means A equals its conjugate transpose. Then the rank-trace inequality says

$\operatorname{rank}(A)\ge\frac{(\operatorname{tr} A)^2}{\operatorname{tr}(A^2)}$
The quantity on the right hand side is known as the stable rank of A. It’s not a rank in any algebraic sense, but it gives a lower bound on rank. And it solves the three problems listed above. See the next post for a proof of the rank-trace theorem.

Stability

First of all, trace is a continuous function of a matrix, and so stable rank is also a continuous function of a matrix, provided the denominator isn’t zero. A small change to a matrix only makes a small change to its stable rank. That’s why stable rank is called stable.

Efficiency

Second, although computing rank takes O(n³) operations, computing stable rank takes only O(n²) operations, though this isn’t immediately obvious.

The trace of A takes n operations: simply sum the elements on the diagonal of A. But how do you take the trace of A²? Squaring A takes n³ operations, and so if you had to square A to find the trace of A² the rank-trace inequality would have no efficiency advantage over finding the rank of A. But you can compute the trace of A² via

$\operatorname{tr}(A^2) = \sum_{i=1}^n \sum_{j=1}^n |a_{ij}|^2$

Formation

Now suppose you don’t have the matrix A per se but you do have a way of probing A, computing the product of vectors with A. Maybe A is too large to fit into memory, or explicitly computing the elements of A would take too long.

There are Monte Carlo algorithms for estimating the traces of A and A² that could be used together to estimate the stable rank of A.

Demonstration

The following Python code illustrates the discussion above.

import numpy as np

np.random.seed(20260904)
n = 5
B = np.random.randn(n, n)
A = B.T @ B + 1e-8 * np.eye(n)  # Gram matrix plus a tiny shift => SPD

rank_A = np.linalg.matrix_rank(A)
tr_A = np.trace(A)
tr_A2 = np.trace(A @ A) # matrix product 
sum_sq = np.sum(A * A) # element-by-element product
stable_rank = (tr_A ** 2) / tr_A2

print(f"A =\n{A}\n")
print(f"rank(A)              = {rank_A}")
print(f"tr(A)                = {tr_A:.12f}")
print(f"tr(A^2) direct       = {tr_A2:.12f}")
print(f"tr(A^2) indirect     = {sum_sq:.12f}")
print(f"stable rank          = {stable_rank:.12f}")

The code above produces the output below.

A =
[[ 1.09945682  0.4899665   0.98901845  0.66983113 -1.35006341]
 [ 0.4899665   0.98531254  0.35067791  0.89757603 -0.72037507]
 [ 0.98901845  0.35067791  4.31233926  0.94556225 -0.54819048]
 [ 0.66983113  0.89757603  0.94556225  1.3494295  -1.33840786]
 [-1.35006341 -0.72037507 -0.54819048 -1.33840786  3.54858332]]

rank(A)              = 5
tr(A)                = 11.295121449420
tr(A^2) direct       = 51.035447533673
tr(A^2) indirect     = 51.035447533673
stable rank          = 2.499826585688

[1] Topological argument: A map from a connected space (such as ℝ^n×n) onto a discrete space (such as ℤ) cannot be continuous, otherwise the inverse images of the points in the range would partition the connected space into disjoint open sets, violating the definition of a connected space.

Second solutions

Posted on 27 August 2026 by John

This post provides a couple examples to go along with two earlier posts.

The pattern we’re illustrating is families of polynomials p_n(x) that each satisfy a differential equation and a three-term recurrence. The differential equations have a second solution q_n(x) that is the larger solution with respect to x but the smaller solution with respect to n.

In both the examples below p_n(x) is a polynomial, and so bounded on the interval [−1, 1], and q_n(x) is not a polynomial, with singularities at ±1. This is analogous to the previous examples with Bessel functions J_n(x) and Q_n(x) that satisfy the same differential equation but have contrasting behavior with respect to x versus n.

Legendre polynomials

The differential equation

$(1-x^2)\,y^{\prime\prime} - 2x\,y^\prime + n(n+1)\,y = 0$

has two solutions for each n, P_n(x) and Q_n(x).

The solutions P_n(x) are the Legendre polynomials. The solutions Q_n(x) are not polynomials but involve a term log((1 + x)/(1 − x)) that blows up at 1 and −1. But for fixed x and increasing n, P_n(x) grows exponentially and Q_n(x) decays exponentially, provided |x| > 1.

Chebyshev polynomials

The differential equation

$(1-x^2)\,y^{\prime\prime} - x\,y^\prime + n^2\,y = 0$

has two solutions for each n, T_n(x) and V_n(x).

The solutions T_n(x) are the Chebyshev polynomials. The solutions V_n(x) are not polynomials but involve a term √(x² — 1) that become vertical up at 1 and −1. But for fixed x with |x| > 1 and increasing n, T_n(x) grows exponentially and V_n(x) decays exponentially.

Junk solutions

Posted on 26 August 2026 by John

When you’re interested in studying a family of functions, it can be useful to look at a differential equation that the functions solve. This is a theme I’ve written about several times, most recently here and here, but also three years ago here.

Orthogonal polynomials are mathematically elegant as well as very useful in applications [1]. Various families of orthogonal polynomials satisfy various differential equations. These equations have a polynomial and non-polynomial solutions. What use are the latter?

If the differential equation modeled something physical, then the second solution would be necessary to have a complete basis of solutions. But if the differential equation is only instrumental in studying the orthogonal polynomials, what use is a non-polynomial solution?

These non-polynomial solutions turn out to be useful. Just as “junk” DNA turned out not to be junk, these “junk” solutions are important. Junk DNA doesn’t directly code for proteins, but it regulates DNA that does code for proteins and serves other purposes. Similarly, these non-polynomial solutions carry information related to the polynomial solutions.

For example, orthogonal polynomials are used to construct numerical integration methods, such as Gaussian quadrature, and the associated non-polynomial solutions describe the error in these integration methods. Incidentally, Gaussian quadrature is based on Legendre polynomials, mentioned in the previous post. For every family of orthogonal polynomials there is a corresponding integration method. See these notes.

Another tie-in to recent posts is that these non-polynomial solutions are the minimal solution to the polynomial family’s three-term recurrence, the solution that takes extra care to compute numerically.

This post has been very high-level, alluding to ideas without going into details. I’d like to write future posts that go into more depth regarding the ideas introduced here.

[1] “Real analysts cannot do without Fourier, complex analysts cannot do without Laurent, and numerical analysts cannot do without Chebyshev [polynomials].” — Lloyd N. Trefethen”

Math