Solving hard problems

We help companies solve hard problems in mathematics, statistics, and computing. Let’s explore how we might work together.

Effective graph resistance

Posted on 21 May 2025 by John

I’ve mentioned the Moore-Penrose pseudoinverse of a matrix a few times, most recently last week. This post will give an application of the pseudoinverse: computing effective graph resistance.

Given a graph G, imagine replacing each edge with a one Ohm resistor. The effective resistance between two nodes in G is the electrical resistance between those the two nodes.

Calculating graph resistance requires inverting the graph Laplacian, but the graph Laplacian isn’t invertible. We’ll resolve this paradox shortly, but in the mean time this sounds like a job for the pseudoinverse, and indeed it is.

The graph Laplacian of a graph G is the matrix L = D − A where D is the diagonal matrix whose entries are the degrees of each node and A is the adjacency matrix. The vector of all 1’s

1 = (1, 1, 1, …, 1)

is in the null space of L but if G is connected then the null space is one dimensional. More generally, the dimension of the null space equals the number of components in G.

Calculating graph resistance requires solving the linear system

Lw = v

where v and w are orthogonal to 1. If the graph G is connected then the system has a unique solution. The graph Laplacian L is not invertable over the entire space, but it is invertable in the orthogonal complement to the null space.

The Moore-Penrose pseudoinverse L⁺ gives the graph resistance matrix. The graph resistance between two nodes a and b is given by

R_ab = (e_a − e_b)^T L⁺ (e_a − e_b)

where e_i is is the vector that has all zero entries except for a 1 in the ith slot.

For more on graph resistance see Effective graph resistance by Wendy Ellens et al.

Multiplying a matrix by its transpose

Posted on 20 May 2025 by John

An earlier post claimed that there practical advantages to partitioning a matrix, thinking of the matrix as a matrix of matrices. This post will give an example.

Let M be a square matrix and suppose we need to multiply M by its transpose M^T. We can compute this product faster than multiplying two arbitrary matrices of the same size by exploiting the fact that MM^T will be a symmetric matrix.

We start by partitioning M into four blocks

$M = \begin{bmatrix}A & B \\ C & D \end{bmatrix}$

Then

$M^\intercal = \begin{bmatrix}A^\intercal & C^\intercal \\ B^\intercal & D^\intercal \end{bmatrix}$

and

$MM^\intercal = \begin{bmatrix} AA^\intercal + BB^\intercal & AC^\intercal + BD^\intercal \\CA^\intercal + DB^\intercal & CC^\intercal + DD^\intercal \end{bmatrix}$

Now for the first clever part: We don’t have to compute both of the off-diagonal blocks because each is the transpose of the other. So we can reduce our calculation by 25% by not calculating one of the blocks, say the lower left block.

And now for the second clever part: apply the same procedure recursively. The diagonal blocks in MM^T involve a matrix times its transpose. That is, we can partition A and use the same idea to compute AA^T and do the same for BB^T, CC^T, and DD^T. The off diagonal blocks require general matrix multiplications.

The net result is that we can compute MM^T in about 2/3 the time it would take to multiply two arbitrary matrices of the same size.

Recently a group of researchers found a way to take this idea even further, partitioning a matrix into a 4 by 4 matrix of 16 blocks and doing some clever tricks. The RXTX algorithm can compute MM^T in about 36/41 the time required to multiply arbitrary matrices, a savings of about 5%. A 5% improvement may be significant if it appears in the inner loop of a heavy computation. According to the authors, “The algorithm was discovered by combining Machine Learning-based search methods with Combinatorial Optimization.”

A bit-twiddling marvel

Posted on 20 May 2025 by John

Pop quiz: what does the following code do?

bool is_leap_year_fast(uint32_t y) {
    return ((y * 1073750999) & 3221352463) <= 126976;
}

It determines whether the year y is a leap year in the Gregorian calendar, of course. :)

It’s very efficient, though I don’t image that would ever matter. But it’s very clever.

Gregorian vs Julian calendar

A year is a leap year in the Julian calendar if and only if it is a multiple of 4. But the Julian year is a bit too long on average to match the earth’s orbit around the sun. You get a much better fit if you add the complication that years divisible by 100 but not by 400 are not leap years.

Presumably no one reading this recalls 1900 not being a leap year, but some of you may need to know some day that 2100 is not a leap year either.

C code

The cryptic function above comes from the recent post A leap year check in three instructions by Falk Hüffner. The function is correct for the next 100 thousand years (more precisely, through 103499) and correct if we anachronistically extent the Gregorian calendar back to the year 0.

The following C code shows empirically that Hüffner’s code is correct, but you’ll have to read his post to understand why it is correct.

