\documentclass[a4paper]{amsart}
\usepackage{graphicx}
\usepackage{amsmath, amsthm}
\usepackage{amssymb}
\usepackage{amsfonts}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Setup

\newtheorem{thm}{Theorem}[section]
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{defn}{Definition}[section]
\newtheorem{example}{Example}[section]
\newtheorem{exercise}{Exercise}[section]
\newtheorem{algorithm}{Algorithm}[section]


\newcommand\be{\begin{equation}}
\newcommand\ee{\end{equation}}
\newcommand\half{\frac{1}{2}}

\newcommand\CC{{\mathbb C}}
\newcommand\EE{{\mathbb E}}
\newcommand\NN{{\mathbb N}}
\newcommand\OO{{\cal O}}
\newcommand\PP{{\mathbb P}}
\newcommand\PPbar{{\overline{\mathbb P}}}
\newcommand\RR{{\mathbb R}}
\newcommand\TT{{\mathbb T}}
\newcommand\ZZ{{\mathbb Z}}

\renewcommand{\theenumi}{\roman{enumi}}    

\DeclareMathOperator{\conv}{co}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\dist}{dist}
\DeclareMathOperator{\GL}{GL}
\DeclareMathOperator{\Span}{Span}
\DeclareMathOperator{\var}{var}
\DeclareMathOperator{\Vol}{Vol}

\newcommand\fhat{{\widehat{f}}}
\newcommand\ghat{{\widehat{g}}}
\newcommand\Lhat{{\widehat{L}}}
\newcommand\LAhat{{\widehat{L^A}}}
\newcommand\Lhhat{{\widehat{L^h}}}
\newcommand\phihat{{\widehat{\phi}}}
\newcommand\phiAhat{{\widehat{\phi^A}}}
\newcommand\phihhat{{\widehat{\phi^h}}}


\newcommand{\va}{{\bf a}}
\newcommand{\vb}{{\bf b}}
\newcommand{\vc}{{\bf c}}
\newcommand{\vd}{{\bf d}}
\newcommand{\vf}{{\bf f}}
\newcommand{\vh}{{\bf h}}
\newcommand{\ve}{{\bf e}}
\newcommand{\vg}{{\bf g}}
\newcommand{\vp}{{\bf p}}
\newcommand{\vq}{{\bf q}}
\newcommand{\vs}{{\bf s}}
\newcommand{\vu}{{\bf u}}
\newcommand{\vU}{{\bf U}}
\newcommand{\vUhat}{\widehat{\bf U}}
\newcommand{\vv}{{\bf v}}
\newcommand{\vw}{{\bf w}}
\newcommand{\vx}{{\bf x}}
\newcommand{\vy}{{\bf y}}
\newcommand{\vz}{{\bf z}}
\newcommand{\vzhat}{\widehat{\bf z}}
\newcommand{\vO}{{\bf 0}}

\newcommand{\vL}{{\bf L}}
\newcommand{\vT}{{\bf T}}

\newcommand{\Balpha}{{\mbox{\boldmath $\alpha$}}}
\newcommand{\Bbeta}{{\mbox{\boldmath $\beta$}}}
\newcommand{\Bgamma}{{\mbox{\boldmath $\gamma$}}}
\newcommand{\Bell}{{\mbox{\boldmath $\ell$}}}
\newcommand{\Blambda}{{\mbox{\boldmath $\lambda$}}}
\newcommand{\Bmu}{{\mbox{\boldmath $\mu$}}}

%
% Uncomment EXACTLY ONE of the next two commands - 
% Uncomment 1st line to avoid text, 2nd line to get text
%
\newcommand{\ignore}[1]{}
%\newcommand{\ignore}{}
%
%
%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


 

\begin{document}

\title[Miscellany]
{Miscellany}
\author[Brad Baxter]{Brad Baxter}


\maketitle

\tableofcontents

\newpage
\section{Introduction}

This is a short collection of miscellanies intended for teaching.

\bigskip
\noindent
{\bf Version: 201912181519}


\newpage
\section{Roots of Unity}

\begin{example}
Let $\omega = e^{2\pi i/3}$. Thus $1, \omega, \omega^2$ are the three
cube roots of unity. Then it is easily checked that $|1 - \omega| = |1
- \omega^2| = \sqrt{3}$, so that
\[
|1 - \omega| |1 - \omega^2| = 3.
\]
\end{example}

\begin{example}
Suppose we take the $4$ points $\pm 1$ and $\pm i$. Then
$|1-i|=\sqrt{2}$, so that
\[
|1-i| |1-(-1)| |1-(-i)| = 4.
\]
\end{example}

These examples lead to a conjecture:
\be
\prod_{k=1}^{n-1} |1 - \omega^k| = n,
\label{u1}
\ee
where $\omega = e^{2\pi i/n}$ and $n \ge 2$, and here is the Matlab
code to check this.

\bigskip
\begin{verbatim}
I=sqrt(-1);
n=5;omega=exp(2*pi*I/n); P=1; for k=1:n-1, P=P*abs(1-omega^k); end; P
\end{verbatim}
\bigskip

In fact, we shall see that a stronger
statement is true:

Let $n > 1$ be an integer and let $\omega=e^{2\pi i/n}$. Thus the
complex numbers $\{\omega^k : k = 0, 1, \ldots, n-1\}$ are the $n$th
roots of unity. Thus
\[
z^n - 1 = \left(z-1\right) \prod_{k=1}^{n-1} \left(z - \omega^k\right).
\]
Hence
\[
\prod_{k=1}^{n-1} \left(1 - \omega^k\right)
= \lim_{z \to 1} \frac{z^n - 1}{z-1} = n,
\]
by de L'H\^opital's rule.


\newpage
\section{The Length of the Day}

\subsection{The Length of the Day at the Solstices}
We shall compute the length of the day at the Summer solstice in the
northern hemisphere. The origin of our coordinate system will be at
the centre of the Earth, the $x$-axis will point directly towards the
Sun, and the $z$-axis will be perpendicular to the Earth's orbital
plane and will be directed into the northern hemisphere. We need the
following orthonormal vectors to describe the motion of a point on the
Earth's surface:
\be
\vu_1 = 
\left(
\begin{array}{c}
\cos \alpha\\
0\\
-\sin \alpha
\end{array}
\right),
\vu_2 = 
\left(
\begin{array}{c}
0\\
1\\
0\\
\end{array}
\right),
\hbox{ and }
\vu_3 = 
\left(
\begin{array}{c}
\sin \alpha\\
0\\
\cos \alpha
\end{array}
\right)
\label{s1},
\ee
where $\alpha = 23.5$ degrees approximately for the Earth.

