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.

Hugging Face Easter Egg

Posted on 3 September 2026 by John

NVIDIA has offered to buy Hugging Face for $12,930,300,000.

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'

Hugging Face emoji

New RSA number factored

Posted on 3 September 2026 by John

Eric Lu announced on X today that he has factored RSA-260, a number N with 260 digits (862 bits) that is the product of two large primes [1].

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

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

Patented application of linear algebra

Posted on 31 August 2026 by John

I just found out Brian Beckman and I got a patent on work we did for GSI Technology [1]. Nearly all the work I do is under an NDA, so I don’t often get a chance to talk about my projects. This work is public now that it’s in a patent; I suppose it has been public since the application was published.

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 patent

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.

Making the unnecessary easier

Posted on 28 August 2026 by John

I watched a few videos this morning, looking for ideas of what I could use AI to do.

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.

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”

Ultraspherical

Posted on 26 August 2026 by John

When I hear the term ultraspherical I think of something extremely spherical. For example, a baseball is spherical, but a billiard ball is more spherical. Maybe a highly polished billiard ball is ultraspherical.

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.

Numerical (in)stability of recurrence relations

Posted on 24 August 2026 by John

The previous post gave several examples of three-term recurrence relations for special functions. These relations can be computationally useful, but they have to be applied carefully.

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.

Author: John

Proof of the rank-trace theorem

Terse proof

Detailed proof

Cyclic trace property

Computing a lower bound on matrix rank

Difficulties in computing rank

Rank-trace inequality

Stability

Efficiency

Formation

Demonstration

Hugging Face Easter Egg

Related posts

New RSA number factored

Related posts

Patented application of linear algebra

Making the unnecessary easier

Related posts

Second solutions

Legendre polynomials

Chebyshev polynomials

Junk solutions

Ultraspherical

Numerical (in)stability of recurrence relations