#include <stdio.h>
#include <stdbool.h>
#include <stdint.h>

bool is_leap_year_fast(uint32_t y) {
    return ((y * 1073750999) & 3221352463) <= 126976;
}

bool is_leap_year(uint32_t y) {
    if ((y % 4) != 0) return false;
    if ((y % 100) != 0) return true;
    if ((y % 400) == 0) return true;
    return false;
}

int main() {
    for (uint32_t y = 0; y <= 102499; y++)
        if (is_leap_year_fast(y) != is_leap_year(y))
            printf("Exception: %d\n", y);
    return 0;
}

Matrices of Matrices

Posted on 20 May 2025 by John

It’s often useful to partition a matrix, thinking of it as a matrix whose entries are matrices. For example, you could look at the matrix 6 × 6

$M = \begin{bmatrix} 3 & 1 & 4 & 0 & 0 & 0 \\ 2 & 7 & 1 & 0 & 0 & 0 \\ 5 & 0 & 6 & 0 & 0 & 0 \\ 0 & 0 & 0 & 2 & 4 & 8 \\ 0 & 0 & 0 & 1 & 7 & 6\\ 0 & 0 & 0 & 3 & 5 & 1 \\ \end{bmatrix}$

as a 2 × 2 matrix whose entries are 3 × 3 matrices.

$\begin{bmatrix} \begin{bmatrix} 3 & 1 & 4 \\ 2 & 7 & 1 \\ 5 & 0 & 6 \\ \end{bmatrix} & \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix} \\ & \\ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \\ \end{bmatrix} & \begin{bmatrix} 2 & 4 & 8 \\ 1 & 7 & 6\\ 3 & 5 & 1 \\ \end{bmatrix} \end{bmatrix}$

M is not a diagonal matrix of real numbers, but its block structure is a diagonal matrix of blocks. It might be possible in some application to take advantage of the fact that the off-diagonal blocks are zero matrices.

You can multiply two partitioned matrices by multiplying the blocks as if they were numbers. The only requirement is that the blocks be of compatible sizes for multiplication: the block widths on the left have to match the corresponding block heights on the right.

For example, if you have two partitioned matrices

$A = \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{bmatrix}$
and
$B = \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix}$

then their product is the partitioned matrix

$C = \begin{bmatrix} A_{11}B_{11} + A_{12}B_{21} & A_{11}B_{12} + A_{12}B_{22} \\ A_{21}B_{11} + A_{22}B_{21} & A_{21}B_{12} + A_{22}B_{22} \end{bmatrix}$

The only requirement is that the matrix products above make sense. The number of columns in A11 must equal the number of rows in B11, and the number of columns in A12 must equal the number of rows in B21. (The rest of the compatibility conditions follow because the blocks in a column of a partitioned matrix have the same number of columns, and the blocks in a row of a partitioned matrix have the same number of rows.)

You can partition matrices into any number of rows and columns of blocks; the examples above have been 2 × 2 arrays of blocks only for convenience.

It’s natural to oscillate a bit between thinking “Of course you can do that” and “Really? That’s incredible.” If you haven’t gone through at least one cycle perhaps you haven’t thought about it long enough.

What could go wrong?

Manipulation of partitioned matrices works so smoothly that you have to wonder where things could go wrong.

For the most part, if a calculation makes logical sense then it’s likely to be correct. If you add two blocks, they have to be compatible for addition. If you “divide” by a block, i.e, multiply by its inverse, then the block must have an inverse. Etc.

One subtle pitfall is that matrix multiplication does not commute in general, and so you have to watch out for implicitly assuming that AB = BA where A and B are blocks. For example, theorems about the determinant of a partitioned matrix are not as simple as you might hope. The determinant of the partitioned matrix

$\begin{vmatrix} A & B \\ C & D \end{vmatrix}$

is not simply det(AD − BC) in general, though it is true if all the blocks are square matrices of the same size and the blocks in one row or one column commute, i.e. if AB = BA, CD = DC, AC = CA, or BD = DB.

Probability of rolling a Yahtzee

Posted on 19 May 2025 by John

IJK (@iconjack) has calculated the probability of rolling a Yahtzee in no more than n rolls.

$\begin{align*} p(n) &= 1 + \frac{53}{13}\left(\frac{5}{6}\right)^{2n+1} + \frac{11\cdot 5^n}{13\cdot 2^{n+5}\cdot3^{3n+1}} \\ &\phantom{=} - \frac{5^n}{8\cdot 3^{2n-2}} - \frac{7\cdot 5^{n+1}}{2^{n+5}\cdot 3^{n-1}} \end{align*}$

The first few numerical values are:

p(1) = 0.0007716
p(2) = 0.0126316
p(3) = 0.0460286

