Complex Matrices; Fast Fourier Transform

Why complex numbers, and the one rule to remember

Complex numbers have slipped into this course already: even a real matrix can have complex eigenvalues, and with them complex eigenvectors. Most of what remains about eigenvalues will be real, but at some point we have to face the question squarely — what changes when the vectors are complex, when the matrices are complex? What is the dot product of two complex vectors?

The answer to every one of those questions is the same small change, and it’s the message of the first half of this lecture: whenever we transpose, we also take complex conjugates.

The second half is the payoff — the most important complex matrix there is, the Fourier matrix, and the idea everyone calls the FFT. The fast Fourier transform is running “in a thousand places” as we speak; it transformed whole industries. And it turns out to be nothing more than a simple matrix factorization — one that cuts the cost of multiplying by an $n \times n$ matrix from $n^2$ down to $n \log_2 n$.

Length and inner product in complex space

Take a vector $\mathbf{z}$ with components $z_1, z_2, \ldots, z_n$ — complex numbers now, so $\mathbf{z}$ lives in $\mathbb{C}^n$, complex $n$-dimensional space, instead of $\mathbb{R}^n$.

What’s its length? The old recipe, $\mathbf{z}^T\mathbf{z}$, is no good. Length squared should be positive, and Strang’s two-component example shows the old recipe failing: take $\mathbf{z} = (1, i)$ in $\mathbb{C}^2$. Then

$$\mathbf{z}^T\mathbf{z} = 1 \cdot 1 + i \cdot i = 1 - 1 = 0,$$

but that vector certainly doesn’t have length zero. What we really want in each slot is $\bar{z}_1 z_1$ — conjugate times the number — because $\bar{z}_1 z_1 = |z_1|^2$, the modulus squared, which is what a length ought to be built from. If $z_1 = i$, then $(-i)(i) = +1$: exactly right. So the correct length squared is

$$\|\mathbf{z}\|^2 = \bar{\mathbf{z}}^T\mathbf{z} = |z_1|^2 + |z_2|^2 + \cdots + |z_n|^2 .$$

For $(1, i)$: conjugate to get $(1, -i)$, multiply, and get $1 + 1 = 2$. Length $\sqrt{2}$ — not the zero the unconjugated formula gave.

The inner product definition is forced on us: once the length is $\mathbf{z}^{\mathrm{H}}\mathbf{z}$, the inner product has to match it.

Hermitian and unitary: the complex dictionary

Two more words need updating — the book gives a little dictionary translating between the real case, $\mathbb{R}^n$, and the complex case, $\mathbb{C}^n$.

Symmetric becomes Hermitian. For real matrices, symmetry means $A^T = A$. That’s no good if $A$ is complex; the right complex version is: flip across the diagonal and take conjugates. A matrix with $A^{\mathrm{H}} = A$ is called a Hermitian matrix. On the diagonal, the entries had better be real — flipping leaves them in place, so conjugating must leave them alone. Off the diagonal, entries pair up as conjugates:

$$A = \begin{bmatrix} 2 & 3+i \\ 3-i & 5 \end{bmatrix}$$

is Hermitian: the $3+i$ above the diagonal must meet a $3 - i$ below it. These matrices are the complex heirs of symmetric matrices, and they inherit the good properties: real eigenvalues and perpendicular eigenvectors.

Orthogonal becomes unitary. But what does perpendicular mean now? For an orthonormal basis $\mathbf{q}_1, \ldots, \mathbf{q}_n$, the test used to be $\mathbf{q}_i^T\mathbf{q}_j$; now we must conjugate too:

$$\mathbf{q}_i^{\mathrm{H}}\mathbf{q}_j = \begin{cases} 0 & i \neq j \\ 1 & i = j \end{cases}$$

— right angles measured in complex $n$-dimensional space, unit lengths measured with the conjugate. Put those columns into a matrix $Q$ and the matrix statement, which used to be $Q^TQ = I$, becomes $Q^{\mathrm{H}}Q = I$. A square matrix with orthonormal columns in this complex sense gets a new name: a unitary matrix — just like an orthogonal matrix, with the word signaling that the entries may be complex.

The Fourier matrix

Now the substance: the most famous complex matrix, one of these unitary-type matrices with orthogonal columns, named after Fourier because it is the Fourier transform. For this half hour Strang moves “from the math department to EE” — electrical engineers start counting at zero, and for this matrix they’re probably right. So rows and columns are numbered $0$ through $n-1$, and the entries could not be simpler:

