eig_circulant.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8cb1b4cf",
   "metadata": {},
   "source": [
    "# Circulant Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fdc98fc",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture describes circulant matrices and some of their properties.\n",
    "\n",
    "Circulant matrices have a special structure that connects them to  useful concepts\n",
    "including\n",
    "\n",
    "- convolution  \n",
    "- Fourier transforms  \n",
    "- permutation matrices  \n",
    "\n",
    "\n",
    "Because of these connections, circulant matrices are widely used  in machine learning, for example, in image processing.\n",
    "\n",
    "We begin by importing some Python packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f45c1ff0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numba import njit\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e695529",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91879243",
   "metadata": {},
   "source": [
    "## Constructing a Circulant Matrix\n",
    "\n",
    "To construct an $ N \\times N $ circulant matrix, we  need only the first row, say,\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} c_{0} & c_{1} & c_{2} & c_{3} & c_{4} & \\cdots & c_{N-1} \\end{bmatrix} .\n",
    "$$\n",
    "\n",
    "After setting entries in the first row, the remaining rows of a circulant matrix are determined as\n",
    "follows:\n",
    "\n",
    "\n",
    "<a id='equation-eqn-circulant'></a>\n",
    "$$\n",
    "C=\\left[\\begin{array}{ccccccc}\n",
    "c_{0} & c_{1} & c_{2} & c_{3} & c_{4} & \\cdots & c_{N-1}\\\\\n",
    "c_{N-1} & c_{0} & c_{1} & c_{2} & c_{3} & \\cdots & c_{N-2}\\\\\n",
    "c_{N-2} & c_{N-1} & c_{0} & c_{1} & c_{2} & \\cdots & c_{N-3}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "c_{3} & c_{4} & c_{5} & c_{6} & c_{7} & \\cdots & c_{2}\\\\\n",
    "c_{2} & c_{3} & c_{4} & c_{5} & c_{6} & \\cdots & c_{1}\\\\\n",
    "c_{1} & c_{2} & c_{3} & c_{4} & c_{5} & \\cdots & c_{0}\n",
    "\\end{array}\\right] \\tag{3.1}\n",
    "$$\n",
    "\n",
    "It is also possible to construct a circulant matrix by creating the transpose of the above matrix, in which case only the\n",
    "first column needs to be specified.\n",
    "\n",
    "Let’s write some Python code to generate a circulant matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d08beca",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def construct_cirlulant(row):\n",
    "\n",
    "    N = row.size\n",
    "\n",
    "    C = np.empty((N, N))\n",
    "\n",
    "    for i in range(N):\n",
    "\n",
    "        C[i, i:] = row[:N-i]\n",
    "        C[i, :i] = row[N-i:]\n",
    "\n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "953e87ea",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# a simple case when N = 3\n",
    "construct_cirlulant(np.array([1., 2., 3.]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f71dd5c",
   "metadata": {},
   "source": [
    "### Some Properties of Circulant Matrices\n",
    "\n",
    "Here are some useful properties:\n",
    "\n",
    "Suppose that $ A $ and $ B $ are both circulant matrices. Then it can be verified that\n",
    "\n",
    "- The transpose of a circulant matrix is a circulant matrix.  \n",
    "- $ A + B $ is a circulant matrix  \n",
    "- $ A B $ is a circulant matrix  \n",
    "- $ A B = B A $  \n",
    "\n",
    "\n",
    "Now consider a circulant matrix with first row\n",
    "\n",
    "$$\n",
    "c = \\begin{bmatrix} c_0 & c_1 & \\cdots & c_{N-1} \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and consider a vector\n",
    "\n",
    "$$\n",
    "a = \\begin{bmatrix} a_0 & a_1 & \\cdots  &  a_{N-1} \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The **convolution** of  vectors $ c $ and $ a $ is defined   as the vector $ b = c * a $  with components\n",
    "\n",
    "\n",
    "<a id='equation-eqn-conv'></a>\n",
    "$$\n",
    "b_k = \\sum_{i=0}^{n-1} c_{k-i} a_i \\tag{3.2}\n",
    "$$\n",
    "\n",
    "We use $ * $ to denote **convolution** via the calculation described in equation [(3.2)](#equation-eqn-conv).\n",
    "\n",
    "It can be verified that the vector $ b $ satisfies\n",
    "\n",
    "$$\n",
    "b = C^T a\n",
    "$$\n",
    "\n",
    "where $ C^T $ is the transpose of the circulant matrix  defined in equation [(3.1)](#equation-eqn-circulant)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75f90572",
   "metadata": {},
   "source": [
    "## Connection to Permutation Matrix\n",
    "\n",
    "A good way to construct a circulant matrix is to use a **permutation matrix**.\n",
    "\n",
    "Before defining a permutation **matrix**, we’ll define a **permutation**.\n",
    "\n",
    "A **permutation** of a set of the set of non-negative integers $ \\{0, 1, 2, \\ldots \\} $ is a one-to-one mapping of the set into itself.\n",
    "\n",
    "A permutation of a set $ \\{1, 2, \\ldots, n\\} $ rearranges the $ n $ integers in the set.\n",
    "\n",
    "A [permutation matrix](https://mathworld.wolfram.com/PermutationMatrix.html) is obtained by permuting the rows of an $ n \\times n $ identity matrix according to a permutation of the numbers $ 1 $ to $ n $.\n",
    "\n",
    "Thus, every row and every column contain precisely a single $ 1 $ with $ 0 $ everywhere else.\n",
    "\n",
    "Every permutation corresponds to a unique permutation matrix.\n",
    "\n",
    "For example, the $ N \\times N $ matrix\n",
    "\n",
    "\n",
    "<a id='equation-eqn-examplep'></a>\n",
    "$$\n",
    "P=\\left[\\begin{array}{cccccc}\n",
    "0 & 1 & 0 & 0 & \\cdots & 0\\\\\n",
    "0 & 0 & 1 & 0 & \\cdots & 0\\\\\n",
    "0 & 0 & 0 & 1 & \\cdots & 0\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "0 & 0 & 0 & 0 & \\cdots & 1\\\\\n",
    "1 & 0 & 0 & 0 & \\cdots & 0\n",
    "\\end{array}\\right] \\tag{3.3}\n",
    "$$\n",
    "\n",
    "serves as  a **cyclic shift**  operator that, when applied to an $ N \\times 1 $ vector $ h $, shifts entries in rows $ 2 $ through $ N $ up one row and shifts the entry in row $ 1 $ to row $ N $.\n",
    "\n",
    "Eigenvalues of  the cyclic shift permutation matrix $ P $ defined in equation [(3.3)](#equation-eqn-examplep) can be computed  by constructing\n",
    "\n",
    "$$\n",
    "P-\\lambda I=\\left[\\begin{array}{cccccc}\n",
    "-\\lambda & 1 & 0 & 0 & \\cdots & 0\\\\\n",
    "0 & -\\lambda & 1 & 0 & \\cdots & 0\\\\\n",
    "0 & 0 & -\\lambda & 1 & \\cdots & 0\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "0 & 0 & 0 & 0 & \\cdots & 1\\\\\n",
    "1 & 0 & 0 & 0 & \\cdots & -\\lambda\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "and solving\n",
    "\n",
    "$$\n",
    "\\textrm{det}(P - \\lambda I) = (-1)^N \\lambda^{N}-1=0\n",
    "$$\n",
    "\n",
    "Eigenvalues $ \\lambda_i $  can be complex.\n",
    "\n",
    "Magnitudes $ \\mid \\lambda_i \\mid $  of these  eigenvalues $ \\lambda_i $ all equal  $ 1 $.\n",
    "\n",
    "Thus, **singular values** of the  permutation matrix $ P $ defined in equation [(3.3)](#equation-eqn-examplep) all equal $ 1 $.\n",
    "\n",
    "It can be verified that permutation matrices are orthogonal matrices:\n",
    "\n",
    "$$\n",
    "P P' = I\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab643157",
   "metadata": {},
   "source": [
    "## Examples with Python\n",
    "\n",
    "Let’s write some Python code to illustrate these ideas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "928a17bf",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit\n",
    "def construct_P(N):\n",
    "\n",
    "    P = np.zeros((N, N))\n",
    "\n",
    "    for i in range(N-1):\n",
    "        P[i, i+1] = 1\n",
    "    P[-1, 0] = 1\n",
    "\n",
    "    return P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62d753d1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P4 = construct_P(4)\n",
    "P4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9c3c1f4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# compute the eigenvalues and eigenvectors\n",
    "𝜆, Q = np.linalg.eig(P4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "569108c0",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "for i in range(4):\n",
    "    print(f'𝜆{i} = {𝜆[i]:.1f} \\nvec{i} = {Q[i, :]}\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "004b21dd",
   "metadata": {},
   "source": [
    "In graphs  below, we shall portray eigenvalues of a shift  permutation matrix   in the complex plane.\n",
    "\n",
    "These eigenvalues are uniformly distributed along the unit circle.\n",
    "\n",
    "They are the **$ n $ roots of unity**, meaning they are the $ n $  numbers  $ z $  that solve $ z^n =1 $, where $ z $ is a complex number.\n",
    "\n",
    "In particular, the $ n $ roots of unity are\n",
    "\n",
    "$$\n",
    "z = \\exp\\left(\\frac{2 \\pi j k }{N} \\right) , \\quad k = 0, \\ldots, N-1\n",
    "$$\n",
    "\n",
    "where $ j $ denotes the purely imaginary unit number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50d9f0fd",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2, 2, figsize=(10, 10))\n",
    "\n",
    "for i, N in enumerate([3, 4, 6, 8]):\n",
    "\n",
    "    row_i = i // 2\n",
    "    col_i = i % 2\n",
    "\n",
    "    P = construct_P(N)\n",
    "    𝜆, Q = np.linalg.eig(P)\n",
    "\n",
    "    circ = plt.Circle((0, 0), radius=1, edgecolor='b', facecolor='None')\n",
    "    ax[row_i, col_i].add_patch(circ)\n",
    "\n",
    "    for j in range(N):\n",
    "        ax[row_i, col_i].scatter(𝜆[j].real, 𝜆[j].imag, c='b')\n",
    "\n",
    "    ax[row_i, col_i].set_title(f'N = {N}')\n",
    "    ax[row_i, col_i].set_xlabel('real')\n",
    "    ax[row_i, col_i].set_ylabel('imaginary')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "875c5012",
   "metadata": {},
   "source": [
    "For a vector of  coefficients $ \\{c_i\\}_{i=0}^{n-1} $, eigenvectors of $ P $ are also  eigenvectors of\n",
    "\n",
    "$$\n",
    "C = c_{0} I + c_{1} P + c_{2} P^{2} +\\cdots + c_{N-1} P^{N-1}.\n",
    "$$\n",
    "\n",
    "Consider an example in which  $ N=8 $ and let $ w = e^{-2 \\pi j / N} $.\n",
    "\n",
    "It can be verified that the matrix $ F_8 $ of eigenvectors of $ P_{8} $  is\n",
    "\n",
    "$$\n",
    "F_{8}=\\left[\\begin{array}{ccccc}\n",
    "1 & 1 & 1 & \\cdots & 1\\\\\n",
    "1 & w & w^{2} & \\cdots & w^{7}\\\\\n",
    "1 & w^{2} & w^{4} & \\cdots & w^{14}\\\\\n",
    "1 & w^{3} & w^{6} & \\cdots & w^{21}\\\\\n",
    "1 & w^{4} & w^{8} & \\cdots & w^{28}\\\\\n",
    "1 & w^{5} & w^{10} & \\cdots & w^{35}\\\\\n",
    "1 & w^{6} & w^{12} & \\cdots & w^{42}\\\\\n",
    "1 & w^{7} & w^{14} & \\cdots & w^{49}\n",
    "\\end{array}\\right]\n",
    "$$\n",
    "\n",
    "The matrix $ F_8 $ defines a  [Discete Fourier Transform](https://en.wikipedia.org/wiki/Discrete_Fourier_transform).\n",
    "\n",
    "To convert it into an orthogonal eigenvector matrix, we can simply normalize it by dividing every entry  by $ \\sqrt{8} $.\n",
    "\n",
    "- stare at the first column of $ F_8 $ above to convince yourself of this fact  \n",
    "\n",
    "\n",
    "The eigenvalues corresponding to each eigenvector are $ \\{w^{j}\\}_{j=0}^{7} $ in order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e227219",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def construct_F(N):\n",
    "\n",
    "    w = np.e ** (-complex(0, 2*np.pi/N))\n",
    "\n",
    "    F = np.ones((N, N), dtype=complex)\n",
    "    for i in range(1, N):\n",
    "        F[i, 1:] = w ** (i * np.arange(1, N))\n",
    "\n",
    "    return F, w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa5813ba",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F8, w = construct_F(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c149c6d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fc58ac3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60f43aff",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# normalize\n",
    "Q8 = F8 / np.sqrt(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d736ac17",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# verify the orthogonality (unitarity)\n",
    "Q8 @ np.conjugate(Q8)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2faa71cb",
   "metadata": {},
   "source": [
    "Let’s verify that $ k $th column of $ Q_{8} $ is an eigenvector of $ P_{8} $ with an eigenvalue $ w^{k} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3151b72f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P8 = construct_P(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7035c229",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "diff_arr = np.empty(8, dtype=complex)\n",
    "for j in range(8):\n",
    "    diff = P8 @ Q8[:, j] - w ** j * Q8[:, j]\n",
    "    diff_arr[j] = diff @ diff.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7cef87c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "diff_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5717925b",
   "metadata": {},
   "source": [
    "## Associated Permutation Matrix\n",
    "\n",
    "Next, we execute calculations to verify that the circulant matrix $ C $ defined  in equation [(3.1)](#equation-eqn-circulant) can be written as\n",
    "\n",
    "$$\n",
    "C = c_{0} I + c_{1} P + \\cdots + c_{n-1} P^{n-1}\n",
    "$$\n",
    "\n",
    "and that every eigenvector of $ P $ is also an eigenvector of $ C $.\n",
    "\n",
    "We illustrate this for $ N=8 $ case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a42f5555",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "c = np.random.random(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b30dfae8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "123ffd22",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "C8 = construct_cirlulant(c)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4645739f",
   "metadata": {},
   "source": [
    "Compute $ c_{0} I + c_{1} P + \\cdots + c_{n-1} P^{n-1} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8ecf854",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N = 8\n",
    "\n",
    "C = np.zeros((N, N))\n",
    "P = np.eye(N)\n",
    "\n",
    "for i in range(N):\n",
    "    C += c[i] * P\n",
    "    P = P8 @ P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f20b4bb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bae6740",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "C8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c03fa699",
   "metadata": {},
   "source": [
    "Now let’s compute the difference between two circulant matrices that we have  constructed in two different ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df823138",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.abs(C - C8).max()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a18ad713",
   "metadata": {},
   "source": [
    "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} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "692e325a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "𝜆_C8 = np.zeros(8, dtype=complex)\n",
    "\n",
    "for j in range(8):\n",
    "    for k in range(8):\n",
    "        𝜆_C8[j] += c[k] * w ** (j * k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be6c9d92",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "𝜆_C8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42f2691d",
   "metadata": {},
   "source": [
    "We can verify this by comparing `C8 @ Q8[:, j]` with `𝜆_C8[j] * Q8[:, j]`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59653474",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# verify\n",
    "for j in range(8):\n",
    "    diff = C8 @ Q8[:, j] - 𝜆_C8[j] * Q8[:, j]\n",
    "    print(diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2777ffad",
   "metadata": {},
   "source": [
    "## Discrete Fourier Transform\n",
    "\n",
    "The **Discrete Fourier Transform** (DFT) allows us to  represent a  discrete time sequence as a weighted sum of complex sinusoids.\n",
    "\n",
    "Consider a sequence of $ N $ real number $ \\{x_j\\}_{j=0}^{N-1} $.\n",
    "\n",
    "The **Discrete Fourier Transform** maps $ \\{x_j\\}_{j=0}^{N-1} $ into a sequence of complex numbers $ \\{X_k\\}_{k=0}^{N-1} $\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "X_{k}=\\sum_{n=0}^{N-1}x_{n}e^{-2\\pi\\frac{kn}{N}i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2eed53b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def DFT(x):\n",
    "    \"The discrete Fourier transform.\"\n",
    "\n",
    "    N = len(x)\n",
    "    w = np.e ** (-complex(0, 2*np.pi/N))\n",
    "\n",
    "    X = np.zeros(N, dtype=complex)\n",
    "    for k in range(N):\n",
    "        for n in range(N):\n",
    "            X[k] += x[n] * w ** (k * n)\n",
    "\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f36d3ea2",
   "metadata": {},
   "source": [
    "Consider the following example.\n",
    "\n",
    "$$\n",
    "x_{n}=\\begin{cases}\n",
    "1/2 & n=0,1\\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79dd47bc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x = np.zeros(10)\n",
    "x[0:2] = 1/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9a7f5ae",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d722d0b1",
   "metadata": {},
   "source": [
    "Apply a discrete Fourier transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d64579ad",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = DFT(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bc3a1a6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aee7a158",
   "metadata": {},
   "source": [
    "We can plot  magnitudes of a sequence of numbers and the  associated discrete Fourier transform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32344277",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_magnitude(x=None, X=None):\n",
    "\n",
    "    data = []\n",
    "    names = []\n",
    "    xs = []\n",
    "    if (x is not None):\n",
    "        data.append(x)\n",
    "        names.append('x')\n",
    "        xs.append('n')\n",
    "    if (X is not None):\n",
    "        data.append(X)\n",
    "        names.append('X')\n",
    "        xs.append('j')\n",
    "\n",
    "    num = len(data)\n",
    "    for i in range(num):\n",
    "        n = data[i].size\n",
    "        plt.figure(figsize=(8, 3))\n",
    "        plt.scatter(range(n), np.abs(data[i]))\n",
    "        plt.vlines(range(n), 0, np.abs(data[i]), color='b')\n",
    "\n",
    "        plt.xlabel(xs[i])\n",
    "        plt.ylabel('magnitude')\n",
    "        plt.title(names[i])\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5febc35",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_magnitude(x=x, X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57ac2cd0",
   "metadata": {},
   "source": [
    "The **inverse Fourier transform**  transforms a Fourier transform  $ X $ of $ x $  back to $ x $.\n",
    "\n",
    "The inverse Fourier transform is defined as\n",
    "\n",
    "$$\n",
    "x_{n} = \\sum_{k=0}^{N-1} \\frac{1}{N} X_{k} e^{2\\pi\\left(\\frac{kn}{N}\\right)i}, \\quad n=0, 1, \\ldots, N-1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79866591",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def inverse_transform(X):\n",
    "\n",
    "    N = len(X)\n",
    "    w = np.e ** (complex(0, 2*np.pi/N))\n",
    "\n",
    "    x = np.zeros(N, dtype=complex)\n",
    "    for n in range(N):\n",
    "        for k in range(N):\n",
    "            x[n] += X[k] * w ** (k * n) / N\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3411831",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "inverse_transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6d1f68f",
   "metadata": {},
   "source": [
    "Another example is\n",
    "\n",
    "$$\n",
    "x_{n}=2\\cos\\left(2\\pi\\frac{11}{40}n\\right),\\ n=0,1,2,\\cdots19\n",
    "$$\n",
    "\n",
    "Since $ N=20 $, we cannot use an integer multiple of $ \\frac{1}{20} $ to represent a frequency $ \\frac{11}{40} $.\n",
    "\n",
    "To handle this,  we shall end up using all $ N $ of the availble   frequencies in the DFT.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be2bb68b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N = 20\n",
    "x = np.empty(N)\n",
    "\n",
    "for j in range(N):\n",
    "    x[j] = 2 * np.cos(2 * np.pi * 11 * j / 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7403a0d9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = DFT(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df7db29b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_magnitude(x=x, X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98945fee",
   "metadata": {},
   "source": [
    "What happens if we change the last example to $ x_{n}=2\\cos\\left(2\\pi\\frac{10}{40}n\\right) $?\n",
    "\n",
    "Note that $ \\frac{10}{40} $ is an integer multiple of $ \\frac{1}{20} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c21f065",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N = 20\n",
    "x = np.empty(N)\n",
    "\n",
    "for j in range(N):\n",
    "    x[j] = 2 * np.cos(2 * np.pi * 10 * j / 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "683c29d7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = DFT(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a723a91c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_magnitude(x=x, X=X)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd1e6e93",
   "metadata": {},
   "source": [
    "If we represent the discrete Fourier transform as a matrix, we discover that it equals the  matrix $ F_{N} $ of eigenvectors  of the permutation matrix $ P_{N} $.\n",
    "\n",
    "We can use the example where $ x_{n}=2\\cos\\left(2\\pi\\frac{11}{40}n\\right),\\ n=0,1,2,\\cdots19 $ to illustrate this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49a120df",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N = 20\n",
    "x = np.empty(N)\n",
    "\n",
    "for j in range(N):\n",
    "    x[j] = 2 * np.cos(2 * np.pi * 11 * j / 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb7d08b9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a2684bd",
   "metadata": {},
   "source": [
    "First use the summation formula to transform $ x $ to $ X $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07240c43",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "X = DFT(x)\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cde58d2",
   "metadata": {},
   "source": [
    "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} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13c6d8d9",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F20, _ = construct_F(20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e568a7f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F20 @ x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88ea05d2",
   "metadata": {},
   "source": [
    "Similarly, the inverse DFT can be expressed as a inverse DFT matrix $ F^{-1}_{20} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef7f9931",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F20_inv = np.linalg.inv(F20)\n",
    "F20_inv @ X"
   ]
  }
 ],
 "metadata": {
  "date": 1718249790.5533202,
  "filename": "eig_circulant.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Circulant Matrices"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}