Lewis Fry Richardson (1881–1953) noticed that the length of the coast of Scotland depended on a 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.”

Related posts

• Minimalist Mandelbrot set
• The fractal nature of Brownian motion
• Randomly generated dragon

[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

One of the cards in the book is for Eugene Frank, President of B & F Instruments. His card lists his phone number as

(215) MErcury 9-7100

At first glance I thought the “E” in “MErcury” had been accidentally capitalized. Then I realized the intention was that someone would dial ME (i.e. 63) and ingore “rcury”. So the phone number would be (215) 639-7100.

This card was from 1968, the height of the space race. Presumably the card was alluding to the Project Mercury or the planet Mercury, or both.

The telephone keypad mapping (ITU E.161 standard) is a poor attempt at making phone numbers more memorable. For starters, there’s no way to encode 0 or 1 [1]. It’s unlikely a phone number will correspond to anything memorable unless you come up with the word first and then try to obtain the phone number, such as 800 FLOWERS.

Inserting extra letters, as Mr. Frank did, greatly increases the chances of encoding a phone number as a word. But then you need to denote which letters count and which ones are filler, so there’s not much advantage. Still, I wanted to play around with it for fun. I found 109 words [1] containing the letters from a telephone encoding of 4228646. (I’m using the file /usr/share/dict/words on my laptop as my list of words.)

Here are some of the more interesting hits.

• semicatholicism
• heartburning
• gladiatorism
• diabetogenic
• galactogenetic
• xanthocreatinine

There are over 30,000 words containing an encoding of the area code 832. One of these is traditional, and so I could write my phone number as

TraDitionAl semICAThOlIcisM.

Another choice for 832 is intercosmic, so

inTErCosmic GAlaCTOGeNetic

is another possibility.

Galactogentic can refer to the production of milk by the mammary glands or to the formation of galaxies (e.g. the Milky Way). Here intercosmic fits with the later sense.

I got greedy and tried to find a word containing the full phone number, 8324228646, but didn’t find anything.

Here’s my business card in the style of the Los Alamos Rolodex cards, created by Grok, using (832) GlAdiATOrIsM as the phone number.

Now suppose you remembered “gladiatorism” but not which letters were capitalized. Then you’d have to try up to 792, i.e. 12 choose 7, possible numbers, so this really isn’t a practical mnemonic. If you remembered “traditional semicatholicism” without capitalization it would be worse, with over a million possibilities (11 choose 3 times 15 choose 7). Some possibilities are counted twice, since different ways of selecting letters can lead to the same phone number, but still there are too many possibilities to try.

Related posts

• Major memory system telephone keypad
• What can you learn from a phone number?
• Phone tones inn musical notation

[1] Not only are there no letters for 0 and 1, the letters O and I represent digits. At one point in time the first digit of an exchange (the middle three digits) could not be a 0 or 1, but these digits could appear anywhere else.

[2] I initially found a list of 185 words, but some of these were duplicates: a word can represent a phone number in more than one way.

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

“king” − “man” + “woman” ≈ “queen”

This should be interpreted as saying that the word vector for king, minus the word vector for man, plus the word vector for woman, is in some sense close to the word vector for queen.

This post will look at another example. Is the expression

“coffee” + “milk” ≈ “latte”

true in some sense?

Notation

In this post I will use “foo” to mean the vector embedding of the word foo.

Coffee + milk

The cosine similarity between “coffee” + “milk” and “latte” is about 0.63. And for reasons given in the previous post, this is a large value of cosine similarity. But there are 11 words that are more similar to “milk” + “coffee” than “latte”. Here are the top 12 matches in order.

1. coffee
2. milk
3. tea
4. drink
5. chocolate
6. cream
7. breakfast
8. ice
9. beer
10. vanilla
11. starbucks
12. latte

There are two questions to resolve. First, why isn’t latte one of the closest words? Second, why is the cosine similarity large even though latte is not one of the best matches?

Concept arithmetic

When word vector arithmetic works, as in the king and queen example, the vectors combine concepts. If you replace the male gender component of king with a female component, you get a vector close to the vector for queen.

But when you add the vectors for milk and coffee, you’re not adding concepts, you’re adding ingredients.

The concepts of milk and coffee are similar in that they’re both common beverages, as are tea and even beer. A latte is a beverage, but it’s not as common as milk, coffee, tea, or beer.

Extremely uneven distribution

If you divide word vectors by their norm, you get a point on a high-dimensional sphere. In the case of the glove-twitter-200 vector embedding, you get a point on a sphere in 200 dimensions. As explained in the earlier post, a fairly large cosine similarity corresponds to a tiny portion of the sphere’s surface area.

In the example of “king” − “man” + “woman”, the vector “queen” is the closest match (except for “king” itself).

But there are a lot of words whose vectors are within a tiny region around “coffee” + “milk”. And by tiny, I mean a region that accounts for a proportion of the sphere on the order of 10⁻²³.

The glove-twitter-200 vector list contains vectors for 1.2 million words. If these vectors were roughly evenly distributed on the sphere when normalized, you’d expect each patch representing 10⁻⁶ of the sphere to contain about a word or two. You wouldn’t expect a patch taking up 10⁻¹² of the sphere to contain more than one word, and you certainly wouldn’t expect a patch taking up 10⁻²³ of the sphere to contain 12 words [1].

Rank order

Rank order based on cosine similarity is more robust than cosine similarity itself. This is an example of a phenomenon that occurs regularly: a metric whose values are dubious might still rank things well. Naive Bayes is another example. It naively computes probabilities in a way that is blatantly wrong, and yet ranking things by these spurious probabilities works well in some cases.

The cosine similarity between “king” − “man” + “woman” and “queen” is roughly the same as the cosine similarity between “coffee” + “milk” and “latte.” But in the former example, rank order picks out queen as the best match; rank order works like you’d expect, because you’re working with attributes that can be decomposed.

Dog + infant = puppy?

I wouldn’t be surprised if the Anglo-Saxon word for puppy was something like dogchild. The language was full of colorful compound words, such as hronrad (“whale-road”) for the sea and nosethyrl (“nose-hole”) for nostril.

Here are the top ten matches for “dog” + “infant” along with their cosine similarities.

1. dog, 0.819
2. infant, 0.809
3. toddler, 0.734
4. dogs, 0.697
5. puppy, 0.688
6. cat, 0.682
7. pet, 0.676
8. child, 0.671
9. newborn, 0.670
10. baby, 0.650

This shows that “puppy” is close to “dog” + “infant”, both in terms of cosine similarity and rank order, though it’s not the closet.

This also shows that you have to take the addition of word vectors with a grain of salt. It’s no surprise that puppy was a good match, but it’s surprising that cat is nearly as good.

[1] I poked around a little to get an idea just how unevenly words are distributed. The closest pair of words is jajaja and jajajaja with a cosine similarity of 0.993. The most isolated word, meaning the word whose nearest neighbor is furthest away, the the Thai word เคยไหม. It’s nearest neighbor is the Russian word боль with a cosine similarity of 0.283.

The glove-twitter-200 vectors were created from a corpus that is about half English and about other languages and strings of symbols that are not words in any language. Presumably เคยไหม would have a much closer neighbor in a corpus containing more Thai words.

I didn’t search the entire corpus, only the 50,000 most frequently occurring vectors, because a full search would require running an O(N²) search with N = 1,200,000.

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.

More Fibonacci posts

• Fibonacci meets Pythagoras
• Certified Fibonacci numbers
• Turning trig identities into Fibonacci identities

[1] K. B. Subramaniam. On a link between Triangular and Fibonacci numbers. The Mathematical Gazette, Vol. 103, No. 558 (November 2019), p. 489.

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.

For vectors x and y,

where θ is the angle between the vectors, and so

This is the cosine similarity between the words represented by x and y.

Small angles have large cosines, and so words with larger cosine similarities are closer together than words with smaller cosine similarities. The cosine similarity between a word and itself equals 1, and we’d expect unrelated words to have a cosine similarity near 0.

You can do a sort of arithmetic with vector embeddings of words. The canonical example is that

“king” − “man” + “woman” ≈ “queen”

Implicit in this equation is that we’re really adding vector representations of the words. Let a, b, c, and d be the vector embeddings of the words king, man, woman, and queen. What we’re really asserting is that

a − b + c ≈ d,

except that’s not true! Or at least it’s not true unless you view it in the right context.

The angle between a − b + c and d is about 49°, which corresponds to a cosine similarity of 0.656. Here I’m using the gensim glove-twitter-200 embedding that represents words as 200-dimensional vectors.

The way to interpret the equation above is not that a 49° degree angle is approximately 0, or that a similarity of 0.656 is approximately 1.

In high dimensions, such as 200-dimensional word embeddings, nearly all vectors are nearly perpendicular. I wrote a post about this here. So the angle between randomly selected words will usually be close to 90°, and so in that context an angle of 49° is relatively small. For example, the angle between the vector representations of king and fireplace is 89.25°.

If you divide word vectors by their norm, you can think of each vector as a point on a high-dimensional sphere, in our case a sphere in 200 dimensions. The proportion of vectors within 49° of a given point is surprisingly small in high dimensions.

Let’s say our point of interest is the north pole of an n-dimensional sphere. We’d like to calculate the proportion of the area of the sphere that is within an angle θ of the pole. I go through the calculations here. (Update: I give an approximation here that’s easier to work with than the exact formula.)

When n = 3, 17% of the area is with 49 degrees of the pole. But when n = 200, the proportion is on the order of 10⁻²⁶, essentially zero.

The vector d above representing queen is within a relatively tiny region around the vector a − b + c.

In terms of cosine similarity, 0.656 is a large similarity. Words with a cosine similarity in this range are quite close, even though we wouldn’t normally think of 0.656 being close to 1. In this context, 0.656 is close to 1.

Related posts

• Angles between words
• Area and volume of a hypersphere cap
• Cosine similarity does not satisfy the triangle inequality

Apps on a phone may automatically interpret a 10-digit number as a phone number or a 16-digit number as a package tracking number. And very often these interpretations are correct, given the kinds of things most people use their phones for.

It’s not surprising that a 10-digit number on a phone is a phone number. It’s more interesting that a 16-digit number is likely a tracking number. It could be other things, such as a credit card number. But people don’t usually write out credit card numbers in a text note; credit card numbers likely saved in some more opaque way.

I run into a variation of this problem routinely, trying to infer what a number represents inside medical notes.

A five-digit number could be a US postal code, or it could be a medical procedure code.

A six-digit number could be a date in MMDDYY format, or it could be a medical record number.

A ten-digit number could be a phone number, or it could be an NPI (National Provider Identifier) number.

It’s interesting that it’s possible make a good guess at what a number means inside unstructured text. Context has been lost, but not all context: you know you’re looking at medical notes. And that meager bit of context can be surprisingly useful.

Now suppose an OCR program sees some character that could be a beta or could be an eszett. It could calculate some kind of distance between between the pixel pattern of the character and the pixel patterns of beta and eszett. But that would be discarding context.

If you’re scanning a Greek document and run into a beta-like symbol, it’s very likely a beta. If you’re scanning a German document and run into a beta-like symbol, it could be a beta. For example, it could be a scientific paper that mentions beta particles or beta carotene. But most likely the symbol is an eszett.

The previous paragraph is saying you should compute the conditional probability of a set of pixels representing a character given the language of the document. You could be more sophisticated and look at the position of the symbol in a word as well. For example, if you see a symbol at the end of a Greek word that could either be ο (omicron) or σ (sigma), it’s likely an omicron because Greek has a different symbol ς for final sigma.

This post is a follow-on to my earlier post on the error rate in Google’s Ngram database. OCR errors are fairly common in that database, so why don’t they “just” fix the errors by using some sort of Bayesian method? OCR software probably does use some sort of Bayesian method, but it’s not that simple.

In that post I looked at the use of the word grok in English. The Ngram database shows the word being used before it was coined in 1961 due to OCR errors. Why didn’t Google compute the probability of a word being “grok” conditional on the publication date? That would be circular. We happen to know exactly when grok was coined, but in general we might try to determine when a word was coined by looking at a large set of scanned books, like the Ngram database!

Now we could compute the probable value of an ambiguously scanned word by conditioning on the language of the surrounding text. That would be a reasonable thing to do in general, but it could also lead to exactly the kind of errors we see in the Ngram data for grok.

Suppose you see an ambiguously scanned word in a book written in English. There is a higher prior probability that the word is an English word than a German word. Now suppose you see “gro?” where ? could be β, ß, or k. Without any context, perhaps the probability of the symbol being a k is small. But grok is an English word and groß is a German word which may lead you to conclude “?” is a k and the ambiguous word is grok.

Assigning higher prior probability to English words in English texts is the best thing to do on average, but in particular instances it will lead to errors. That’s life.

The Ngram database includes millions of scanned books. Google had to use OCR algorithms that work well on average. A linguist with a special interest in a particular word can be more careful and create a more sophisticated probability model (explicit or implicit) customized for their interests. Google did what they could operating at such a large scale.

Related posts

• Uppercase eszett
• A Bayesian view of Amazon resellers

The proof was simplified over the years, and formalized in Coq in 2005. Everyone is satisfied that the theorem is true, but there has never been a satisfying proof, one that a human could read and say “I see now why any map can be colored using only four colors.” And there may never be one, but see this post for a contrary prediction.

The IBM mainframe that ran the calculations completing the proof of the four color theorem did not generate the proof. It simply executed the FORTRAN program that Haken and Appel (and Koch [1]) gave it.

I don’t see the recent proof of finite-time blowup for solutions to the Navier-Stokes equations as entirely different. Computers did higher-level tasks for the OpenAI team than the mainframe did for Haken and Appel, and these tasks were not as directly programmed as the tasks that were given to the mainframe, but still machines do what they are told to do.

Related posts

• The seven color map theorem
• Detecting errors with the four color theorem

[1] John A. Koch was a programmer who worked on the four color proof with Haken and Appel. I don’t know how much credit he deserves, but I suspect it may be more than he was given.

By all accounts, the best programmers have had the biggest productivity boost from AI. And top tier mathematicians are using AI to settle long-standing mathematical conjectures. AI is a powerful tool, but tools don’t come to life and make things on their own.

I routinely have naive amateurs [1] send me proofs of open conjectures, and naturally more recent such proofs involve AI. I’ll get an email saying something like “I’ve solved the Collatz conjecture using ChatGPT, but I’m not a mathematician so I need some help verifying the proof.” And of course the supposed proof is rubbish.

The recent Navier-Stokes proof is impressive, but AI didn’t initiate the proof any more than LaTeX did. Nor did a child steer AI into proving the conjecture. Professional mathematicians were able to use AI to pursue their ideas at superhuman speed. But someone without an understanding of the Navier-Stokes problem, and familiarity with recent ideas for approaching the problem, could not have directed AI to produce a proof.

Computer scientists have been saying “garbage in, garbage out” from the beginning. A variation on this aphorism for the age of AI would be “mediocrity in, mediocrity out.”

***

[1] Amateurs can and do make contributions to mathematics. For example, in 2022 David Smith, a retired print technician, discovered a single shape that can be used to create an aperiodic tiling of the plane. By “naive amateurs” I mean people who literally do not know what they are talking about.

Quite a few other mathematical conjectures have been settled recently using AI, and these have also been accompanied with formal proofs, using Lean 4 in particular.

Until very recently, generating machine-verifiable formal proofs has been excruciatingly tedious. In 2005, Henk Barendregt and Freek Wiedijk wrote

To give an indication of how much work is needed for formalisation, we estimate that it takes approximately one work-week (five work-days of eight work-hours) to formalise one page from an undergraduate mathematics textbook.

That was the rule of thumb: forty hours per page. And this in the context of undergraduate textbooks. Research publications are much denser than textbooks. Furthermore, page 100 of a textbook probably depends mostly on material on pages 1 through 99. A sentence in a research article could cite anything that has been published before.

Say a research article takes 20 times more effort to formalize than page in an undergraduate textbook. Then formalizing the 166-page paper from OpenAI would take 132,800 person-hours. It took OpenAI 17 hours to verify their proof in Lean. I hesitate to use the word “revolutionary,” but lowering the cost of anything by four orders of magnitude is revolutionary.

I’ve used AI to generate formal proofs to check my work just for a little blog post. I wouldn’t dream of doing that if I had to pay someone a week’s salary to check my work.

Formal verification doesn’t just apply to mathematics. You could, for example, formally verify that a set of security policies are consistent and that, given certain assumptions, they accomplish their purpose. You could formally verify that a smart contract imposes a certain maximum liability. You could verify the correctness of mission-critical algorithms. These problems are easier than formalizing mathematics research, and it is easier to quantify the return on investment.

Related posts

• Automation and validation
• Formal methods let you explore the corners
• When are formal methods worth the effort?

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.

Related posts

• Engineering a waterpark
• Euler-Lagrange equations
• Data privacy

[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.

grok (v.)

“understand empathically,” 1961, an arbitrary formation by U.S. science fiction writer Robert A. Heinlein (1907-1988) in his book “Stranger in a Strange Land.” In the book it is a transliteration of a Martian word and is said to mean etymologically “to drink.” It attained popular use in 1960s-70s counterculture but is perhaps obsolete now except in internet technology circles.

I don’t believe anything in the statement above is disputed. And yet Google’s Ngram Viewer tells a very different story.

The plot implies that use of the word grok had been increasing before Heinlein’s book came out and is now much more common than it was in the 1970s. Note that the plot ends before the Grok AI came out in late 2023.

Apparently the Ngram data is unreliable, mainly for two reasons: OCR errors and inaccurate date attribution. Presumably the blip around 1900 was due to the former, OCR causing words like crok or grog to be cataloged as grok. And presumably the rise in usage before 1961 was due to the latter, misattributing the date of sources published after 1961.

The supposed rise in usage before 1961 is interesting. You’d expect some lag between the time a word circulates in conversation and when it appears in books, but apparently this lag can be smaller than the effect of date misattribution.

Etymonline speculates that grok is “perhaps obsolete now except in internet technology circles.” That matches my experience. Even in technological circles, the word was uncommon before Grok was released. Maybe it was more common in print than in conversation.

Related posts

Previous posts with Ngram stats. The effects are so large that they’re probably directionally correct after adjusting for a substantial error rate.

• Duodecimal vs. Hexadecimal
• Orwellian vs. Huxleyian

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

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.

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

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

Then by the Cauchy-Schwarz inequality we have

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.

Therefore

More generally, trace has the cyclic property

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

Then

but

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

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

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.

129303 is the Unicode code point for the Hugging Face emoj (U+1F917), which you can verify with the following Python code.

>>> import unicodedata
>>> 129303 == 0x1F917
True
>>> unicodedata.name(chr(0x1F917))
'HUGGING FACE'

Related posts

• Prevent characters from displaying as emoji
• Unicode, Tolkien, and Privacy
• Unicode surrogates
• Making flags in Unicode

RSA numbers are challenge problems posed to gauge the security of RSA encryption, which rests on the difficulty of factoring large numbers [2]. The naming scheme is confusing because RSA-n might have n digits or n bits. For example, RSA-768 is smaller than RSA-260 because the former has 768 bits and the latter has 260 digits.

RSA-260 is the largest RSA number factored so far. What does the news of its factorization say about the security of RSA?

Based on equations here, an RSA key with 862 bits would have a security level of 74 bits, i.e. the same security level as symmetric encryption with a 74-bit key. The minimum recommended RSA key size now is 2048 bits, which has a security level of 107 bits.

Security levels are on a logarithmic scale: each additional bit of security doubles the effort required to break the encryption by brute force. So breaking a 2048-bit RSA key would take 2³⁴, roughly 10¹⁰, times more effort than factoring RSA-260. All this depends on numerous assumptions, such as the state of factorization algorithms and the non-existence of CRQC [3].

Update: A week after the initial announcement, Eric Lu released an article explaining how he factored RSA-260. He says “In total, I estimate that this factorization cost about 4,900 GPU-days, or 13.5 GPU-years, which is about $400k at current market prices.”

Related posts

• RSA implementation flaws
• Generating and inspecting an RSA key
• Martin Gardner’s RSA article
• RSA munitions T-shirt

[1] N = pq = 22112825529529666435281085255026230927612089502470015394413748319128822941402001986512729726569746599085900330031400051170742204560859276357953757185954298838958709229238491006703034124620545784566413664540684214361293017694020846391065875914794251435144458199

p = 4397328654844826923795068102505872571721883526553349659561256924505973939597593482272505698004801207988043088656411102133523080581

q = 5028695206842569864686141618253083416610081090075366674776775706538324961364412200138116378509733307971876652984898985905923678379

[2] The ability to efficiently factor large primes would break RSA. It’s possible that there’s a way to break RSA without being able to factor large numbers. More on that here.

[3] Cryptographically-relevant quantum computer. Quantum computers exist, but so far they’re cryptographically irrelevant. So far quantum computers cannot factor 21 without cheating.

Brian did most of the work on the project. My contribution was to mathematically formalize low-level operations on sheets of bits using linear algebra over a binary field. Lots of Hadamard products and outer products, if I remember correctly. When you can reduce computations to algebra, you can prove that a sequence of operations is correct, and you can find optimizations by simplifying expressions.

The patent mentions a programming language called Tartan. I suggested calling it plaid because it used matrices with mask patterns that reminded me of a plaid pattern, and Brian countered saying we should call it Tartan. I like that name better.

Figure 5B from the patent.

[1] Brian Beckman and John D. Cook. Compiler for a parallel processor. U.S. Patent 12,717,871 B2. Applicant/Assignee: GSI Technology Inc., Sunnyvale, CA.

In one video, someone had an agent monitor tech news sites every 30 minutes to notify him of a variety of developments. No doubt that’s less effort than visiting a bunch of sites every half hour, but not monitoring tech news in real time takes even less effort.

Another video mentioned having Grok Bot order a sandwich through DoorDash. If I want a sandwich, I make a sandwich.

One video showed how to manage dozens messaging services. Maybe you could just not use dozens of messaging services.

And of course there are videos on creating agents to monitor the agents that monitor your news, order your sandwiches, and manage your messages.

All the use cases I saw were ways to using technology mitigate problems caused by technology, making it easier to do things that don’t need to be done, or at least things that I don’t need to do.

Of course different people have different needs. Some people have a professional need to monitor news in real time, for example, and having an agent help with that could be big win. I suspect, however, that the use cases that you’ll see most often in YouTube videos have been made up to appeal to a wide audience rather than to scratch the author’s itch.

Productivity is deeply personal. As I said in an earlier post, the scripts I’ve found most useful are of zero interest to anyone else because they are so specific to my work. I’ve mostly automated tasks with Python and bash, not with AI.

I’m not trying to avoid using AI. As I said at the top of the post I’m looking for more ways to take advantage of it. But I don’t want to fall into the trap of doing more easily what doesn’t need to be done. Or, to put a finer point on it, I don’t want to find ways for my business to do things that we do not need to do, things that other businesses may need to do.

Related posts

• Productive productivity
• Maybe you only need it because you have it

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

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

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.

Some say this doesn’t matter, we are heading toward dark software factories with source code never inspected, and maybe we should eliminate source code altogether—“source code is the new assembly code”.

Others have a different view. Humans sometimes need to debug source code. Source code may need to be audited by humans for compliance. Code should be clear enough for a human to inspect and reason about the algorithms and behavior. Also, good code quality can make the code more legible to agents and reduce unnecessary context.

Developer teams can have many different ideas of what constitutes high-quality software and good coding style. Though there are many valid ways to write code, there is also wide consensus on general principles of code quality. For example, avoiding very large single functions or modules, avoiding code duplication, avoiding undisciplined feature creep or patchy code and avoiding unnecessarily deep class hierarchies or function call chains. Some coding style choices are testable empirically for impact on developer productivity. Furthermore, some code complexity measures can be computed objectively and programmatically.

My experiences are with GPT 5.5 (Extra High reasoning) and 5.6 (Extra High, Max and occasionally Ultra). Much of my experience is “out of the box” usage of Codex, with simple AGENTS.md file, though I am working on improving the engineering of guidance files, and it is helping. My source code is mostly Python. Unfortunately it is difficult to generalize any one set of experiences universally, since developers have different code bases, languages, models, harnesses and AGENTS.md files. One-shotting a simple computer game or website would be very different from developing a complex research code in a new domain.

At first glance, the AI-written code is not incomprehensible. It does not use odd variable names like “iiii” or “a87275,” and it does not look like it came out of an obfuscated code competition. But, in my experience, still the generated code has deficiencies:

• The agent has a tendency to write much more code than is necessary (commonly 2-3X more—see also related findings here). Though it is capable of deleting code, its primary impulse seems to be to write more code.
• You can work with the agent to shorten the code, but it takes work. The agent it seems is not fluent in finding structural simplifications and then extracting commonalities. In one session I spent 1/2 hour having the agent write a few hundred lines of code, and 4 hours to get it to shorten and simplify the code. You can imagine the kind of technical debt this would accumulate.
• It behaves as though code simplification is much more out of its reach than code generation. At times it just completely fails to do some simplification task I ask it to do.
• It has no instinct for when to break a file into multiple files for conceptual clarity, even if a file becomes over 10,000 lines long.
• It can reinvent a similar but different helper function in different code modules rather than designing a simple reusable function once.
• It can make massive function argument lists with 10-20 arguments rather than recognizing that the parameters may form a coherent concept representable as an abstraction or parameter object.
• Importantly, it can define functions based on abstractions that do not model the underlying domain well and are hard to decipher. When I called it on this, it said: “You’re right. The code is naming implementation mechanics instead of stating intent … it forces the reader to mentally execute several layers of infrastructure just to discover that it means.”
• It often invents terminology that cannot instantly be understood by the reader (the source code analogy of Don’t Make Me Think).
• It can repeat the same expression multiple times instead of defining a variable with a meaningful name to represent the quantity.
• It can hardwire unexplained “magic constants” into the code instead of defining them with meaningful names.

Indeed, when pressed, the models are sometimes capable of doing better. For example, for a hard design problem, 5.6 Ultra was capable of creating a good object design that was a good match to the problem domain, when I asked it to look hard at the problem—better than the less sophisticated models.

I would certainly expect that with more engineering of the agent guidance files, many or most of these problems would get better. However, it should not require extreme measures to get coding agents to write good code.

I have not compared other coding agents, but it would not surprise me if they had similar issues. Rightly, the coding models have been optimized for their software development utility, and this has undoubtedly succeeded in a revolutionary way.

It seems there is not yet a widely accepted code-quality benchmark playing the role for frontier coding models that SWE-bench has played for software engineering capability (though there are efforts). It’s especially interesting because many aspects of code quality are verifiable, making the problem seemingly quite amenable to treatment in post-training. I am hoping that someone can put together a good benchmark for this problem, and that the frontier labs will embrace these kinds of evaluations in model development.

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”

Using this line of thought, the term ultraspherical polynomial is inexplicable. This is an example of the arcane terminology I wrote about recently. In this post I’ll explain what it conveys.

A spherical polynomial is a polynomial that naturally falls out of solving Laplace’s equation in spherical coordinates, using separation of variables. Legendre polynomials are spherical polynomials.

Gegenbauer polynomials are so called because a man named Gegenbauer studied them, just as Legendre polynomials take their name from Legendre. Gegenbauer polynomials are also called ultraspherical polynomials. Why is that?

There are two possible reasons. I’m not sure which is the historical reason, but both are plausible and are useful mnemonics.

Ultraspherical polynomials are not extremely spherical, they’re beyond spherical in some sense. More modern terminology uses the hyper- prefix rather than ultra-, which helps a bit.

Ultraspherical polynomials are beyond spherical in two ways. Gegenbauer polynomials are a generalization of Legendre polynomials, so they’re beyond Legendre polynomials in this sense.

More importantly, Gegenbauer polynomials fall out of solving Laplace’s equation on a hypersphere, i.e. a sphere in ℝⁿ for n > 3, just as Legendre polynomials fall out of the case n = 3. It makes sense to call these polynomials hyperspherical because they fall out of solving an equation on a hypersphere. Unfortunately the classical term is ultraspherical rather than hyperspherical.

I think, but I’m not sure, that at one time higher dimensional spheres were called hyperspheres, but the the higher dimensional analog of spherical coordinates was called ultraspherical coordinates. If so, it would be understandable that the adjective modifying coordinates would be applied to the polynomials that result from solving equations in these coordinates.

Several years ago I wrote a post on stable and unstable recurrences. In that post I show that the stability of the recurrence relation for Bessel functions produces depends on which kind of Bessel function and which direction the recurrence is applied.

In the forward direction, computing higher order values from lower order values, works well for Bessel functions of the second kind Y_n but not for Bessel functions of the first kind J_n. In the reverse direction, the recurrence is stable for J_n but not for Y_n.

I didn’t explain in that post why this is. In this post I will.

Second order linear difference equations have two independent solutions, just like second order linear differential equations. For both kinds of equations, all solutions are linear combinations of the two solutions. Suppose one solution grows with n and the other decays. You may want to compute the decaying solution, but in doing so you might pick up a small component of the growing solution due to rounding error. This post illustrates this phenomena for differential equations, and this post illustrates it for difference equations.

When you look at a plot of Bessel functions in a text book, you’ll probably see a few plots of J_n(x) andY_n(x) for a few small values of n. The functions seem to behave roughly the same way, like sine and cosine. And that’s true, as functions of x.

But it’s not true for J_n(x) andY_n(x) as functions of n for fixed x. As n increases, J_n(x) decays to zero and Y_n(x) goes off to −∞.

That’s the source of numerical instability. And there will be similar instability problems for other recurrences where the ratios of the two independent solutions goes to zero or infinity as a function of n.

There are techniques for computing the solution that does not diverge, the so-called minimal solution, such as Miller’s algorithm mentioned here.

where a and b are functions of x and possibly n. This is called a three-term recurrence formula.

It’s amazing how often you can run into three-term recurrence formulas. There are theorems that give conditions for such recurrences to hold, but I haven’t reached the bottom of that rabbit hole [1].

For this post I just want to give examples.

NB: before using any of the recurrences below, see the next post for a numerical pitfall to avoid.

Bessel functions of the first and second kind:

Modified Bessel functions of the first and second kind:

Chebyshev polynomials of the first and second kind:

Hermite polynomials (physicists’ convention):

Legendre polynomials:

[1] See Bochner’s theorem for orthogonal polynomials, the Nikiforov–Uvarov method, and Infeld-Hull factorization.

When you’re looking at a probability density f(x) for the first time, it helps to ignore the normalizing constant. Concentrate on the part of the function involving x and know that the normalizing constant is whatever it has to be. For example, about half of the ink that it takes to write down a beta or chi-squared density is devoted to the normalization constant; the rest of the expression is easier to understand.

This post will do the opposite of the advice above and focus on normalization constants because this ties into the previous post on modified Bessel functions.

The von Mises probability distribution on a circle has two parameters, μ and κ, and its density function is

The normalizing constant is 2π I₀(κ). The factor of 2π is unsurprising for anything defined on a circle. The more interesting part is I₀, the modified Bessel function of order 0.

The von Mises-Fisher distribution is the generalization of the von Mises distribution to a sphere in p dimensions. The density function is

where the normalization constant C_p(κ) is

where I_{p/2 − 1} is the modified Bessel function of order p/2 − 1. The values of x and μ are in bold face because they are now vectors, points on the unit sphere.

When p = 2, we have the “sphere” in two dimensions, i.e. the circle, and the von Mises-Fisher distribution reduces to the von Mises distribution. But where did the cosine go? The inner product of x and μ is the cosine of the angle between the two vectors.

When p = 3, obviously an important special case, the von Mises-Fisher distribution is known as the Fisher distribution. In that case the normalizing constant C₃(κ) can be written without using modified Bessel functions because when ν = ½ + n for an integer n, I_ν(x) is an elementary function.

Things are given their names for reasons. Discovering those reasons may help you understand their motivation and use.

Pure math perspective

For each integer n, the modified Bessel function I_n is essentially the Bessel function J_n evaluated along the imaginary axis. Specifically,

From a certain shallow perspective, that’s the end of the story: modified Bessel functions are modified in the sense that the argument is multiplied by i. And there’s a fiddly constant term up front for no apparent reason.

But of course that’s not the end of the story or else this wouldn’t be worth an entire post.

The equation above is analogous to the relationships between circular and hyperbolic functions

These relationships are interesting because the circular and hyperbolic functions are independently meaningful. If you view these equations merely as definitions you lose their significance. Circular and hyperbolic functions were widely used before Euler discovered the connection between them.

Similarly, there’s a reason the modified Bessel functions were given a name their own. If you were led to Bessel functions and modified Bessel functions separately by different applications, you would regard the equation

as a discovery rather than just a definition. The following section explains why someone would be interested in modified Bessel functions.

Before we move on, I’d like to explain the reason for the term i⁻ⁿ term. In general

for all real ν. The reason for the exp(νπi/2) term is that it makes I_ν(x) real for all real x.

Applied math perspective

Bessel functions often arise from solving problems with radial symmetry. Solving the wave equation in cylindrical coordinates using separation of variables leads to Bessel’s differential equation

and its solutions J_n and Y_n, Bessel functions of the first and second kind.

Solving the heat equation in cylindrical coordinates with separation of variables leads to the modified Bessel equation

and its solutions I_n and K_n, the modified Bessel functions of the first and second kind.

This is the reason behind the complex analysis perspective above: the change of variables sending x to ix changes the sign of the x² term in Bessel’s equation.

Bessel functions describe radially symmetric oscillations, such as the vibrations of a drum head. Modified Bessel functions describe radially symmetric exponential growth or decay [1], such as in the heat in a cylinder.

Other modified functions

Struve functions are closely related to Bessel functions. The (modified) Struve functions also satisfy Bessel’s (modified) differential equation, but with a non-zero right hand side. The modified Struve functions are proportional to the unmodified Struve functions evaluated along the imaginary axis, with a proportionality constant that makes the modified Struve functions real for real arguments.

There’s a similar relationship between the Mathieu functions and modified Mathieu functions. The general pattern is that “modified” in the context of special functions means “evaluated at ix and multiplied by a constant to make the function real for real arguments.”

[1] The functions I_n grow exponentially and the functions K_n decay exponentially. For this reason, A&S didn’t tabulate I_n and K_n per se. Instead it tabulated e^−xI_n and e^xK_n because these functions varied less over their range.

The more widely useful a function is, the more likely it is that the function will be discovered independently multiple times. Independent discoveries lead to varying definitions and notations. For example, there are two widely used definitions of Hermite polynomials, one used in probability and another used in physics, that only differ by a scaling factor. This also explains why there are so many variations on the definitions of the Fourier transform and spherical coordinates.

Special functions were discovered and applied before they were studied systematically. As with most mathematics, practice preceded theory. In hindsight, some names and conventions were less than ideal, at least from the perspective of someone seeking to organize a theory.

Functions can have arcane names for several reasons, one being that their usefulness became apparent long ago. If your instinct is that things with strange names are no longer important, your instinct might be backward. The strange name may be an indication that something is so important that its usefulness became apparent long ago.

Sometimes special functions have bland, uninformative names because the names stuck before anybody could think of something better. Bob looks into an interesting family of functions [1], then later he finds another interesting family of functions. These become known as “Bob’s functions of the first kind” and “Bob’s functions of the second kind.” These names are quite understandable at the time, though in the future people will want to know what distinguishes the functions, other than the fact that Bob discovered them, and what the groupings have in common other than the order in which Bob found them.

I started this post intending to discuss modified Bessel functions and explain what exactly is modified about them, but my preface became its own post. “Modified” is an example of the bland terminology mentioned above. There are Bessel functions and modified Bessel functions. Without more context, the “modified” term isn’t very informative. But it does provide a clue that there’s some kind of close relationship between the modified and unmodified functions. That’ll be the topic of my next post.

[1] Math education doesn’t place much emphasis on history and motivation. You may have to do some digging to find out why Bob was interested in his functions. What else was Bob known for? Maybe they’re related.

If you wanted to be more accurate, you’d need to take into account the fact that the orbit of Mars is tilted about 1.85° relative to the Earth’s orbit. How much difference does that make?

To simplify things, let’s assume the Earth orbits the sun in a circle of radius 1 and Mars orbits the sun in a circle of radius 1.5. The distance between Earth and Mars over time would be basically sinusoidal.

How much does inclination contribute to this distance? In other words, what is the difference between the distance accounting for the inclination of Mars’ orbit and the distance if we assume the two orbits are in the same plane?

This plot gives the answer.

The effect is not large, about three orders of magnitude smaller than the main effect, but it’s interesting how erratic it is.

The plots were made with the following code.

from numpy import *

R = 1.5
T = R**1.5 # Kepler's third law

def f(t, theta):
    return sqrt(
        (cos(t) - R*cos(t/T)*cos(theta))**2 +
        (sin(t) - R*sin(t/T))**2 +
        (R*sin(theta)*cos(t/T))**2
    )

The first plot graphs f(t, θ) and the second graphs f(t, θ) − f(t, 0).