Stop using inheritance. Use composition instead.

These are notes from a tech talk I gave at work in Jan 2024. I've removed some details specific to the company.

Introduction

• Inheritance and composition are two popular methods of code reuse.

• Inheritance is when you reuse a class by inheriting from it and extending/modifying it.

• Composition is when you reuse a class by creating and storing an instance as a variable or attribute, and interacting with it via its public interface.

• Aim: to convince you that inheritance is fundamentally toxic and that composition is almost always the right choice.

The scenario

• Scenario: implement "OrderedDict" ("ODict"): a dict type that maintains key insertion ordering, and exposes methods like getNthInsertion(n), popfirst(), poplast(), …

• How it might look:
  
  

Composition	Inheritance
class ODict:   def init():     self.base_dict = dict()     <initialise ordering data here>   def insert(k,v):     self.base_dict.insert(k,v)     <update ordering data here>	class ODict:   def init():     super.init()     <initialise ordering data here>   def insert(k,v):     superclass.insert(k,v)     <update ordering data here>

• By OOP theory, this is a valid use of inheritance! Satisfies "is-a" principle and Liskov substitution.

• They are superficially very similar approaches, but with big differences:

Composition	Inheritance
User knows nothing about relationship between ODict and dict.	New subtype relationship between ODict and dict, affecting all ambient code
Public interface of ODict is isolated to class definition.	Method delegation means that public interface of ODict is the superposition of public interfaces of ODict, dict, superclass of dict, …, spread across different classes and files.
One-way communication. Base dict doesn't know about ODict. All interaction is by ODict invoking methods of dict.	Bi-directional communication due to "open recursion" feature of inheritance. dict can invoke methods on ODict without the user initiating it.

• These differences cause consequences that make dealing with inheritance-based code a living nightmare.

Consequence 1: Loss of encapsulation ⇒ refactoring is harder

• In general, reducing the interconnectedness of systems improves maintainability and comprehensibility

• Encapsulation definition 1: hiding the information internal to a system, e.g. by making the implementation details private.

• Encapsulation reduces and controls the size of the surface at which the external world can connect and interact with a system, which protects against over-interconnectedness.

• Encapsulation applies to structs, classes, modules, services, systems, software teams, …

• In composition, ODict connects to base dict only via dict's public interface. And base dict doesn't connect to ODict.

• In inheritance, ODict connects to base dict via dict's internal implementation details. And base dict can interact with ODict in unclear/unpredictable ways.

• Encapsulation is the mortal enemy of refactoring: the difficulty of refactoring is proportional to the amount of interconnectedness of a system, and inheritance leads to more interconnectedness.

• "Fragile base class" problem: if you inherit from a class, then that base class becomes fragile because its descendants will depend on its implementation details

Consequence 2: Loss of abstraction ⇒ testing is hard

• Abstraction = exposing only the essential characteristics of a system, so that users focus on what's important.

• Abstraction is the most important concept in software design.

• Encapsulation is an implementation-level concern. Abstraction is a design-level concern.

• But encapsulation can be a tool to achieve abstraction, by forcing the user to depend only on a public interface that may not represent the true implementation details hidden beneath.

• Example: a spreadsheet in Google Sheets isn't implemented as a 2D array under the hood. But the 2D array representation is what matters to the user.

• Well-abstracted code is easy to test, because the underlying implementation can be replaced with a simplified or faked version, so that we isolate our testing to the code we wish to test.

• This is an example of the design principle of inversion of control. Composition adapts easily to inversion of control, inheritance doesn't.

• In practice, composition can take advantage of language/library features like monkey-patching or dependency injection.

• Due to the subclass's potential dependence on implementation details of the superclass, it's hard to replace or mock the functionality of a superclass.

• Even if the subclass interacted with the superclass in a clear way, most languages do not provide support for mocking/faking a superclass.

• This is why <SubclassExample> tools are difficult to test! They all inherit from <SuperclassExample>, which has an enormous interface which <SubclassExample> tools depend on. If we could rewrite the <SubclassExample> interface now, we would use a composition-based approach.

Consequence 3: The yo-yo problem

• Encapsulation definition 2: bundling of related data and functionality.

• Encapsulation reduces cognitive load for the user/developer.

• A class written with composition groups all its data and functionality in the class definition, and all interactions with underlying classes takes place via their public interfaces.

• A class in an inheritance hierarchy has its data and functionality spread over its entire inheritance tree.

• Thus, the experience of working with the code involves yo-yoing up and down the inheritance tree, usually across different files, just to build a mental map of the logic or to make any changes.

• A key aspect of clean design is small focused public interfaces. But when we overuse inheritance, our class inherits the interface of all parent classes. This obscures the public interface of the class and makes good design harder to practice.

• You might insist that you're a good developer, and you use inheritance just enough but never overdo it, and that your inheritance trees will never be more than 2-3 deep. But usually you're not the only one working on a codebase, and inheritance for developers like a drug: it's easy and enticing and it solves the problem you have right now, so devs reach for it without considering the long-term harm.

• The <xxx> codebase is an example: it was built by 4 talented senior developers in a small span of time, but these developers didn't have an aversion to inheritance. Within a couple of years, the codebase had inheritance trees of hundreds of classes, including trees eight levels deep. This means you need to have eight files open at once just to understand what the class does.

• Just say no to inheritance!

Consequence 4: Multiple inheritance

• The ODict implemented with inheritance achieves code reuse by inheriting from one underlying class: dict.

• But eventually, we might want to implement its ordering logic in terms of some other class. Say, a binary tree ("bintree").

• At this point we must decide between:

• we can inherit from both dict and bintree. But then we deal with the many many serious pitfalls of multiple inheritance.

