2048 in Python

Let's write 2048 in Python. Disclaimer: this is 2am bored-on-the-couch code and is not representative of my professional coding style ;)

First let's implement the action of the left arrowkey. This isn't the shortest or most elegant code, but it is easy to reason about. We treat each row individually, ignore the blank elements in it, and compress the remaining elements into the left. For the element-combining logic, we first decide if anything needs to be combined, and if anything does, we consider cases. There are actually very few cases that provoke combination, and these can be enumerated as: AA, AAB, ABB, AABB, AABC, ABBC, ABCC. We can handle AAAA, ABAA, etc. through these cases too.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

def left(x): news = [] score = 0 for row in x: c = [v for v in row if v] if all(a != b for a, b in zip(c, c[1:])): new = c[:] elif len(c) == 2: new = [c[0]*2]; score = c[0]*2 elif len(c) == 3: if c[0] == c[1]: new = [c[0]*2, c[2]]; score = c[0]*2 elif c[1] == c[2]: new = [c[0], c[1]*2]; score = c[1]*2 else: if c[0] == c[1] and c[2] == c[3]: new = [c[0]*2, c[1]*2]; score = c[0]*2 + c[1]*2 elif c[0] == c[1]: new = [c[0]*2, c[2], c[3]]; score = c[0]*2 elif c[1] == c[2]: new = [c[0], c[1]*2, c[3]]; score = c[1]*2 elif c[2] == c[3]: new = [c[0], c[1], c[2]*2]; score = c[2]*2 new += [None] * (4 - len(new)) news.append(new) return news, score

The remaining arrowkey actions can be implemented by their symmetry with the left arrowkey action, via rotation. The do method will carry out the arrowkey actions. We only place a new tile on the board if the action resulted in a state change.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

def cw(x): return [[x[3-j][i] for j in range(4)] for i in range(4)] def do(x, d): if d == 'l': a, s = left(x) if d == 'd': a, s = left(cw(x)) a = cw(cw(cw(a))) if d == 'r': a, s = left(cw(cw(x))) a = cw(cw(a)) if d == 'u': a, s = left(cw(cw(cw(x)))) a = cw(a) return (placenew(a) if a != x else a), s

The placenew function puts a new 2 or 4 into an empty square, following the probabilities of the standard web version of the game.

1 2 3 4 5 6

def placenew(x): es = [(r,c) for r in range(4) for c in range(4) if x[r][c] is None] r, c = random.choice(es) y = deepcopy(x) y[r][c] = 4 if random.random() < 0.1 else 2 return y

This function determines if our position can be escaped. If it can't be escaped with any actions, the game is over.

1 2	def escapable(x, dirs='udlr'): return any(do(x,d)[0] != x for d in dirs)

Let's code up a play strategy where we pick actions at random.

1 2 3 4 5 6 7 8 9 10 11 12 13

x = [ [None, None, 2, None], [None, None, None, None], [None, None, None, None], [None, 2, None, None], ] def playrandom(x): ts = 0 while escapable(x): x, s = do(x, random.choice('udlr')) ts += s return ts

How successful is this strategy? Let's play lots of random games and see the average score.

1 2 3 4 5 6 7 8 9 10 11 12 13

def mean(v): return sum(v) / len(v) def stdev(v): u = mean(v) return (sum([(k-u)**2 for k in v]) / (len(v)-1))**0.5 scores = [] while True: s = playrandom(x) scores.append(s) if len(scores) > 1: print(int(mean(scores)), int(stdev(scores)), max(scores))

This shows that the mean score at endgame with the random strategy is about 810. Can we do better? Anyone who's played the game a couple times will know that ideally you wanna keep your big tiles in one corner, say the bottom left corner. So alternating down and left actions might prove more effective. We can add in right actions if we get stuck with downs and lefts. In a true disaster we can try go up and back down.

1 2 3 4 5 6 7 8 9 10 11 12

def playwell(x): ts = 0 while escapable(x): start = ts x, s = do(x, 'd'); ts += s x, s = do(x, 'l'); ts += s if not escapable(x, 'dl'): x, s = do(x, 'r'); ts += s if not escapable(x, 'dlr'): x, s = do(x, 'u'); ts += s x, s = do(x, 'd'); ts += s return ts

This gives a mean endgame score of about 2400. Nice!

Intuition for Lagrange multipliers with multiple constraints

This is an attempt to explain the method of Lagrange multipliers without resorting to arguments about tangency of level curves, which I don't find particularly helpful and which also don't work in the case of multiple constraints.

Task: maximise $f : \mathbb R^n \to \mathbb R$ subject to a number of constraints $g_i = 0$ for each $g_i : \mathbb R^n \to \mathbb R$.

For intuition, take $n=2, i=1$ so that we are maximising $f(x,y)$ subject to $g(x,y) = 0$. The first observation is that $g(x,y)=0$ describes a curve $S$ (specifically a level curve) on the XY-plane. We are walking along this curve and looking for the maximum value of $f$. By the fact of being a level curve of $g(x,y)$, the gradient $\nabla g$ is orthogonal to $S$ at all points on $S$. This is obvious: if the gradient wasn't orthogonal, then we could nudge ourselves along $S$ in the direction of the gradient and get to a higher value of $g$ than where we were before; but the value of $g$ does not change along $S$ by construction.

Now, consider the gradient $\nabla f$ along $S$. By the reasoning we just used, if $\nabla f$ is not orthogonal to $S$, then we can nudge ourselves along $S$ in the direction of the gradient (by "in the direction of the gradient", I mean we could project $\nabla f$ onto the tangent line at our current position on $S$, giving us a vector on that tangent line, then take a step along $S$ in that direction) and reach a higher value of $f$. At the point on $S$ where $f$ is maximal, there is no direction along $S$ where we can move in order to increase $f$, and so $\nabla f$ must be orthogonal to $S$.

At this point, both $\nabla f$ and $\nabla g$ are orthogonal to $S$, meaning they are both in its normal space. $S$ is $(n-1)$-dimensional, therefore the normal space is 1-dimensional and therefore $\nabla f$ and $\nabla g$ are parallel. We can't say that $\nabla f = \nabla g$ because the gradients might have different magnitudes, or one might point in the opposite direction; but we don't care about the exact direction or the magnitudes of the gradients, we just care about the orthogonality condition. So the sought point $(x_0,y_0)$ satisfies $\nabla f(x_0,y_0) = \lambda \nabla g(x_0,y_0)$ for some constant $\lambda$, and it also satisfies $g(x_0,y_0) = 0$.

Now we will do a trick: we will artificially construct a function $L(x,y,\lambda)$ which, when we equate its gradient to 0, gives us the two conditions above. The first condition is $\lambda f - \lambda g = 0$, so one function that spits out this condition when its gradient is equated to 0 is $L(x,y, \lambda) = f(x,y) - \lambda g(x,y)$. And as a freebie, the third element of the $\nabla L$ vector, when equated to 0, gives the third condition $g(x,y) = 0$ – magic!

So the maximal point of $f(x,y)$ along $g(x,y)=0$ is one of the points where $\nabla L(x,y, \lambda) = 0$. There might be several such points, and some will be local minimums and some will be local maximums, but we can easily evaluate $f$ at each one and pick the maximum.

Multiple constraints

Now suppose we are maximising $f(x,y,z)$ subject to $g_1(x,y,z) = g_2(x,y,z) = 0$.

Things are more interesting now. For one, instead of one level curve we now have two level surfaces for $g_1,g_2$ respectively. Their gradients are still vectors, but they won't generally coincide. The maximum of $f$ doesn't necessarily lie at a point where the gradient of either $g_1$ or $g_2$ is 0.

It's difficult to reason about this setup in the abstract, so let's make this concrete. Picture the level surfaces $g_1=0, g_2=0$ as spheres in space, and imagine the two spheres intersect at the circle $g_1=g_2=0$, like two overlapping soap bubbles. Label the circle $S$.

Consider the gradients of $g_1,g_2$ at a point $p$ on $S$. The gradient of $g_1$ is a vector based at $p$ orthogonal to the surface of the first sphere, and the gradient of $g_2$ is a vector based at $p$ orthogonal to the second sphere. Both gradients are orthogonal to $S$ at $p$, and moreover, the hyperplane that passes through $p$ and intersects both vectors is orthogonal to $S$ at $p$. In visual terms, if we are an ant walking along $S$, the hyperplane looks like as a flat wall at $p$ that impedes our progress. Think hard about this before proceeding.

