\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}