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

```

\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

```{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.

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`.

and solving

$$

\textrm{det}(P - \lambda I) = (-1)^N \lambda^{N}-1=0

\textrm{det}(P - \lambda I) = (-1)^N(\lambda^N - 1)=0

$$

It can be verified that permutation matrices are orthogonal matrices:

$$

P P' = I

P P^\top = I

$$

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))

```

```{code-cell} ipython3

P4 = construct_P(4)

P4 = construct_cyclic_shift_matrix(4)

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.

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')

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

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

$$

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

```

```{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

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)

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

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):

$$

```{code-cell} ipython3

def inverse_transform(X):

def inverse_DFT(X):

N = len(X)

w = np.e ** (complex(0, 2*np.pi/N))

```

```{code-cell} ipython3

inverse_transform(X)

inverse_DFT(X)

```

Another example is

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.

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)

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

```