The probability is 0.95 after 23 rolls, and 0.99 after 32 rolls.

Here’s a plot.

Trinomial Coefficients and Kings

Posted on 16 May 2025 by John

The term trinomial coefficient is used a couple different ways. The most natural use of the term is a generalization of bionomial coefficients to three variables, i.e. numbers of the form

$\frac{n!}{i! \, j! \, k!}$

where i + j + k = n. These numbers are the coefficients of y^j z^k in the expansion of (x + y + z)ⁿ.

But there’s another use of the term trinomial coefficient, the one that we are interest in for this post, and that is the coefficients of x^k in the expansion of (1 + x + 1/x)ⁿ. These numbers can be visualized in a triangle analogous to Pascal’s triangle.

In Pascal’s triangle, each entry is the sum of the two entries above it. In the trinomial triangle, each entry is the sum of the three entries above it.

If you’ve been reading this blog regularly you know I’ve been interested in chess puzzles lately. The thing that make me interested in trinomial coefficients is that they relate to moves of a king on a chessboard.

If you note the number paths a king can take to a square using the minimal number of moves, you get the trinomial triangle. The reason is simple: the number of ways to get to a square equals the number of ways to get there via the square up and to the left, plus the number of ways to get their via the square above, plus the number of ways to get their via the square up and to the right.

Riff on an integration bee integral

Posted on 16 May 2025 by John

I saw an integral posted online that came from this year’s MIT Integration Bee.

$\int_0^1 x^{2024} (1 - x^{2025})^{2025}\, dx$

My thoughts on seeing this were, in order:

It looks like a beta function.
The answer is a small number.
You can evaluate the integral using the substitution u = 1 − x²⁰²⁵.

I imagine most students’ reactions would be roughly the opposite. They’d first see the substitution. Some might think of the value being small, and none would think of beta functions.

Size estimate

The value of the integral is small because you have numbers between 0 and 1 being raised to large powers. In real applications, being able to estimate the size of integral may be more valuable than being able to compute the integral.

But on a quiz, like the Integration Bee, knowing that the result is small would be worthless, except possibly to check a result. You know from the context of the problem being on a timed quiz that there must be a trick [1].

Incidentally, if you changed the problem slightly, say by replacing one of the instances of 2025 with 2025.03, the integration problem becomes much harder, but you could still estimate that the value is small.

Beta functions

The first time I saw someone evaluate an integral by recognizing a beta function it seemed like he was pulling a rabbit out of a hat. I was completely surprised that he thought something was common knowledge that I’d never seen.

Years later I went to work in biostatistics at MDACC and worked with the beta function and beta random variables all the time. If you look through my publications, you’ll see the word “beta” appears several times.

The beta function is defined as follows. You could take the first equality as the definition and the second as a theorem.

$B(x, y) = \int_0^1 t^{x-1} (1 - t)^{y-1} \,dx = \frac{\Gamma(x) \, \Gamma(y)}{\Gamma(x+y)}$

The integrand in the definition of the beta function is the probability density function for a beta(x, y) random variable, and so B(x, y) is the normalizing constant for a beta(x, y) probability distribution.

The integral in the Integration Bee question is not a beta function. I thought maybe it could be turned into a beta function, but the u substitution is simpler.

If you change the integrand above to

$\int_0^1 x^{2024} (1 - x)^{2025}\, dx$

then you do have a beta function. The value of this integral is B(2025, 2026) which equals 2024! 2025! / 4050!, which is an extremely small number.

Numerical estimation

The integral in the previous section is roughly the same as that of x²⁰²⁴ (1 − x)²⁰²⁴. This function has its maximum at ½, and so the integrand is less than 2⁻⁴⁰⁴⁸ ≈ 10⁻¹²¹⁹.

If we want to evaluate 2024! 2025! / 4050! numerically we’ll need to be indirect about it because 2024! would overflow and the final result would underflow. We could calculate the log base 10 of the value using Python as follows.

>>> from scipy.special import gammaln
>>> from math import log
>>> (gammaln(2025) + gammaln(2026) - gammaln(4051))/log(10)
-1220.576093208249

So our upper bound got us to within an order of magnitude or two of the result.

[1] When you see the current year used in a math problem, the quiz writer is winking at you. The year is often meant to suggest that a number is large but somewhat arbitrary. Here you suspect that the key to the question would be the same if 2024 and 2025 were replaced with 1997 and 1998. Also, the year is often a hint that you don’t want to use brute force. In this case, it’s a warning not to expand the integrand using the binomial theorem.

Another little chess puzzle

Posted on 14 May 2025 by John

Here’s another little chess puzzle by Martin Gardner, taken from the same paper as earlier.

