\documentclass{article}
\input{../macros.tex}

\begin{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Lecture 26: Complex matrices; fast Fourier transform}

Matrices with all real entries can have complex eigenvalues! 
So we can't avoid working with complex numbers.  In this lecture we
learn to work with complex vectors and matrices.

The most important complex matrix is the Fourier
matrix $F_n$, which is used for Fourier transforms.  
Normally, multiplication by $F_n$ would require $n^2$ multiplications.
The fast Fourier transform (FFT) reduces this to roughly $n \log_2n$
multiplications, a revolutionary improvement.

\subsection*{Complex vectors}

\subsubsection*{Length}

Given a vector $\vb z = \small \fourbyone {z_1}{z_2}\vdots{z_n}\in \C^n$ with complex entries, how do we find its length?  Our old definition:
$${\small \vb z^T\vb z = \onebyfour {z_1}{z_2}\cdots{z_n} \fourbyone
{z_1}{z_2}\vdots{z_n}}$$ is no good; this quantity isn't
always positive.  For example:
$$\onebytwo 1i \twobyone 1i = 0.$$ We don't want to define the length
of $\twobyone 1i$ to be $0$.  The correct definition is: $|\vb z|^2 =
\overline{\vb z}^T \vb z = |z_1|^2 + |z_2|^2 + \cdots + |z_n|^2$.  Then we have:
$$\left(\text{length} \twobyone 1i\right)^2 = \onebytwo 1{-i} \twobyone 1i = 2.$$
To simplify our notation we write $|\vb z|^2 = \vb z^H \vb z$, where
$\vb z^H = \overline{\vb z}^T$.  The $H$ comes from the name Hermite, and
$\vb z^H \vb z$ is read ``$\vb z$ Hermitian $\vb z$''.

\subsubsection*{Inner product}

Similarly, the inner or dot product of two complex vectors is not just $\vb y^T \vb x$.  We must also take the complex conjugate of $\vb y$: $$\vb y^H \vb x = \overline {\vb y}^T \vb x = \overline y_1 x_1 + \overline y_2 x_2 + \cdots + \overline y_n x_n.$$

\subsection*{Complex matrices}

\subsubsection*{Hermitian matrices}

Symmetric matrices are real valued matrices for which $\mb A^T = \mb
A$.  If $\mb A$ is complex, a nicer property is $\overline{\mb A}^T = \mb
A$; such a matrix is called {\em Hermitian} and we abbreviate
$\overline{\mb A}^T$ as $\mb A^H$.  Note that the diagonal entries of a
Hermitian matrix must be real.  For example,
$$\overline{\mb A}^T = \mb A = \twobytwoc 2{3+i}{3-i}5.$$
Similar to symmetric matrices, Hermitian matrices have real
eigenvalues and perpendicular eigenvectors.

\subsubsection*{Unitary matrices}