• or we can inherit from dict, and keep bintree as an attribute. But now we are inheriting and composing, and suffer the downsides of both. And how do we decide whether we inherit from dict and compose with bintree, or vice versa?

• You might say that an ordered dict "is not" a binary tree, so it shouldn't inherit from bintree. But in general, as a design develops, you will eventually encounter classes that satisfy the "is-a" relationship for two different things.

• Composition simply doesn't have this problem. It is more flexible than inheritance by its very nature. And we get this flexibility "for free", without any disadvantages.

Consequence 5: The single-instance problem

• When ODict inherits from dict, the inheritance mechanism creates one instance of the underlying dict (the one that you get via super() in Python).

• But as a design develops, you will often find that it naturally generalises to be implemented in terms of more than one instance.

• E.g. maybe you want to store the reverse dictionary mapping (value->key) as well, in which case ODict is naturally expressed in terms of two dictionary instances.

• With composition, you just add another attribute! dict1, dict2. Free flexibility once again.

• With inheritance, you can add dict2 as an attribute, but once again you now have a Frankenstein combination of inheritance and composition. Your code ends up interacting with one of its underlying objects via super(), and the others via attribute access.

Conclusion

• Each of these consequences cause severe and long-term maintenance impacts on the code:

• Loss of encapsulation ⇒ code becomes harder to refactor

• Loss of abstraction ⇒ code becomes harder to test

• Yo-yo ⇒ code becomes harder to understand

• Multiple inheritance ⇒ code becomes harder to extend

• Single-instance problem ⇒ code doesn't generalise to >1 instance

• Implementation inheritance is a fundamentally bad feature. Even when used "correctly" – i.e. obeying "is-a" and Liskov substitution etc. –  it still incurs all these problems.

• It follows that you should only use inheritance over composition if its advantages outweigh all of these disadvantages.

• It's very rare that this is the case. It generally only happens when you're interacting with external classes that you don't control, where you might be forced to inherit from those classes to work with some interface, or where the consequence of not inheriting might be that you have to implement dozens of delegation methods.

• Interface inheritance (i.e. inheriting from a pure interface class that has no data or functionality) is great! It enables polymorphism, which is a powerful tool for abstraction. Use interface inheritance to your heart's content – it does not suffer from any of the problems we've mentioned.

• So, next time you see a problem, and you think "I'll use inheritance" – stop. Take your hands off the keyboard. Think really really hard about whether that's what you want to do to yourself and to the people around you. Consider all alternatives. Reach for inheritance as a last resort.

Farmer (IOI 2004) solution and proof

Easy- and medium-difficulty competitive programming problems generally involve a single core technique, for example one graph algorithm, or one computational geometry tool, or one greedy method. Harder problems generally involve more than one technique.

Farmer is probably the easiest problem on the IOI 2004 paper. The task goes: fences consist of fenceposts and the segments in between the posts. Some fences are closed, and these have equal numbers of fenceposts and segments. Other fences are open , and these consist of one less segment than the number of fenceposts. You are given $M$ closed fences (loops) with $m_1, \dots, m_M$ fenceposts each and $K$ open fences (strips) of $k_1, \dots, k_K$ fenceposts each. You must choose any $Q$ fenceposts from among these fences. Whenever you pick two consecutive fenceposts, you earn the the segment in between. What is the maximum number of segments you can earn?

When presented with a non-trivial problem, it's a good idea to try make some observations.

• We obviously want to always choose sequences of consecutive posts where possible.
• Loops are "high-value", in the sense that if we pick all the posts in a loop with $n$ posts we earn $n$ segments. But if we pick a strip of $n$ posts, or a portion of a loop with $n$ posts, we earn just $n-1$ segments.
• Loops are equal value regardless of size. Two loops of size 3 each has the same value as one loop of size 6.
• When presented with a long strip and a short strip, we can always take a portion of the long strip of equal size to the short strip, earning the same number of segments. So there is never a situation where picking a short strip is advantageous to picking a portion of a long strip.
• We earn more segments by taking $n$ posts from a long strip than by taking $n$ posts from two or more short strips. In combination with the previous observation, this means we should always take posts from long strips over short strips.

At this point we can start forming an algorithm. The best case scenario would clearly be to take only full loops, in which case we earn exactly $Q$ segments. To determine if this is possible, we must determine if some subset $m_1, \dots, m_M$ sums to $Q$. This is the subset sum problem which can be solved by dynamic programming. So this is what we check first, and if it's possible than we're done. If not, then we cannot achieve $Q$, since any other selection of posts must select at least one strip or at least one portion of a loop.

If we can't achieve $Q$, then there are two cases to consider:

1. Suppose $\sum_i m_i > Q$. Loops are higher-value than strips, and all loops are equal value, so suppose we just keep taking full loops in any order while we haven't depleted the quota. For each of these loops, we get a number of segments equal to the number of posts. Eventually we will have taken a total of say $u$ posts, earning $u$ segments, and we will reach a loop which has more posts than our remaining quota $Q-u$. If we take a strip of $q$ posts from this loop, we earn $Q-u-1$ segments, for a total of $Q-1$ segments earned. We know that $Q$ segments are unattainable, so if we earn $Q-1$ we have achieved the best possible result. Hence this greedy algorithm solves this case.
2. Suppose instead $\sum_i m_i < Q$. Start by taking all the loops, because they are highest-value. Now only strips remain. By our earlier observations on strips of differing lengths, we can just greedily consume strips in descending order of length while quota remains.

Even in this relatively simple problem, there are two intertwined components that require different techniques: a dynamic programming component followed by a greedy component. One of the most important skills demanded of a competitive programmer is to be able to disentangle a seemingly intractable problem into its components in this way.

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.