Lossless compression, invertible change of variables and entropy

In the last article, we introduced the notion of information, and defined a bit as a unit of information that halved the probability mass on the sample space. We used that to more generally define the information associated with an observation, in intuitive terms of a "number of bits" it possesses.

$$I(X=x)=\log_{1/2}(P(X=x))$$
All this makes perfect sense of course, when talking about a sample space of IID uniform binary sequences (e.g. fair coin tosses), then discovering each term in the sequence (i.e. a bit) literally gives you one bit of information.

But e.g. consider the sample space $\{HH,HT,TH,TT\}$ and consider receiving the information that the result is either $HH$ or $TT$. You haven't actually received a physical bit, but you have halved your sample space, so you have received an information-theoretic bit. 

That's uncomfortable. Maybe we should start looking for the source of our discomfort.

Well, out of all the words in italic sequence above, the only phrase I'm uncomfortable is "physical bit". That's a term that makes sense when we're talking about physical computers, but what does it really mean in an abstract, mathematical, information-theoretic sense? What exactly does it mean to receive the "physical bit" $H$, and how does it differ from receiving information like "the two tosses are the same"?

And then you realize, it doesn't differ. "Receiving the physical bit $H$" just means receiving the information that the first bit is an $H$. That's all the bit means. There is no such thing as a physical bit. (I mean, there is -- we'll use "physical bit" interchangeably with "character" from now on, but the notion of a physical bit particularly having to be something like "$H$" or "$T$", rather than some arbitrary encoding, is something arbitrary.)

Instead of encoding our sample space in terms of "(what was the first coin toss?)(what was the second coin toss?)", we could encode our sample space in terms of "(were the two coin tosses the same?)(what was the first coin toss?)". There is absolutely nothing fundamental about either encoding. And the notion of a physical bit doesn't really exist (think of unfair coins), all we have is the abstract idea of a bit.

This is of course completely expected -- most of our data (not on computers) isn't even in the form of binary sequences.

What we have is an invertible change of variables between two encodings.

---

Often what determines our choice of encoding are cost considerations.

What do I mean? Suppose we have two coin tosses, as before, but instead of being fair, the probability of a heads is 90%. Then HH will appear 81% of the time, HT and TH will both appear 9% of the time, and TT will appear 1% of the time.

Then consider the following two encodings:

• (first toss is an H?)(second toss is an H?)
• (result is HH?)(else: result is HT?)(else: result is TH?)

In the first encoding, the data will always occupy 2 characters. 

In the second encoding, the data occupies 1 character 81% of the time, 2 character 9% of the time and 3 character 10% of the time, i.e. an expected value of 1.29 characters.

This sort of knowledge about probability distributions is used for data compression all the time. For example, images are very likely to have long, continued strings of the same colour, so rather than writing an extra character for each pixel, you just store information like "5 reds, 6 greens..." (this is the idea behind run-length encoding, by the way).

This use of change of encoding is known as data compression -- and the basic idea is to get more common strings to have shorter encodings (the probability of a string's occurrence captures the probabilities of each individual character and the correlations between the characters). 

(and we can have a non-invertible change of variables too, of course, which leads to lossy compression, where we project our data by discarding dimensions we don't care about.)

---

The natural question to ask is what the optimal data encoding should be. This is a very fundamental question -- about how much information is actually present in an object, how many Yes/No questions are actually needed, on average, to specify this particular object in the sample space.

(And if you don't understand why these are the same question: you should re-read the first article in this course: Information and bits.)

Well, the formulation of the problem basically spells out "binary tree". We need to organize the elements of our sample space into a binary tree so that the expected number of yes/no questions we need to ask/bits we need to provide are minimized, so these physical bits are actual bits of information. We organize our data space into a binary tree, with each split (bit) dividing the data into two halves of the probability mass -- this way each physical bit (character) represents an actual bit of information.

Such an encoding -- an optimal binary tree -- is known as entropy encoding; here, the optimal length of code for an object is precisely equal to its information, $\log_{1/2}(P(X))$. As a result, the expected number of bits necessary to specify an object is $\sum P(X)\log_{1/2}(P(X))$, which is known as the entropy of the system.

Simple entropy encodings include the Huffman coding and the Arithmetic coding.

---

In machine learning, we are often more interested in lossy compression -- the problem of determining the optimal representation for a given expected value of "distortion" (calculated via a loss function). Often practically, we have a given size (rate) of representation and calculate the minimum expected loss (distortion) that can be achieved through a compression algorithm.

In general, the distortion is represented by a loss function $d(x,\hat{x})$ -- this loss function is not an information theoretic construct, but derives from the practical motivation to reduce distortion, i.e. the preferences/utility of the decision-making agent choosing the compression algorithm. Typical loss functions include the Hamming (discrete/Dirac) distortion $\delta(x\ne\hat{x})$, appropriate for symbols that are essentially categorical, the mean-squared-error, appropriate for continuous variables, the cross-entropy, appropriate for probability distributions etc.

One may imagine that in the encoder-decoder setting as above, $Y$ should be a lossless encoding  for $\hat{X}$ and thus its length should be the entropy of $\hat{X}$. However, in general, we are only concerned with the information that $\hat{X}$ provides on $X$, i.e. the mutual entropy $I(X;\hat{X})$. Thus the rate-distortion function -- representing the minimal rate for a given distortion -- ought to be given by $R(D)=\inf\limits_{\hat{X}\mid E(d(\hat{X},X)) < D}I(X;\hat{X})$, where the minimization is over probability distributions $p(\hat{X}\mid X)$.

Fourier series and Hilbert spaces

The idea behind Fourier series is to try and express some function on a domain $[-L,L]$ into a sum of complex exponentials of the form $\frac{1}{\sqrt{2L}}e^{2\pi i \ nx/L}$. One of the reasons this is interesting is that the complex exponentials are orthonormal system under the dot product $\int f(x)\overline{g(x)}\ dx$.

One can start by considering the vector space $V$ of all square-integrable functions on $[-L,L]$ -- this gives us a vector space with an inner product. Specifically, we're interested in the subspace $V_n$ that is the span of complex exponentials upto $n$ and $-n$.Then given a vector $f$ in $V$, we can ask for its projection $f_n$ onto $V_n$.

As the complex exponentials are already orthonormal, it is easy to calculate this projection in their basis: 

\[\begin{gathered}
  {a_k} = \left\langle {f,\frac{1}{\sqrt{2L}}{e^{2\pi i\;nx/L}}} \right\rangle  = \int\limits_{ - L}^L {f(x)\frac{e^{ - 2\pi i\;kx/L}}{\sqrt{2L}}dx}  \hfill \\
  {f_n}(x) = \sum\limits_{|k|\le n} {{a_k}\frac{e^{2\pi i\;kx/L}}{\sqrt{2L}}}  \hfill \\
\end{gathered} \]
Notably this implies by Cauchy-Schwarz that:

\[{\left| f \right|^2} \geqslant \sum\limits_{|k| \leqslant n} {{{\left| {{a_k}} \right|}^2}} \]
This really just is Cauchy-Schwarz, and is known as Bessel's inequality. If we can show that the Fourier series approaches $f$, i.e. that $\left\|f-f_n\right\|\to 0$, then it would be obvious that

\[{\left| f \right|^2} = \sum\limits_{|k| \in \mathbb{Z}} {{{\left| {{a_k}} \right|}^2}} \]
Which is just the Pythagoras theorem, and is known as Parseval's theorem. Obviously, these theorems exist in the general theory of Hilbert spaces.

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.

Fancy polynomials

Orthogonal polynomials

Given a dot product on functions of the form $\langle f, g \rangle = \iint_{\mathbb{R}} w(x)f(x)g(x)\, dx$, we can consider orthogonal polynomials, i.e. of the form $\langle f, g\rangle = 0$. A question is if we can generate sequences $p_n$ of these polynomials of consecutive degree -- such a set would then be a basis for polynomials up to that degree (do you see why?).

What we want is to pin down what $p_n$ must be given $p_1,\dots p_{n-1}$. What this means is trying to express the "non-leading part" of $p_n$ in terms of these former terms. We can access a "non-leading part" of the polynomial by normalizing the polynomials to monic and considering the $n-1$-degree polynomial $p_n-xp_{n-1}$. What are the components of this guy in the basis of $p_1,\dots p_{n-1}$? Well, writing

\[{p_n} - x{p_{n - 1}} = \sum\limits_{m < n} {{a_m}{p_m}} \]
Then for all $m<n$,

\[\begin{align}
  {a_m} &= \frac{{\left\langle {{p_n} - x{p_{n - 1}},{p_m}} \right\rangle }}{{\left\langle {{p_m},{p_m}} \right\rangle }} \\
   &=  - \frac{{\left\langle {x{p_{n - 1}},{p_m}} \right\rangle }}{{\left\langle {{p_m},{p_m}} \right\rangle }} \\
   &=  - \frac{{\left\langle {{p_{n - 1}},x{p_m}} \right\rangle }}{{\left\langle {{p_m},{p_m}} \right\rangle }} \\
\end{align} \]
Now, $xp_m$ has degree $m+1$. So if $m+1<n-1$, ${\left\langle {{p_{n - 1}},x{p_m}} \right\rangle }=0$. So the only $m$s that we need to bother about are $n-2$ and $n-1$. Therefore:

\[{p_n} - x{p_{n - 1}} = {a_{n - 2}}{p_{n - 2}} + {a_{n - 1}}{p_{n - 1}}\]
So:

\[\begin{align}
  {p_n} &=  - \frac{{\left\langle {{p_{n - 1}},x{p_{n - 2}}} \right\rangle }}{{\left\langle {{p_{n - 2}},{p_{n - 2}}} \right\rangle }}{p_{n - 2}} + \left[ {x - \frac{{\left\langle {{p_{n - 1}},x{p_{n - 1}}} \right\rangle }}{{\left\langle {{p_{n - 1}},{p_{n - 1}}} \right\rangle }}} \right]{p_{n - 1}} \\
   &=  - \frac{{\left\langle {{p_{n - 1}},x{p_{n - 2}} - {p_{n - 1}}} \right\rangle  + \left\langle {{p_{n - 1}},{p_{n - 1}}} \right\rangle }}{{\left\langle {{p_{n - 2}},{p_{n - 2}}} \right\rangle }}{p_{n - 2}} + \left[ {x - \frac{{\left\langle {{p_{n - 1}},x{p_{n - 1}}} \right\rangle }}{{\left\langle {{p_{n - 1}},{p_{n - 1}}} \right\rangle }}} \right]{p_{n - 1}} \\
   &= - \frac{{\left\langle {{p_{n - 1}},{p_{n - 1}}} \right\rangle }}{{\left\langle {{p_{n - 2}},{p_{n - 2}}} \right\rangle }}{p_{n - 2}} + \left[ {x - \frac{{\left\langle {{p_{n - 1}},x{p_{n - 1}}} \right\rangle }}{{\left\langle {{p_{n - 1}},{p_{n - 1}}} \right\rangle }}} \right]{p_{n - 1}}
\end{align} \]

Examples:

• Legendre polynomials: $w$ is the indicator for $[-1,1]$. Sequence: $1, x, x^2-\frac13, x^3-\frac35x,\dots$ 
• Chebyshev polynomials: $w$ is $(1-x^2)^{-1/2}$ on $[-1,1]$. Sequence: $\cos(n\arccos x)$
• Laguere polynomials: $w$ is $e^{-x}$ on $\mathbb{R}^{\ge 0}$. Sequence: $1, x-1, x^2-4x+2, x^3-9x^2+18x-6$
• Hermite polynomials: $w$ is $e^{-x^2}$ everywhere. Sequence: $1, x, x^2-\frac12, x^3-\frac32 x$

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

Covariance matrix and Mahalanobis distance

In the article Random variables as vectors, we discussed that random variables were vectors, and their covariance was their dot product.

The basic motivation for coming up with this idea was from contrasting $\mathrm{Var}(X+X)=4\mathrm{Var}(X)$ to the formula for variables with zero covariance $\mathrm{Var}(X+Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)$. These correspond to the geometric cases of adding two parallel and perpendicular vectors -- a more general addition is expressed through the cosine rule. What's the "cosine rule for random variables"?

Well, it's $\mathrm{Var}(X+Y)=\mathrm{Var}(X)+\mathrm{Var}(Y)+2\mathrm{Cov}(X,Y)$. To me, this -- like the dot product form of the cosine rule -- is highly suggestive of a bilinear form, specifically the Gram matrix, called the covariance matrix, of the random vector $\mathbf{X}=\left[ {\begin{array}{*{20}{c}}X\\Y\end{array}} \right]$ (which is really to be seen as a "matrix", because the random variables are to be understood as row vectors).

$$\Sigma ({X_1}, \ldots {X_n}) = \left[ {{\rm{Cov}}({X_i},{X_j})} \right]$$
One may compare this Gram matrix interpretation -- $\Sigma=\mathbf{X}\mathbf{X}^T$ (note: not $\mathbf{X}^T\mathbf{X}$, the way we've defined $X$ -- this is important!) -- to the variance formula $\sigma^2=XX^T$, and realise that the covariance matrix is the "right" measure of variance of a random vector (note how if we made random variables column vectors, this would all become $X^TX$, etc.).

(yeah, yeah, you need to subtract the mean, etc.)

Analogously, one may define a cross-covariance matrix $K_{\mathbf{X}\mathbf{Y}}=\mathrm{E}((\mathbf{X}-\mu_{\mathbf{X}})(\mathbf{Y}-\mu_{\mathbf{Y}})^T)$ measuring the covariance between two random vectors.

---

It is rather natural to see this, being a bilinear form, as related to some notion of distance -- the standard deviation, after all, can be seen as a "natural distance unit" in one dimension (in the sense that the "unlikeliness" of a data point depends on its distance from the mean in units of standard deviation).

Suppose we wish to find the variance across some direction, i.e. the variance of some random variable $u_1X+u_2Y=\mathbf{X}\hat{u}$ with $|\hat{u}|=1$ -- this is clearly just $\hat{u}^T\Sigma\hat{u}$. So this defines a natural distance scale in the direction of $\hat{u}$, so that the norm of a vector $\vec{v}$ is defined as:

$$\|\vec{v}\|=\frac{\vec{v}^T\vec{v}}{\hat{v}^T\Sigma\hat{v}}$$
It is not too hard to show -- from the bilinearity of the expression -- that this is equivalent to:

$$\|\vec{v}\|=\vec{v}^T\Sigma^{-1}\vec{v}$$
Another way to interpret is this $\Sigma^{-1}$ maps the distribution into a spherical one (one with identity covariance matrix), and this norm is just the norm of the data point in this spherical distribution, which is adjusted for variances and covariances. This measure of distance is called the Mahalanobis distance.

Mixed states I: density matrix, partial trace, the most general Born rule

In the last article, we saw that sub-systems entangled with other sub-systems did not have well-defined pure states themselves -- just like correlated random variables don't have their own probability distributions. Since pretty much everything you see in the real world is entangled with something -- has correlations with some other thing -- this is a problem. One can't just consider the "state of the entire universe" when you just want to study a single electron or something.

Wait -- why can't we just consider the marginal distributions, like we do in statistics? OK, suppose we start with the system -- with $|\phi\rangle$ and $|\varphi\rangle$ an orthonormal basis:

$$|\psi\rangle = \frac1{\sqrt2}|\phi\rangle\otimes|\varphi\rangle+\frac1{\sqrt2}|\varphi\rangle\otimes|\phi\rangle$$
Naively, you may think that the state of the first sub-system $|\psi_1\rangle$ may be given by $|\psi'_1\rangle=\frac1{\sqrt2}|\phi\rangle+\frac1{\sqrt2}|\varphi\rangle$. Certainly, if we're measuring the subsystem with an operator with eigenvalues $|\phi\rangle$ and $|\varphi\rangle$, you have 50% probabilities of each. But to say that two things are in the same state requires that they produce the same outcome for any measurement, not just that one. Does our sub-system behave exactly like $|\psi'_1\rangle$ for all observables? Recall that in the last article, we showed that collapsing the first sub-system onto $|\chi\rangle$ collapses the entire system into the state:

$$|\chi\rangle\otimes\left(
\langle\chi|\varphi\rangle|\phi\rangle+\langle\chi|\phi\rangle|\varphi\rangle\right)$$
To calculate the probability amplitude of this collapse, we may take the inner product of this with the original state -- you can compute this, and see the answer comes down to $1/\sqrt2$, i.e. there's a probability of $1/2$ of the first subsystem collapsing to any such eigenstate $|\chi\rangle$. You can use any observable in this two-dimensional state space, and the sub-system would collapse into either eigenstate with probability exactly $1/2$.

This is a completely different situation from if the state of the first subsystem were simply a pure state like $|\psi'_1\rangle$.

The situation we're dealing with is called a mixed state -- an example of a mixed state, in line with the motivating examples we had at the beginning of the course -- is unpolarised light. In fact, the state we described above models precisely unpolarised light involving two photons (is it obvious why?).

The basic idea behind mixed states is that we have some uncertainty as to what the state of a particle is -- we don't know if the particle has state $|\phi\rangle$ or $|\varphi\rangle$ -- it has a 50% chance of either. This is a classical probability, rather than a quantum one, and different from the state being a superposition of these states, as we just saw above.

Does this sort of thing occur with multivariate distributions in statistics? Suppose we have a multivariate distribution $\psi(x,y)$ and extract the marginal $x$-distribution $\phi(x) = \int_y \psi(x,y) dy$. Certainly this $\phi(x)$ gives us the right probability densities of each $x$-value. But the analog of considering general states like $|\chi\rangle$ is to make a transformation of the domain -- like a Fourier transform -- and consider probability densities in the transformed domain.

As an exercise, write down a multivariate Fourier transform expression for $\hat{\psi}(\omega_1,\omega_2)$ and use it to compute the $\omega_1$-marginal probabilities $\hat{\phi}(\omega_1)$ -- compare this to what you would get if you were to Fourier-transform the $x$-marginal $\phi(x)$ directly.

---

But saying "it has 1/2 probability of being in $|\phi\rangle$ and 1/2 probability of being in $|\varphi\rangle$" is clearly an overdetermination. As we saw above, this resulting state has a 1/2 probability of collapsing onto any state -- this is a statement that doesn't depend on $|\phi\rangle$ and $|\varphi\rangle$, the behaviour of the state is the same if you describe it instead as having "1/2 probability of being in $\frac1{\sqrt2}(|\phi\rangle+|\varphi\rangle)$ and 1/2 probability of being in $\frac1{\sqrt2}(|\phi\rangle-|\varphi\rangle)$". These two are different statistical ensembles but in the same mixed state.

You can see similarly that "50% left-polarised + 50% right-polarised" is the same mixed state as "50% left-circular + 50% right-circular" -- they're both just unpolarised light.

What's the general condition for two statistical ensembles to produce the same observations?

Given statistical ensemble $\left(\left(p_i,|\psi_i\rangle\right)\right)$ and $\left(\left(p_i,|\psi'_i\rangle\right)\right)$, they are the same mixed state if for all $|\chi\rangle$, the probabilities of collapsing onto $|\chi\rangle$ is the same, i.e.

$$\sum_i p_i|\langle\psi_i|\chi\rangle|^2=\sum_i p_i|\langle\psi'_i|\chi\rangle|^2$$
Well, each side of this equation is just the evaluation of a quadratic form for the vector $|\chi\rangle$ -- and two quadratic forms are identically equal on all vectors if and only if their matrix representations are the same. Well, what's the matrix representation? In the basis of $|\psi_i\rangle$, it's just the matrix of probabilities $p_i$. The way to write this in Bra-ket notation is to factor out the $|\chi\rangle$s:

$$\left\langle \chi  \right|\left( {\sum\limits_i {{p_i}\left| {{\psi _i}} \right\rangle \left\langle {{\psi _i}} \right|} } \right)\left| \chi  \right\rangle  = \left\langle \chi  \right|\left( {\sum\limits_i {{p_i}\left| {{{\psi '}_i}} \right\rangle \left\langle {{{\psi '}_i}} \right|} } \right)\left| \chi  \right\rangle
$$
This quadratic form in between, representing a mixed state, is called the density matrix and can be used to completely specify mixed states. In this sense, it is a generalisation of the state vector, which can only be used to represent pure states.

$$\rho={\sum\limits_i {{p_i}\left| {{\psi _i}} \right\rangle \left\langle {{\psi _i}} \right|} }
$$
You may confirm that indeed:

$$\frac12|\phi\rangle\langle\phi|+\frac12|\varphi\rangle\langle\varphi|=\frac12\left(\frac{|\phi\rangle+|\varphi\rangle}{\sqrt2}\frac{\langle\phi|+\langle\varphi|}{\sqrt2}\right)+\frac12\left(\frac{|\phi\rangle-|\varphi\rangle}{\sqrt2}\frac{\langle\phi|-\langle\varphi|}{\sqrt2}\right)$$
In fact, there is a simpler way to see that those two ensembles are the same: the density matrix is simply the Gram matrix of the ensemble -- you take the states in the ensemble, weighted by $\sqrt{p_i}$ in a matrix $Y$, and $\rho=Y^*Y$. Well, $Y^*Y=Y'^*Y' \iff Y'=UY$  for some unitary $Y$, i.e. the ensembles are rotations of each other.

---

Properties of the density matrix, generalised Born's rule, etc.

Here's something that's obvious: the density matrix is nonnegative-definite ("positive-semidefinite") Hermitian and unit-trace -- and all such matrices can represent density matrices.

Well, so it's a Hermitian operator -- does it represent any interesting observable? Not really. It's an observable, sure, but not an interesting one (you might say it measures something's being in one of the ensemble states -- written in an orthonormal basis -- and whose eigenvalues are the mixing ratios, etc. -- but what if two mixing ratios are the same? Its behaviour is just bizarre and useless, really).

We saw earlier that the probability of a density matrix collapsing into a state $|\chi\rangle$ is given by $\langle\chi|\rho|\chi\rangle$.

This is completely different from the generalised Born's rule we saw earlier which took the form $\langle\psi|L|\psi\rangle$! There, the state was the vector and the information on the projection space was the quadratic form in between. Here, the state is the quadratic form in between while the state being projected onto is the vector. This is just a generalisation of the simple Born rule $\langle\chi|\psi\rangle\langle\psi|\chi\rangle$, as far as I can see. If anyone comes up with a connection between it and the generalised Born rule for pure states, tell me.

This brings the question, though -- what's the most generalised Born's rule we can come up with? What is the probability of a mixed state collapsing into some eigenspace of a Hermitian projection operator?

Well, given the ensemble $((p_i,|\psi_i\rangle))$ (you can start writing ensembles with their density matrices now if you like, like $\sum p_i|\psi_i\rangle\langle\psi_i|$ -- but I just want to reaffirm that our result will indeed be in terms of the density matrix), the probability is:

$$\sum_i p_i\langle\psi_i|L|\psi_i\rangle$$
This is hardly useful -- it's not in terms of the density matrix at all. But look at each term -- what's $p_i\langle\psi_i|L|\psi_i\rangle$? $L|\psi_i\rangle$ is the $i$th column of $L$ in the $(|\psi_i\rangle)$-basis -- the inner product $\langle\psi_i|L|\psi_i\rangle$ is the $i$th entry of this column. Multiplying this by $p_i$ gives us the dot product of the $i$th row of $\rho$ with the $i$th column of $L$. The sum of these for all $i$ gives us the trace of $\rho L$:

$$\ldots = \mathrm{tr}(L\rho)$$
This is the most general form of Born's rule. Note that our derivation could have also applied to finding the expectation value of a general operator $A$ under the density matrix $\rho$ (recall that Hermitian projection operators are basically "indicator variables" whose expectation values represent probabilities), indeed generally:

$$\langle A\rangle_\rho=\mathrm{tr}(L\rho)$$
(note that $\mathrm{tr}(V)=\sum_{i} \langle i | V|i \rangle $ for any basis $(|i\rangle)$, which you should show.)

It is also trivial to show that upon collapse given by Hermitian projection operator $L$, the density matrix collapses to:

$$\rho'=\frac1{\mathrm{tr}\left(L\rho L\right)}{L\rho L}=\frac1{\mathrm{tr}\left(L\rho\right)}{L\rho L}$$
Generalising the pure state collapse to $|\psi'\rangle=\frac{1}{\langle \psi | L | \psi\rangle}L|\psi\rangle$. One may check that the above expression reduces to $|\chi\rangle\langle\chi|$ in the case where $L=|\chi\rangle\langle \chi|$.

---

Partial trace, trace

We started our discussion considering the pure state $\frac1{\sqrt2}|\phi\rangle\otimes|\varphi\rangle+\frac1{\sqrt2}|\varphi\rangle\otimes|\phi\rangle$ and asking for the mixed state of the first sub-system. We computed the inner product of this state with its projection under the operator $|\chi\rangle\langle\chi|\otimes1$ -- this tells us the evaluation of the quadratic form $\langle\chi|\rho_A|\chi\rangle$ at all vectors $|\chi\rangle$, which determines the quadratic form $\rho_A$ of the first state.

So what exactly did we do -- in general? Starting with a density matrix $\rho$ on $H_1\otimes H_2$, we compute the probability of the first sub-system appearing in state $|\chi\rangle$: it's $\mathrm{tr}((|\chi\rangle\langle\chi|\otimes 1)\rho)$. So we try to find a density matrix $\rho_1$ satisfying, for all states $|\chi\rangle$:

$$\mathrm{tr}((|\chi\rangle\langle\chi|\otimes 1)\rho)=\mathrm{tr}(|\chi\rangle\langle\chi|\rho_1)$$

Exercise: Let $V$ be an operator on $H_1\otimes H_2$. We define its partial trace on $H_2$ as $\mathrm{tr}_2(V)=\sum_{j}\langle j|V|j\rangle $ for basis $(|j\rangle)$ of $H_2$ (where the inner product is done by extending operators by tensoring them with the identity). Show that the density matrix $\rho_1$ is given by:

$$\rho_1=\mathrm{tr}_2(\rho)$$
I.e. show that for operators of the form $A\otimes I$: $\mathrm{tr}[(A\otimes I)\rho]=\mathrm{tr}_1(A\,\mathrm{tr}_2\rho)$.

Projection operators, generalised Born's rule, position basis, wavefunction

At the end of the last article, I asked you to investigate Born's rule for continuous variables like position and momentum.

Well, the problem is that if $x$ is continuously distributed (i.e. we have an operator $X$ whose eigenvalues form a continuous spectrum $\Sigma_X$), typically $P(x=\lambda)=0$ -- and this gives us very little information about the actual probability distribution. What we're really interested in is $P(x\in B)$ for $B$ some subset of $\Sigma_X$.

Technically, we need $B$ to be a "Borel subset", or "measurable subset". We will be omitting several such technicalities in the article, such as the need for the spectral theorem to define a "projection-valued measure" or "spectral measure" on an operator with a continuous spectrum -- this is something that will be covered in the MAO1103: Linear Algebra course.

First, let's think about $P(x\in B)$ in the countable case. One can write $B=\{\lambda_1,\ldots\lambda_n\}$, and then simply say that

$$P(x\in B)=\sum |\langle\psi|\phi_k\rangle|^2$$
But the term on the right is a Pythagorean sum -- specifically, it is the length-squared of the vector formed by summing all the projections of $|\psi\rangle$ onto the eigenstates $|\phi_1\rangle\ldots|\phi_k\rangle$. But this is the same as the length of the projection of $|\psi\rangle$ onto the span of these eigenstates.

(Note on notations: From here onwards, we will use the notation $|\lambda\rangle$ to refer to the eigenvector corresponding to the eigenvalue $\lambda$ (if the eigenspace has dimension more than 1, we'll figure something out). We will use the notation $\{|B\rangle\}$ to refer to the span of the eigenvectors corresponding to the eigenvalues in $B$.)

So we could just define a Hermitian projection operator $L_X(B)$ for any subset $B$ of the spectrum of $X$ -- it is an easy exercise to write down an explicit form for $L_X(B)$ in terms of the eigenvectors of $X$.

Then the probability $P(x\in B)$ is simply $|L_X(B)|\psi\rangle|^2$. Recalling that a Hermitian projection operator satisfies $L^*=L=L^2$, we can write the generalised Born's rule as:

$$\begin{align}P(x\in B) &= |L_X(B)|\psi\rangle|^2\\ &= \langle\psi|L_X(B)|\psi\rangle\end{align}$$
Well, this is interesting! In the last article, you proved that the expected value of an observable $X$ given a state $|\psi\rangle$ is given by $\langle\psi|X|\psi\rangle$. But here we have a probability given by the same expression. So we want to interpret our projection operators as some sort of "observable" -- we can omit the "Hermitian", since all observables are Hermitian.

There's another place you might've seen something like this, and that is with indicator variables in probability and statistics -- the expected value of an indicator variable for an event is the probability that the event occurs.

Try to interpret these projection operators as observables that are analogous to "indicators" in some sense. If you think a little about it, you might see exactly what these observables represent: the eigenvalues of $L_X(B)$ are all 1 and 0 -- if the value "1" is realised, the state has been projected into the $\{|B\rangle\}$ -- and if the value "0" is realised, it hasn't.

So projection operators are a special type of observable, measuring the answer to "Yes/No questions" -- if the answer to "is the system in one of the states $\{|B\rangle\}$?" is yes, the observable $L_X(B)$ takes the value 1 -- if the answer is no, then it takes the value 0. So it is precisely an "indicator variable" for $\{|B\rangle\}$.

We have seen such projection operators, of course, in the context of polarisation -- where the operator represented whether or not the photon has passed through. Indeed, one may formulate quantum mechanics entirely in terms of projection operators, as any question can be formulated with some number of Yes/No questions (the key reason why this can be done, as we will see -- is that these "yes/no questions" all commute, i.e. the corresponding projection operators share an eigenbasis). Let's not.

---

Well, this can be generalised in the straightforward way to an operator with a continuous spectrum, resulting in the same expression. We can also calculate probability densities using this result. Let $X$ be an operator with continuous spectrum $\Sigma_X$ -- then we can write the state $|\psi\rangle$ in the eigenbasis of $X$:

$$ |\psi\rangle = \int_{\Sigma_X} |x\rangle\, \Psi(x)\, dx $$
Where $\Psi(x)\, dx=\langle\psi|x\rangle$ are the coefficients of the state in the eigenbasis, i.e. the probability amplitudes -- we call $\Psi(x)$ the wavefunction, and it represents probability amplitude densities. Then for some set $M\subseteq \Sigma_X$ of eigenvalues $L_X(B)|\psi\rangle$ is the projection:

$$ L(M)|\psi\rangle = \int_B |x\rangle\, \Psi(x)\, dx $$
And one may calculate the dot product, noting that complex dot products require taking the complex conjugate:

$$ \langle \psi | L_X(B) | \psi \rangle = \int_B \Psi^*(x)\, \Psi(x)\, dx $$
Which gives us an expression for the probability density function on $\Sigma_X$ as:

$$\begin{align}\rho(x) &=\Psi^*(x)\,\Psi(x) \\
&=|\Psi(x)|^2\end{align}$$
And this applies to any operator with a continuous spectrum, like position and momentum.

---

Some texts define the eigenvectors $|x\rangle$ of a continuous-spectrum observable differently from us -- it is often conventional to let $|x\rangle$ be infinitely large so that $\langle x_1|x_2\rangle = \delta(x_1-x_2)$. This is so that the amplitudes $\langle\psi|x\rangle$ are not infinitesimal, but instead $\langle\psi|x\rangle=\Psi(x)$ (without multiplication by $dx$). For consistency with discrete spectra, we do not use this convention.

From polarisation to quantum mechanics: states, observables, Born's law

Like most texts on the theory, I will motivate the mathematics of quantum mechanics from the example of polarisation -- mostly because it's a very accessible example of stuff being wavelike. From this example, we will be able to motivate: the state vector (generalising the polarisation), state vector collapse (the event of polarisation), observables and their eigenvalues (stuff like energy, number of photons, etc.), eigenstates and their orthogonality (polarisation basis), noncommuting operators and uncertainty (the noncommuting of lenses).

The key feature of quantum mechanics -- the fundamentally probabilistic nature -- comes from the following two facts, confirmed by experiments (the famous experiments here are the double-slit experiment and photoelectric effect respectively):

• Everything is a wave -- objects behave as waves, following the superposition principle and the waves represent densities of observations at large scales.
• Everything is a particle -- which manifests itself in the form of some stuff, like energy and momentum, coming in little quanta.

This is the principle of wave-particle duality. You may realise how this implies a probabilistic description, but the following example should make it quite clear: consider a wave of light, with energy $hf$ (so it's a single photon) polarised at angle $\theta$ to the horizontal -- and it passes through a horizontal polarising filter. Well, then the wave that passes through would be a horizontally polarised wave with energy $hf\cos^2\theta$, right?

But this is impossible, since energy levels in quantum mechanics are quantised -- you can't have $\cos^2\theta$ of a photon, you can only have integer multiples of a photon. But the fact that energy drops as $E\cos^2\theta$ is something that you can verify at your home, using sunglasses -- what the heck?

The key point is that the empirical verification of the $\cos^2\theta$ business that you can do at home is on a macroscopic level, when you have a large number of photons $E=Nhf$. So something occurs with the photons on a microscopic level such that when you try it with a large number of photons, $\cos^2\theta$ of the photons pass through.

Well, this is essentially the "definition" of probability! A single photon passes through the filter with a probability of $\cos^2\theta$ so that for a large number of photons, $\cos^2\theta$ of the photons pass through. This is a non-trivial result -- wave-particle duality makes no mentions of probability as such, it just tells us that stuff is both a particle and a wave, but this simple condition in itself implies a probabilistic, non-deterministic reality.

---

Similar thought experiments can illustrate the probabilistic nature of other things (the "things" in question here will soon be called "eigenstates"): position is easy -- consider a standing wave photon in a box (this can easily be constructed). This is uniformly distributed throughout the box -- so how much of the energy is in some chunk of the box?

Momentum is trickier, but shouldn't be too hard if you're familiar with Fourier transforms -- what's the analog of a "box" in momentum-space? Well, consider a concentrated pulse of light -- this can be written, via a Fourier transform, as the sum of several light waves of different momenta (i.e. frequencies), each wave with some lower energy. Taking "some chunk" of this "box" amounts to filtering some specific frequencies of the light. This can be done easily, e.g. with a colour filter -- so how much of the energy is contained in the waves with these specific momenta?

In both cases, the key point is that you can't have a fraction of the energy of the photon at these positions/momenta, so you must have a probability of measuring the photon to be in a specific range of positions or a specific range of positions -- to be in a specific region or in a specific region of momentum-space.

---

The fundamental point here can be made for any quantity $X$: if you can filter out the "part" of a collection of particles that has $X$ in a certain subset of its range, then on a microscopic level, is probabilistic. The act of "filtering out the parts with a certain $X$", applied to a single particle, is just the act of checking if a particle is in a certain $X$-interval, and is called measurement. Any quantity that you can measure is called an observable. 

Something like polarisation is really a form of measurement -- you're finding out whether or not the photon is in a certain polarisation $|\phi_{\parallel}\rangle$. You may have another observable, corresponding to a different polarisation -- even one that is orthogonal to the first polarisation -- $|\phi_{\perp}\rangle$ and still get that the photon is in $|\phi_\perp\rangle$. There is nothing wrong with this, as we just know beforehand that the photon is in $|\phi_\parallel\rangle$ or $|\phi_\perp\rangle$. If you perform the polarisation with $|\phi_\perp\rangle$ after the polarisation with $|\phi_\parallel\rangle$, you will find that the photon doesn't pass through, as you know for sure that the photon is not in both $|\phi_\parallel\rangle$ and $|\phi_{\perp}\rangle$.

Now, you may have certain psychological issues with this, as have many in history -- however, you might want to note that the aim of quantum mechanics is not to fix your psychological problems but to explain nature. You need to accept logical positivism and learn to shut up and calculate to be comfortable with quantum mechanics.

So whatever calculus we invent to describe these probabilistic phenomena, it is going to apply to all observables.

In our first example, the polarisation of the photon can be represented by a unit vector which we will denote as $|\psi\rangle$. The polarising filter has two special axes, represented by unit vectors $|\phi_{\parallel}\rangle$ and $|\phi_\perp\rangle$ -- these are special in the sense that an incoming photon polarised as $|\phi_{\parallel}\rangle$ or $|\phi_\perp\rangle$ will simply be scaled, by factors of 1 and 0 respectively -- so these form an eigenbasis for a certain operator.

Well, we said that the photon passes through (with polarisation $|\phi_{\parallel}\rangle$) with probability $\cos^2\theta$ -- this arises simply from considering the amplitude of $|\psi\rangle$ in the direction of $|\phi_\parallel\rangle$. So we can write the probability that the photon ends up in a state $|\phi\rangle$ as $|\langle\psi|\phi\rangle|^2$ where $\langle\psi|\phi\rangle$ is called the corresponding "probability amplitude".

This expression, $P(x=\lambda)=|\langle\psi|\phi_\lambda\rangle|^2$ is called Born's rule.

Let's get back to the eigenbasis -- what exactly is this an eigenbasis of? We said that the corresponding eigenvalues are 1 and 0, so this gives us a complete description of the operator. Note that this operator depends only on the observable (namely "number of photons in the $|\phi_{\parallel}\rangle$ direction), not on the state or any other feature of the observation. So we decide to call this operator/matrix the "observable", and its eigenvalues are the values of the observable that can be measured.

To find properties of these observables, the natural way is to note that the only feature we've really required of them is Born's rule, i.e. the probabilistic interpretation -- so we can apply the axioms of probability and see what they apply in the context of these observables.

• $P(E)\ge 0$ -- imply that the observables are over either the reals or complexes, so that $|\langle\psi|\phi\rangle|^2\in \mathbb{R}$ in the first place. The nonnegativity then follows.
• $P(\Omega)=1$ and $P\left(\bigcup_i E_i\right) = \sum_i P(E_i)$ for disjoint $E_i$ -- this, along with the second axiom, implies that $\sum |\langle\phi|\psi\rangle|^2 = 1=|\langle\psi|\psi\rangle|^2$ where the sum is taken over all eigenstates $|\phi\rangle$ of the operator. As this must be true for all states $|\psi\rangle$, the thing on the left must be a Pythagorean sum, so the $|\phi\rangle$s must form an orthogonal basis. This implies that all observables are normal operators.

The latter fact is very important, and can also be seen in the following way -- if you a system is in one eigenstate, it cannot possibly collapse onto another eigenstate (the probabilistic interpretation is: if you know for sure the value of the symbol is a thing, it's that thing) -- so we must have $|\langle \phi_1|\phi_2\rangle|^2=0$ for all eigenstates $|\phi_1\rangle$ and $|\phi_2\rangle$.

Another restriction we add is that the observables be not only normal, but Hermitian operators in particular, so they have real eigenvalues. This may seem an odd choice, but it makes sense, as any normal operator may be uniquely written as $X_H+iX_{AH}$ where $X_H$ and $X_{AH}$ are Hermitian, and $X_H$ and $X_{AH}$ commute, so any complex observation can be done unambiguously as two real observations. So we stick to real eigenvalues.

This also makes it essential that we allow complex operators rather than just real ones (the two choices were given to us from the first probability axiom), so that this decomposition is possible. Later, we will see concrete examples of this with commutators $[X,Y]$, which must be multiplied by $i$ to turn Hermitian. We will also see more fundamental reasons to choose complex numbers in QM.

---

Exercise: Show that the expected value of an observable $X$ given a state $\psi$ can be given as $\langle \psi|X|\psi \rangle$ (i.e. $\psi^*X\psi$ in conventional notation).

Exercise: Explain Born's rule with other observables, like position and momentum. Explain why it holds in general.

The Cauchy Riemann Equations: what do they really mean?

Question: Geometrical Interpretation of Cauchy Riemann equations?

One might think that being differentiable on $\mathbb{R}^2$ is sufficient for differentiability on $\mathbb{C}$. But the Jacobian of an arbitrary such function doesn't have a natural complex number representation.

$$
\left[ {\begin{array}{*{20}{c}}
{\partial u/\partial x} & {\partial u/\partial y} \\
{\partial v/\partial x} & {\partial v/\partial y}
\end{array}} \right]
$$
Another way of putting this is that no complex-valued derivative (see below for an example, known as the Wirtinger derivative) you can define for an arbitrary function fully captures the local behaviour of the function that is represented by the Jacobian.

$$
\frac{df}{dz} = \left(\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right) + i\left(\frac{\partial v}{\partial x}-\frac{\partial v}{\partial y}\right)
$$
The idea is that we should be able to define a complex-valued derivative "purely" for the value $z$, without considering directions, i.e. we want to consider $\mathbb{C}$ one-dimensional in some sense (the sense being "as a vector space"). More precisely, the derivative in some direction in $\mathbb{C}$ should determine the derivative in all other directions in a natural manner -- whereas on $\mathbb{R}^2$, the derivatives in *two* directions (i.e. the gradient) determines the directional derivatives in all directions.

If you think about it, this is quite a reasonable idea -- it's analogous to how not every linear transformation on $\mathbb{R}^2$ is a linear transformation on $\mathbb{C}$ -- only spiral transformations are.

$$
\left[ {\begin{array}{*{20}{c}}
{a} & {-b} \\
{b} & {a}
\end{array}} \right]
$$
How would we generalise differentiability to an arbitrary manifold? Here's an idea: a function is differentiable if it is locally a linear transformation. So on $\mathbb{R}^2$, any Jacobian matrix is a linear transformation. But on $\mathbb{C}$, only Jacobians of the above form are linear transformations -- i.e. the only linear transformation on $\mathbb{C}$ is multiplication by a complex number, i.e. a spiral/amplitwist. So a complex differentiable function is one that is locally an amplitwist (geometrically), which can be stated in terms of the components of the Jacobian as:

$$
\begin{align}
\frac{\partial u}{\partial x} & = \frac{\partial v}{\partial y} \\
\frac{\partial u}{\partial y} & = - \frac{\partial v}{\partial x} \\
\end{align}
$$
This is precisely why you shouldn't (and can't) view complex differentiability as some basic first-degree smoothness -- there is a much richer structure to these functions, and it's better to think of them via the transformations they have on grids.

One might observe that the Cauchy-Riemann equations can also be restated as: $\frac{\partial f(z)}{\partial\bar{z}}=0$, where the derivative is the Wirtinger derivative defined above. This seems completely bizarre to me, but apparently it makes sense with some differential geometry, and I'm given the keyword "complexifying the tangent bundle".

Trace, Laplacian, the Heat equation, divergence theorem

The aim of this article is to help build an intuition for the trace of a matrix, "the sum of the elements on the diagonal" -- the basic idea is that the trace is an "average" of some sort, an average of the action of an operator or a quadratic form. We'll make this idea clearer with an example from classical physics: the heat equation.

---

Consider an $n$-dimensional space with some temperature distribution $T(\vec{x},t)$. We wish to set up a differential equation for this function.

In the case that $n = 1$, this differential equation is exceedingly easy to write down, considering the difference $(T(x+dx)-T(x))-(T(x)-T(x-dx))$ as the double-derivative upon division by $dx^2$. More rigorously, what we're doing here is applying a localised version of the fundamental theorem of calculus. I.e. we're writing down:

$$\begin{align}
\lim_{\Delta x \to 0} \frac{1}{\Delta x}(T'(x + \Delta x) - T'(x)) &= \lim_{\Delta x \to 0} \frac{1}{{\Delta x}}\int_x^{\Delta x} {T''(x)dx}  \\
& = T''(x)
\end{align}
$$
More generally, we may consider the $n$-dimensional case.

Analogously to before, one may try to look at temperature flows in each direction -- here, we have an integral, done on the boundary of an infinitesimal region $V$ (this symbol will also represent the volume of the region):

$$ \frac{{\partial T}}{{\partial t}} = \lim_{V \to 0} \frac{\alpha }{V}\int_{\partial V} {\hat u\,dS \cdot \vec \nabla T} $$
At this point, one may apply the divergence theorem, converting this to:

$$\frac{{\partial T}}{{\partial t}} = \mathop {\lim }\limits_{V \to 0} \frac{\alpha }{V}\int\limits_V {\vec \nabla  \cdot \vec \nabla T\;dV}  = \alpha{\left| {\vec \nabla } \right|^2}T$$
In this sense, the divergence theorem is analogous to the fundamental theorem of calculus for manifolds with boundaries that are more than one-dimensional (see the bottom of the page for a link to a formalisation/an abstraction based on this analogy). But there are more ways to intuitively understand this. Note how the Laplacian is the trace of the Hessian matrix (note: we use $\vec{\nabla}^2$ to refer to the Hessian and $\left|\vec\nabla\right|^2$ to refer to the Laplacian):

$${\left| {\vec \nabla } \right|^2}T = {\mathop{\rm tr}} \left({\vec{\nabla} ^2}T\right)$$
The trace of a matrix is fundamentally linked to some notion of averaging -- the simplest interpretation of this is that it is the mean of the eigenvalues. But more relevant to our situation, it can be shown that the trace of a matrix is the expected value of the quadratic form defined by the matrix on the unit sphere -- or on a general sphere $S$:

$${\mathop{\rm tr}} A = \frac{1}{S}\int_S {\frac{{\Delta {x^T}A\,\Delta x}}{{\Delta {x^T}\Delta x}}\,dS} $$
One may check that taking the limit as $\Delta x \to 0$, substituting $\nabla^2$ for the operator and writing ${\overrightarrow \nabla ^2}f\,d\vec x = \overrightarrow \nabla  f$, one gets the original "average of directional derivatives" expression.

Can you interpret the other coefficients of the characteristic polynomial in terms of statistical ideas?

---

Further reading:

• Using the "infinitesimal region" idea to define divergence, curl and Laplacian rigorously: Khan Academy
• An abstraction based on the "analogy" between FTC, Divergence Theorem, Navier-Stokes Theorem, etc. Stokes' theorem (Wikipedia)

SVD, polar decomposition, normal matrices; a re-look at transposes and FTLA

Back in Null, row spaces, transpose, fundamental theorem of algebra, we first introduced some hand-wavy intuition for the transpose and the orthogonality of the row space and the null space (and the following fundamental theorem of linear algebra). Here, we solidify this intuition a bit more clearly.

Consider the "symmetric collapse" discussed in the above article. Our study of the transformation relied specifically on looking at it in a specific basis -- an orthogonal basis -- comprised of the column space and the null space. In this basis, the transformation is a scaling on both axes. In the more general case of an asymmetric collapse -- in which we rotated our space before collapsing, we looked at a basis formed by the row space and the null space -- the basis got rotated and scaled into the new basis, that was the column space and an arbitrary other vector (that could be perpendicular to the column space).

A sensible question to ask is if any transformation can be written in this fashion -- as a transformation of an orthogonal basis into another orthogonal basis. Analogous to how an eigenvalue decomposition of a matrix writes it as scalings in some basis, we're looking to represent the matrix as a spiral (i.e. a scaling combined with a rotation) in some basis. But let's stick with the first formulation of our question -- for any linear transformation $A: \mathbb{R}^n \to \mathbb{R}^m$, can we find an orthogonal basis on $\mathbb{R}^n$ that is mapped to an orthogonal basis $\mathbb{R}^m$?

One could, e.g. consider the images under $A$ of the angle of each orthonormal basis in $\mathbb{R}^2$ (i.e. look at the function $AU(\theta)\vec{e}_1 \cdot AU(\theta)\vec {e}_2$ for varying $\theta$ where $U(\theta)$ is the rotation matrix by angle $\theta$) and apply the intermediate value theorem, etc. And such a proof could in principle be extended to $\mathbb{R}^n$.

(See here for a thorough explanation.)

Here's another, more insightful way you might come to prove this -- we've been visualising linear transformations so far by looking at the image of the basis vectors, but another way to do visualise these transformations is by looking at the image of the unit circle under the transformation.

Why does this make sense? Well, one can find an ellipse passing through any two vectors centered at the origin. Elaborate on this argument. Is the resulting ellipse unique? (Hint: no, unless you mark points on the circumference)

Specifically, consider the axes of the image ellipse $\sigma_1 u_1$, $\sigma_2 u_2$ where $u_1$, $u_2$ are unit vectors. In the original unit circle, any pair of orthogonal vectors on the circle can be axes, so consider the pre-image of the axes of the image ellipse $v_1$, $v_2$.  So we have:
$$ A v_1 = \sigma_1 u_1\\
A v_2 = \sigma_2 u_2 $$
Or in general:
$$AV = U\Sigma\\
A = U\Sigma V^*$$
Where $\Sigma$ is diagonal and positive-definite, while $U$ and $V$ are orthogonal/unitary. This is called the Singular-Value Decomposition (SVD) of $A$.

In a sense, one can view this as an alternative to the eigen-decomposition. In the eigendecomposition, one looks for a single basis in which the transformation is a scaling. In the singular value decomposition, one looks at scaling and then "re-interpreting" in another basis, but requires that the bases be orthogonal, and that the diagonal matrix be positive and real-valued.

The entries $\Sigma$ are called the singular values of $A$, the columns of $U$ and $V$ respectively are the left-singular vectors and the right-singular vectors of $A$ respectively.

(Exercise: You know that $\Sigma$ is the scaling of the orthogonal basis, i.e. of the right-singular basis. Convince yourself that the rotation of the basis is given by $UV^*$.)

This gives us a much better intuition for the transpose. The SVD of $A^*$ is clearly:

$$A^* = V\Sigma U^*$$
I.e. the transpose has precisely the opposite rotational effect as $A$ and the same scaling. This is as opposed to the inverse matrix, which has both the opposite rotational and scaling effect as the matrix. For a rotation matrix (or generally an orthogonal matrix $A^*A=I$), the transpose equals the inverse, analogous to how the conjugate equals the inverse for a unit complex number $\bar{z}z=1$. A Hermitian matrix, $A=A^*$, by contrast, is one for which is irrotational, $UV^*=1$, i.e. for which the SVD equals the eigendecomposition.

Use the SVD to get some intuition for transpose identities like $(AB)^*=B^*A^*$

---

It's instructive to consider the SVD in the case of our original motivating example -- an asymmetric matrix representing a collapse of $\mathbb{R}^2$ into a line. What are the singular bases of this transformation? Well, it maps the orthogonal basis formed by row space and the null space into the orthogonal basis formed by the column space and the left-nullspace.

Think about its transpose.

Arranging the singular values from largest to smallest, we then have the following relation between the SVD and the items in the fundamental theorem of linear algebra, where $n$ is the dimension of the domain, $m$ is the dimension of the codomain, and $r$ is the dimension of the image/column space:

• The last $n - r$ singular values are zeroes. 
• The first $r$ singular values are positive.
• The last $n - r$ right-singular vectors span the null space.
• The first $r$ right-singular vectors span the row space.
• The last $m-r$ left-singular vectors span the left-null space.
• The first $r$ left-singular vectors span the column space.

Note that the terms kernel, coimage, cokernel and image are also used for the null space, row space, left-null space and column space, sometimes in a more general setting.

This is the full form of the Fundamental Theorem of Linear Algebra, which generalizes the rank-nullity theorem, row rank equals column rank and row space being perpendicular to null space.

Spend some time thinking about the SVD of non-square matrices, relating them to square collapse matrices. Think about their transposes.

Show that the right-singular vectors $V$ of a matrix $A$ are given by the eigenvectors of $A^*A$ (hint: start by considering the two-dimensional case, relating the right-singular vectors to a maximisation/minimisation problem, and extend the idea to more dimensions).

From this, it is clear that the left-singular eigenvectors $U$ (which are the right-singular eigenvectors of $A^*$) are given by the eigenvectors of $AA^*$ and the singular values are the square roots of the singular values of $A^*A$. Well, some of them are (which ones?).

We observed earlier that the rotational effect of the matrix $A = U \Sigma V^*$ can be given by $UV^*$. The scaling effect is given by $\Sigma$ on the basis of $V$. Hence we can write:

$$A = (UV^*)(V\Sigma V^*)$$
Letting $W=UV^*$ and $R = V\Sigma V^*$, this gives us a representation of $A$ as:

$$A=WR$$
Where $W$ is orthogonal and $R$ is positive-semidefinite. This is known as the right-polar decomposition of $A$. Analogously, one may consider the left-polar decomposition:

$$A = (U\Sigma U^*)(UV^*) = R'W$$

Interpret the above decomposition like we did the right decomposition. Note how $R' = WRW^*$, and how a right-polar decomposition leads to a left-polar decomposition of the transpose, and vice versa.

When is the polar decomposition unique? Compare the situation to the polar decomposition of complex numbers.

Here's a question: when is $R=R'$? One way of putting it is that $WRW^* = R$, i.e. $R$ commutes with $W$. This is a bit difficult to work with. Instead, one may show that the matrices $R$ and $R'$ can be given by $(A^*A)^{1/2}$ and $(AA^*)^{1/2}$ respectively (prove it!) -- noting that a principal square root can uniquely be defined for positive semidefinite matrices. So $R=R'\Leftrightarrow A^*A=AA^*$. These are known as normal matrices. Below is an exercise that provides more intuition into the behaviour of normal matrices.

---

What does it mean for a matrix to commute with its transpose? Earlier, when discussing commuting matrices, we referred to it as "matrices that do not disturb each other" -- specifically, they preserve each other's (generalised) eigenspaces. In the case of commuting with its transpose, it's easy to show that this means having the same eigenvectors (prove this!).

Here's a fact about the eigenvectors of the Hermitian transpose: the eigenvector of $A$ corresponding to the eigenvalue $\lambda$ is orthogonal to all eigenvectors of $A^*$ corresponding to any eigenvalue other than $\lambda^*$ (prove this!).

From these two facts, it follows that the following are equivalent (write down the proofs clearly!):

• $A$ is normal.
• $A$ commutes with $A^*$.
• $R$ commutes with $W$.
• $A$ is unitarily diagonalisable.

This is known as the spectral theorem.

---

(Anyone ever thought about how weird the word "normal" is? Sometimes, it means "perpendicular", sometimes -- as in "orthonormal" and "normalisation" -- it means "unit length", because "norm". What does it mean in this context? It's probably just referring to its eigenvectors being normal/orthogonal, but I like to think it's referring to the fact that the two alternative Hermitian-valued "norms" of the matrix, $A^*A$ and $AA^*$, are equal, so the matrix has a single "norm".)

Fundamentally, normal matrices are "analogous" to complex numbers, or they "generalise" complex numbers, in the sense that each of its eigenvalues acts as a complex number, transforming, acting as spirals on each of the orthogonal eigenvectors within its own copy of $\mathbb{C}$. One may construct the following table of analogies between normal matrices and complex numbers:

Complex numbers	Normal matrices
Zero (sorta)	Singular
Non-zero	Invertible
Real	Hermitian
Positive real	Positive-definite
Nonnegative real	Positive-semidefinite
Imaginary	Anti-Hermitian
Unit	Unitary
Conjugate	Hermitian transpose
Norm-squared	Gram matrix $A^*A$
Magnitude	$(A^*A)^{1/2}=R=V\Sigma V^*$
Argument	$AR^{-1}=W=UV^*$
Real Part	$\frac12(A+A^*)$
Imaginary Part times $i$	$\frac12(A-A^*)$

Exercise: an EP matrix or range-Hermitian matrix is a weakened version of a Hermitian matrix -- the row space of the matrix equals the column space. Although this was a bit hard to understand in the first article and was only briefly mentioned towards the end, we now have the intuition to comprehend them. Explain why a matrix is range-Hermitian if and only if it is unitarily similar to a matrix of the form of a block matrix:

$$\left[ {\begin{array}{*{20}{c}}C&0\\0&0\end{array}} \right]$$
Where $C$ is a non-singular square matrix and the zeroes are zero block matrices. This decomposition is called the core-nilpotent decomposition. Hence, show that being range-Hermitian is a weakened form of being normal.

Invariant and generalised eigenspaces; Jordan normal form

We defined "eigenvectors" -- or really "eigenlines" -- in order to understand the behaviour of linear transformations as scalings across certain axes (which may be complex, and the scalings may be complex too). But simply thinking of eigenlines as 1-dimensional spaces that a transformation leaves invariant (the fancy phrase here is: "the 1-dimensional subspaces on which it is an endomorphism"), it is natural to wonder about higher-dimensional invariant spaces -- subspaces on which some transformation $A$ acts as an endomorphism.

The problem is that any transformation has all sorts of useless invariant subspaces -- for instance any transformation $A:F^n\to F^m$ (where $m \le n$) has the entirety of $F^n$ as an invariant subspace (for any $F$), and rotations -- although fully described by their eigenvectors and eigenvalues -- have a bunch of real and complex planes as unnecessary invariant subspaces. And if $A$ has an eigenvalue with geometric multiplicity $>1$, there are an infinite number of useless invariant subspaces.

Specifically, if the goal is to find useful representations of defective matrices (non-diagonalisable), invariant subspaces seem completely useless -- they certainly have no hope of giving us any sort of unique representation. Perhaps more on the point, our "eigenlines" have corresponding eigenvalues that tell us how the transformation behaves within an eigenline. Our invariant subspaces currently have nothing of the sort -- the transformation can have any sort of behaviour on the invariant subspace -- rotation, skewing/scaling, shearing, skewering -- and we'd have no idea. We need a convenient way to write down the behaviour of the transformation on an invariant subspace.

Here's something we can start to think about: ordinary eigenvectors satisfy $(A-\lambda I)v=0$, which gives us a one-dimensional solution space. In analogy with solutions to linear differential equations (linear homogenous if you use the conventional terminology, but I reserve "affine" for linear non-homogenous), an equation like

$$(A-\lambda_1 I)(A-\lambda_2 I)v=0$$
(where ${\lambda _1},{\lambda _2}$ are both eigenvalues of $A$) would have a 2-dimensional solution space, etc.

Note that when $\lambda_1, \lambda_2$ are not eigenvalues, we don't have a 2-dimensional solution space (what does the solution space look like then?). Why does it work with differential equations for any $\lambda_1, \lambda_2$? (hint: what do the eigenvalues look like?)

It's sensible to ask: are these solution spaces the same as our invariant subspaces? I.e. is every member of a $k$-dimensional invariant subspace a solution to an equation of the form

$$(A - {\lambda _1}I)...(A - {\lambda _k}I)v = 0$$
for eigenvalues $\lambda_1,...\lambda_k$?

The answer is yes. I encourage you to try and prove it for yourself -- it is instructive to first consider special cases: (i) rotation in a plane, where indeed $(A-iI)(A+iI)=A^2+1=0$ is the minimal polynomial of $A$ (ii) more generally, $F^n$ is an invariant subspace for all isomorphisms, and indeed for all $v$ in this subspace (i.e. $\forall v \in F^n$), $p(A)v=0$ where $p$ is the characteristic polynomial of $A$ by the Cayley-Hamilton theorem.

The key to proving that every invariant subspace is given by solutions to an equation of the form $(A - {\lambda _1}I)...(A - {\lambda _k}I)v = 0$ (and vice versa) lies in recognising that on any $k$-dimensional invariant subspace, $A$ is acts as an endomorphism, and therefore the Cayley-Hamilton theorem applies to it, with a $k$-order characteristic polynomial.

I encourage you to spend some time thinking about this -- try relating it to differential equations. Come up with another proof of the statement -- an inductive one. See if this results in a better intuition for the Cayley-Hamilton theorem.

Is it true that there $2^n$ invariant subspaces of any transformation on an $n$-dimensional linear space? What about the identity transformation?

---

Now, let's discard the invariant subspaces we don't want. We already know how to handle cases with distinct eigenvalues -- i.e. we have distinct eigenvalues in $\lambda_1...\lambda_k$ -- we just get an eigenvector for each eigenvalue. So we're really just concerned with subspaces of the form ${(A - \lambda I)^k}v = 0$. This is analogous to linear differential equations with repeated roots being weirder than ones with distinct roots.

Note that we still know how to handle ${(A - \lambda I)^k}v = 0$-like equations when the algebraic multiplicity is accounted for by geometric multiplicity -- when this is the case, you can reduce the power from $k$ (by subtracting from it the geometric multiplicity). This nuance doesn't exist with differential equations, because distinct eigenvectors have distinct eigenvalues.

Vectors satisfying such an equation are called generalised eigenvectors of order $k$ where $k$ is the minimum value for which it does satisfy the equation, and the invariant subspaces formed by generalised eigenvectors of the same eigenvalue are called generalised eigenspaces. The dimension of the generalised eigenspace always equals the algebraic multiplicity, unlike the eigenspace, whose dimension equals the geometric multiplicity.

Check this, and that $k$ is the difference between algebraic and geometric multiplicity.

What kind of transformations precisely do generalised eigenvectors with degree greater than 1 correspond to? Clearly, skews and rotations are out of the question. But some insight can be gained from looking at the nature of skews and rotations on a plane.

In two dimensions, a characteristic polynomial with a positive discriminant yields a skew along some axis, a negative discriminant yields a rotation, and the case we're interested -- the presence of repeated roots -- corresponds to the point "between" skews and rotations (speaking hand-wavily), shears.

In a more general setting, if one has: ${(A - \lambda I)^k}v_k = 0$  for generalised eigenvector  $v_k$of degree $k$, then one can extract generalised eigenvectors of each degree lower:

$$\begin{array}{l}{(A - \lambda )^i}{(A - \lambda )^{k - i}}{v_k} = 0
\\ \Rightarrow {v_i} = {(A - \lambda I)^{k - i}}{v_k}\end{array}$$
Implying that the generalised eigenvectors with the same eigenvalue (can) form a basis for the corresponding generalised eigenspace.

$$\begin{array}{*{20}{r}}{(A - \lambda I){v_1} = 0 \Leftrightarrow A{v_1} = \lambda {v_1}}\\{{{(A - \lambda I)}^2}{v_2} = 0 \Leftarrow (A - \lambda I){v_2} = {c_1}{v_1} \Leftrightarrow A{v_2} = {c_1}{v_1} + \lambda {v_2}}\\{{{(A - \lambda I)}^3}{v_3} = 0 \Leftarrow (A - \lambda I){v_3} = {c_2}{v_2} \Leftrightarrow A{v_3} = {c_2}{v_2} + \lambda {v_3}}\\{ \vdots \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,}\\{{{(A - \lambda I)}^k}{v_k} = 0 \Leftarrow (A - \lambda I){v_k} = {c_{k - 1}}{v_{k - 1}} \Leftrightarrow A{v_k} = {c_{k - 1}}{v_{k - 1}} + \lambda {v_k}}\end{array}$$
This gives us a very clear picture of how these general "sheary" transformations look in the basis of the generalised eigenvectors -- there is a shear in each plane $\left\langle {{v_{k - 1}},{v_k}} \right\rangle$. Suitable scalings could of course be chosen to make all the $c_i=1$.

It's hard to overstate the significance of this -- what we've just found is that any defective matrix can be decomposed into shears on each of its generalised eigenspaces. This completely classifies the diagonalisability of matrices. If it's defective, it's a shear.

Draw out some of these transformations in three or more dimensions.

Notice the directions of the implication signs -- can we make them double-sided? What if we have some geometric multiplicity? How would our shears look like then? How many dimensions do you need to visualise this?

Think about why this characterisation of defective matrices makes sense. What effect does adding a 1 to the subdiagonal have on the determinant? Why? (hint: area of a parallelogram) What about the other coefficients of the characteristic polynomial? (hint: think of these in terms of traces of some matrix). So these matrices are precisely those which have the same characteristic polynomial as a diagonalisable matrix without actually being similar to one.

---

Clearly, the generalised eigenspaces of a transformation are pairwise disjoint (i.e. intersect only at the origin). Since all eigenvalues are being considered, their union is all of $F^n$. Thus the union of their bases forms a basis for $F^n$. This gives us a representation of the transformation $A$ in this basis.

From the last section, it is clear that the effective transformation on a $k$-dimensional eigenspace with eigenvalue $\lambda$ (called a "Jordan block") is given by:

$$\left[ {\begin{array}{*{20}{c}}\lambda &0&0& \cdots &0\\1&\lambda &0& \cdots &0\\0&1&\lambda & \cdots &0\\ \vdots & \vdots & \vdots & \ddots & \vdots \\0&0&0& \cdots &\lambda \end{array}} \right]$$
The Jordan normal form of a transformation is then the matrix formed by putting all Jordan normal blocks along the diagonal -- i.e. the representation of $A$ in its generalised eigenbasis.

Some writers define the Jordan normal form with the 1's on the superdiagonal. It should be clear to you from the work we've done that this form is obtained by taking the basis vectors in reverse order (i.e. changing the basis to $\langle e_n,...,e_1 \rangle$. If this isn't clear to you, go back and work through the previous section once more. Loop.

(Stay tuned for the next article to see how we can -- instead of defining generalised eigenspaces whose dimension is defined by the algebraic, rather than geometric multiplicity -- "force" algebraic multiplicity to equal geometric multiplicity, i.e. diagonalise any matrix -- with a ring extension.)

All matrices can be diagonalised over R[X]/(X^n)

This post follows from my answer to the math stackexchange question What kind of matrices are non-diagonalisable?

---

Non-diagonalisable 2 by 2 matrices can be diagonalised over the dual numbers -- and the "weird cases" like the Galilean transformation are not fundamentally different from the nilpotent matrices.

The intuition here is that the Galilean transformation is sort of a "boundary case" between real-diagonalisability (skews) and complex-diagonalisability (rotations) (which you can sort of think in terms of discriminants). In the case of the Galilean transformation $\left[\begin{array}{*{20}{c}}{1}&{v}\\{0}&{1}\end{array}\right]$, it's a small perturbation away from being diagonalisable, i.e. it sort of has "repeated eigenvectors" (you can visualise this with MatVis). So one may imagine that the two eigenvectors are only an "epsilon" away, where $\varepsilon$ is the unit dual satisfying $\varepsilon^2=0$ (called the "soul"). Indeed, its characteristic polynomial is:

$$(\lambda-1)^2=0$$
Whose solutions among the dual numbers are $\lambda=1+k\varepsilon$ for real $k$. So one may "diagonalise" the Galilean transformation over the dual numbers as e.g.:

$$\left[\begin{array}{*{20}{c}}{1}&{0}\\{0}&{1+v\varepsilon}\end{array}\right]$$
Granted this is not unique, this is formed from the change-of-basis matrix $\left[\begin{array}{*{20}{c}}{1}&{1}\\{0}&{\epsilon}\end{array}\right]$, but any vector of the form $(1,k\varepsilon)$ is a valid eigenvector. You could, if you like, consider this a canonical or "principal value" of the diagonalisation, and in general each diagonalisation corresponds to a limit you can take of real/complex-diagonalisable transformations. Another way of thinking about this is that there is an entire eigenspace spanned by $(1,0)$ and $(1,\varepsilon)$ in that little gap of multiplicity. In this sense, the geometric multiplicity is forced to be equal to the algebraic multiplicity*.

Then a nilpotent matrix with characteristic polynomial $\lambda^2=0$ has solutions $\lambda=k\varepsilon$, and is simply diagonalised as:

$$\left[\begin{array}{*{20}{c}}{0}&{0}\\{0}&{\varepsilon}\end{array}\right]$$
(Think about this.) Indeed, the resulting matrix has minimal polynomial $\lambda^2=0$, and the eigenvectors are as before.

---

What about higher dimensional matrices? Consider:

$$\left[ {\begin{array}{*{20}{c}}0&v&0\\0&0&w\\0&0&0\end{array}} \right]$$
This is a nilpotent matrix $A$ satisfying $A^3=0$ (but not $A^2=0$). The characteristic polynomial is $\lambda^3=0$. Although $\varepsilon$ might seem like a sensible choice, it doesn't really do the trick -- if you try a diagonalisation of the form $\mathrm{diag}(0,v\varepsilon,w\varepsilon)$, it has minimal polynomial $A^2=0$, which is wrong. Indeed, you won't be able to find three linearly independent eigenvectors to diagonalise the matrix this way -- they'll all take the form $(a+b\varepsilon,0,0)$.

Instead, you need to consider a generalisation of the dual numbers, sometimes called (in computing mathematics and non-standard analysis) the "hyperdual numbers", with the soul satisfying $\epsilon^n=0$. Then the diagonalisation takes for instance the form:

$$\left[ {\begin{array}{*{20}{c}}0&0&0\\0&{v\epsilon}&0\\0&0&{w\epsilon}\end{array}} \right]$$

---

*Over the reals and complexes, when one defines algebraic multiplicity (as "the multiplicity of the corresponding factor in the characteristic polynomial"), there is a single eigenvalue corresponding to that factor. This is of course no longer true over the hyperdual numbers, because they are not a field, and $ab=0$ no longer implies "$a=0$ or $b=0$".

In general, if you want to prove things about these numbers, the way to formalise them is by constructing them as the quotient $\mathbb{R}[X]/(X^n)$, so you actually have something clear to work with.

(Perhaps relevant: Grassmann numbers as eigenvalues of nilpotent operators -- the Hyperdual numbers are not the same as the Grassmann numbers, and the algebra of the Grassmann numbers is definitely different from that of nilpotent and shear matrices, but go see if you can make sense of it.)

Something important to note is that the diagonalisation is not of the form $D=P^{-1}AP$, as the eigenvector matrices are not invertible. However, it is still true that $PD=AP$ -- nonetheless, this limitation prevents this formalism for being any good for e.g. dealing with polynomial-ish differential equations with repeated roots, for instance, as far as I can see. The infinitesimal-perturbation/"take a limit" approach we talked about in Limiting Cases II: repeated roots of a differential equation are still the right approach for that.

Intuition, analogies and abstraction

$$-1=\sqrt{-1}\sqrt{-1}=\sqrt{(-1)(-1)}=\sqrt{1}=1$$
I bet you've seen the fake "proof" above that minus one and one are equal. And the standard explanation as to why it's wrong is that the statement $\sqrt{ab}=\sqrt{a}\sqrt{b}$ only applies when $\sqrt{a}$ and $\sqrt{b}$ are real, or something like that (maybe only one of them needs to be real -- something like that -- who cares?).

But if you're like me, that isn't a very satisfactory proof. Why does the identity not hold for complex numbers? For that matter, why does it hold for real numbers? Well, that is a good question, and one way of answering it would be to try and prove the identity for real numbers, and see what properties of the real numbers (or of the real square root, in particular) you use. And if this article were being filed under "MAR1104: Introduction to formal mathematics", that's how I might explain things -- but that doesn't give us too much insight -- not about square roots and complex numbers, anyway.

Let's think about what $\sqrt{ab}=\sqrt{a}\sqrt{b}$ means.

What does the square root of a real number mean, anyway? It's some property related to multiplying a real number by itself. What does multiplication mean? What does a real number mean? The picture I have in my head of the real numbers is of a line. But what exactly is this line? -- the real numbers are just a set. Why did you put them on this line in this specific way? In doing so, you gave the real numbers a structure, a specific type of structure called an "order", defined by the operation $<$.

But there are other ways to think about/structure the real numbers. One way is to think of real numbers as (one-dimensional) scalings. You can scale things like mass, and volume, using real numbers, representing the scalings as real numbers. Scaling a mass by 2 is equivalent to multiplication by 2. So this gives the real numbers a multiplicative structure, defined by the operation $\times$ (or whatever notation -- or lack thereof -- you prefer). And the "real line" then just represents the image of "1" under all scalings.

So the way to think about square roots is to think of numbers as linear transformations called scalings, and think about the scaling that when done twice, gives you the number you're taking the square root of. So what's $\sqrt{-1}$? What's $-1$? $-1$, multiplicative, is a reflection. What's its square root? Try to think of a (linear!) transformation that when done twice gives you a reflection. It can't be done in one dimension. And can you think of another such transformation? Can you prove these are the only two? Are you sure -- what about if you add a dimension?

So the natural way to think about square roots of numbers that may or may not be complex, is with so-called "Argand diagrams", on the complex plane, the image of "1" under all complex numbers multiplicative.

Click "edit graph" to play with a and b!

To simplify things, consider only unit complex numbers (this is okay, because all complex numbers can be written as a real multiple of a unit complex number and a real number). The product of complex numbers $a$ and $b$ involves rotating by $a$, then rotating by $b$. The square roots of $a$ and $b$ involve going halfway around the circle as $a$ and $b$, and the square root of $ab$ goes halfway around the circle as $ab$.

So it seems like the identity should hold, doesn't it? $\sqrt{ab}$ goes half as much as $a$ and $b$ put together -- this seems to be exactly what $\sqrt{a}\sqrt{b}$ does -- go around half as much as $a$, then half as much as $b$. Isn't $\frac{\theta+\phi}2=\frac{\theta}2+\frac{\phi}2$?

The problem is that $\sqrt{ab}$ doesn't really go $\frac{\theta+\phi}2$ around the circle, if $\theta+\phi$ is greater than $2\pi$. You can see this in the diagram courtesy of Desmos above -- $ab$ has gone a full circle, and its square root is defined to halve the argument of $ab$, but the argument isn't $\arg (ab)=\arg (a) + \arg (b)$, rather:

$$\arg (ab) \equiv \arg (a) + \arg (b) \pmod{2\pi}$$
But halving is not an operation that the $\bmod$ equivalence relation respects -- not in general, anyway. It is not true that

$$\arg (ab)/2 \equiv (\arg (a) + \arg (b))/2 \pmod{2\pi}$$
Instead:

$$\arg (ab)/2 \equiv (\arg (a) + \arg (b))/2 \pmod{\pi}$$
Let's recall from basic number theory -- on integers, the general result regarding multiplication on mods. If $a\equiv b\pmod{m}$, then $na\equiv nb \pmod{nm}$, certainly, and also $na\equiv nb \pmod{m}$ iff $n$ is an integer*. But $1/2$ isn't an integer, which is why only the former result is relevant.

This is also why $(ab)^2=a^2b^2$ does hold for complex numbers.

*when $n$ isn't an integer, we need $na$, $nb$ to be integers for the statement to even be well-defined in standard number theory, and then you have a result for division on mods involving $\gcd(d,m)$, etc. This isn't a concern for us here because we're dealing with divisibility over the reals -- if you want to be formal, a real number is divisible by another real number if the former can be written as an integer multiple of the latter.

So there you have it -- I just demonstrated a very fundamental analogy between two seemingly incredibly unrelated ideas: complex numbers modular arithmetic -- square roots of complex numbers don't multiply naturally, because mod doesn't respect division. It's almost as if somehow, somewhere, somehow magically, exactly the same kind of math was used to derive results, to prove things, about these unrelated objects.

As if they're just two instances of the same thing.

I wonder what that thing could be.

---

Let's talk about something completely unrelated (no, genuinely -- completely unrelated -- I won't tell you this is an instance of the "same thing" too). Let's talk about logical operators, specifically: do $\forall$ and $\exists$ commute? I.e. is $\forall t, \exists s, P(s,t)$ equivalent to $\exists s, \forall t, P(s,t)$?

You just need to read the statements aloud to realise they don't. To use a classical example, "all men have wives" and "there is a woman who is the wife of all men" are two very different statements (okay, in this case both statements are false, so they're equivalent in that sense, so you get my point).

But let's think more deeply about why they don't commute. What do $\forall t, \exists s, P(s,t)$ and $\exists s, \forall t, P(s,t)$ mean, anyway? $\forall$ and $\exists$ are just infinite $\land$ and $\lor$ statements , i.e. $\forall t$ is just an $\land$ statement ranging over all possible values that $t$ can take and $\exists s$ is just an $\lor$ statement ranging over all possible values $s$ can take.

So $\forall t, \exists s, P_{st}$ just means (letting $s$ and $t$ be natural numbers for simplicity, but they don't have to):

$$({P_{11}} \lor {P_{21}} \lor ...) \land ({P_{12}} \lor {P_{22}} \lor ...) \land ...$$
And $\exists s, \forall t, P(s,t)$ means:

$$({P_{11}} \land {P_{12}} \land ...) \lor ({P_{21}} \land {P_{22}} \land ...) \lor ...$$
This is a bit complicated, so let's instead look at the simpler case where you have only 2 by 2 statements -- i.e. just construct the analogy between $\forall,\exists$ and actual $\land,\lor$ statements.

So the question is if:

$$({P_{11}} \lor {P_{21}}) \land ({P_{12}} \lor {P_{22}}) \Leftrightarrow ({P_{11}} \land {P_{12}}) \lor ({P_{21}} \lor {P_{22}})$$
This is interesting. Maybe you see where this is going. Let me just do a notation change -- I'll use "$\times$" for $\land$, "$+$" for $\lor$, "$=$" for $\Leftrightarrow$" and some new letters for the propositions. Under this new notation, where $\times$ is invisible as always, we're asking if:

$$(a + b)(c + d) = ac + bd$$

Aha! This is Freshman's dream, isn't it? And we know it's not true -- it's a dream, after all, don't be delusional -- and we know why it's not true too.

But wait -- we aren't talking about elementary algebra here. I just gave you some silly notation and made it look like Freshman's dream. But here's the thing: the proof (or algebraic proof -- a counter-example is also a proof, but that isn't so interesting... not here, anyway) that these propositions aren't equivalent is exactly the same as in algebra. We expand out the brackets (because we know that $\land$ distributes over $\lor$ -- we also know that $\lor$ distributes over $\land$, incidentally, something that is not true in standard algebra) and point out that there are extra terms, and point out that these extra terms change the value of the expression (they aren't zero).

So there's some kind of relationship between the boolean algebra and an elementary algebra. A lot of proofs that can be done in one of these algebras can be written almost identically in the other. Not all these proofs, mind you -- then the algebras would just be isomorphic to each other -- but some of them can. Maybe a lot of important ones can.

An abstraction that produces such proofs simultaneously for both elementary algebra and boolean algebra may be more complicated than you think -- there's no real sense in which a statement is "always zero" in boolean algebra. Take for instance, distributivity of $\lor$ over $\land$ -- $a+bc=(a+b)(a+c)$. This is not true in elementary algebra, because the extra term $ab+ac$ is not always equal to zero ($a^2\ne a$ is not really an example, because $a^2=a$ for $a\in\{0,1\}$ -- but $a(b+c)=0$ is not true for all $a,b,c\in\{0,1\}$). It's just that it leaves the value of the existing terms unchanged in this specific instance.

---

I've just illustrated two examples here -- the first one is a type of group, by the way, but you've probably seen dozens of other such "connections between different areas of mathematics" yourself. I've made these sorts of analogies fundamental to a lot of the articles I've written here (I think). You might've just thought of them as interesting insights, but in reality, abstract mathematics/abstract algebra -- or really just mathematics in general -- is all about these analogies.

In a sense, mathematics is largely about abstraction. I mean, that's not what mathematics fundamentally is -- fundamentally, math is just logic -- but it's how mathematics largely functions. Whenever one talks of axioms, you could think of them as fundamental defining ideas of mathematical objects, and you can also think of them as "interfaces" between mathematics and reality (see my introduction to linear transformations). There are a massive number of different physical phenomena that we can study, and rather than prove everything from scratch for each one of them, it is much better -- and more insightful in terms of understanding the connections between things -- to show that they satisfy a certain set of axioms that apply to a whole range of things, and then deduce that all the logical consequences of these axioms -- all theorems -- are satisfied by the objects.

If we can do that with physical phenomena, we can sure as well do it with mathematical phenomena too -- instead of proving something from scratch for every new mathematical object, we prove that it is a group, or a ring, or a field, or a module, or an algebra, or a topology, or a geometry of some sort, by verifying it matches the axioms -- and then use all the abstract knowledge we have about these things and deduce they must necessarily apply to our new object, because they are logical consequences of our axioms.

Abstract mathematics is, in this sense, all about generalising things by finding the "smallest set of axioms" the thing requires.

(Well, not really -- the most general statement is "true", and everything else is just a logical deduction from this statement. So in that sense mathematics is all about finding special cases. But in order to know what to take a special case of, and what special case that "what" is of "true", you need to generalise.)

List some weird analogies you've seen before in math. Something about divisibility sound familiar?