$$(F_n)_{ij} = W^{ij}, \qquad i, j = 0, 1, \ldots, n-1,$$

the exponent is the row number times the column number. Column zero is all ones ($W^0$); column one holds the powers $1, W, W^2, \ldots, W^{n-1}$; column two holds the powers of $W^2$; and so on. Every entry is a power of one special number — none are zero, this is a full matrix, and it’s symmetric (in the old, not-quite-perfect $A^T = A$ sense, since these entries are complex).

What is $W$? It’s the number whose $n$-th power is one:

$$W = e^{2\pi i/n} = \cos\frac{2\pi}{n} + i\sin\frac{2\pi}{n}.$$

But forget the rectangular coordinates — “it’s never good to work with the real and imaginary parts when we’re taking powers.” In the exponential form, the tenth power is instantly $e^{20\pi i/n}$. Geometrically $W$ sits on the unit circle, one $n$-th of the full way around. For $n = 6$ it’s at $60°$, and its powers march around the circle — $W^2$ at $120°$, then $W^3$, $W^4$, $W^5$, and $W^6$ coming back to $1$: the six sixth roots of one, with $W$ the primitive one.

Orthogonal columns make the inverse easy

Why is this matrix so remarkable? First: its columns are orthogonal — but only if you remember to conjugate. Try columns one and three of $F_4$ (in EE counting): $(1, i, -1, -i)$ and $(1, -i, -1, i)$. Without conjugating, the products are $1, -i^2, 1, -i^2$ — four ones, not zero, and the columns look anything but perpendicular. Nevertheless they are orthogonal. Conjugate the first vector — the $i$ becomes $-i$, the $-i$ becomes $i$ — and the inner product is $1 - 1 + 1 - 1 = 0$. Right angles in the complex sense.

The columns aren’t quite orthonormal — each has length squared $1+1+1+1 = 4$, so length $2$ — but divide by $2$ and $\frac{1}{2}F_4$ has genuinely orthonormal columns. Then, exactly like an orthogonal matrix, the inverse is immediate:

$$\left(\tfrac{1}{2}F_4\right)^{\mathrm{H}}\left(\tfrac{1}{2}F_4\right) = I \quad\Longrightarrow\quad \text{the inverse is the conjugate transpose.}$$

So anything good we learn about $F_4$, we know a similar fact about its inverse. Fourier knew that. What Fourier didn’t notice — Gauss noticed it but didn’t make a point of it — is that this matrix is so special it can be broken into pieces with lots of zeros, so that multiplying by it, or by its inverse, becomes very, very fast.

The fast Fourier transform

Here is the idea. $F_{64}$ has a neat connection to $F_{32}$, half as big ($F_6$ to $F_3$, $F_8$ to $F_4$ — same story). The link is in the $W$‘s: the $W$ for $F_{64}$ is $e^{2\pi i/64}$, one 64th of the way around the circle; square it — double the angle — and you get $e^{2\pi i/32}$, exactly the $W$ for $F_{32}$. That’s the hope, and here’s the factorization that realizes it:

$$F_{64} \;=\; \begin{bmatrix} I & D \\ I & -D \end{bmatrix} \begin{bmatrix} F_{32} & 0 \\ 0 & F_{32} \end{bmatrix} P$$

Read it right to left. $P$ is a simple odd–even permutation matrix: it takes the even-numbered components of the vector first — $x_0, x_2, x_4, \ldots$ — and then the odds — $x_1, x_3, x_5, \ldots$ That’s hard-wired in the computer, essentially instantaneous. The middle matrix does a size-32 Fourier transform separately on the even part and the odd part; its two big zero blocks are the key, because half the matrix costs nothing. The left matrix puts the pieces back together, and its only real content is $D$, the diagonal matrix

$$D = \mathrm{diag}\left(1,\, W,\, W^2,\, \ldots,\, W^{31}\right),$$

powers of $W$. Multiplying by $I$ costs nothing; multiplying by the diagonal $D$ costs just 32 multiplications. So the count drops from $64^2$ to

$$2 \times 32^2 + 32 \quad\text{(two half-size transforms, plus the fix-up)}.$$