The motion of a point at latitude $\theta$ is then
\be
\left(
\begin{array}{c}
x(t)\\
y(t)\\
z(t)
\end{array}
\right)
= \left(\vu_1 \cos t + \vu_2 \sin t\right) \cos \theta
+ \sin \theta \vu_3.
\label{s2}
\ee
In particular, we have
\begin{align}
x(t) &= \cos \theta \cos \alpha \cos t + \sin \theta \sin \alpha\nonumber\\
y(t) &= \cos \theta \sin t.
\label{s3}
\end{align}
At the Summer solstice, day corresponds to $x(t) > 0$. Thus, solving
$x(t) = 0$, we obtain the length of day as a function of the latitude $\theta$:
\be
L(\theta) = 2 \cos^{-1} \left(-\tan \theta \tan \alpha\right), \quad
|\theta| \le 90 - \alpha,
\label{s4}
\ee
and this gives the length of the day in {\em degrees}. Thus the length
of the day in hours is given by$L_h(\theta) = (24/360) L(\theta) =
(1/15) L(\theta)$, i.e.
\be
L_h(\theta) = \frac{2}{15} \cos^{-1} \left(-\tan \theta \tan \alpha\right), \quad
|\theta| \le 90 - \alpha,
\label{s5}
\ee
For $\theta \in (90 - \alpha, 90)$, $L_h(\theta) = 24$;
similarly $L(\theta) = 0$ for $\theta \in (-90, -90 + \alpha)$.

The following Matlab code generates the ratio of the longest day to
the shortest day.

\bigskip
\begin{verbatim}
alpha= 23.5*pi/180;
theta=0:pi/100:(pi/2) - alpha;
y = acos(-tan(alpha)*tan(theta));
R = y ./ (pi - y);
plot(theta,R)
plot((180/pi)*theta,R)
grid
\end{verbatim}
\bigskip

\begin{figure}[tb]
\vspace{6pt}
\includegraphics[scale=0.6, bb=0 200 600 600,clip]{d1.pdf}
\parbox{28.5pc}{\caption{Annual variation in day-length at 51 degrees North}\label{d1}}
%\end{center}
\end{figure}

\subsection{The Variation in the Length of the Day during the Year}
We solve the equation
\be
\left(
\begin{array}{c}
x(t)\\
y(t)
\end{array}
\right)^T
\left(
\begin{array}{c}
\cos u\\
\sin u
\end{array}
\right)
= 0,
\label{s6}
\ee
where $0 \le u \le 360$ measures orbital time in degrees, i.e. one
year corresponds to $360$ degrees. Expanding \eqref{s6}, we obtain
\be
\cos \theta \cos \alpha \cos u \cos t + \cos \theta \sin u \sin t
=
- \sin \theta \sin \alpha \cos u,
\label{s7}
\ee
or
\be
\cos \left(t - \beta\right)
= \frac{\sin \theta \sin \alpha \cos u}{\gamma},
\label{s8}
\ee
where
\be
\gamma^2 = \cos^2 \theta \cos^2 \alpha \cos^2 u + \cos^2 \theta \sin^2 u.
\label{s9}
\ee
Hence the sunrise and sunset times are given by
\be
t_{\pm} - \beta = \pm \cos^{-1} \left(\frac{-\sin\theta \sin\alpha
  \cos u}{\gamma}\right),
\label{s10}
\ee
and the length of the day is then
\be
t_{+} - t_{-} = 2\cos^{-1} \left(\frac{-\sin\theta \sin\alpha
  \cos u}{\gamma}\right),
\label{s11}
\ee

\bigskip
\begin{verbatim}
%
% Displays the yearly variation in the length of the day
% (in hours) at latitude theta, where |theta| < pi/2 - alpha.
%
alpha= 23.5*pi/180;
theta=51*pi/180;
u=0:pi/1000:2*pi;
A = -sin(theta)*sin(alpha)*cos(u);
B = cos(theta)*sqrt( (cos(alpha)^2)*(cos(u).^2) + (sin(u).^2) );
D = 2*acos(A ./ B)*12/pi;
plot(u,D)
%
% D is quite close to sinusoidal
%
%hold on
%plot(u, 12+(max(D)-12)*cos(u),'r')
%hold off
\end{verbatim}
\bigskip


\newpage
\section{Distance Seen and Height}

If we take the Earth to be a perfect sphere of radius $R$, then the
distance seen $D$ at height $H$ is given by
\be
D = R\theta
\label{h1}
\ee
where
\be
\left(R + H\right)\cos \theta = R.
\label{h2}
\ee
It's useful to nondimensionalize these by introducing
\be
d = \frac{D}{R} \quad\hbox{ and }\quad
h = \frac{H}{R}.
\label{h3}
\ee
Thus \eqref{h1} and \eqref{h2} become
\be
d = \theta \quad\hbox{ and }\quad \left(1 + h\right) \cos \theta = 1.
\label{h4}
\ee
Eliminating $\theta$ from \eqref{h4} we obtain
\be
\cos d = \left[\frac{1}{1 + h}\right].
\label{h5}
\ee
If $h$ is small, then $d$ must also be small, so we have
\be
1 - d^2/2 + \cdots = 1 - h + \cdots,
\label{h6}
\ee
or 
\be
d^2 = 2h.
\label{h7}
\ee
Returning to our original variables, we find
\be
D^2 = 2 H R.
\label{h8}
\ee

\begin{example}
Taking $R = 6.4 \times 10^6$ m and $H = 100$ m. Then $\sqrt{2HR}
\approx 36$ km.
\end{example}

\begin{example}
Here is some MATLAB code to illustrate the approximation's worth.
\begin{verbatim}
R=6.4e6; h = 0:100:100000;
dtrue = R*acos((1 + h/R).^(-1));
dapprox = (2*R*h).^(1/2);
\end{verbatim}
\end{example}


