Two ways to understand overfitting (and don't peek!)

Assorted comments on overfitting

The very first time I was introduced to the notion of overfitting -- by three diagrams of regression curves on scatter plots labelled as you may guess -- I became very uncomfortable. I was probably eleven years old, and did not understand Bayesian statistics, did not realize that seeing 51 heads out of a hundred didn't imply that the coin genuinely had a heads rate of 0.51.

Or rather: I didn't realize that I did realize that. If you had made me bet on the number of heads that would come up in the next hundred flips, I would not provide odds that would indicate an honest belief of seeing 51 heads.

(Because deep down, I had a non-uniform prior.)

And the same principle applies to drawing regression curves. You may insist that your curve with all its squiggles is "unbiased" or give any one of the terms (one for every squiggle) used to describe non-Bayesian estimators, but ultimately -- that's not the curve you'll bet on. You know it's just far more likely that those squiggles are the result of noise which will not be the same (or knowably correlated) in the next sample, than for them to actually be a determining feature.

Yes, maybe the 79th coin toss will always be a tail because of a tiny AI hidden in the coin that counts, or maybe it was the result of factors that affected the 79th coin toss that you just didn't measure. While your data should affect your beliefs, they shouldn't completely overrule your priors.

And that's the key idea behind overfitting (and really the Bayesian notion of probability in general) -- how will your model, based on your data, of perform if you exposed it to data it hasn't yet seen. Because using it on data it hasn't seen is the point of your model -- that is your purpose in building it.

So overfitting occurs when a model learns features specific to your particular data set that don't generalize well. There are two ways this can occur:

• The training set is a biased sample: E.g. MNIST digits are all centered (but suppose you're testing on non-centered digits), or your medical database is all from a particular country. 
• There is noise: I mean, of course there's noise -- it's a statistical problem. Even if the world were deterministic, you still don't have all the information in the world. And you should avoid your model using this noise to make predictions, since noise is, by definition, unpredictable.

The first is a more tractable problem -- it can often be solved by data augmentation (if you have a very good picture of exactly how the data set is biased) or transfer learning (if there's a good chance the model is picking up on important features so you can just train it on the ones it missed). It's important to develop interesting transfer learning algorithms to solve this anyway, as the way humans learn often involves biased samples (e.g. personal experience) and reasoning capacity to unbias their knowledge. 

(Not that humans do always do this -- people often do form beliefs based on mere personal experience, but humans are capable of reasoning more clearly.)

The second is a problem that requires algorithmic solutions that hint our neural network towards Bayesian solutions. Solutions like "well, parameters are just a priori unlikely to be very large, so let's penalize that" (Lasso/Ridge regression), or more complicated (to explain in a simple Bayesian way) regularization algorithms like cross-validation and early stopping.

Another unclassified comment on overfitting: it's very easy to mistakenly "peek" at the test data. Simply in the act of saying "hey, this model works well on the test data, let's choose it" you are already performing a simple algorithm that checks several models and chooses them based on their performance on the test data -- i.e. you are kinda training the model on the test data, even if it is not seen in your code, just in your choice of hyperparameters.

That's why developing some kind of "theory" of hyperparameter optimization and regularization techniques is of importance, so you actually have a theoretical justification for picking your models. 

(Of course, this is hard. We've always been picking models, haven't we? For example when we decide to model something as belonging to a particular family of distributions so we only have to optimize in a 1-dimensional parameter space instead of the literal theory-space. In a sense, machine learning is the way to avoid modeling, due to the universal approximation theorem -- and the hope is that we can eventually make the priors as human-like as possible, which is the eventual goal of hyperparameter optimization and regularization algorithms.) 

Polynomial interpolation and Vandermonde

Suppose you want to find the minimum-degree polynomial $p(x) = \sum {a_jx^j}$ passing through some points $(x_i, y_i)$. This amounts to solving the system:

$$\sum {a_j{x_i}^j}=y_i$$
A first bit of intuition: it seems completely reasonable that any set of $n+1$ points with $x_i$s distinct (so you're actually describing a function), can be interpolated with a polynomial of $n$ degree. What this means is that the matrix $X_{ij}={x_i}^j$, called the Vandermonde matrix, should be square for it to be invertible. Indeed, it's easy to show by considering the degree of the determinant polynomial of the matrix that:

$$\det X = \prod_{1\le i < j \le n} {(x_i-x_j)}$$

---

So the question is of course if there's a simple general expression for the inverse of the Vandermonde matrix.

Here's an idea: if all the $y_i$s were zero, then an $n+1$-degree (NOT $n$-degree! this is a different category of problem, which is not uniquely determined) polynomial would be a constant times $(x-x_0)\dots(x-x_n) $. If we didn't want it to be zero at some particular $x_{i_0}$, we could exclude $x-x_{i_0}$ from the product. Then we can play with the constants so it takes the value we really want ($y_{i_0}$).

In other words, the polynomial

$$L_n^{i_0}(x) = \frac{1}{\prod_{i\ne {i_0}}(x_{i_0}-x_i)}\prod_{i\ne i_0}(x-x_i)$$
Called the Lagrange polynomial gives zeroes at all $x_{i}$ except $x_{i_0}$, where it gives 1. Therefore the polynomial:

$$\sum_i{y_i L_n^i(x)}$$
Which is conveniently of $n$-degree, is the true interpolating polynomial.

---

Conveniently, this approach also tells you what to do when you get a new data point: just add a polynomial that is zero at the existing points and the right adjustment at the added data point. I.e. given $P_n$ is the $n$-degree interpolating polynomial for $(x_0,y_0)\dots(x_n, y_n)$, we want to add $p_{n+1}$, the $n+1$-degree polynomial that is zero at all these points but $y_{n+1}-P_n(x_{n+1})$ at $x_{n+1}$. I.e.

$$P_{n+1}(x)=P_n(x)+\left(y_{n+1}-P_n(x_{n+1})\right)L_{n+1}^{n+1}$$
Alternatively we can also see the added term as a polynomial proportional to $\prod_{i<n+1}{(x-x_i)}$, with the coefficient given by $(y_{n+1}-P_n(x_{n+1})/\prod_{i<n+1} {(x_{n+1}-x_i)}$.

It is easy to show that this coefficient, which we will write as $f[x_0,\dots x_{n+1}]$, can be written recursively as:

$$f[x_0,\dots x_n]=\frac{f[x_1,\dots x_n] - f[x_0,\dots x_{n-1}]  }{x_n-x_0}$$
This is known as Newton's divided difference, and is the coefficient on $x^n$ in the interpolating polynomial -- one may observe that this is a discrete analog of the higher-order derivative (with some attention given to the denominators). It should be perfectly natural that this occurs of course.

Matrix decompositions

Some decompositions

The algorithms for the QR and LU decompositions are fairly self-explanatory from their definition.

The Cholesky decomposition can be computed for a nonnegative-definite matrix through a simultaneous LU decomposition from both sides. Starting with a matrix in the block form:

\[A = \left[ {\begin{array}{*{20}{c}}
  {{a_1}}&{{v^T}} \\
  v&{{A_2}}
\end{array}} \right]\]
We perform an LU decomposition first on the columns, from the left, then on the rows, from the right:

\[\begin{gathered}
  \left[ {\begin{array}{*{20}{c}}
  {{a_1}}&{{v^T}} \\
  v&{{A_2}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&0 \\
  {v/\sqrt {{a_1}} }&1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&{{v^T}/\sqrt {{a_1}} } \\
  0&{{A_2} - v{v^T}/{a_1}}
\end{array}} \right] \\
   = \left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&0 \\
  {v/\sqrt {{a_1}} }&1
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  1&0 \\
  0&{{A_2} - v{v^T}/{a_1}}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
  {\sqrt {{a_1}} }&{{v^T}/\sqrt a } \\
  0&1
\end{array}} \right] \\
\end{gathered} \]
The process can then be continued to give an LU decomposition where the $L$ and $U$ are transpose, thus the Cholesky decomposition.

(Why does the matrix have to be nonnegative-definite?)

Eigenvalue algorithm

A lot of factorisations are fundamentally related to the "eigen"-stuff: this includes the change-of-basis factorizations: diagonalization, SVD, Schur -- and also the polar decomposition. Eigenstuff is fundamentally analytic (whatever that means), so terminating algorithms won't work, we need an asymptotic algorithm.

A simple such algorithm to calculate the eigenvectors and eigenvalues is power iteration. It is related to the flows of systems of linear differential equations, in which the eigenspace with the largest eigenvalue provides the only stable equilibrium of the system. Unless you choose an initial vector $v$ that happens to itself be an eigenvector, $A^nv$ will approach the projection of $v$ onto the principal eigenspace of $A$.

So that gives you the largest eigenvalue. There are many ways to get the remaining eigenvalues -- one inefficient way I thought of was to subtract off the part corresponding to the stretching of the primary eigenvector, but this is incredibly inefficient and numerically unstable.

An actually sensible approach is as follows: consider the matrix $(A-\mu I)^{-1}$, which has the same eigenspaces as $A$. The largest eigenvalue of this is the eigenvalue of $A$ closest to $\mu$. By varying $\mu$ across the real line, one can discover all the eigenvalues of $A$.

Schur algorithm

You can use the eigenvalue algorithm to calculate the Schur decomposition directly via its definition -- however, there is a more efficient method.

A triangular matrix is invariant under conjugation by the $Q$ in its own $QR$ decomposition (because the $Q$ is the identity matrix). So similar to power iteration, one may iteratively QR-factorize $A$ and replace it with $RQ$. This is called the QR algorithm.

Apparently this is expensive in its crude form, so you instead bring into a form called "upper Hessenberg form" (upper-triangular but with a subdiagonal) and that apparently makes the computations less expensive. The mechanism to bring a matrix into upper Hessenberg form is based on householder reflections, and involves doing householder reflections on the parts of the columns below the subdiagonal, so the corresponding householder reflections on the rows (since the Schur is a change-of-basis transformation) do not interfere with the column of interest and screw up all your zeroes. It's not very interesting.

Least-squares

Suppose you wanted to solve $Ax\approx b$, i.e. find the $x$ that minimizes $\|Ax-b\|$. This is relevant when $A$ is not surjective.

The idea is that $\|Ax-b\|$ is minimised when $Ax$ is the projection $b_A$ of $b$ onto the image of $A$. One can equivalently write $A^{T}(Ax-b)=0$, as $Ax-b$ is perpendicular to the column space of $A$. Thus it suffices to solve $A^TAx=A^Tb$.

If $A$ is full-rank (this does not mean surjective), we can actually have a unique $x$, given by $x=(A^TA)^{-1}A^Tb$, and $A^+=(A^TA)^{-1}A^T$ is called the "Moore-Penrose pseudoinverse" of $A$.

An alternative, more efficient algorithm is to use the QR factorization (the one with the non-square $Q$) to construct the projector as $QQ^T$ (which works, as $Q$ definitionally provides a basis for the image, and $QQ^T$ satisfies the defining property of a Hermitian projector). Then one can solve $QRx=QQ^Tb$, or equivalently $Rx=Q^Tb$. In the full-rank case, $R^{-1}Q^T$ is the Moore-Penrose pseudoinverse.

To add:

• SVD algorithm
• Algorithm complexity for each thing

Machine learning as function approximation; statistical and cognitive motivations

The most general formulation of machine learning is that it is the attempt to approximate (or "learn") some mathematical object. For example:

• Regression attempts to approximate a function.
• Classification attempts to approximate an equivalence relation.
• Dimensionality reduction attempts to approximate a parameterization.
• Generative neural networks attempt to approximate a distribution (i.e. a sampler).
• Representation attempts to approximate a compression algorithm.
• Game-theoretic learning (e.g. reinforcement learning, no-regret learning) attempts to approximate a distribution of strategies. 
• Neural networks in NLP attempt to approximate a vector space that effectively represents desired relationships between words.

The Data Science course addresses several "toy algorithms" that discover low-order or simplified approximations. There are two -- different -- explanations of what more "AI-esque" machine learning is, which also lead to two different motivations for machine learning, and two different ways to understand each neural network architecture and each of some ideas in machine learning.

1. The statistical motivation -- machine learning is a non-linear generalisation of "linear" statistical techniques (data mining) like linear regression, PCA and linear decision boundaries. One can always do these linear techniques with some transformation of the data that makes relationships linear (while making sure the transformation is not absurd), but you need a way to "train" what the right such function is. In this sense, machine learning acts as a function approximator.
2. The cognitive science motivation -- a computer should be able to do whatever a brain can, but how exactly does a brain do the stuff it does? To take a simple example, the brain can recognise digits -- well, whatever the brain does, it takes an image as input and outputs a digit, i.e. it's a function. So once again, we need a function approximator.

Great. So machine learning is about making function (or object) approximators. The basic idea is that we're looking for a function that minimises the overall error for a population of data -- it's basically a calculus of variations problem, isn't it? Well, except it isn't, because we don't have access to the entire population, so we need to avoid overfitting (i.e. we need to consider a Bayesian prior). This is also what we meant by "making sure the transformation is not absurd" as we mentioned.

As a general rule, I'd say that if a feature of human brains are present at birth, we should expect to have to hard-code it, while if a feature is learned by humans, we should definitely get our AI to learn it, too (this is the "converse" of a general rule I have when trying to organize knowledge about biology in my head, which is we pretend that evolution is not an algorithm but a hard-coding, simply because actually trying to simulate an evolutionary algorithm in a highly complicated environment is hard). So for example, the transformations that a spatial transformer network finds acceptable are hard-coded, because we didn't actually need to stand on our head or squish our eyes with a truncheon to learn how to read squished-up text. 

Anyway, we want a universal function approximator -- a system that can generate a function arbitrarily close to any given function given sufficiently many parameters. A basic example of such a system is polynomial regression, but for most applications this has the wrong Bayesian prior (it gives zero prior probabilities to high-order polynomials, but most machine learning applications require functions with heavy non-local effects).

Another function approximator is a neural network. That a neural network (even of single layer) is a universal approximator is called the universal approximation theorem, i.e. functions can be written as linear combinations of some scaled and translated sigmoid functions.

(Exercise: explain why the universal approximation theorem is true for the sigmoid function. What other kinds of functions is it true for? It's actually not that hard at all. If you do get stuck, check out the visuals in Michael Nielson's e-book. A rigorous proof can be found here.)

In fact, the universal approximation theorem is not actually particularly important at all to the success of neural networks -- like we said, plenty of systems are universal approximators, but they don't have the right Bayesian prior (and this matters when you have limited data). The single-layer neural network actually has a bad implied prior for many tasks, which is why we usually study "deep" neural networks, which have a surprisingly good prior.

I have seen many explanations as to why this is so: deep learning means doing things in steps, deep learning corresponds to decomposing "hierarchy" or "structure" in the world and our world is "inherently hierarchial", etc. But honestly, these all seem like terrible rationalizations -- I don't even see how these claims are testable. You would need to create a system that isn't "inherently hierarchial" and demonstrate that deep learning doesn't do very well on it (which, as far as I can see, makes no sense). The right way to "explain why deep learning works so well" is to compute the implied prior and find out how it scales with depth.

Someone do this.