What does it mean for complex vectors $\vb q_1, \vb q_2, ..., \vb q_n$ to be
perpendicular (or orthonormal)?  We must use our new definition of the
inner product.  For a collection of $\vb q_j$ in complex space to be
orthonormal, we require:
$$\overline{\vb q}_j \vb q_k = \left\{ \begin{array}{cl} 0 & j \neq
  k\\ 1 & j = k\end{array}\right.$$ We can again define $\mb Q =
  \onebyfour {\vb q_1} {\vb q_2} \cdots {\vb q_n}$, and then $\mb Q^H
  \mb Q = I$.  Just as ``Hermitian'' is the complex equivalent of
  ``symmetric'', the term ``{\em unitary}'' is analogous to
  ``orthogonal''.  A {\em unitary matrix} is a square matrix with
  perpendicular columns of unit length.

\subsection*{Discrete Fourier transform}

A {\em Fourier series} is a way of writing a periodic function or {\em signal} as a sum of functions of different frequencies:
$$f(x) = a_0 + a_1 \cos x + b_1 \sin x + a_2 \cos 2x + b_2 \sin 2x +
\cdots.$$  When working with finite data sets, the {\em discrete Fourier transform} is the key to this decomposition.

\smallskip

In electrical engineering and computer science, the rows and columns
of a matrix are numbered starting with $0$, not $1$ (and ending with
$n-1$, not $n$).  We'll follow this convention when discussing the Fourier matrix:
$$F_n =  \left[ \begin{array}{ccccc}
1 & 1 & 1 & \cdots & 1\\
1 & w & w^2 & & w^{n-1}\\
1 & w^2 & w^4 && w^{2(n-1)}\\
\vdots & & & \ddots & \vdots\\ 
1 & w^{n-1} & w^{2(n-1)} & \cdots & w^{(n-1)^2}
\end{array}\right].$$
Notice that $F_n = F_n^T$ and $(F_n)_{jk} = w^{jk}$, where $j,k = 0, 1, ...,
n-1$ and the complex number $w$ is $w = e^{i\cdot2\pi/n}$ (so $w^n = 1$).  The columns of this matrix are
orthogonal.

All the entries of $F_n$ are on the unit circle in the
complex plane, and raising each one to the $n$th power gives $1$.  We could write $w =
\cos (2 \pi/n) + i \sin(2 \pi /n)$, but that would just make it harder
to compute $w^{jk}$.

\bigskip

Because
$w^4 = 1$ and $w = e^{2\pi i/4} = i$, our best example of a Fourier matrix is:
$$F_4 = \left[ \begin{array}{cccc}
1 & 1 & 1 & 1\\
1 & i & i^2 & i^3 \\
1 & i^2 & i^4 & i^6 \\
1 & i^3 & i^6 & i^9
\end{array}\right]
 = \left[ \begin{array}{rrrr}
1 & 1 & 1 & 1\\
1 & i & -1 & -i \\
1 & -1 & 1 & -1 \\
1 & -i & -1 & i
\end{array}\right].$$To find the Fourier transform of a vector with four
components (four data points) we multiply by $\mb F_4$.  

It's easy to check that the columns of $\mb F_4$ are orthogonal, as long as we remember to conjugate when
computing the inner product. However,
$\mb F_4$ is not quite unitary because each column has length
$2$.  We could divide each entry by $2$ to get a matrix whose
columns are orthonormal: $$\frac{1}{4} \mb F_4^H \mb F_4 = \mb I.$$

\subsubsection*{An example}

The signal corresponding to a single impulse at time zero is (roughly) described by $\small \fourbyone 1000$.  To find the Fourier transform of this signal we compute:
$$ \left[ \begin{array}{rrrr}
1 & 1 & 1 & 1\\
1 & i & -1 & -i \\
1 & -1 & 1 & -1 \\
1 & -i & -1 & i
\end{array}\right]
 \fourbyone 1000
= \fourbyone 1111.$$
A single impulse has all frequencies in equal amounts.

If we multiply by $\mb F_4$ again we almost get back to $(1,0,0,0)$:
$$ \left[ \begin{array}{rrrr}
1 & 1 & 1 & 1\\
1 & i & -1 & -i \\
1 & -1 & 1 & -1 \\
1 & -i & -1 & i
\end{array}\right]
 \fourbyone 1111
= \fourbyone 4000 = 4 \fourbyone 1000.$$
Because $\frac{1}{\sqrt{n}}\mb F_n$ is unitary, multiplying by $\mb F_n$ and dividing by the scalar $n$ inverts the transform.

\subsection*{Fast Fourier transform}

Fourier matrices can be broken down into chunks with lots of zero entries;
Fourier probably didn't notice this.  Gauss did, but didn't realize how significant a discovery this was.

There's a nice relationship between $\mb F_n$ and $\mb F_{2n}$ related
to the fact that $w_{2n}^2 = w_n$:
$$\mb F_{2n} = 
\twobytwo {\mb I} {\mb D} {\mb I} {-\mb D} 
\twobytwo {\mb F_n}{\mb 0}{\mb 0}{\mb F_n}
\mb P,$$
where $\mb D$ is a diagonal matrix and $\mb P$ is a $2n$ by $2n$ permutation matrix: $$\mb P=
 \small \left[\begin{array}{ccccccc}
1 & 0 & 0 & 0 & \cdots & 0 & 0\\
0 & 0 & 1 & 0 & \cdots & 0 & 0\\
 &  & & \vdots &  & &\\
0 & 0 & 0 & 0 & \cdots & 1 & 0\\
0 & 1 & 0 & 0 & \cdots & 0 & 0\\
0 & 0 & 0 & 1 & \cdots & 0 & 0\\
 &  && \vdots &  &  &\\
0 & 0 & 0 & 0 & \cdots & 0 & 1\\
\end{array}\right].
$$ So, a $2n$ sized Fourier transform $\mb F$ times $\vb x$ which we might think would
require $(2n)^2 = 4n^2$ operations can instead be performed using two
size $n$ Fourier transforms ($2n^2$ operations) plus two very simple
matrix multiplications which require on the order of $n$
multiplications.  The matrix $\mb P$ picks out the even
components $x_0$, $x_2$, $x_4$, ... of a vector first, and then the
odd ones -- this calculation can be done very quickly.

Thus we can do a Fourier transform of size 64 on a vector by
separating the vector into its odd and even components, performing a
size 32 Fourier transform on each half of its components, then
recombining the two halves through a process which involves
multiplication by the diagonal matrix $\mb D$.
$$D =  \small \left[\begin{array}{cccccc}
1\\
&w\\
&&w^2\\
&&&\ddots\\
&&&&w^{n-1}
\end{array}\right].$$

Of course we can break each of those copies of $\mb F_{32}$ down into two
copies of $\mb F_{16}$, and so on.  In the end, instead of using $n^2$
operations to multiply by $F_n$ we get the same result using about
$\frac{1}{2}n \log_2 n$ operations.

A typical case is $n = 1024 = 2^{10}$.  Simply multiplying by $\mb
F_n$ requires over a million 
calculations.  The fast Fourier transform can be completed with only
$\frac{1}{2}n\log_2n = 5\cdot1024$ calculations.  This is $200$ times faster!

\smallskip

This is only possible because Fourier matrices are special matrices
with orthogonal columns.  In the next lecture we'll return to dealing
exclusively with real numbers and will learn about positive definite
matrices, which are the matrices most often seen in applications.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}