\newpage
\section{The Railway problem}

This is a very old chestnut indeed. We imagine a straight piece of
rail of unit length which, under thermal expansion, becomes a circular
arc of length $1+\delta$, where $1 \gg \delta > 0$. The rail will bow
upwards, attaining a maximum height $h$ at its centre, and the problem
is to determine $h$, which is surprisingly large. 

If we let $R$ denote the radius of the circular arc after expansion,
and $\theta$ denote the half-angle subtended at the centre of the
circle, then we have the equations
\begin{align}
2R\theta &= 1+\delta, \label{r1}\\
R - h = R \cos \theta, \label{r2}\\
\frac{1}{2} = R \sin \theta. \label{r3}
\end{align}
Of course
\be
h = R\left(1 - \cos \theta\right).
\label{r4}
\ee
Eliminating $R$ from \eqref{r1} and \eqref{r3}, we obtain
\[
\frac{\sin \theta}{2\theta} = \frac{1/2}{1+\delta}
\]
or
\be
\frac{\sin \theta}{\theta} = \frac{1}{1+\delta}.
\label{r5}
\ee
Now $0 < \delta \ll 1$ implies that $\theta$ is also small, so that
\be
\frac{\sin \theta}{\theta} = 1 - \frac{1}{6}\theta^2 + \cdots = 1 -
\delta + \cdots,
\label{r6}
\ee
or
\be
\theta^2 \approx 6 \delta.
\label{r7}
\ee
Substituting this approximation in \eqref{r1} yields
\be
R = \frac{1+\delta}{2\theta} \approx \frac{1 + \delta}{2\sqrt{6\delta}}.
\label{r8}
\ee
Substituing \eqref{r8} in \eqref{r4} then provides
\[
h \approx R \theta^2 /2 \approx
\left(\frac{1+\delta}{2\sqrt{6\delta}}\right) \frac{6\delta}{2} =
  \frac{\sqrt{6\delta} \left(1+\delta\right)}{4} 
\]
i.e.
\be
h \approx \sqrt{3\delta/8}.
\label{r9}
\ee
This is at the root of the surprising size of $h$: $\sqrt{\delta}$
dominates $\delta$ for small $\delta$. 

\begin{example}
Suppose $\delta = 10^{-4}$, which corresponds to expansion of $10$ cm
for a rail of length one kilometre. In this case $h = \sqrt{3 \times
  10^{-4}/8} = 6.1m$.
\end{example}

\newpage
\section{A Derivation of the FFT}
We choose $n = 2^M$ and illustrate the Fast Fourier Transform
algorithm, which computes the DFT in $O(M 2^M)$ operations.

Our primary data are the values $\{f(2\pi j/2^M): j = 01, 2, \ldots,
2^M - 1\}$ of our function evaluated at the $2^M$-th roots of one. For
each $m \in \{0, 1, \ldots, M-1\}$, we define
\be
F^{(m)}_{jk}
= \sum_{p=0}^{2^m - 1} f\left(e^{2\pi i\left(\frac{p}{2^m} +
  \frac{k}{2^M}\right)}\right) e^{-2\pi ijp/2^m}.
\label{fft1}\ee
for $j = 0, 1, \ldots, 2^m - 1$ and $k = 0, 1, \ldots,
2^{M-m}-1$. Thus $F^{(m)} \in \CC^{2^m \times 2^{M-m}}$. In other
words, each $F^{(m)}$ contains $2^M$ numbers, but their sizes are as
follows:

\begin{align*}
F^{(0)} &\hbox{\ is\ } 1 \times 2^M;\\
F^{(1)} &\hbox{\ is\ } 2 \times 2^{M-1};\\
F^{(2)} &\hbox{\ is\ } 2^2 \times 2^{M-2};\\
&\vdots \\
F^{(M-1)} &\hbox{\ is\ } 2^{M-1} \times 2;\\
F^{(M)} &\hbox{\ is\ } 2^{M} \times 1.\\
\end{align*}

\noindent
In other words, $F^{(0)}$ is a row vector, $F^{(m)}$ has twice the
number of rows as $F^{(m-1)}$, but half the number of columns,
and $F^{(M)}$ is a column vector.

\begin{example}
When $M = 3$ and $m = 2$, there are $2$ $4$-transforms.
\end{example}

\begin{example}
When $M = 3$ and $m=1$, there are $4$ $2$-transforms.
\end{example}