Since the hyperplane is orthogonal to $S$ at $p$, any vector $v$ based at $p$ lying on this hyperplane is also orthogonal to $S$ at $p$ (but note that it probably won't be orthogonal to either sphere!). The hyperplane is precisely the span of the two vectors, so $v$ is a linear combination of $\nabla g_1, \nabla g_2$.

Now what can we say about the maximality condition on $f$? If we are at a point on $S$ where the gradient of $f$ is not orthogonal to $S$, then once again we can nudge ourselves in the direction of the gradient to increase $f$. But if we are at a point $p$ where $f$ is maximal, then the gradient of $f$ will be orthogonal to $S$ at $p$. In other words, the gradient of $f$ will be some linear combination of $\nabla g_1, \nabla g_2$, since these two vectors span the hyperplane of orthogonal vectors to $S$ at $p$. Concretely, $\nabla f = \sum_i \lambda_i \nabla g_i$ for constants $\lambda_i$

Constructing our artificial function $L$ as before, but satisfying the new conditions we easily see that the right form is $L(x,y,z, \lambda_1, \lambda_2) = f(x,y,z) - \sum_i \lambda_i g_i(x,y,z)$. Equating the gradient $\nabla L$ to 0 gives us exactly the conditions we want: $\nabla f = \sum_i \lambda_i \nabla g_i$, and $g_i=0$ for all $i$.

Notes on ordinary least squares

Vector and matrix derivatives

These are definitions. They are all intuitive. (I've picked one convention for the definitions. The other convention gives all results transposed.)

Taking the derivative of a scalar by a column vector is just taking the grad:

$$\frac{\mathrm{d}a}{\mathrm{d} \bf x} = \nabla a = \left[\frac{\partial a}{\partial x_1}\ \frac{\partial a}{\partial x_2} \dots \ \frac{\partial a }{\partial x_n}\right]^\top$$

Taking the derivative of a vector by a scalar is componentwise differentiation. The result is a row vector by convention.

$$\frac{\mathrm{d}\bf{x}}{\mathrm{d}a} = \left[\frac{\mathrm{d}x_1}{\mathrm{d}a}\ \frac{\mathrm{d}x_2}{\mathrm{d}a}\ \dots\frac{\mathrm{d}x_n}{\mathrm{d}a}\right]$$

These basic derivatives are easily obtained from the definitions:

$$\frac{\mathrm{d}}{\mathrm{d}\bf x}\mathbf{x}^\top b = b$$

$$\frac{\mathrm{d}}{\mathrm{d}\bf x}\mathbf{x}^\top \mathbf x = 2\mathbf x$$

$$\frac{\mathrm{d}}{\mathrm{d}\bf x}\mathbf{x}^\top B \mathbf x = 2B\mathbf x$$

Note: taking the derivative of a row vector by a column vector gives a Jacobian matrix. The denominator of the derivative ranges across rows and the numerator ranges across columns. This is consistent with the two definitions above.

$$\frac{\mathrm{d} \mathbf y}{\mathrm{d}\bf x} = \begin{pmatrix}\frac{\mathrm d y_1}{\mathrm d x_1} & \dots & \frac{\mathrm d y_n}{\mathrm d x_1} \\ \vdots & \ddots  & \vdots\\ \frac{\mathrm d y_1}{\mathrm d x_n} & \dots & \frac{\mathrm d y_n}{\mathrm d x_n} \end{pmatrix}$$

This yields an additional useful derivative, but we won't use it:

$$\frac{\mathrm{d}}{\mathrm{d}\bf x}\mathbf{x}^\top B = B$$

Also note that $\bf x^\top y = y^\top x$ and these are both scalars. Same goes for $\mathbf{x}^\top A \mathbf{y} = \mathbf{y}^\top A^\top \mathbf{x}$. But $\bf xy^\top \neq yx^\top$; both sides are matrices and the rhs is the transpose of the lhs.

Algebra

We have a matrix of features $X$ whose first column is 1s. We have a column vector of responses $y$. We wish to find a column vector of coefficients $\beta$ (the first entry of which represents the constant term) such that

$$X\beta \approx y$$

Specifically we wish to minimise the sum of squares of elements in the vector $y-X\beta$. Fortunately, sum of squares is simply

$$\begin{align*}(y-X\beta)^\top (y-X\beta) &= y^\top y - y^\top X\beta - \beta^\top X^\top y + \beta^\top X^\top X \beta\\ &= y^\top y - 2\beta^T X^\top y + \beta^\top X^\top X \beta \end{align*}$$

and we wish to minimise this quantity over $\beta$. It is quadratic in $\beta$ so it has one critical point, and the coefficient of $\beta^T\beta$ is positive so the critical point is a minimum. Taking the derivative w.r.t. $\beta$ using the rules above and setting equal to 0 gives

$$-2X^\top y + 2X^\top X  \beta= 0$$

Rearrange to give the normal equation

$$\beta = (X^\top X)^{-1}X^\top y$$

Geometric derivation

Minimising the sum of squares $\|X \beta - y\|^2$ corresponds geometrically to choosing $\beta$ to minimise the Euclidean distance between two vectors. The $y$ vector is fixed. $\beta$ allows us to choose a linear combination of the column of $X$ as our other vector. This other vector lies in the column space of $X$ (i.e. its image when interpreted as a linear map), and this column space won't contain $y$ (if it did, that would mean $X$'s features can perfectly recreate the response for all samples).

The minimum distance between a vector and a plane is located at the orthogonal projection of the vector onto the plane. This point is $X\beta$ (for the optimal $\beta$), and so the residual vector $y-X\beta$ is orthogonal to the plane. This implies that $(y-X\beta) \cdot Xv$ is 0 for all vectors $v$ (since all such vectors lie on the plane), i.e. $(y-X\beta)^\top Xv = 0$ for all $v$. Since this holds for all $v$, then indeed $(y-X\beta)^\top X = 0$. This rearranges to give the normal equation.

$\hat y = X\beta = X(X^\top X)^{-1}X^\top y$ is the projection of $y$ onto the column space of $X$, and more generally the matrix $H = X(X^\top X)^{-1}X^\top$ projects a response vector to its corresponding prediction vector. It is known as the "hat matrix" since it "puts the hat on y", and it is a projection (i.e. idempotent) operator: $H^2 = H$.

Since $Hy=\hat y$, the diagonal entry $H_{ii}$ gives the weighting of the response $y_i$ in the prediction $\hat y_i$, compared to the weighting of other responses $y_j$. This is known as the leverage of the $i$th observation and is a number between 0 and 1 (follows from the idempotency and symmetry conditions of the hat matrix, see Wikipedia). It tells us: for this prediction on this sample, how much did the observed response of the sample play a role, compared to the responses on other samples? In a linear regression model, points at extreme $x$ values (which are not necessarily outliers!) have high leverage. If these points are outliers, it can cause problems. See https://online.stat.psu.edu/stat462/node/170/ 

Statistics

First, the variance of a random vector $x$ is defined to be the variance-covariance matrix $\text{Var}(x) = E[(x-E[x])(x-E[x])^\top]$.

Fix our dataset $X$ (i.e. it is a constant, not a random matrix). We suppose our dataset was generated by the model $y = X\beta + \epsilon$ where $y$ is a random vector, $\beta$ is an unknown parameter and $\epsilon \sim N(0,\sigma^2)$ is a random vector with parameter $\sigma$. This explains the usual regression assumptions of normal iid (implying homoscedastic) errors.

Now let $\hat \beta = (X^\top X)^{-1}X^\top y$ be an estimator of $\beta$. We would like to know its bias and its variance.

First, substitute $y$ into our expression for $\hat \beta$ to give (after expanding and simplifying) 

$$\hat \beta = \beta + (X^\top X)^{-1}X^\top \epsilon$$

The only random variable here is $\epsilon$ which has mean 0, so

$$E[\hat \beta] = \beta$$

meaning our estimator is unbiased. We now compute the variance-covariance matrix of $\hat \beta$, which will be 

$$\begin{align*}\text{Var}(\hat \beta) &=E[(\hat \beta - \beta)(\hat \beta - \beta)^\top]  \\&= E[\left((X^\top X)^{-1}X^\top \epsilon\right)\left((X^\top X)^{-1} X^\top \epsilon\right)^\top]\\&= E[(X^\top X)^{-1}X^\top \epsilon \epsilon^\top X (X^\top X)^{-1}]\\&= (X^\top X)^{-1}X^\top E[\epsilon \epsilon^\top] X (X^\top X)^{-1}\\\end{align*}$$

Now $E[\epsilon \epsilon^\top]$ is a component-wise expected value of the matrix $\epsilon \epsilon^\top$. For all off-diagonal entries we have $E[\epsilon_i \epsilon_j] = E[\epsilon_i]E[\epsilon_j] = 0$ by independence. On diagonals we have entries $E[\epsilon_i^2]$. Now $\epsilon_i \sim N(0, \sigma^2)$ therefore $\sigma^2 = \text{Var}(\epsilon_i) = E[\epsilon_i^2] - E[\epsilon_i]^2$, hence $E[\epsilon_i^2] = \sigma^2$ (well isn't that a neat little trick!). So $E[\epsilon \epsilon^\top] = \sigma^2 I_n$, a diagonal matrix which hence commutes in multiplication. [N.B. you can also simply observe that $\epsilon \epsilon^\top$ is the variance-covariance matrix of $\epsilon$, simplifying the argument above.] Continuing, 

$$\begin{align*} \text{Var}(\hat \beta) &= (X^\top X)^{-1}X^\top X \sigma^2 (X^\top X)^{-1}\\ &= \sigma^2 (X^\top X)^{-1} \end{align*}$$

What is this saying? The variance of a parameter estimate reflects our uncertainty about the value of that parameter; the parameter itself is a constant. Our estimator is a random vector and this is its covariance structure. All the covariances in the matrix are proportional to $\sigma^2$, which makes intuitive sense: the variance of the slope of the line of best fit is naturally proportional to the variance of the errors $\sigma^2$. But the covariances are also based on the inverse of an expression in $X$. The more samples contained in $X$, the larger the values of $X^\top X$, and so the smaller the inverse term in the covariance. This corresponds to the fact that the more samples we have, the more information (and hence the more certainty) we have about the parameters $\beta$. All the uncertainty in $\beta$ comes from us having to infer its value from our limited data. In the extreme case where the number of samples goes to infinity, there is no uncertainty left, and our $\beta$ can be read straight out of the data.

On the diagonal of this matrix are the variances of the regression coefficients $\hat \beta_i$. These are variances of the sampling distribution of the $\hat \beta_i$. We can take the square root of these to get the standard errors for the coefs.

We might also ask about the consistency of $\hat \beta$, i.e. whether $\hat \beta \to_p \beta$ as the dataset grows. This is the formal posing of the extreme case above. Answering the question is a bit tricky to do in the general case, and it involves reasoning formally about the distribution of $X$ which we have thus far avoided. $\hat \beta$ is indeed consistent, and the argument for this in short is that in the expression $\hat \beta = \beta + (X^\top X)^{-1}X^\top \epsilon$, the size of the $X^\top \epsilon$ can be cleverly bounded by the law of large numbers using the independence of the samples and the errors as the randomly sampled $X$ grows large. See 9.1 of http://cameron.econ.ucdavis.edu/e240a/asymptotic.pdf .

Simple OLS

Let $$X = \begin{pmatrix}1 & x_1 \\ \vdots & \vdots\\ 1 & x_n\end{pmatrix}, \quad y=(y_1\ \dots \ y_n)^\top, \quad \beta = (\beta_0 \ \beta_1)^\top$$

and write $\bar x = \frac 1 n \sum_i x_i$ and likewise for other vectors. (the $\beta$s should have hats on them but we are lazy now)

Then $$X^\top X = \begin{pmatrix}n & \sum_i x_i \\ \sum_i x_i &\sum_i x_i^2\end{pmatrix} = \begin{pmatrix}n & n\bar x \\ n \bar x & n\overline{x^2}\end{pmatrix}$$

$$(X^\top X)^{-1} = \frac{1}{n^2 \overline{x^2} - n^2 \bar x ^2} \begin{pmatrix}n\overline{x^2} & -n \bar x \\ -n \bar x & n\end{pmatrix} = \frac{1}{n(\overline{x^2} - \bar x^2)} \begin{pmatrix}\overline{x^2} & -\bar x \\ -\bar x & 1\end{pmatrix}$$

$$X^\top y = \begin{pmatrix}n\bar y \\ \sum_i x_i y_i\end{pmatrix} = \begin{pmatrix}\bar ny \\ n \overline{xy}\end{pmatrix}$$

Thus

$$\beta = \frac{1}{\overline{x^2} - \bar x^2} \begin{pmatrix}\overline{x^2}\bar y - \bar x \overline{xy} \\ -\bar x \bar y + \overline{xy} \end{pmatrix}$$

So the slope is $$\beta_1 = \frac{\overline{xy} - \bar x \bar y}{\overline{x^2} - \bar{x}^2} = \frac{\text{Cov}(x,y)}{\text{Var}(x)}$$

Recalling that $\text{Corr}(x,y) = \frac{\text{Cov}(x,y)}{\sqrt{\text{Var}(X)\text{Var}(Y)}}$, and letting $\sigma_x,\sigma_y$ be the standard deviations of $x,y$, we have

$$\beta_1 = \text{Corr}(x,y) \frac{\sigma_y}{\sigma_x}$$

which is an intuitive way of understanding the line of best fit in simple regression. The slope shows the correlation between the independent and the dependent variable, scaled appropriately by each one's standard deviation so it's in the right units. $\sigma_y, \sigma_x$ go in the numerator, denominator respectively to match the rise/run on the plot.

Notes on entropy

A measure of the uncertainty of a distribution. What does this mean? When we sample from the distribution, how certain can we be of what we'll get?

Take a discrete probability distribution over outcomes A, B, C. Suppose it's 0.9, 0.05, 0.05. Then we are pretty certain that when we sample from it we'll get A. But if the distribution is 0.33, 0.33, 0.33, then we have very little certainty.

How to quantify this? Suppose we were encoding a message whose components were sampled from the distribution. We'd want to use short encodings for common components and longer encodings for rare components. The entropy of the distribution is the expected length of a component in the most efficient possible binary prefix-free encoding. Since "expected length" is a bit hard to conceptualise (since it involves two concepts, the probability distribution and the encoding), it might be easier to think of the entropy as the inverse concept to the *density* or the *compressibility* of a distribution. Low entropy means can transmit lots of samples in few bits. High entropy means takes lots of bits to transmit samples.

Indeed it can be shown that if a component is sampled with probability p, then its optimal length is $\log_2(1/p) = -\log_2(p)$. Let's check intuition:

• if a component is sampled w/ prob 1, then we don't encode it at all, because there's no uncertainty in the distribution. We know what the message will be without even receiving it.
• if a component is sampled w/ prob 0.5, it's pretty likely; we should probably spend $\log_2(2)$ = 1 bit on it.
• if a component is very unlikely, say probability 0.01, we are happy to spend a lot of bits on it because we'll almost never have to transmit it. The number of bits is inversely proportional to the prob. But we don't want to spend 1/0.01 = 100 bits on it, that would be absurd, and misses the tree structure that we have. We do want to spend $\log_2(1/0.01) \approx 7$ bits.

This is the role played by the negative log. Note that negative log just comes from a log law applied to $\log_2(1/p)$ which is the more intuitive formula, because it captures the inverse proportionality of length to probability, as well as the $\log_2$ that comes from the binary tree structure of prefix codes. It transforms probabilities into suggested encoding lengths. 

To compute the average length of a code, we use the usual expected value formula over samples $x$ from the distribution, to give the formula for entropy of $p$:

$$H(p) = -\sum_x p(x) \log_2 p(x)$$

Cross entropy, KL divergence, mutual information

$$H(p,q) = -\sum_x q(x) \log_2 p(x)$$

What could this correspond to? Suppose we have two probability distributions over the same space. We come up with an optimal encoding for one space - but then we need to use that encoding for the other space. How many bits on average will we need for the expected component? I.e. what we are doing is randomly sampling a component from distribution $q$ and determining its length under the encoding that we came up with for distribution $p$. If it costs us more to use longer encodings, then cross-entropy is the cost of encoding a message from $Q$ under the optimal encoding scheme of $P$. If we are encoding $P$, then there is no additional cost, so $H(p,p) = H(p)$. N.B. Cross-entropy is not symmetric!

Why do we care? Cross-entropy can tell us how different two distributions are. If $p, q$ are very different, then we expect $H(p)$ to be very different to $H(p,q)$. Note that $H(p)$ and $H(q)$ might be exactly the same (e.g. because they assign the same actual probability values to their outcomes, even if they do so in a different way. p might assign 0.5, 0.25, 0.25 whereas q assigns 0.25, 0.5, 0.25), but still we'll see that the cross-entropy $H(p,q)$ is higher. Indeed, the difference $H(q,p) - H(p)$ is the Kullback-Leibler divergence, expressing the distance of a probability distribution $q$ compared to a reference distribution $p$. KL can be interpreted as the *additional cost* of encoding the distribution $Q$ under the encoding of $P$, compared to encoding the distribution of $P$ under $P$.

$$D_{KL}(P \| Q) = H(Q,P) - H(P) \\=- \sum_x P(x) \log_2(Q(x))+\sum_x P(x)\log_2(P(x)) - \\ = -\sum_x P(x) \log\frac{Q(x)}{P(x)}$$

Why do we care? As seen previously, we often want to compare two distributions in ML. For e.g. to evaluate how a multilabel classification performed on a single instance, we can do the cross-entropy of the ground truth distribution of labels for that instance (a one-hot encoding consisting of 0's and 1's) with the output of the activation function, which is a fuzzier distribution. We precisely want to quantify the distance between the network's output and the labels.

There is also a concept of *mutual information* between r.v.s $X,Y$, where $I(X;Y) = D_{KL}(P_{X,Y} \| P_X \otimes P_Y)$, the additional cost of encoding $X,Y$ as independent random variables when in fact their joint distribution is the product of marginals.

Summary: cross-entropy is a natural way to compare distributions. It measures in appropriate units, i.e. "multiplicative" units (or log of multiplicative units, really) - e.g. if one distribution gives p=1 and the other gives p=0.01, the difference in probability is huge, even if the absolute difference is only 0.99. And in fact this is penalised as log(100) = 6. If one distribution gave p=1 and the other gave p=0.001, this is penalised as -log(1000). Note that this is quite different to something like integrating the difference between the curves, which doesn't account for how much worse a p=0.0001 prediction is compared to a p=0.01 prediction when the label is 1.

Entropy of a probability distribution measures the length (or "cost", if we are paying per bit) of a message encoded in the optimal encoding for that distribution.

Cross-entropy measures the length (cost) of encoding a message from one distribution using the optimal encoding from another.

Kullback-Liebler divergence measures the *additional cost* of encoding a message from one distribution using the optimal encoding from another, compared to encoding a message from the second distribution.

Mutual information (MI) measures the information shared between two r.v.s by looking at the additional cost of encoding a message from their joint distribution as if they were independent variables, compared to encoding a message assuming that their joint distribution is the product of marginals.

Extra info: $h(x) = \log_2(1/p(x)) = -\log_2(p(x))$ is the *information* of event x, which we can think of as the "surprise" of x. If an unlikely event happens, this carries lots of information. So entropy is just the expected information of a distribution. "Expected surprise", or "expected uncertainty".

Notes on rejection sampling and importance sampling

Scenario: we know a distribution in terms of its pdf $f(x)$ and we wish to sample from it. 

Sampling from a pdf can be difficult even when we know the pdf formula. Given the closed form expression of the normal distribution pdf for example, there is no obvious way to generate samples from the distribution.

Approach 1: inverse transform sampling

From $f$ get the cdf $F$ and then invert it to $F^{-1}$. Then sample $U(0,1)$ and pipe into $F^{-1}$ to get values on the x-axis that are distributed as $f$.

Intuition: think of the cdf in your head. Pick values uniformly 0-1 on the y-axis, walk along a horizontal line until you hit the cdf, and take the x-coord.

This works well if you can compute $F^{-1}$. But this is often difficult, e.g. for the normal distribution.

Approach 2: rejection sampling

MC algorithm. Let $M$ be the maximum of $f$. Sample x uniformly in some finite interval on the x-axis. For each sample, sample $y \sim U(0,M)$. If $y < f(x)$, keep the x sample, otherwise reject. Pretty intuitive.

The obvious downside is when the desired distribution is very far from uniform, e.g. all its mass lies in a small interval. Then almost all our uniform samples will get rejected.

We can repair this by noting that instead of sampling $x$ uniformly, we can sample according to any *proposal distribution*, and we should choose one whose pdf $g$ looks similar to $f$. Choose $Q$ a scaling factor such that the subgraph of $Qg$ fully contains $f$ ($Q$ fulfils a similar role to $M$). Now repeatedly sample $x$ values from $g$, for each one sample $y \sim U(0,Qg(x))$ and keep the sample if $y < f(x)$.

Think hard about why this works for any proposed distribution -- it's not obvious. For intuition, in the uniform case it is as if we are sampling $(x,y)$ pairs uniformly from a 2D rectangle of the desired pdf, and keeping the samples that fall in the subgraph. In the non-uniform case, instead of a 2D rectangle, we are instead sampling uniformly from a wonky shape which contains the desired pdf. And sampling uniformly from the wonky shape is equivalent to sampling $x$ coords from the wonky shape itself (since short bits of the wonky shape are less likely to be uniformly sampled than tall bits), and then sampling along the vertical line from the bottom to the top of the wonky shape.

An easier way to think about this is: think about throwing darts at a rectangular board containing the pdf, and accepting samples that fall below the pdf. This certainly works, and is equivalent to sampling uniformly on the x and y. But the board does not have to be rectangular; it can be any shape as long as it contains the whole pdf, and if we sample from the shape uniformly and discard the points outside the pdf subgraph, that's equivalent to sampling the pdf subgraph uniformly. So say our wonky shape is itself a pdf (or a positive multiple of a pdf); then sampling uniformly from the wonky shape is equivalent to sampling an x coord according to the pdf of that wonky shape, and then sampling uniformly along the vertical to the top of that wonky shape.

See https://en.wikipedia.org/wiki/Rejection_sampling#Description

Importance sampling

Importance sampling is not actually a sampling method. It is a method of estimating a parameter of a distribution when it is difficult to sample from that distribution, or when we seek a lower-variance estimator than the naive mean estimator.

Suppose we have a rv $X$ distributed according to pdf $p$ and we wish to estimate $E[X]$. If we can sample from $p$, we can estimate $E[X] = \frac1n \sum_{i=1}^n x_i$ where $x_i$ are iid samples, and this is usually a pretty good estimate. But if we *cannot* sample from $p$, it is a bit trickier. One thing we can do is to estimate the integral using MC integration. Sample iid $u_1, \dots, u_n \sim U(a,b)$ where $[a,b]$ covers all (or at least most) of the support of $X$, then use

$$E[X] = \int x p(x) dx \approx \frac{b-a}{n} \sum_{i=1}^n u_i p(u_i)$$

Here we are estimating the area under a curve by sampling $x$ coords uniformly under the curve, taking the average of the heights of the curve at those $x$s to get an average height of the curve, and then multiplying that height by $b-a$ to get the area under the curve.

This is an unbiased estimator, but it doesn't have particularly low variance. Consider the case of a distribution where almost all the mass is around $[-0.01, 0.01]$. If we sample uniformly from $[-1,1]$, then most of our samples are taken in areas which don't have much impact on the estimate of the parameter, i.e. which aren't *important* to the estimation. This estimator would have very high variance. We'd be better off mostly sampling from a tighter interval around 0, by encouraging the sampling of more important values. This is the key.

In the estimate above, each $u_i$ is distributed according to a uniform pdf $q(x) = \frac{1}{b-a}$ supported on $[a,b]$. Rewrite it as

$$E_p[X] = \int xp(x) dx = \int x\frac{p(x)}{q(x)}q(x)dx \\= \frac{1}{q(x)} \int x p(x) q(x) dx = (b-a)\int xp(x)q(x) dx \\= (b-a)E_q[xp(x)] = \frac{b-a}{n} \sum_{i=1}^n u_i p(u_i)$$

This is mathematically how the MC integration works. We are evaluating the expectation of $X$ w.r.t. distribution $p$ by evaluating the expectation of $p(x)$ w.r.t. $q$, and then using the usual mean estimator for that expectation. But we could replace the uniform distribution with any other distribution! If the pdf of that distribution is $q(x)$ and its samples are $y_i$, our estimator is

$$E_p[X] = \frac1n \sum_{i=1}^n y_i \frac{p(y_i)}{q(y_i)}$$

So we can choose the distribution of $q$ to be close to the distribution of $p$. As long as $q$ is non-zero on the support of $p$, this works and is unbiased.

Intuitively, we are sampling from $p$ by instead sampling from $q$ and then *reweighting* those samples according to what we know about the distributions of $p$ and $q$. The ratio $p(x)/q(x)$ is called the *sampling ratio* and gives the importance of each sample of $q$ to the distribution $p$.

Even in cases where we can directly sample from the target distribution, importance sampling can give us a lower-variance estimate. For example to sample from $N(0,1)$ we can do $\hat \theta = \frac 1 n \sum_{i=1}^n x_i$ where the $x_i \sim N(0,1)$ are iid, satisfying $\text{var}(\hat \theta) = \frac{1}{n}$. But if we instead sampled from $N(0,\frac {1}{2})$, more of our samples would come from around the mean of $N(0,1)$, so the estimator would have lower variance. But to use importance sampling in this way, we need an a priori idea of the parameter value to know where to focus our search.

Effective sample size

We sample 30 elements independently from a population distributed with variance $\sigma^2$. Then the sample size $n$ is 30, and the estimator of the mean $\mu$ has variance $\sigma^2/n$. So the more elements we sample, the more likely it is that the estimator is close to the true mean.

Now suppose 20 of those elements were correlated to each other. Then the estimator of the mean will now have variance bigger than $\sigma^2/n$, meaning we are less confident that the estimator gives us an estimate close to the mean. Another way to think about this is that although our actual sample size was 30, the effective sample size $n_e$ will be somewhat smaller. The ESS $n_e$ is defined as the quantity such that $\sigma^2/n_e$ gives the variance of the estimator.

For instance if all elements of the sample are 100% correlated, then we only really have one datapoint from which to compute the mean, so $n_e=1$.

Rephrasing: given a sample, if all its elements are independent, then we can use the sample mean as a good estimate of the mean. But if the elements are not independent, then we lose confidence in the precision of that mean. The ESS tells us: "when it comes to calculating the variance of the mean estimator, your sample size might be $n$ but *effectively* your sample size is $n_e$".

Kish's ESS

Suppose our sample is independent but weighted. Take a sample {a,b} with weights {1,2}. Then it is as if we actually had three samples {a,b,b} where the last two are 100% correlated. If we were to take $(a+2b)/3$ as our estimator of weighted mean, the variance of this mean is $\frac 5 9 \sigma^2$, which is quite a bit smaller than the $\frac 1 3 \sigma ^2$ that we would get if we had three 1-weighted elements in the sample. So the effective sample size is not 3, in fact it's less than 2 (since we only have two distinct elements and we've unbalanced them), it is 9/5 = 1.8.

This value is given by Kish's formula for ESS, sum(weights)^2 / sum(square of each weight).

Notes on covariance, contravariance, invariance

These are words used to describe data types (“elevateds”) composed of other data types (“constituents”).

• Covariance: the subtype relationship of constituents is conserved.
• Contravariance: the subtype relationship of constituents is reversed.
• Invariance: no subtype relationship conserved.

Covariance applies when the elevated only exposes instances of the constituent. e.g. when the elevated is a read-only list of constituents, or the elevated returns a constituent. e.g. a ConstList<Cat> is a subtype of ConstList<Animal> because anything that wants to consume an Animal is happy to receive Cats.

Contravariance applies when the elevated only consumes instances of the constituent. e.g. when the elevated is a function of the constituent. e.g. an Action<Animal> is a subtype of an Action<Cat> because anything that wants to send off a Cat is happy to send it to a service receiving Animals.

Invariance applies when the elevated consumes and exposes instances of the constituent. e.g. when the elevated is a mutable list of a constituent. e.g. List<Cat> is not a subtype of List<Animal> because the user might want to take the list and store a Dog in it, and List<Animal> is not a subtype of List<Cat> because the user expects the list to contain cats.

Then function types are contravariant in the inputs and covariant in the outputs. So Cat -> Animal has non-trivial subtypes Animal -> Cat, Cat -> Cat, Animal -> Animal.

Don’t try to think of this in terms of calls to the function itself. This is just a statement about the function type. When we say Animal -> Animal is a subtype of Cat -> Animal we don’t mean that we can substitute an Animal as the input when a Cat is asked for. We mean that we can substitute the function itself with an instance of the subtype. We are saying that in the expression Animal a = f(Cat), we can replace f with a Animal -> Animal (or indeed an Animal -> Cat or Cat -> Cat) and it’ll still work.

Notes on quadratic forms, the Hessian, and positive definiteness

A quadratic form in two variables is a polynomial in two variables with all degrees 2: $f(x,y) = ax^2 + bxy + cy^2$. We can associate a two-variable quadratic form given in this form with the matrix $A = [[a, \frac12 b], [\frac12 b, c]]$ since $f(x,y) = (x,y)^\top A(x,y)$. This generalises naturally: an $n$-variable quadratic form corresponds to a unique symmetric $n\times n$ matrix.

$n$-dimensional quadratic forms are useful because they act as quadratic approximation terms in $n$ dimensions. To understand what this means, consider the Taylor expansion around $f(x)$ given by $f(x+c) =f(x) + f'(x)c + \frac12 f''(x)c^2 + \dots$ Each term here captures a property of the local behaviour of $f$. The $f(x)$  is a constant. The $f'(x)c$ term is a linear approximation without a constant term, since we've already included the constant term in the expansion. The $\frac12 f''(x)c^2$ is a quadratic approximation without a constant or linear term, since those have already been captured. And so on. The constant term captures only the constant behaviour. The linear term captures only the linear behaviour i.e. the behaviour of the slope. The quadratic term captures only the quadratic behaviour, i.e. the concavity.

How to visualise 2-variable quadratic forms? They are paraboloids, which can be elliptic (imagine a parabola rotated about its vertical, but not necessarily symmetric), hyperbolic (a saddle), or a parabolic cylinder (think of a half pipe).

We are given a function of two variables $f(x,y)$ and wish to approximate it at $f(x+a,y+b)$. The approximation is 

$$f(x+a,y+b) = f(x,y) + f_x(x,y)a + f_y(x,y)b + \tfrac12 f_{xx}(x,y)a^2 + f_{xy}(x,y)ab + \tfrac12 f_{yy}(x,y)b^2 + \dots$$

This is intuitive. We capture the constant behaviour first. Then for the linear behaviour, we capture the x and y behaviours separately and combine them linearly to describe the plane approximation. For the quadratic behaviour, we can't just capture x and y and combine them linearly, since the quadratic surface is more complex than a linear plane, and can't just be captured by combining two independent dimensions. Instead we capture the behaviour of the x and y directions separately and also of combining the x and y together.

Rewriting in vector form, we approximate $f(v+c)$ around $f(v)$ by

$$f(v+c) = f(v) + \nabla_f(v)\cdot c + \tfrac12 c^\top H_f(v)c + ...$$

where the Hessian is $H_f = [[f_{xx}, f_{xy}], [f_{yx}, f_{yy}]]$ and we evaluate this at $v$. The Hessian therefore represents the component of our approximation corresponding to "concavity information", in the same way that the grad captures "gradient/slope information". Neat.

Since partials commute for $f$ satisfying basic differentiability conditions, the Hessian is symmetric. Suppose the evaluated Hessian matrix is $[[1,-2], [-2,1]]$. What information about the concavity can we read off this matrix? $f_{xx}(v)$ and $f_{yy}(v)$ are both 1, so the double derivatives in the x and y directions are both positive. Is this enough to conclude that we are at a local minimum, and that whatever direction we head in will take us "upwards" (i.e. bending away from the tangent plane)?

Real symmetric matrices have real orthogonal eigenvectors. For our example Hessian, there is eigenvector $[1,1]$ with eigenvalue $-1$ and eigenvector $[-1,1]$ with eigenvalue $3$. Indeed we can diagonalise this Hessian using these eigenvectors/values. What it tells us is that if we head in the direction of the first eigenvector, the double derivative in this direction is -1, and if we head in the direction of the second eigenvector the double derivative is 3. So although the x and y directions both have positive double derivatives, we have found a direction $[1,1]$ where we don't have positive concavity. So we are at a saddle point, not a minimum.

Is there any other information about the concavity that's not present in its eigenvectors and eigenvalues? No! The eigenvectors and eigenvalues will fully describe the quadratic approximation at the point. Specifically the eigenvectors are the axes and the eigenvalues are the scalings. If we take level sets of this surface, we get hyperbolas with orthogonal axes, and the eigenvectors and eigenvalues give the axes positions and scalings.

In this example the Hessian at $v$ had a positive and a negative eigenvalue, so we concluded that $v$ is a saddle. Generalising: if the Hessian has two positive eigenvalues, then we are at a local minimum (i.e. the quadratic term in the Taylor series is an elliptic paraboloid), and any direction in which we move will take us upward. And if the Hessian has two negative eigenvalues, then we are at a local maximum.

And if we are working in more than two variables? Then we are at a local minimum if all the Hessian's eigenvalues are positive, and a local maximum if all the Hessian's eigenvalues are negative. But these are precisely the conditions for positive/negative definiteness!

So: positive definiteness is a condition on a matrix which means that its corresponding quadratic form is convex. Since the Hessian of a function at a point captures all the information about the function's convexity at that point, then the Hessian being positive definite at a point means the function is convex at that point. And if the Hessian is not positive or negative definite, then we can read off its eigenvectors and values to see the directions in which the function moves up or down.

And how does positive definiteness relate to quadratic optimisation? Take the function $f(X) = \tfrac 12 x^\top A x + x^\top b + c$. The second derivative of this function is $A$ everywhere. If $A$ is positive definite, then the second derivative is "positive" everywhere (i.e. from any starting point, whatever direction we go to will lead to the surface bending away from its tangent plane), and we know this from the eigenvalues all being positive. But this precisely means that $f$ is convex! This is the relationship between convex optimisation and positive semidefiniteness.

Bounded Counting Knapsack

Counting knapsack problems look like this: given $n$ items with respective integer weights $w_1,\dots,w_n$, how many ways can we form combinations of these items that have total weight $s$?

There are three basic variations of counting knapsack.

1. In the “unbounded” variation we allow combinations to contain repeated items. For instance with 3 items $a,b,c$ with weights 1,3,5 we can form a combination of $3a, 2b, 1c$ with total weight 16. The partition function $p(s)$ in number theory is given by solving this problem for $n=s$ and $w_i=i$ for all $i$.
2. In the “binary” or “0/1” variation, we can only use each item once. The weights need not be distinct, so if the $w_i$ are 1,1,3,5 respectively then we can form a combination 1,1,5 that sums to 7, but we cannot form a combination 1,3,3.
3. In the "bounded" variation, we are provided with integers $c_1, \dots, c_n$ and we are allowed to use up to $c_i$ instances of item $i$ in each combination. This is a generalisation of the two variations above, where in the unbounded variation we have $c_i = \infty$ for all $i$ and in the 0/1 variation we have $c_i = 1$ for all $i$.

The unbounded and 0/1 variations are classic combinatorics examples that can be solved in $O(ns)$ time and $O(n)$ space using basic dynamic programming. Denote $[i,k]$ to be the number of combinations of items $1, \dots, i$ with total weight $k$, so that the solution to the problem is given by $[n,s]$. Then the unbounded variation has recurrence $$[i,k] = [i-1,k] + [i,k-w_i]$$ This recurrence has the interpretation "to form a combination of weight $k$ using items $1, \dots, i$ we can either add the $i$th item to a combination of weight $k-w_i$ formed by items $1, \dots, i$, or we can just take any existing combination of weight $k$ formed by items $1, \dots, i-1$. The 0/1 variation has the very similar recurrence $$[i,k] = [i-1,k] + [i-1,k-w_i]$$ with an analogous interpretation; the only difference is that we are prevented from using item $i$ more than once in a combination. The $O(n)$ space bound is given by observing that in both cases, we can evaluate $[i, k]$ by referring only to $[j,k]$ for $j \geq i-1$, so any $[j,k]$ with $j<i-1$ can be discarded. In all variations we have base cases $[i,k] = 1$ if $i=k=0$ and $[i,k] = 0$ otherwise.

For the bounded variation things get more interesting. We might expect a solution to look like the solution to the solution for the unbounded variation, except that the unbounded solution has no "memory" for how many times it has used an item, so it has no inbuilt way of preventing us from exceeding our allowance for an item. The 0/1 variation, on the other hand, generalises naturally to the bounded case. We just replace the binary decision in its recurrence – "do we use item $i$ or not?" – by a more general decision "how many instances of item $i$ do we use?". The resulting recurrence is

$$[i,k] = \sum_{q=0}^{c_i} [i-1, k-qw_i]$$

When $c_i = 1$, this is exactly equivalent to the 0/1 variation recurrence. However, when $c_i = \infty$, this looks quite different to the recurrence for the unbounded case. Why? It is instructive to think hard about the answer to this question.

What is the time complexity of this solution? It is not $O(ns)$ because clearly the recurrence is no longer $O(1)$. Overall, for each item $i$ we sum over $c_i$ values $s$ times, so we could say it is $O(s\sum_i c_i)$. In the worst case our items have small $w_i$ and $c_i \approx s$, so $\sum_i c_i$ starts to look like $ns$ and we're looking at $O(s^2n)$.

At this point we observe that all $w_i$ are distinct: if we had $w_i = w_j$ for $i \neq j$, we could replace items $i,j$ with a new item $i+j$ with weight $w_{i+j} = w_i + w_j$ and $c_{i+j} = c_i + c_j$. For any $i,k$ we have $[i,k-w_i] = 0$ whenever $k-w_i \leq 0$, so for any fixed $w_i$ we are summing over at most $k/w_i$ terms. Using the distinctness of $w_i$ we have

$$\sum_{i=1}^n \frac{k}{w_i} \leq \sum_{i=1}^n \frac{k}{i} \leq k \log n$$

What are we saying here? If we fix a value $k$ we can evaluate $[i,k]$ across all $i$ in $O(k \log n)$. $k$ is bounded by $s$ and we must evaluate this for each $k$, so the algorithm in fact runs in $O(s^2 \log n)$. Of course, if the $c_i$ are known a priori to be small, this is a very loose bound and it's more useful to think of the algorithm as running in $O(n \sum_i c_i)$. Either way it's a far cry from the $O(ns)$ solutions for the other variations.

Generally we speed up an algorithm by noticing and getting rid of redundant computation, and the case is no different here. The critical observation is best conveyed through example. Fix $w_i = 3, c_i=4$. We have

$$\begin{aligned} [i,k] &= [i-1,k] + [i-1,k-3] + [i-1,k-6] + [i-1,k-9] + [i-1,k-12]\\ [i,k+1] &= [i-1, k+1] + [i-1, k-2] + [i-1, k-5] + [i-1, k-8] + [i-1, k-11] \\ [i,k+2] &= [i-1, k+2] + [i-1, k-1] + [i-1, k-4] + [i-1, k-7] + [i-1, k-10] \end{aligned}$$

No repeated work yet. But see what happens for $[i,k+3]$:

$$[i,k+3] = [i-1, k+3] + [i-1, k] + [i-1, k-3] + [i-1, k-6] + [i-1, k-9]$$

All except the first element of this summation have been seen before in the solution to $[i,k]$! In fact we can write

$$[i,k+3] = [i-1,k+3] + [i,k] - [i,k-12]$$

What is happening here is that to solve $[i,k]$ we are use $[i,k-w_i]$ as a cumulative array that sums values across $[i-1, \dots]$. Unfixing $w_i, c_i$ we arrive at the $O(1)$ recurrence

$$[i,k] = [i-1, k] + [i,k-w_i] - [i,k-(c_i+1)w_i]$$

To me this is a remarkable result. In some informal sense the recurrence appears to "lie between" the recurrences for the 0/1 and unbounded variations, in that it starts with $[i-1, k] + [i, k-w_i]$ like the unbounded variation, but then it corrects itself using $[i-1,k-(c_i+1)w_i]$ to avoid using item $i$ too many times. It makes sense if I squint, but I don't think I could have come up with that correction procedure through combinatorial reasoning. The recurrence seems to say "we make a combination of items up to $i$ with weight $k$ either by just taking a combination of weight $k$ using items up to $i-1$, or by adding item $i$ to a combination $k-w_i$ of items up to $i$, but if it's possible that we used too many of item $i$, uncount the number of combinations that contain too many uses of that item." This gives us a $O(ns)$ time, $O(n)$ space solution to the bounded counting knapsack problem.

In the next entry we will consider another counting knapsack variation which generalises the unbounded counting knapsack problem into another "dimension", with beautiful results.

Some reasons for studying mathematics

A different sort of post today.

I've never had much interest in doing maths. I'm fundamentally incurious about mathematical objects, I don't feel joy or excitement when I prove a theorem, and I rarely enjoy the process of solving mathematical problems. The typical questions that I encounter in my course of mathematical study don't appeal to my intellectual hunger in the same way that, say, a computer science problem can take over a corner of my brain for weeks. This is how I know that mathematics is not for me. So why did I choose to study it?

When I studied algorithms in high school I learned a bunch of graph search techniques and dynamic programming and computational geometry etc, but more importantly I developed a sort of mental toolbox for reasoning about systems and processes. Items in the toolbox could also be described as "problem-solving techniques" or "approaches" or "paradigms" or "meta-algorithms"; they can feel abstract and fuzzy and difficult to describe, but certainly seem to exist as distinct concepts in my brain so that I can manifest them at will. They're obviously useful in the context of algorithmic problem solving, but they're also applicable to problems in totally unrelated fields. There are hundreds of tools in my toolbox of various levels of abstraction and usefulness. I'll have a go at identifying some of the tools in my algorithms toolbox:

• Thinking in terms of growth rates
• The notion of "greediness"
• "Sliding window" methods: we can solve a sequence of overlapping problems by solving the first problem and efficiently transforming that solution into a solution for the subsequent problem, and so on.
• Divide and conquer
• Laziness vs eagerness
• Invariants
• "X reduces to Y", meaning in a precise sense "X is at least as hard as Y, and we can't solve X without solving Y"

This mental toolbox represents the greatest benefit I received from studying algorithmic problem solving.

The study of mathematics has given me a whole new set of tools to add to my toolbox. As in the examples above, these tools tools are used constantly in the study of maths, but they have little to do with mathematical problems themselves. These mathematical tools are often far more insightful and valuable to me than the problems that they help to solve, and I can apply them to problems outside of mathematics. On its own, this is a good enough reason for me to study maths (it's not the only reason – see below). In this entry I will list and briefly describe some of the tools I've added to my toolbox in the past few years of maths education. Most of these tools can be expressed as statements about systems or structures. I may update this list over time as new tools occur to me. All tools are stated informally.

• Action to structure: we can study actions on a structure by treating the actions as first-class objects themselves. (for instance studying symmetries via the symmetry group. Reverse of "Structure to action".)
• Basis: in some cases, the behaviour of a system is wholly described by its behaviour on a small set of representative inputs.
• Continuity: the property of being able to induce arbitrarily small changes in the output of a system via correspondingly small changes in the input.
• Contradiction (proof by): one can assert a proposition by supposing that it doesn't hold and subsequently taking a logical sequence of steps to arrive at an impossibility.
• Convexity: the property of "not having dents"; not necessarily in a geometric setting.
• Density arguments: if a property holds for a subset $A$ of a set $X$ and any element of $X$ is the limit of elements in $A$, then we may be able to extend the property to all elements of $X$. (This is useful when $A$ consists of "simple" examples where the property is easily explored. The obvious example is the treatment of simple functions in integration theory.)
• Diagonalisation: if we have a countable number of instances from some class, we can generate a new instance of that class simply by ensuring that it differs from each previous instance at exactly one point.
• Differentiation: we can learn about a system by observing its instantaneous change.
• Duality: Two seemingly disparate domains sometimes pair up conveniently, so that we can study items in one domain by looking at the corresponding item in the other domain.
• Eigen-: we can learn about a system by studying the inputs on which it is particularly well-behaved.
• Exchange arguments: we can show that a state is optimal by taking any other state and showing that it doesn't get any worse if we bring it a step closer to the proposed optimal state.
• Graphs: complex systems and structures often have natural representations as networks.
• Integration: we can learn about a system by observing its cumulative effect.
• Isomorphism: rather than speaking about the equality of two objects, it's often more meaningful to consider equivalence with respect to some structural property.
• Lifts: methods that transform solutions for a problem in one domain to solutions in a broader domain, hence "lifting" them.
• Limits: properties that hold as we get close to something may or may not hold once we get there, depending on the attributes of the surrounding space.
• Locality vs globality: properties that hold at a micro level often don't hold at the macro level and vice versa.
• Orthogonality: a stronger or more precise form of independence. If A, B are orthogonal then we can change A in any way without affecting B and vice versa.
• Pigeonhole principle: if you're placing balls into boxes with more balls than boxes, at least one box will have more than one ball. (And in reverse: if there are less balls than boxes, at least one box will be empty.)
• Probabilistic method: informally, if a property has non-zero chance of holding for a randomly-selected element of a set, then the property holds for at least one element of the set. If the mean of the values in a set is $x$ then at least one of the values is at least $x$ and at least one of the values is at most $x$.
• Projection: we can learn about a complex structure by studying what it looks from a fixed perspective. I'm stunned at how often the notion of projection appears in mathematics and how abstract it can become. We can project in the literal geometric sense – as in projective geometry, or projecting a vector onto another – but we can also project the integers onto the set {0,1} via the equivalence relation of parity, or we can project a vector space to another vector space via a linear map, or we can project a sigma-algebra to a sub-sigma-algebra via conditional expectation. Projection is closely linked to the notion of homomorphism – so much so that I wish our study of homomorphisms had begun with an explicit acknowledgement that homomorphisms are precisely the maps that behave "projectively".
• Representatives: we can sometimes study the interactions of complex structures (often equivalence classes) by instead considering the interactions of individual representatives from those complex structures. (e.g. dealing with cosets via coset representatives.)
• Structure to action: we can learn about a structure by constructing an action of that structure on some class of objects and studying the action instead. (for instance studying a group via its group action on a set. The reverse of "Action to structure".)
• Substructure: we can learn about a structure by studying the smaller structures embedded within it. (see subgroup, subspace, submodule, subgraph, etc.)
• Symmetry: we can learn about a structure by studying the ways in which we can transform the structure without changing the way it looks.
• Two-player game: some notions can be understood in terms of a two-player game oppositional game. For instance epsilon-delta proofs interpreted as a game between an epsilon-picker and a delta-picker.

Why else do I study maths if I fundamentally don't care for it? I find maths to be an exceptional form of brain training. Most of this training happens via two mechanisms:

1. Structuring and summarising information: as part of my maths education I'm expected to be able to regurgitate in exam conditions something like fifty theorems and proofs per course. Each proof can consist of pages of dense equation manipulation or tricky argumentation, which is far too much to learn by rote. Instead, I need to memorise the core ideas of each theorem and proof in summarised form, but in a way where I can elaborate on those summaries as much as required until I arrive at their full unsummarised forms. I also need to develop a rich mental map of the relationships between the objects of study. Both of these processes develop my reasoning abilities.
2. Coping with abstraction: the objects of study in pure mathematics become increasingly abstract throughout a course of study, but we still need to reason about them. So we develop skills and techniques for coping with abstraction. For instance, we can fall back to geometric intuition (vector space = "space with arrows that can be added to each other and scaled"; orbit-stabiliser theorem = "sending a fixed vertex of a polyhedron to all other vertices then counting symmetries which fix that vertex in place"), as long as we keep in mind the limitations of that geometric intuition. Or we can comprehend abstraction in terms of sufficiently distinct representative examples whose generalisation "spans" the full abstraction (for instance it's difficult to understand abstract measure theory in full generality, so we develop most of our intuition by comparing and contrasting canonical examples: Lebesgue measure, counting measure, Hausdorff measure, Dirac measure, ...). If we think of an abstract concept as a space of instances of examples ("measure theory" as the study of the set of measures), this technique corresponds to choosing a set of examples at different points on the boundary of the space that help us to reason about any instance inside the boundary. I suspect that this skill is exercised more in mathematical study than in any other area of study.

Studying maths for me feels a lot like going to the gym. My interest in number theory or complex analysis is about as deep as my interest in push ups or squats. To me their value is as means of self-improvement rather than as ends in themselves. I deeply respect my lecturers and mathematically-minded peers, but as my study has progressed and the mathematical objects in consideration have become more contrived, I increasingly struggle to understand the motivation behind a career in research mathematics. But of course many academic mathematicians say the exact same thing about careers in tech or finance. Different strokes for different folks...

Developing intuition for Noether's isomorphism theorems

In abstract algebra, Noether's isomorphism theorems are extremely general statements about algebraic objects and the relationships between them. I've studied the theorems in separate courses on groups, rings, and modules over commutative rings (including vector spaces), so I've spent plenty of time trying to grasp their essence. In this post I share some of the intuition that I've developed. If nothing else, this intuition helps to remember the statement of the theorems. I will present the theorems in the context of abelian groups, where the statements are especially straight-forward, in order to avoid clouding the intuition with details about normal subgroups and so on. I will number the theorems 1, 2, 3, 4.

Theorem 1. If $f$ is a homomorphism between abelian groups $A, B$ then $$\frac{A}{\text{Ker } f} \cong \text{Im } f$$

This is the core theorem from which theorems 2 and 3 easily follow. By taking the quotient of $A$ with the kernel of $f$, we are forming a group that looks like $A$ except where we've combined the elements that together map to the same element in the codomain. I think of this as forcing $f$ to become injective by taking all the sets of elements where it fails to be injective and combining those elements into a single element. Since our version of $f$ is now injective, and it maps onto its image (by definition of image), then obviously it's a bijection onto its image. Moreover, the structure-preserving nature of the homomorphism means that structure is preserved through the bijection – i.e. there is an isomorphism between the LHS and RHS.

To me, the statement becomes almost tautological under this interpretation. It seems to say "if we take a homomorphism and disregard all the points where it fails to be injective, then it becomes injective".

Theorem 2. If $A, B$ are subgroups of an abelian group $G$, then $$\frac{A+B}{A} \cong \frac{B}{A \cap B}$$

I think of this statement as a broad generalisation of a basic result from number theory: if $a, b$ are positive integers then $a/\gcd(a,b) = \mathrm{lcm}(a,b)/b$ (this rearranges into the more familiar form $ab/\gcd(a,b) = \mathrm{lcm}(a,b)$. If we draw the diamond-shaped partial order of the "divides" relation with the positive integers $a, b, \gcd(a,b), \mathrm{lcm}(a,b)$ and then we draw the partial order of the "subgroup" relation with groups $A, B, A \cap B, A + B$, we see that the number-theoretic statement on the "divides" partial order looks very similar to what the second isomorphism theorem says about the "subgroup" partial order.

Indeed, $A \cap B$ as the "greatest common subgroup" of $A, B$ and $A + B$ is the "smallest common supergroup" of $A, B$. We know that the quotient operation feels a lot like division, so it seems natural to me to relate the quotient operation with greatest common subgroups and smallest common supergroups to the division operation with greatest common divisors and least common multiples.

The number theory statement can in fact be derived from the second isomorphism theorem applied to the $(\mathbb Z, +)$, so we're not saying anything groundbreaking. The benefit of this view of things is that you probably already have a good intuition for the properties of GCDs and LCMs, so it's just a matter of translating that intuition to the world of abstract algebra.

This isn't a complete intuition. In particular it's not clear under this intuition why we get an isomorphism between the LHS and RHS rather than just a bijection. But it's a useful starting point, or at the very least a mnemonic device.

Theorem 3. If $A \leq B \leq G$ are abelian groups, then $$\frac{G/A}{B/A} \cong \frac{G}{B}$$

This is probably the easiest theorem to remember out of the three. Quotients feel like fractions, and indeed we can "multiply both sides" of the fraction to cancel out quotients of quotients. Unfortunately I don't yet have a good intuition for the preservation of structure.

It's worth keeping in mind that the proofs of the three theorems are all in the "follow your nose" style. For the first theorem we want to establish an isomorphism between the LHS and RHS; it turns out that the most natural possible mapping from the LHS to the RHS turns out to be this isomorphism. For the remaining theorems we again pick the most natural mappings between the LHS and RHS (it'll be easiest to take the mapping from the RHS to the LHS). Again this mapping will happen to be a homomorphism, and the second and third theorems fall out immediately from applying the first theorem to this homomorphism. This is quite remarkable and it means that even if you haven't developed a fully-formed intuition for the theorems, it's easy to convince yourself rigorously that they're correct.

Theorem 4. Let $P$ be a subgroup of an abelian group G. Then there exists a bijective correspondence$$\left\{ \text{subgroups of } \frac{G}{P} \right\} \longleftrightarrow \left\{ \text{subgroups of G that include P}\right\}$$

The right-to-left map looks like this: if $Q \leq P$ is a subgroup satisfying $P \subseteq Q$, then $P$ is also a subgroup of $Q$, so we can form the quotient $Q/P$. Since $Q \leq G$ it would seem natural that $Q/P \leq G/P$. Indeed this is the case, so $Q/P$ is in the LHS. In other words, we map from the right to the left simply by taking the quotient of the subgroup with $P$.

For the left-to-right map, notice that any subgroup of $G/P$ has elements of the form $g+P$ where $g$ is a coset representative. Denote the set of all such representatives by $H$. We can show that $H$ is a subgroup of $G$ (it inherits its structure from $G/P$). Moreover, the representatives of the identity element of the subgroup of $G/P$ are the elements of $P$ itself, so the elements of $P$ are contained within $H$. Thus $P \subseteq H \leq G$ so $H$ lands in the RHS. I visualise this map as "unpacking" a subgroup of the quotient into its representatives, and it makes intuitive sense to me that (1) these representatives form a subgroup of the quotient numerator, and (2) the identity of the quotient (i.e. the coset $0+P$) is in the subgroup and unpacks to the set $P$.

Some work has to be done to prove that these maps are inverses of each other, but if we think of the maps above as packing elements into a quotient and then unpacking the quotient back to elements, the fact that the maps are inverses of each other shouldn't come as a surprise.

Permutations as Sets of Cycle Graphs (and why we care)

This post describes how to transform difficult problems involving permutations into easier problems involving graphs, then lists some problems that require this transformation and gives their solutions. In the process it demonstrates general techniques that may be useful in the face of such problems.

Introduction to permutations

A permutation of a set is a particular ordering of all the elements of the set. If we take the letters a, b, c, d, e, then c, d, b, e, a is a valid permutation, as is a, b, c, d, e or e, d, c, a, b or any other arrangement of these letters. a, c, e is not a permutation as it doesn't use all the elements, and a, a, b, c, d, e is also not a permutation because it uses an element more than once.

A permutation is a bijection from a set to itself, and permuting a sequence corresponds to applying the bijection element-wise to the sequence. If the set is $S = \{a, b, c, d, e\}$ and the permutation $f ~\colon~ S \to S$ is defined by $f(a) = b, f(b) = a, f(c) = e, f(d) = c, f(e) = d$, then applying this permutation to $a, b, c, d, e$ gives the sequence $f(a), f(b), \dots, f(e) = b, a, e, c, d$. The image of $f$ is also its domain so if $x \in S$, then $f(f(x)) \in S$, $f(f(f(x))) \in S$, and so on. In our permutation, $f(a) = b$ and $f(b) = a$, so $f(f(a)) = a$ and $f(f(b)) = b$. We say that $a$ is in a cycle of length 2 because it takes two applications of the permutation to $a$ to get back to $a$. $b$ is also in a cycle of length 2, and in fact we say that $a$ and $b$ are in the same cycle because to get back to $a$ by repeated applications of the permutation, we encounter $b$ and vice versa. In the same vein $c$ is in a cycle of length 3 along with $d$ and $e$. Every element of the set appears in exactly one cycle (note that it is possible to have a cycle of length 1) and cycles don't interact with each other. We say that any permutation can be written as a product of disjoint cycles.

Permutations can be represented as the union of a set of cycle graphs

Let's leave the world of algebra aside for a second and talk about graph theory. A directed graph is a set of vertices $V$ and edges $E$ where each edge maps from one vertex to another . A functional graph is a directed graph in which there is exactly one outgoing edge from each vertex. A permutation has a natural representation as a functional graph: let $V$ be the set on which the permutation is defined, and let there be an edge from $v_1$ to $v_2$ iff $f(v_1) = v_2$. The resulting graph is a specific class of functional graph where each vertex has exactly one incoming edge; in fact, the graph is the union of a set of cycle graphs, with each cycle graph corresponding to a disjoint cycle in the permutation.

Why do we care? This transformation allows us to transform problems about permutations into problems about graphs, and we know how to solve problems about graphs. In the rest of the post, I'll describe three different problems of increasing difficulty and show how to use this transformation to solve them. Without this transformation, these problems would be extremely difficult or impossible. The first problem describes a functional graph directly and the following two problems describe permutations.

Problem 1: a game involving a tennis ball

You are playing a game with $n-1$ other people and one tennis ball. You all sit in a circle and you hold the tennis ball in your hand. On the count of three, two things happen simultaneously:

You yell out a number $k$ between 2 and $n$ inclusive.

Everyone in the circle points at exactly one other person in the circle. You can assume that all participants except for you choose someone in a uniformly random fashion. (By a simple symmetry argument, this means that it doesn't matter who you point to, so we can assume you also point uniformly randomly)

You then throw the ball to the person you're pointing to, and they pass it to the person they're pointing to, and so on, so that the ball is passed a total of $k$ times. Whoever holds the ball at the end of $k$ passes loses the game.

What value of $k$ should you yell to minimise your chance of losing?

The underlying graph structure is immediately evident from the problem description: each player is a vertex and they point to each other to form a functional graph. The graph doesn't necessarily correspond to a permutation because it's possible for a vertex to have no incoming edges.

The solution follows from this line of reasoning: suppose you yell "2" and you lose. This means the ball was passed to whoever you pointed to ("Jane") and Jane returned it to you: you are in a cycle of size 2. Indeed, if you yell "2" you lose iff are in a cycle of size 2. Since each pointing arrangement (i.e. "valid graph" i.e. $n$-vertex functional graph with no cycles of length 1) is equally likely (proof left as exercise to the reader, but it's intuitive that if each person points uniformly randomly at someone else then all valid graphs are equally likely), we conclude that the probability of you losing is exactly the probability of randomly pulling a graph where you are in a cycle of size 2 from the uniform probability distribution over all valid graphs. This probability is equal to the number of valid graphs where you are in a cycle of size 2 divided by the total number of valid graphs.

Suppose you yelled "4" instead. If your graph happens to include you in a cycle of size 4, you will lose; but if you are in a cycle of size 2, you will also lose! And the number of valid graphs where you are in a cycle of size 2 OR a cycle of size 4 is obviously bigger than the number of valid graphs where you are only in a cycle of size 2. It follows that yelling "4" is always a worse choice than yelling 2. The same argument can be made for yelling 6, 8, and so on. In fact, the same argument can be made for yelling any composite number; it's always better to yell one of its prime factors instead.

In fact, this restricts us to yelling only primes. We just have to choose which prime. Some back-of-the-envelope combinatorial calculations show us the probability of being in a cycle of size 2 is $1/(n-1)$ which is larger than the probability of being in a cycle of size 3 which is larger than the probability of being in a cycle of length 5 and so on. The probability of being a member of a large cycle is smaller than the probability of being in a small cycle. We conclude that the solution to the problem is to choose $k$ to be the largest prime no greater than $n$.

The most interesting observation from this solution actually comes from the last few calculations that I skimmed: it turns out that large cycles in random functional graphs or large disjoint cycles in random permutations are extremely unlikely. In fact, the chance that you are in a cycle of length $n$ is $(n-1)!/(n-1)^{n-1}$. When $n=100$, this value in the order of $10^{-44}$.

Problem 2: 100 Prisoners Mk2

This is a variation on the classic 100 Prisoners problem.

The director of a prison offers 100 death row prisoners, who are numbered from 1 to 100, a last chance. A room contains a cupboard with 100 boxes in a row. The director randomly puts one prisoner's number in each closed box. The prisoners enter the room, one after another. Each prisoner may open and look into 50 boxes in any order. The boxes are closed again afterwards. If, during this search, every prisoner finds his number in one of the boxes, all prisoners are pardoned. If just one prisoner does not find his number, all prisoners die. Before the first prisoner enters the room, the prisoners may discuss strategy — but may not communicate once the first prisoner enters to look in the drawers.

The prisoners have also bribed a guard such that, before the prisoners commence the exercise, the guard goes into the room and is allowed to look into all the boxs and then swap the numbers in two boxes if the prisoners' strategy calls for it. The guard has no further communication with the prisoners after entering the room.

There is a strategy allowing all prisoners to survive. Find it.

The solution is to observe that the boxes can be numbered themselves from 1 to 100, in which case box $a$ containing the number of prisoner $b$ represents a permutation $f$ for which $f(a) = b$. Now suppose a prisoner $a$ enters the room and opens box $a$. $a$ points to $b$, which in turn points to $c$, and so on; since he is effectively re-applying the permutation (or walking the cycle graph, however you want to interpret it), he will eventually return to box $a$. But returning to box $a$ means opening a drawer containing prisoner number $a$ inside, which is the terminating condition.

The number of boxes the prisoner opens is exactly the size of the cycle in the graph containing his box. But what if this cycle is larger than size 50? Well, there can only be one such cycle in a graph of 100 vertices, and this is where the bribed guard plays his role! The guard first enters the room and determines if the graph corresponding to the arrangement of numbers contains a cycle of size 51 or more. If so, he simply breaks it into two smaller cycles by swapping two numbers in boxes. The resultant graph has no cycles larger than 50 vertices. Now all prisoners can follow this same strategy and all will survive. Brilliant!

The observation to take away from this problem is the same underlying graph/permutation can be interacted with in different ways (in this case it's walked along in different ways) but we take advantage of the fact that the graph/permutation itself is constant.

Problem 3: repeated number in array

Although the description for this problem is simple, it's an extremely difficult problem with a beautiful solution. It supposedly took Don Knuth 24 hours to solve; I discussed it with my team of Google engineers for over an hour and we concluded that it didn't have a solution. We were wrong.

Unlike the last problem, this problem benefits immensely from reducing it to a graph problem. However, now that you've seen the general strategy to solve problems of this type, I recommend you to attempt to solve it yourself before reading the solution.

Given an array of $n$ elements ranging from $1$ to $n-1$ in which exactly one element appears more than once, find the element in linear time and constant space without modifying the array.

Examples of arrays: [1, 2, 3, 4, 5, 5] or [4, 4, 4, 4, 1] or [6, 5, 2, 3, 5, 4, 1].

There's a trivial linear time and space solution using a flag array, and a fairly simple $O(n \log n)$ constant space solution using a modified binary search, but neither fulfils the requirements of the problem.

To solve it, notice that the first $n-1$ elements of the array are a permutation of the integers $1\dots n-1$, but where zero or the elements have been 'overwritten' with the integer $x$ in this range. If we try to construct the graph representing this almost-permutation, we get a functional graph where the vertex $V_x$ representing $x$ potentially has several incoming edges, and other vertices (the vertices representing overwritten array elements) now have no incoming edges. All components but one in the graph are cycle graphs, and the remaining component is rather messy: it comprises a cycle graph containing $V_x$ with a bunch of path graphs that enter the cycle at $V_x$. Note that the answer to the problem is $x$.

So far, we have ignored the final array element $k$ (represented by vertex $V_n$). $V_n$ has no incoming edges (since an incoming edge would mean there is an element $n$ somewhere in the array, which isn't possible because elements only range up to $n-1$) and one outgoing edge to vertex $V_k$ in the messy component of the functional graph described above. If we follow the edge from $V_n$ we go to $V_k$, and $V_k$ is either on a path graph to the messy component's cycle or it's on the cycle itself. So $V_n$ is the start of a path graph that ends at the messy component's cycle at vertex $V_x$.

Ignoring the vertices not in this path graph or the cycle it leads to, we end up with a rho-shaped graph, also sometimes (misleadingly) called a "linked list with a cycle", and if we find the length of the cycle, we can derive the length of the path graph preceding the cycle and hence find the vertex $V_x$ giving us the value $x$.

But we have efficient algorithms for finding cycles and their lengths! In particular, the hare and tortoise algorithm solves it in linear time and constant space.

In summary, the solution to the problem is to run a linear-time constant-space cycle detection algorithm on the implicit functional graph formed by the array beginning at its last element, and output the index in the array that corresponds to the vertex which starts the cycle.

The take-away here is that even problems that involve not-quite-permutations can be solved by relating the properties of the not-quite-permutation to the properties of permutations.

Conclusion

Problems that require interpreting an array as a permutation, or a permutation as a graph, are fairly rare. However, they often don't have alternative solutions, and if you haven't seen these tricks and techniques before, they can prove to be unsolvable.

Furthest 'Close' Pair in a Set

Another great algorithms question I heard recently:

Given a set $S$ of real numbers, call a pair of numbers $x,y \in S$ close if $x$ and $y$ appear consecutively in a sorted representation of $S$ (i.e. $x,y$ are close iff $S \cap (x,y) = \emptyset$). Find a close pair $x,y \in S$ that maximises $|x-y|$.

The trivial $O(n \log n)$ solution is to sort $S$ using any $O(n \log n)$ comparison sort, iterate through each pair of elements in $S$, and pick the pair whose absolute difference is greatest. However, there is a simple and elegant $O(n)$ solution.

If $S$ contained integers bounded in a reasonably-sized interval, we could modify the trivial solution described above to use a $O(n)$ sort such as radix sort. Unfortunately the elements of $S$ are arbitrary real numbers so no such property holds. However, we will try to use a bucket sort approach anyway. The first step is to decide on a bucketing scheme, and this is where the first "aha!" observation in the solution comes into play:

Claim: a solution $x,y$ to this problem satisfies $|x-y| \geq \frac{\max(S) - \min(S)}{|S|-1}$.

Where does that right-hand-side expression come from? To me it represents a sort of discrete mean value theorem. There are $|S|-1$ close pairs in $S$: the number of consecutive pair of elements in the sorted representation of $S$. Each of these pairs represents a disjoint interval in $\left[\min(S), \max(S)\right]$, and together they partition this entire interval. Say we divide the entire contiguous interval covered by $S$ into $|S|-1$ disjoint uniform partitions of length $L$. In any other partitioning of this interval (including the partitioning denoted by the close pairs in $S$), at least one partition must have length at least $L$ (similarly, at least one partition must have length at most $L$). This is intuitive and easy to visualise, but for the sticklers there's an easy proof by contradiction: suppose that the close pair $x,y \in S$ maximising $|x-y|$ has $|x-y| < \frac{\max(S) - \min(S)}{|S|-1}$. Then the size of the union of all the $|S|-1$ close pair intervals in $S$ is less than $\max(S) - \min(S)$. But the close pair intervals partition $\left[\min(S), \max(S)\right]$ so the size of the union of all close pair intervals must equal $\left|\left[\min(S), \max(S)\right]\right| = \max(S) - \min(S)$.

We will proceed with our bucketing using $|S|-1$ buckets each of size $\frac{\max(S) - \min(S)}{|S|-1} = M$. The first bucket contains elements $s \in S$ where $s < \min(S) + M$, the second bucket contains $s \in S$ where $\min(S) + M \leq s < \min(S) + 2M$, etc. The bucketing process is a linear time operation because for each $s \in S$, a simple mathematical expression tells us which bucket will contain $s$, and we can put it in a bucket in constant time. (Buckets would probably be implemented as linked lists or arrays).

Now what? Well, we've asserted that the elements in the solution pair have absolute distance at least $M$, so they can't both be in the same bucket. We only have to consider pairs of elements in different buckets. Also, if a bucket contains more than 2 elements, we can ignore all but the minimum and maximum in that bucket: no element $x$ in the bucket can form a close pair with an element $y$ of another bucket, because either the minimum or the maximum in the bucket will lie in between $x$ and $y$. So we will proceed as follows: for each bucket, find the minimum in that bucket and the maximum in that bucket, and delete all other elements of the bucket. This is a linear-time operation as every element in $S$ is processed at most twice (once to check if it's a maximum or minimum in its bucket, and potentially once more to erase it from its bucket).

Now we have $|S|-1 = O(n)$ buckets, each with at most two elements. We can simply iterate through each bucket in increasing order and build a list of the elements in all buckets, by adding first the minimum then the maximum in each bucket. This list has $O(n)$ elements, more importantly it's sorted, and most importantly, both the elements in the solution pair are guaranteed to be in the list. To finish off, we simply iterate through each consecutive pair of the list and pick the pair with the greatest difference.

In my opinion this is too much of an "aha!" problem to give a useful signal in an interview (either you see the $M$ trick or you don't), but it's fascinating nonetheless.

Edit: a colleague of mine found a similar solution that follows easier reasoning. Use $|S|+1$ buckets of equal size across the range $\left[\min(S), \max(S)\right]$. Since there are only $|S|$ elements, by [an alternative formulation of] the pigeonhole principle at least one of the buckets must be empty, so the difference between the closest pair is at least one bucket-width. This allows us to ignore all elements that aren't the min or max in a bucket. We now proceed as before.