The task is to place two rooks, two bishops, and two knights on a 4 by 4 chessboard so that no piece attacks any other.

As before, there are two basic solutions, here and here, plus symmetries.

Multiplying by quaternions on the left and right

Posted on 14 May 2025 by John

The map that takes a quaternion x to the quaternion qx is linear, so it can be represented as multiplication by a matrix. The same is true of the map that takes x to xq, but the two matrices are not the same because quaternion multiplication does not commute.

Let q = a + bi + cj + dk and let _qM be the matrix that represents multiplication on the left by q. Then

$_qM = \begin{bmatrix} a & -b & -c & -d \\ b & a & -d & c \\ c & d & a & -b \\ d & -c & b & a \\ \end{bmatrix}$

Now let M_q be the matrix that represents multiplication on the right by q. Then

$M_q = \begin{bmatrix} a & -b & -c & -d \\ b & a & d & -c \\ c & -d & a & b \\ d & c & -b & a \\ \end{bmatrix}$

Can prove both matrix representations are correct by showing that they do the right thing when q = 1, i, j, and k. The rest follows by linearity.

You might speculate that the matrix representation for multiplying on the right by q might be the transpose of the matrix representation for multiplying on the left by q. You can look at the matrices above and see that’s not the case.

In this post I talk about how to represent rotations with quaterions, and in this post I give an equation for the equivalent rotation matrix for a rotation described by a quaternion. You can prove that the matrix representation is correct by multiplying out _qM and M_q* . Keep in mind that q in that case is a unit quaterion, so the sum of the squares of its components add to 1.

Embeddings, Projections, and Inverses

Posted on 14 May 2025 by John

I just revised a post from a week ago about rotations. The revision makes explicit the process of embedding a 3D vector into the quaternions, then pulling it back out.

The 3D vector is embedded in the quaternions by making it the vector part of a quaternion with zero real part:

(p₁, p₂, p₃) → (0, p₁, p₂, p₃)

and the quaternion is returned to 3D by cutting off the real part:

(p₀, p₁, p₂, p₃) → (p₁, p₂, p₃).

To give names to the the process of moving to and from quaterions, we have an embedding E : ℝ³ to ℝ⁴ and a projection P from ℝ⁴ to ℝ³.

We can represent E as a 4 × 3 matrix

$E = \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}$

and P by a 3 × 4 matrix

$P = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}$

We’d like to say E and P are inverses. Surely P undoes what E does, so they’re inverses in that sense. But E cannot undo P because you lose information projecting from ℝ⁴ to ℝ³ and E cannot recover information that was lost.

The rest of this post will look at three generalizations of inverses and how E and P relate to each.

Left and right inverse

Neither matrix is invertible, but PE equals the identity matrix on ℝ³ , and so P is a left inverse of E and E is a right inverse of P.

$PE = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}$

On the other hand, EP is not an identity matrix, and so E is not a left inverse of E, and neither is P a right inverse of E.

$EP = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix}$

Adjoint

P is the transpose of E, which means it is the adjoint of E. Adjoints are another way to generalize the idea of an inverse. More on that here.

Pseudo-inverse

The Moore-Penrose pseudo-inverse acts a lot like an inverse, which is somewhat uncanny because all matrices have a pseudo-inverse, even rectangular matrices.

Pseudo-inverses are symmetric, i.e. if A⁺ is the pseudo-inverse of A, then A is the pseudo-inverse of A⁺.

Given an m by n matrix A, the Moore-Penrose pseudoinverse A⁺ is the unique n by m matrix satisfying four conditions:

A A⁺ A = A
A⁺ A A⁺ = A⁺
(A A⁺)* = A A⁺
(A⁺ A)* = A⁺ A

To show that A⁺ = P we have to establish

EPE = E
PEP = A⁺
(EP)* = EP
(PE)* = PE

We calculated EP and PE above, and both are real and symmetric, so properties 3 and 4 hold.

We can also compute

$EPE = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} = E$

and

$PEP = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} =P$

showing that properties 1 and 2 hold as well.

Solving hard problems

Effective graph resistance

More spectral graph theory posts

Multiplying a matrix by its transpose

Related posts

A bit-twiddling marvel

Gregorian vs Julian calendar

C code

Related posts

Matrices of Matrices

What could go wrong?

More linear algebra posts

Probability of rolling a Yahtzee

Trinomial Coefficients and Kings

Related posts

Riff on an integration bee integral

Size estimate

Beta functions

Numerical estimation

Another little chess puzzle

Multiplying by quaternions on the left and right

Related posts

Embeddings, Projections, and Inverses

Left and right inverse

Adjoint

Pseudo-inverse

Related posts