diff --git a/lectures/eig_circulant.md b/lectures/eig_circulant.md index 270aa6046..f0d021530 100644 --- a/lectures/eig_circulant.md +++ b/lectures/eig_circulant.md @@ -17,21 +17,23 @@ kernelspec: This lecture describes circulant matrices and some of their properties. -Circulant matrices have a special structure that connects them to useful concepts -including +Circulant matrices are useful because multiplying by them is closely connected to convolution, and their eigenvectors can be constructed using the Discrete Fourier Transform. + +We use circulant matrices to connect several useful concepts, including * convolution * Fourier transforms * permutation matrices -Because of these connections, circulant matrices are widely used in machine learning, for example, in image processing. +For background on eigenvalues and eigenvectors, see {doc}`linear_algebra`; for another use of Fourier transforms and convolution, see {doc}`hoist_failure`. + +Circulant matrices are also widely used in machine learning, for example, in image processing. We begin by importing some Python packages ```{code-cell} ipython3 import numpy as np -from numba import jit import matplotlib.pyplot as plt ``` @@ -60,14 +62,28 @@ c_{1} & c_{2} & c_{3} & c_{4} & c_{5} & \cdots & c_{0} \end{array}\right] $$ (eqn:circulant) +This pattern can be formalized as follows. + +```{prf:definition} Circulant matrix +:label: def-circulant-matrix + +An $N \times N$ matrix $C$ is **circulant** if there are numbers $c_0, \ldots, c_{N-1}$ such that + +$$ +C_{ij} = c_{(j-i) \bmod N}, +\qquad 0 \leq i,j \leq N-1. +$$ + +Equivalently, each row is obtained from the previous row by shifting entries one step to the right. +``` + It is also possible to construct a circulant matrix by creating the transpose of the above matrix, in which case only the first column needs to be specified. Let's write some Python code to generate a circulant matrix. ```{code-cell} ipython3 -@jit -def construct_cirlulant(row): +def construct_circulant(row): N = row.size @@ -83,14 +99,16 @@ def construct_cirlulant(row): ```{code-cell} ipython3 # a simple case when N = 3 -construct_cirlulant(np.array([1., 2., 3.])) +construct_circulant(np.array([1., 2., 3.])) ``` ### Some Properties of Circulant Matrices Here are some useful properties: -Suppose that $A$ and $B$ are both circulant matrices. Then it can be verified that +Suppose that $A$ and $B$ are both circulant matrices of the same order and constructed using the same cyclic shift convention. + +Then it can be verified that * The transpose of a circulant matrix is a circulant matrix. @@ -111,16 +129,18 @@ Now consider a circulant matrix with first row The **convolution** of vectors $c$ and $a$ is defined as the vector $b = c * a $ with components $$ - b_k = \sum_{i=0}^{n-1} c_{k-i} a_i + b_k = \sum_{i=0}^{N-1} c_{k-i} a_i $$ (eqn:conv) +Here and below, indices such as $k-i$ are interpreted modulo $N$. + We use $*$ to denote **convolution** via the calculation described in equation {eq}`eqn:conv`. It can be verified that the vector $b$ satisfies -$$ b = C^T a $$ +$$ b = C^\top a $$ -where $C^T$ is the transpose of the circulant matrix defined in equation {eq}`eqn:circulant`. +where $C^\top$ is the transpose of the circulant matrix defined in equation {eq}`eqn:circulant`. @@ -176,7 +196,7 @@ $$ and solving $$ -\textrm{det}(P - \lambda I) = (-1)^N \lambda^{N}-1=0 +\textrm{det}(P - \lambda I) = (-1)^N(\lambda^N - 1)=0 $$ @@ -189,7 +209,7 @@ Thus, **singular values** of the permutation matrix $P$ defined in equation {eq It can be verified that permutation matrices are orthogonal matrices: $$ -P P' = I +P P^\top = I $$ @@ -200,8 +220,7 @@ $$ Let's write some Python code to illustrate these ideas. ```{code-cell} ipython3 -@jit -def construct_P(N): +def construct_cyclic_shift_matrix(N): P = np.zeros((N, N)) @@ -213,7 +232,7 @@ def construct_P(N): ``` ```{code-cell} ipython3 -P4 = construct_P(4) +P4 = construct_cyclic_shift_matrix(4) P4 ``` @@ -224,7 +243,7 @@ P4 ```{code-cell} ipython3 for i in range(4): - print(f'𝜆{i} = {𝜆[i]:.1f} \nvec{i} = {Q[i, :]}\n') + print(f'𝜆{i} = {𝜆[i]:.1f} \nvec{i} = {Q[:, i]}\n') ``` In graphs below, we shall portray eigenvalues of a shift permutation matrix in the complex plane. @@ -249,7 +268,7 @@ for i, N in enumerate([3, 4, 6, 8]): row_i = i // 2 col_i = i % 2 - P = construct_P(N) + P = construct_cyclic_shift_matrix(N) 𝜆, Q = np.linalg.eig(P) circ = plt.Circle((0, 0), radius=1, edgecolor='b', facecolor='None') @@ -287,7 +306,7 @@ F_{8}=\left[\begin{array}{ccccc} \end{array}\right] $$ -The matrix $F_8$ defines a [Discete Fourier Transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform). +The matrix $F_8$ defines a [Discrete Fourier Transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform). To convert it into an orthogonal eigenvector matrix, we can simply normalize it by dividing every entry by $\sqrt{8}$. @@ -332,7 +351,7 @@ Q8 @ np.conjugate(Q8) Let's verify that $k$th column of $Q_{8}$ is an eigenvector of $P_{8}$ with an eigenvalue $w^{k}$. ```{code-cell} ipython3 -P8 = construct_P(8) +P8 = construct_cyclic_shift_matrix(8) ``` ```{code-cell} ipython3 @@ -353,15 +372,11 @@ Next, we execute calculations to verify that the circulant matrix $C$ defined i $$ -C = c_{0} I + c_{1} P + \cdots + c_{n-1} P^{n-1} +C = c_{0} I + c_{1} P + \cdots + c_{N-1} P^{N-1} $$ and that every eigenvector of $P$ is also an eigenvector of $C$. -```{code-cell} ipython3 - -``` - We illustrate this for $N=8$ case. ```{code-cell} ipython3 @@ -373,10 +388,10 @@ c ``` ```{code-cell} ipython3 -C8 = construct_cirlulant(c) +C8 = construct_circulant(c) ``` -Compute $c_{0} I + c_{1} P + \cdots + c_{n-1} P^{n-1}$. +Compute $c_{0} I + c_{1} P + \cdots + c_{N-1} P^{N-1}$. ```{code-cell} ipython3 N = 8 @@ -403,7 +418,7 @@ Now let's compute the difference between two circulant matrices that we have co np.abs(C - C8).max() ``` -The $k$th column of $P_{8}$ associated with eigenvalue $w^{k-1}$ is an eigenvector of $C_{8}$ associated with an eigenvalue $\sum_{h=0}^{7} c_{j} w^{h k}$. +The $j$th column of $Q_{8}$ is an eigenvector of $C_{8}$ associated with eigenvalue $\sum_{k=0}^{7} c_k w^{j k}$. ```{code-cell} ipython3 𝜆_C8 = np.zeros(8, dtype=complex) @@ -430,7 +445,7 @@ for j in range(8): The **Discrete Fourier Transform** (DFT) allows us to represent a discrete time sequence as a weighted sum of complex sinusoids. -Consider a sequence of $N$ real number $\{x_j\}_{j=0}^{N-1}$. +Consider a sequence of $N$ real numbers $\{x_j\}_{j=0}^{N-1}$. The **Discrete Fourier Transform** maps $\{x_j\}_{j=0}^{N-1}$ into a sequence of complex numbers $\{X_k\}_{k=0}^{N-1}$ @@ -498,7 +513,7 @@ def plot_magnitude(x=None, X=None): if (X is not None): data.append(X) names.append('X') - xs.append('j') + xs.append('k') num = len(data) for i in range(num): @@ -526,7 +541,7 @@ x_{n} = \sum_{k=0}^{N-1} \frac{1}{N} X_{k} e^{2\pi\left(\frac{kn}{N}\right)i}, \ $$ ```{code-cell} ipython3 -def inverse_transform(X): +def inverse_DFT(X): N = len(X) w = np.e ** (complex(0, 2*np.pi/N)) @@ -540,7 +555,7 @@ def inverse_transform(X): ``` ```{code-cell} ipython3 -inverse_transform(X) +inverse_DFT(X) ``` Another example is @@ -551,7 +566,7 @@ $$ Since $N=20$, we cannot use an integer multiple of $\frac{1}{20}$ to represent a frequency $\frac{11}{40}$. -To handle this, we shall end up using all $N$ of the availble frequencies in the DFT. +To handle this, we shall end up using all $N$ of the available frequencies in the DFT. Since $\frac{11}{40}$ is in between $\frac{10}{40}$ and $\frac{12}{40}$ (each of which is an integer multiple of $\frac{1}{20}$), the complex coefficients in the DFT have their largest magnitudes at $k=5,6,15,16$, not just at a single frequency. @@ -614,7 +629,7 @@ X = DFT(x) X ``` -Now let's evaluate the outcome of postmultiplying the eigenvector matrix $F_{20}$ by the vector $x$, a product that we claim should equal the Fourier tranform of the sequence $\{x_n\}_{n=0}^{N-1}$. +Now let's evaluate the outcome of postmultiplying the eigenvector matrix $F_{20}$ by the vector $x$, a product that we claim should equal the Fourier transform of the sequence $\{x_n\}_{n=0}^{N-1}$. ```{code-cell} ipython3 F20, _ = construct_F(20) @@ -624,13 +639,9 @@ F20, _ = construct_F(20) F20 @ x ``` -Similarly, the inverse DFT can be expressed as a inverse DFT matrix $F^{-1}_{20}$. +Similarly, the inverse DFT can be expressed as an inverse DFT matrix $F^{-1}_{20}$. ```{code-cell} ipython3 F20_inv = np.linalg.inv(F20) F20_inv @ X ``` - -```{code-cell} ipython3 - -```