We now define a mapping constructing $F^{(m)}$ from
$F^{(m-1)}$. Specifically, we divide the sum over $p$ in 
\eqref{fft1} into even $p$ and odd $p$, as follows
\be
F^{(m)}_{jk} = E_m + O_m,
\label{fft2}
\ee
where
\be
E_m
= 
\sum_{q=0}^{2^{m-1}-1}
f(\exp(2\pi i\left(\frac{q}{2^{m-1}} + \frac{k}{2^N}\right) \exp(-2\pi i j
q/2^{m-1})
\label{fft3}
\ee
and
\be
O_m
= 
\sum_{r=0}^{2^{m-1}-1}
\left(f(\exp(2\pi i\left(\frac{r}{2^{m-1}} + \frac{k+2^{N-m}}{2^N}\right) \exp(-2\pi i j
q/2^{m-1})\right) \exp(-\pi i j/2^{m-1}).
\label{fft4}
\ee
Now $F^{(m-1)}$ is a $2^{m-1} \times 2^{M-m+1}$ matrix, but it is
useful to slightly abuse notation noting that $\ZZ \ni j \mapsto
F^{(m-1)}_{jk}$ is a $2^{m-1}$-periodic sequence. With this abuse of
notation in mind, we obtain
\be
F^{(m)}_{jk}
= F^{(m-1)}_{jk} + e^{-\pi i j/2^{m-1}} F^{(m-1)}_{j,k+2^{M-m}},
\label{fft5}
\ee
for $j = 0, 1, \ldots, 2^m-1$, $k=0, 1, \ldots, 2^{M-m}-1$.



\newpage
\section{Constrained Optimization}
Suppose we are considering investing money in two assets
whose returns are independent random variables $X_1$ and $X_2$. Their
distribution is unknown, but we do know the mean $\mu_k = \EE X_k$ and
the variance $\sigma_k^2 = \var X_k$, for $k = 1, 2$, and we shall
assume that these variances are strictly positive.

Being risk-averse, we want to divide our investment between the two assets
to minimize our risk. More formally, we have
\be
Y = s_1 X_1 + s_2 X_2, \quad\hbox{ where } s_1 + s_2 = 1.
\label{LA1}
\ee
Now, by the independence of $X_1$ and $X_2$, we have
\be
\var Y = f(\vs) = s_1^2 \sigma_1^2 + s_2^2 \sigma_2^2,
\qquad \vs = \left(s_1, s_2\right)^T \in \RR^2.
\label{LA2}
\ee

Thus our problem is as follows:
\be
\begin{array}{ll}
  \hbox{minimize} & f(\vs)\\
  \hbox{subject to} &g(\vs) = 1,
\end{array}
\label{LA3}
\ee
where
\be
g(\vs) = s_1 + s_2, \quad \vs = \left(s_1, s_2\right)^T \in \RR^2.
\label{LA4}
\ee

Now the function $f(\vs)$ satisfies
\[
\nabla f(\vs) =
\left(
\begin{array}{c}
  2\sigma_1^2 s_1\\
  2\sigma_2^2 s_2\\
\end{array}
\right)
\]
and
\[
D^2 f(\vs) =
\left(
\begin{array}{cc}
  2\sigma_1^2 & 0\\
  0 & 2\sigma_2^2\\
\end{array}
\right),
\]
and all higher derivatives vanish. In other words, $f(\vs)$ is a
quadratic and satisfies
\be
f(\vs + \vh)
= f(\vs) + \vh^T \nabla f(\vs) + \frac{1}{2} \vh^T D^2 f(\vs) \vh.
\label{LA5}
\ee
Further, the constraint function $g(\vs)$ is linear and satisfies
\be
g(\vs + \vh) = g(\vs) + \vh^T \nabla g(\vs) = g(\vs) + \vh^T \ve,
\label{LA6}
\ee
where
\be
\nabla g(\vs) \equiv \ve =
\left(
\begin{array}{c}
  1\\
  1\\
\end{array}
\right).
\label{LA6a}
\ee

One way to understand such problems is via {\em line search}: we
choose a point $\vs \in \RR^2$ and a search direction $\vd \in \RR^2$
and consider the univariate function
\be
\phi(t) = f(\vs + t \vd), \qquad t \in \RR.
\label{LA7}
\ee
Thus
\be
\phi(t) = f(\vs) + t \vd^T \nabla f(\vs) + \frac{1}{2} t^2 \vd^T D^2
f(\vs) \vd,
\label{LA8}
\ee
but we also require the search direction to satisfy the linear
constraint:
\be
1 = g(\vs + t \vd) = g(\vs) + t \vd^T \nabla g(\vs) = 1 + t \vd^T
\nabla g(\vs),
\label{LA9}
\ee
or
\be
\vd^T \nabla g(\vs) = 0
\label{LA10}
\ee
When do we know we are at a minimum? In this case, we must have
$\phi^\prime(0) = 0$ for any $\vd$ satisfying \eqref{LA10}. Hence
\be
\vd^T \nabla f(\vs) = \vd^T \nabla g(\vs) = 0,
\label{LA11}
\ee
which implies that
\be
\nabla f(\vs) = \lambda \nabla g(\vs),
\label{LA12}
\ee
for some $\lambda \in \RR$. In other words, we have
\be
\left(
\begin{array}{c}
  2\sigma_1^2 s_1\\
  2\sigma_2^2 s_2\\
\end{array}
\right)
= \lambda
\left(
\begin{array}{c}
  1\\
  1\\
\end{array}
\right),
\label{LA13}
\ee
which imply that
\be
s_k = \frac{1}{2}\lambda \sigma_k^{-2}, \quad k = 1, 2,
\quad\hbox{ and } s_1 + s_2 = 1.
\label{LA14}
\ee
Thus
\[
\lambda = \frac{2}{\sigma_1^{-2} + \sigma_2^{-2}}
\]
and
\be
s_k = \frac{\sigma_k^{-2}}{\sigma_1^{-2} + \sigma_2^{-2}}, \quad k = 1, 2.
\label{LA15}
\ee
The resulting minimal variance is then given by
\begin{align*}
\sigma^2
&\equiv f(\vs)\\
&=
\sigma_1^2 \frac{\sigma_1^{-4}}{\left(\sigma_1^{-2} +  \sigma_2^{-2}\right)^2}
+
\sigma_2^2 \frac{\sigma_2^{-4}}{\left(\sigma_1^{-2} +
  \sigma_2^{-2}\right)^2}\\
&= \frac{1}{\sigma_1^{-2} + \sigma_2^{-2}},\\
\label{LA16}
\end{align*}
or
\be
\sigma^{-2} = \sigma_1^{-2} + \sigma_2^{-2}.
\label{LA17}
\ee

\begin{example}
When $\sigma_1^2 = 1/10$ and $\sigma_2^2 = 1/5$, the minimal variance
is given by $\sigma^{-2} = 10 + 5 = 15$, or $\sigma^2 = 1/15$.
\end{example}

In general, we have
\[
Y = s_1 X_1 + s_2 X_2 = \frac{\sigma_1^{-2} X_1 + \sigma_2^{-2}
  X_2}{\sigma_1^{-2} + \sigma_2^{-2}}
\]
and
\[
\EE Y = s_1 \mu_1 + s_2 \mu_2 = \frac{\sigma_1^{-2} \mu_1 + \sigma_2^{-2}
  \mu_2}{\sigma_1^{-2} + \sigma_2^{-2}}
\]

Thus, if $\sigma_1 \gg \sigma_2$, then $\EE Y \approx \mu_2$, which is
to be expected, whilst $\sigma_1 = \sigma_2$ implies $\EE Y = (\mu_1 +
\mu_2)/2$.

\subsection{Lagrange Multipliers}
The above technique is much more general. Suppose we have a
risk-metric for investments in $n$ assets which is given by
\be
f(\vs) = \vs^T A \vs, \qquad \vs \in \RR^n,
\label{LA18}
\ee
where $A \in \RR^{n \times n}$ is a symmetric, positive definite
matrix. We want to solve the constrained optimization problem
\be
\begin{array}{ll}
  \hbox{minimize} & f(\vs)\\
  \hbox{subject to} &g(\vs) = 1,
\end{array}
\label{LA19}
\ee
where
\be
g(\vs) = \vw^T \vs, \qquad \vs \in \RR^n,
\label{LA20}
\ee
where $\vw \in \RR^n$ is some fixed vector. Then a similar argument
implies that
\be
\nabla f(\vs) = \lambda \nabla g(\vs),
\label{LA21}
\ee
where
\be
\nabla f(\vs) = 2 A \vs
\quad\hbox{ and }\quad \nabla g(\vs) = \vw.
\label{LA22}
\ee
Hence
\be
\vs = \frac{1}{2}\lambda A^{-1} \vw \quad\hbox{ and }\quad 1 =
\vw^T\vs,
\label{LA23}
\ee
which implies
\be
\lambda = \frac{2}{\vw^T A^{-1} \vw}
\label{LA24}
\ee
and
\be
\vs = \frac{A^{-1} \vw}{\vw^T A^{-1} \vw}.
\label{LA25}
\ee

\begin{exercise}
Prove that \eqref{LA25} implies that the corresponding minimal
risk-metric is given by
\be
f(\vs) = \vw^T A^{-1}\vw.
\label{LA26}
\ee
\end{exercise}

\newpage

\section{The Cholesky Factorization}
Let $\PP_n$ denote the set of all non-negative definite symmetric
matrices in $\RR^{n \times n}$. Given any $A_n \in \PP_n$, there is a
unique lower triangular matrix $L_n \in \RR^{n \times n}$, with
positive diagonal elements, for which $A_n = L_n L^T$, and this is
called the Cholesky factorization. This section provides a
constructive proof of this result, the factorization being obvious
when $n=1$.

Let us now consider the problem of computing the Cholesky
factorization $A_{n+1} = L_{n+1} L_{n+1}^T$, where
\be
A_{n+1}
=
\left(
\begin{array}{cc}
  A_n & \va\\
  \va^T & b\\
\end{array}
\right)
\in \RR^{(n+1)\times (n+1)},
\label{ch1}
\ee
where $\va \in \RR^n$, $b \ge 0$ and we assume that we have already
computed the Cholesky factorization $A_n = L_n L_n^T$. We define
\be
L_{n+1}
=
\left(
\begin{array}{cc}
  L_n & \vO\\
  \vp^T & q\\
\end{array}
\right),
\label{ch2}
\ee
where $\vp \in \RR^n$ and $q \ge 0$ are to be determined. Then
\be
A_{n+1}
=
\left(
\begin{array}{cc}
  A_n & \va\\
  \va^T & b\\
\end{array}
\right)
= 
\left(
\begin{array}{cc}
  L_n & \vO\\
  \vp^T & q\\
\end{array}
\right)
\left(
\begin{array}{cc}
  L_n^T & \vp\\
  \vO^T & q\\
\end{array}
\right),
\label{ch4}
\ee
and $\vp$ and $q$ must therefore satisfy the equations
\be
L_n \vp = \va
\label{ch5}
\ee
and
\be
\|\vp\|^2 + q^2 = b.
\label{ch6}
\ee
It is \eqref{ch6} that presents the difficulty: we must prove that
\be
b \ge \|\vp\|^2 = \|L_n^{-1} \va\|^2
\label{ch6a}
\ee
to ensure that $q^2 \ge 0$. To this end, we shall first deal with the
simpler case when $A_n = I_n$.

\begin{lem}
Let $A_n = I_n$. Then $b \ge \|\va\|^2$.
\label{lemmach1}
\end{lem}

\begin{proof}
For any $\vv \in \RR^n$ and $w \in \RR$ we have
\begin{align*}
0
&\le
\left(
\begin{array}{c}
  \vv\\
  w\\
\end{array}
\right)^T
\left(
\begin{array}{cc}
  I_n & \va\\
  \va^T & b\\
\end{array}
\right)
\left(
\begin{array}{c}
  \vv\\
  w\\
\end{array}
\right)\\
&= \vv^T \vv + 2w \vv^T \va + b w^2\\
&= \|\vv + w\va\|^2 + \left(b - \|\va\|^2\right) w^2.
\end{align*}
Setting $w=1$ and $v = -a$, we obtain $0 \le b - \|\va\|^2$, as desired.
\end{proof}

To extend this result to the original case, we use the following trick
to relate the general $A_{n+1}$ to the case where $A_n = I_n$.
\begin{align*}
&\left(
\begin{array}{cc}
  L_n^{-1} & \vO\\
  \vO^T & 1\\
\end{array}
\right)
\left(
\begin{array}{cc}
  A_n & \va\\
  \va^T & b\\
\end{array}
\right)
\left(
\begin{array}{cc}
  L_n^{-T} & \vO\\
  \vO^T & 1\\
\end{array}
\right)\\
&\qquad=
\left(
\begin{array}{cc}
  L_n^{-1}A_nL_n^{-T} & L_n^{-1}\va\\
  \left(L_n^{-1}\va\right)^T & b\\
\end{array}
\right)\\
&\qquad=
\left(
\begin{array}{cc}
  I_n & L_n^{-1}\va\\
  \left(L_n^{-1}\va\right)^T & b\\
\end{array}
\right).\\
\end{align*}

Hence Lemma \ref{lemmach1} implies that
\[
b \ge \|L_n^{-1} \va\|^2,
\]
which is \eqref{ch6a}, as required.

Students are often be more familiar with the square-root defined by
$A^{1/2} = QD^{1/2}A^TQ$, where $A = QDQ^T$ is the spectral
factorization of $A$, rather than the Cholesky factorization $A =
LL^T$. Thus $(A^{1/2})^2 = L L^T$, and it can be shown that $L =
A^{1/2}W$, where $W$ is an orthogonal matrix. [Essentially the
argument is as follows. If we compute the SVD $L = USV^T$, where $U$
and $V$ are orthogonal matrices and $S$ is the diagonal matrix of
singular values of $L$, then $LL^T = (USV^T)(VSU^T)=US^2 U^T = A =
QDQ^T$. Hence $U=Q$ and $S = D^{1/2}$. Thus $L = QD^{1/2}V^T =
A^{1/2} W$, where $W = QV^T$.]

With this in mind, we see that \eqref{ch6a} becomes
\be
b \ge \|L_n^{-1}\va\|^2 = \|A_n^{-1/2}\va\|^2 = \va^T A_n^{-1}\va.
\label{ch7a}
\ee

Once we know condition \eqref{ch7a}, it's possible to remove all of
the scaffolding used above, although I believe most readers will find
the more circuitous route described above useful: it's often good to
leave some scaffolding in place!

\begin{lem}
Let $A_n \in \RR^{n \times n}$  be any symmetric non-negative definite matrix
and define $A_{n+1} \in \RR^{(n+1)\times (n+1)}$ by \eqref{ch1}. Then
$A_{n+1}$ is non-negative definite if and only if
\be
b \ge \va^T A_n^{-1}\va.
\label{ch7b}
\ee
\label{lemmach2}
\end{lem}

\begin{proof}
For any $\vv \in \RR^n$ and $w \in \RR$ we have
\begin{align*}
0
&\le
\left(
\begin{array}{c}
  \vv\\
  w\\
\end{array}
\right)^T
\left(
\begin{array}{cc}
  A_n & \va\\
  \va^T & b\\
\end{array}
\right)
\left(
\begin{array}{c}
  \vv\\
  w\\
\end{array}
\right)\\
&= \vv^T A_n \vv + 2w \vv^T \va + b w^2\\
&= \|A_n^{1/2}\vv + wA_n^{-1/2}\va\|^2 + \left(b - \va^T A_n^{-1} \va\right) w^2.
\end{align*}
[How did I complete the square here? The key point is that
  $\vv^TA_n\vv = (A_n^{1/2}\vv)^T (A_n^{1/2}\vv)$, which implies that
  we must then write the second term as $\vv^T \va = (A_n^{1/2}\vv)^T (A_n^{-1/2}\va)$.]
If $A_{n+1}$ is non-negative definite, then setting $w=1$ and $v = -a$
we obtain $0 \le b - \va^T A_n^{-1}\va$. Conversely, if $b - \va^T
A_n^{-1} \va \ge 0$, then $A_{n+1}$ is non-negative definite.
\end{proof}



\newpage
\section{Conics}

\subsection{The Ellipse and Hyperbola}

Let's begin with the ellipse and the hyperbola, 
which we shall define as contours of the
functions
\be
f_{\pm}(\vx) = \|\vx+\vs\| \pm \|\vx-\vs\|, \qquad \vx \in \RR^n.
\label{c1}
\ee

The key trick is to observe that
\be
4\vx^T\vs
= \|\vx+\vs\|^2 - \|\vx-\vs\|^2 = f_+(\vx) f_{-}(\vx).
\label{c2}
\ee

If $f_+(\vx) = \alpha$, then $f_-(\vx) = 4\vx^T\vs/\alpha$. Adding
these equations gives
\be
f_+(\vx) + f_-(\vx) = 2\|\vx+\vs\| = \alpha + \frac{4\vx^T\vs}{\alpha} = 
\frac{\alpha^2 + 4\vx^T\vs}{\alpha},
\label{c3}
\ee
and squaring both sides yields the quadratic form
\be
4\|\vx+\vs\|^2 = \left(\frac{\alpha^2 + 4\vx^T\vs}{\alpha}\right)^2.
\label{c4}
\ee
The matrix occurring in  this quadratic form is
\be
M = 4\left(I_n - \frac{4}{\alpha^2}\vs \vs^T\right).
\label{c5}
\ee

Similarly, if 
$f_-(\vx) = \alpha$, then $f_+(\vx) = 4\vx^T\vs/\alpha$ and adding
them yields \eqref{c4} and \eqref{c5}. The distinction between the
contours of $f_\pm$ lies in the eigenvalues of $M$, which are $1$
(with multiplicity $n-1$) and
\[
\lambda = 4\left(1 - \frac{4\|\vs\|^2}{\alpha^2}\right).
\]
If $f_+(\vx) = \alpha$, then the triangle inequality implies that
\[
\|\vs\| \le \frac{1}{2}\left( \|\vs+\vx\| + \|\vs-\vx\| \right) 
= \frac{\alpha}{2},
\]
i.e. 
\[
4 \|\vs\|^2 \le \alpha^2,
\]
which implies $\lambda \ge 0$,
with inequality if and only if $\vx$ and $\pm \vs$ are collinear.
Thus $M$ is non-negative definite on
contours of $f_+$ and $M$ is positive definite when the contour
is not the line segment joining $\pm \vs$.

In contrast, the triangle inequality also implies that
\[
\alpha 
= \|\vx+\vs\| - \|\vx-\vs\| 
\le \|\vx+\vs - \left(\vx - \vs\right)\|
= 2 \|\vs\|,
\]
or $4\|\vs\|^2/\alpha^2 \ge 1$ on contours of $f_-$, i.e. $\lambda \le 0$.

\subsection{The Reflector Property}

We have
\be
\nabla f_{\pm}(\vx) = \frac{\vx + \vs}{\|\vx + \vs\|}
\pm \frac{\vx-\vs}{\|\vx-\vs\|}.
\label{cr1}
\ee
Then
\be
\left(\frac{\vx+\vs}{\|\vx+\vs\|}\right)^T\nabla f_{\pm}(\vx) 
= 1
\pm \frac{\left(\vx+\vs\right)^T\left(\vx-\vs\right)}{\|\vx+\vs\|\|\vx-\vs\|}.
\label{cr2}
\ee
and
\be
\left(\frac{\vx-\vs}{\|\vx-\vs\|}\right)^T\nabla f_{\pm}(\vx) 
=
\frac{\left(\vx+\vs\right)^T\left(\vx-\vs\right)}{\|\vx+\vs\|\|\vx-\vs\|}
\pm 1.
\label{cr3}
\ee
Thus
\be
\left(\frac{\vx+\vs}{\|\vx+\vs\|}\right)^T\nabla f_{\pm}(\vx) 
=
\pm \left(\frac{\vx-\vs}{\|\vx-\vs\|}\right)^T\nabla f_{\pm}(\vx),
\label{cr4}
\ee
which is the reflector property for the ellipse and hyperbola.

\subsection{Conic sections really are conic sections}

Let us take a cone of semi-angle $\theta$ in $\RR^n$ whose axis is the line
generated by a unit vector $\vu \in \RR^n$, i.e.
\be
C = \{\vx \in \RR^n : \frac{\vu^T\vx}{\|\vx\|} = \pm \cos \theta\}.
\label{cs1}
\ee
In other words, the equation for the cone is
\be
\left(\vu^T\vx\right)^2 = \|\vx\|^2 \cos^2 \theta,
\label{cs2}
\ee
or
\be
\vx^T\left(I_n \cos^2 \theta - \vu \vu^T\right)\vx = 0.
\label{cs3}
\ee
Now let $\vv_1, \vv_2, \ldots, \vv_n$ be any orthonormal basis for
$\RR^n$ and consider the hyperplane $P$ with
normal vector $\vv_n$ at signed distance $z_n$ from the origin. In
other words,
\be
P = \{\vx = \sum_{k=1}^n z_k \vv_k : z_1, z_2, \ldots, z_{n-1}\in \RR\}.
\label{cs4}
\ee
\noindent
If we let $V \in \RR^{n \times n}$ be the orthogonal matrix with
columns $\vv_1, \ldots, \vv_n$ and substitute $\vx = V \vz$ in
\eqref{cs2}, 
then we obtain
\be
\left(\vz^T V^T \vu\right)^2
=
\|\vz\|^2 \cos^2 \theta.
\label{cs5}
\ee
Setting
\be
\vU = V^T \vu,
\label{cs6}\ee
we see that \eqref{cs5} becomes
\be
\left(\sum_{k=1}^n z_k U_k\right)^2
=
\left(\sum_{k=1}^n z_k^2\right) \cos^2 \theta,
\label{cs6a}
\ee
or
\be
\sum_{k, \ell=1}^{n-1} z_k z_\ell U_k U_\ell + 2 z_n U_n
\sum_{k=1}^{n-1} z_k U_k + z_n^2 U_n^2
= \left(\sum_{k=1}^{n-1} z_k^2\right)\cos^2 \theta + z_n^2 \cos^2 \theta.
\label{cs7}\ee
Hence, writing 
\[
\vzhat = 
\left(
\begin{array}{c}
z_1\\ 
z_2\\
\vdots\\ 
z_{n-1}
\end{array}
\right)
\quad\hbox{ and }\quad
\vUhat = 
\left(
\begin{array}{c}
U_1\\ 
U_2\\
\vdots\\ 
U_{n-1}
\end{array}
\right)
\]
\eqref{cs7} becomes the quadratic form
\be
\vzhat^T M \vzhat - 2z_n U_n \vzhat^T\vUhat + z_n^2\left(\cos^2 \theta -
U_n^2\right) = 0,
\label{cs8}\ee
where the matrix $M \in \RR^{(n-1)\times (n-1)}$ is given by
\be
M = I_{n-1}\cos^2 \theta - \vUhat \vUhat^T.
\label{cs9}\ee

\begin{example}
Let us choose $\vv_n = \vu$, so that $\vUhat = 0$ and $U_n = 1$. Then
\eqref{cs8} becomes
\[
\left(\sum_{k=1}^{n-1} z_k^2\right) \cos^2 \theta - z_n^2 \sin^2
\theta = 0,
\]
or
\[
\sum_{k=1}^{n-1} z_k^2 = z_n^2 \tan^2 \theta.
\]
\label{conicsex1}
\end{example}

\noindent
The eigenvalues of $M$ are $\cos^2 \theta$ (with multiplicity $n-2$)
and 
\[
\mu := \cos^2 \theta - \|\vUhat\|^2.
\]
Now 
\[
1 = \|\vu\|^2 = \sum_{k=1}^n \left(\vu^T\vv_k\right)^2 = \|\vUhat\|^2 + U_n^2,
\]
which implies
\be
\mu = \cos^2 \theta - \left(1 - U_n^2\right) = U_n^2 - \sin^2 \theta.
\label{cs10}\ee

\begin{example}
Let $n=3$ and suppose $\mu = 0$. Then
\[
\|\vUhat\|^2 = \cos^2 \theta
\]
and 
\[
U_n = \pm \sin \theta.
\]
If $\vq_1=\vUhat/\|\vUhat\|$ and $\vq_2 \in
\RR^2$ is orthogonal to $\vUhat$, then the matrix 
\[
Q
=\left(
\begin{array}{cc}
\vq_1 & \vq_2\\
\end{array}
\right) \in \RR^{2 \times 2}
\]
is orthogonal and
\[
D := Q^T M Q =
\left(
\begin{array}{cc}
\mu & 0\\
0 & \cos^2 \theta\\
\end{array}
\right).
\]
If we let $\vzhat = Q \vx$, then 
\be
\cos^2 \theta\  x_2^2 \pm 2z_3 \sin \theta \cos \theta x_1 +z_3^2
\left(\cos^2 \theta - \sin^2 \theta\right) = 0.
\label{cs12}\ee

\label{conicsex2}
\end{example}


\newpage
\section{The Birthday Problem}

This is a traditional probabilistic problem: given $n$ people, whose
birthdays are assumed to be uniformly distributed over the $N=365$
days of the year (ignoring leap years), find the probability that at
least two of them share a birthday. Now
\be
\PP(\hbox{at least two share a birthday}) = 1 - \PP(\hbox{all
  birthdays distinct}) =: 1 - p_n.
\label{bd1}
\ee
Now
\be
p_n = \frac{N(N-1)(N-2)\cdots(N-n+1)}{N^n}
\label{bd2}
\ee
and, dividing numerator and denominator by $N^n$, we obtain
\be
p_n = 
\left(1 - \frac{1}{N}\right)\left(1 - \frac{2}{N}\right)
\cdots\left(1 - \frac{n-1}{N}\right).
\label{bd3}
\ee
In one sense the problem is now solved. The surprise is that $p_n$
tends to zero rather quickly. Indeed, $p_{23} = 0.5073$, by direct
calculation. However, plotting $p_n$ reveals a suspiciously Gaussian
curve, as we see in Figure \ref{bday1}. Why does $p_n$ decay so
quickly and can we understand the seemingly Gaussian behaviour?

\bigskip
\begin{verbatim}
N=365; n=100; p=1; prob=zeros(1,n);
%
% prob(k) = (1 - 1/N)(1 - 2/N)...(1-k/N)
%         = prob(all k+1 bdays different) 
%
for k=1:n
  p=p*(1-k/N);
  prob(k)=p;
end
\end{verbatim} 
\bigskip

\begin{figure}[tb]
\vspace{6pt}
\includegraphics[scale=0.6, bb=0 200 600 600,clip]{bday1.pdf}
\parbox{28.5pc}{\caption{The probability that all $n$ birthdays are different}\label{bday1}}
%\end{center}
\end{figure}

First take logarithms:
\be
\log p_n = \sum_{k=1}^{n-1} \log \left(1 - \frac{k}{N}\right).
\label{bd4}
\ee
Now
\[
\log \left(1 + x\right) = x - \frac{x^2}{2} + \frac{x^3}{3} -
\frac{x^4}{4} + \cdots,
\]
and the series is convergent for $|x| < 1$. Thus
\be
\log \left(1 - x\right) = - x - \frac{x^2}{2} - \frac{x^3}{3} -
\frac{x^4}{4} - \cdots \le -x,
\label{bd5}
\ee
for $0 \le x < 1$. If $|x| \ll 1$, then we also have the approximation
$\log \left(1 + x\right) \approx -x$. Thus
\be
\log p_n \le -\sum_{k=1}^{n-1} \frac{k}{N} = -\frac{n(n-1)}{2N},
\label{bd6}
\ee
which implies
\be
p_n \le e^{-n(n-1)/(2N)}.
\label{bd7}
\ee
This explains the rapid decay and the Gaussian resemblance, as we see
in Figure \ref{bday2}.

\begin{figure}[tb]
\vspace{6pt}
\includegraphics[scale=0.6, bb=0 200 600 600,clip]{bday2.pdf}
\parbox{28.5pc}{\caption{x: $p_n$; line: $\exp(-n(n-1)/(2N))$}\label{bday2}}
%\end{center}
\end{figure}

\newpage
\section{The Bike Problem via the Inclusion--Exclusion Formula}
Suppose $n$ cyclists randomly permute their bikes. What is the
probability that at least one cyclist has the correct bike?

More formally, the sample space $X$ consits of all positive
permutations of the integers $1, 2, \ldots, n$, i.e.
\be
X = \{(i_1, i_2, \ldots, i_n) : i_1, \ldots, i_n \hbox{ a permutation of } 1,
\ldots, n\}.
\label{bikep1}
\ee
We shall assign each of these permutations the same probability
$1/n!$.

Further, define
\be
A_k = \{x \in X : i_k = k\},
\qquad\hbox{ for } k = 1, 2, \ldots, n.
\label{bikep2}
\ee
Thus $A_k$ is the set of outcomes for which cyclist $k$ gets bike
$k$. We want to calculate the probability
\[
\PP\left(A_1 \cup A_2 \cup \cdots \cup A_n\right).
\]

\begin{example}
If $n=3$, then the sample space is
\[
X = \{(1 2 3), (1 3 2), (2 1 3), (2 3 1), (3 1 2), (3 2 1)\}.
\]
Then $A_1 = \{(1 2 3), (1 3 2)\}$, $A_2 = \{(1 2 3), (3 2 1)\}$
and $A_3 = \{(1 2 3), (2 1 3)\}$, whilst
$\PP(A_1 \cup A_2 \cup A_3) = 4/6 = 2/3$.
\label{bpex1}
\end{example}

The solution requires the inclusion--exclusion formula:
\be
\PP(A_1 \cup A_2 \cup \cdots \cup A_k)
= \sum_{\ell=1}^n (-1)^{\ell-1} \sum_{1 \le k_1 < k_2 < \cdots <
  k_\ell \le n} \PP(A_{k_1} \cap A_{k_2} \cap \cdots \cap A_{k_\ell}).
\label{bikep3}
\ee

\begin{exercise}
\[
\PP(A_{k_1} \cap \cdots \cap A_{k_m}
= \frac{(n-m)!}{n!}.
\]
\end{exercise}

Thus
\[
\sum_{1 \le k_1 < k_2 < \cdots <
  k_\ell \le n} \PP(A_{k_1} \cap A_{k_2} \cap \cdots \cap A_{k_\ell})
= \binom{n}{m} \frac{(n-m)!}{n!} = \frac{1}{m!}.
\]
Hence
\be
\PP(A_1 \cup \cdots \cup A_n) = 1 - \frac{1}{2!} + \frac{1}{3!} +
\cdots + \frac{(-1)^{n-1}}{n!} \to 1 - e^{-1}.
\label{bikep4}
\ee

\subsection{The Inclusion--Exclusion Formula}
For each subset $A$ of the sample space $X$, the {\em indicator
function} $I_A : X \to \{0, 1\}$ is defined by $I_A(x) = 1$ if and
only if $x \in A$.

\begin{example}
For Example \ref{bpex1} we have
\[
I_{A_1}(1 2 3) = I_{A_1}(1 3 2) = 1
\]
but
\[
I_{A_1}(2 1 3) = 
I_{A_1}(2 3 1) = 
I_{A_1}(3 1 2) = 
I_{A_1}(3 2 1) = 0.
\]
\end{example}

The indicator function has some crucial properties. Firstly
\[
1 - I_A(x) = I_{A^c}(x)
\]
where $A_c$ is the complement of $A$, i.e. $X \setminus A$. Further,
\[
I_{A \cap B}(x) = I_A(x) I_B(x).
\]
Further, we use de Morgan's Law:
\[
\left(A_1 \cup \cdots \cup A_n\right)^c = A_1^c \cap \cdots \cap A_n^c.
\]
Thus
\begin{align*}
I_{A_1 \cup \cdots \cup A_n}(x)
&= 1 - I_{(A_1 \cup \cdots \cup A_n)^c}(x)\\
&= 1 - \prod_{k=1}^n \left(1 - I_{A_k}(x)\right).
\end{align*}



\newpage
\begin{thebibliography}{99}
\bibitem{afb} Beardon (2005), 
{\em Algebra and Geometry}, CUP.

\end{thebibliography}

\end{document}