But why stop? Recursion is the key word. Each $F_{32}$ breaks into two copies of $F_{16}$ with its own $I$-and-$D$ fix-up and its own size-32 odd–even permutation. The stacked permutations sort the vector into even-evens ($x_0, x_4, x_8, \ldots$), even-odds ($x_2, x_6, x_{10}, \ldots$), odd-evens, and odd-odds — four pieces that get multiplied separately by $F_{16}$. The $2 \times 32^2$ gives way to $2(2 \times 16^2 + 16)$, plus the original 32. Keep going: the transforms in the middle shrink to two-point, then one-point transforms — the count in the middle “goes down to zip” — while cheap factors pile up on both sides: permutations on the right, $I$‘s and $D$‘s on the left. From 64 to 32 to 16 to 8 to 4 to 2 to 1 is six steps — $\log_2 64$ — and each fix-up level costs about 32 multiplications, so the final count is $6 \times 32$: that is, $\log_2 64$ times $\tfrac{1}{2} \cdot 64$.

Feel the difference at a typical size, $n = 2^{10} = 1024$. Then $n^2$ is over a million — $1024 \times 1024$. The new count is $1024 \times \tfrac{10}{2} = 5 \times 1024$. The calculation shrinks by a factor of about 200, just from factoring the matrix properly: 200 Fourier transforms where before we could do one — and in real scientific computing, Fourier transforms are happening all the time.

And notice what the whole thing hinged on: a special matrix with orthonormal columns. Next lecture it’s back to real numbers, and to the matrices you see most in applications — positive definite matrices.

Problems

Work these with the one rule in hand — transpose means conjugate-transpose.

Problem 26.1 Complex length and inner product

Let $\mathbf{z} = \begin{bmatrix} 2 \\ 1 - i \end{bmatrix}$ and $\mathbf{y} = \begin{bmatrix} i \\ 1 + i \end{bmatrix}$ in $\mathbb{C}^2$.

(a) Compute $\mathbf{z}^T\mathbf{z}$ and $\mathbf{z}^{\mathrm{H}}\mathbf{z}$, and explain which one is the true length squared.

(b) Compute the inner product $\mathbf{y}^{\mathrm{H}}\mathbf{z}$.

Show solution

(a) Without conjugating: $\mathbf{z}^T\mathbf{z} = 2^2 + (1-i)^2 = 4 + (1 - 2i + i^2) = 4 - 2i$ — a complex number, which cannot be a length squared. With conjugating:

$$\mathbf{z}^{\mathrm{H}}\mathbf{z} = |2|^2 + |1-i|^2 = 4 + (1^2 + (-1)^2) = 6.$$

The true length squared is $6$, so $\|\mathbf{z}\| = \sqrt{6}$. Length squared must be a positive real number, and only the Hermitian form guarantees that.

(b) Conjugate $\mathbf{y}$ first: $\bar{\mathbf{y}} = (-i,\; 1 - i)$. Then

$$\mathbf{y}^{\mathrm{H}}\mathbf{z} = (-i)(2) + (1-i)(1-i) = -2i + (1 - 2i - 1) = -4i.$$

The inner product of two complex vectors is usually complex — that’s expected.

Problem 26.2 Building a Hermitian matrix

For which values of $a$, $b$, $c$ is

$$A = \begin{bmatrix} a & 4 - 2i \\ b & c \end{bmatrix}$$

a Hermitian matrix? For your answer, verify that $\mathbf{z}^{\mathrm{H}} A \mathbf{z}$ is real for $\mathbf{z} = \begin{bmatrix} 1 \\ i \end{bmatrix}$ when $a = 1$, $c = 3$.

Show solution

Hermitian means $A^{\mathrm{H}} = A$: the diagonal entries must equal their own conjugates, so $a$ and $c$ must be real (any real numbers), and the off-diagonal entries must be conjugates of each other, so $b = \overline{4 - 2i} = 4 + 2i$.

Check with $a = 1$, $c = 3$, $\mathbf{z} = (1, i)$: first

$$A\mathbf{z} = \begin{bmatrix} 1 & 4-2i \\ 4+2i & 3 \end{bmatrix}\begin{bmatrix} 1 \\ i \end{bmatrix} = \begin{bmatrix} 1 + (4-2i)i \\ 4+2i + 3i \end{bmatrix} = \begin{bmatrix} 3 + 4i \\ 4 + 5i \end{bmatrix}.$$

Then $\mathbf{z}^{\mathrm{H}}(A\mathbf{z}) = (1)(3+4i) + (-i)(4+5i) = 3 + 4i - 4i + 5 = 8$ — real, as it must be for a Hermitian matrix.

Problem 26.3 Fourier matrix by hand

Write down $W$ for $n = 4$ and use it to build $F_4$. Then verify that columns $1$ and $2$ (EE counting: the columns $(1, i, -1, -i)$ and $(1, -1, 1, -1)$) are orthogonal, and use orthogonality to write down $F_4^{-1}$.

Show solution

$W = e^{2\pi i/4} = i$, so $(F_4)_{jk} = i^{jk}$ gives

$$F_4 = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{bmatrix}.$$

Orthogonality of columns 1 and 2 — conjugate the first: $(1, -i, -1, i)$, then

$$(1)(1) + (-i)(-1) + (-1)(1) + (i)(-1) = 1 + i - 1 - i = 0. \checkmark$$

Each column has length squared $4$, so $\left(\tfrac{1}{2}F_4\right)^{\mathrm{H}}\left(\tfrac{1}{2}F_4\right) = I$, giving $F_4^{\mathrm{H}} F_4 = 4I$ and

$$F_4^{-1} = \tfrac{1}{4}F_4^{\mathrm{H}} = \tfrac{1}{4}\begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{bmatrix}.$$

The inverse transform is just the conjugate transform, scaled — replace $i$ by $-i$ and divide by $4$.

Problem 26.4 Counting the FFT's savings

(a) Using the factorization of $F_n$ into two copies of $F_{n/2}$ plus a diagonal fix-up of size $n/2$, write the multiplication count for one level of the recursion when $n = 16$.

(b) The full FFT costs $\tfrac{1}{2}n\log_2 n$ multiplications. For $n = 2^{12} = 4096$, compare this with $n^2$ and give the speedup factor.

Show solution

(a) One level replaces the $16^2 = 256$ multiplications by two transforms of size 8 plus a diagonal $D$ with 8 entries:

$$2 \times 8^2 + 8 = 128 + 8 = 136.$$

Already nearly a factor of 2, and the recursion keeps going on each $F_8$.

(b) Direct: $n^2 = 4096^2 = 16{,}777{,}216$. FFT: $\tfrac{1}{2}\, (4096)(\log_2 4096) = \tfrac{1}{2}(4096)(12) = 24{,}576$. The speedup is

$$\frac{n^2}{\tfrac{1}{2}n\log_2 n} = \frac{2n}{\log_2 n} = \frac{2(4096)}{12} \approx 683.$$

About 680 Fourier transforms for the price of one — the same effect as Strang’s factor of 200 at $n = 1024$, but larger because the gap widens as $n$ grows.

Problem 26.5 Roots of unity on the circle

Let $W = e^{2\pi i/6}$, the primitive sixth root of one. (a) Locate $W, W^2, \ldots, W^6$ on the unit circle by their angles. (b) Show that $W^2$ is the number that would serve as the "$W$" for $F_3$. (c) Compute $1 + W + W^2 + W^3 + W^4 + W^5$ and explain what this says about the columns of $F_6$.

Show solution

(a) $W$ is at $60°$; the powers double, triple, … the angle: $W^2$ at $120°$, $W^3$ at $180°$ (that’s $-1$), $W^4$ at $240°$, $W^5$ at $300°$, and $W^6$ back at $0°$, equal to $1$. Six equally spaced points — the six sixth roots of one.

(b) Squaring doubles the angle: $W^2 = e^{2\pi i \cdot 2/6} = e^{2\pi i/3}$, which is one third of the way around — exactly the primitive cube root of one that builds $F_3$. This is the connection ($W_n^2 = W_{n/2}$, here in the $6 \to 3$ form) that powers the FFT recursion.

(c) The sum is a geometric series: $\displaystyle\sum_{k=0}^{5} W^k = \frac{W^6 - 1}{W - 1} = \frac{1 - 1}{W - 1} = 0$. (Geometrically, six equally spaced unit vectors balance out.) But this sum is precisely the inner product of column 0 of $F_6$ (all ones, its own conjugate) with column 1 (the powers of $W$) — so it verifies that those two columns of the Fourier matrix are orthogonal.