{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d32fcbe2-702c-4f06-84d2-17de984db13e",
   "metadata": {},
   "source": [
    "# Eigen-decomposition of a Matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13b9b54c-2ec3-431d-970c-62e0e18a78c7",
   "metadata": {},
   "source": [
    "In many applications, we need to find the eigenvalues and eigenvectors of a matrix. This can be done using the NumPy linear algebra module `numpy.linalg`. For example, let us consider a random matrix $\\mathbf{A}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c23ad13a-3958-442e-a1da-504af3a117ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(suppress=True)\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b8a86a50-6cb8-428a-8ff9-34834207ab42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.33222557 -0.63034335  0.07890825]\n",
      " [ 1.35934868 -0.3405165  -0.49860337]\n",
      " [-0.19721574  0.19602863 -0.83804111]]\n"
     ]
    }
   ],
   "source": [
    "A = np.random.randn(3,3)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8806eb0-bf51-40c4-922f-35b31556da5a",
   "metadata": {},
   "source": [
    "To find its eigenvalues and eigenvectors, use the function `numpy.linalg.eig` as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "d5e1c78c-ec05-46f2-9946-ff86c98cc42e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigenvalues = [-0.03188706+0.90169563j -0.03188706-0.90169563j -0.78255793+0.j        ]\n",
      "eigenvectors = \n",
      "[[-0.19177332-0.50350317j -0.19177332+0.50350317j  0.18268415+0.j        ]\n",
      " [-0.83200937+0.j         -0.83200937-0.j          0.43354694+0.j        ]\n",
      " [-0.00783031+0.13193422j -0.00783031-0.13193422j  0.88241915+0.j        ]]\n"
     ]
    }
   ],
   "source": [
    "w, v = np.linalg.eig(A)    # find eigenvalues and eigenvectors\n",
    "print(f'eigenvalues = {w}')    # each entry is an eigenvalue\n",
    "print(f'eigenvectors = \\n{v}')    # each column is the corresponding eigenvector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce2d31be-6a72-4496-ae07-07986ee6fe5d",
   "metadata": {},
   "source": [
    "Each entry of `w` is an eigenvalue, and each column `v[:,i]` is the eigenvector corresponding to the eigenvalue `w[i]`. Note that even for a real matrix, the eigenvalues and eigenvectors can have complex numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebbfcb24-38f1-40b3-8f0a-4571536597d9",
   "metadata": {},
   "source": [
    "Let $\\mathbf{D}$ be a diagonal matrix formed by the eigenvalues, and $\\mathbf{R}$ be a matrix whose columns are the corresponding eigenvectors. We can check that $\\mathbf{A}$ can be decomposed as $\\mathbf{A} = \\mathbf{R} \\cdot \\mathbf{D} \\cdot \\mathbf{R}^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "df1c3eac-c5ea-4de9-8ab1-518dfc1b9760",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.33222557+0.j -0.63034335+0.j  0.07890825-0.j]\n",
      " [ 1.35934868-0.j -0.3405165 -0.j -0.49860337-0.j]\n",
      " [-0.19721574-0.j  0.19602863-0.j -0.83804111+0.j]]\n"
     ]
    }
   ],
   "source": [
    "D = np.diag(w)    # convert a vector to a diagonal matrix\n",
    "R = v\n",
    "A_decomp = np.dot(R, np.dot(D, np.linalg.inv(R)))    # decompose A into eigenmodes\n",
    "print(A_decomp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50404e62-cdab-4747-9ef3-09da6a9dc7e8",
   "metadata": {},
   "source": [
    "In the special (but common) case where a matrix $\\mathbf{A}$ is Hermitian (or symmetric if $\\mathbf{A}$ is real), the eigenvalues will be all real. In this case, it is convenient to use the function `numpy.linalg.eigh` instead, which returns the eigenvalues in ascending order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "483c48f8-e943-44bd-812a-29506f78b5b3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.33222557  0.36450267 -0.05915374]\n",
      " [ 0.36450267 -0.3405165  -0.15128737]\n",
      " [-0.05915374 -0.15128737 -0.83804111]]\n"
     ]
    }
   ],
   "source": [
    "As = (A + A.T) / 2    # construct a real symmetric matrix\n",
    "print(As)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7b1d6971-398f-48fa-820a-ef23d45df661",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eigenvalues = [-0.88217369 -0.46595103  0.50179268]\n",
      "eigenvectors = \n",
      "[[-0.04201506 -0.41556531  0.90859243]\n",
      " [ 0.29489901  0.86372198  0.40867949]\n",
      " [ 0.95460427 -0.2851137  -0.08626047]]\n"
     ]
    }
   ],
   "source": [
    "w, v = np.linalg.eigh(As)    # find eigenvalues and eigenvectors of a Hermitian matrix\n",
    "print(f'eigenvalues = {w}')    # eigenvalues in ascending order\n",
    "print(f'eigenvectors = \\n{v}')    # each column is the corresponding eigenvector"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a344d3e1-87b0-4dfb-b587-5af32e9da5f1",
   "metadata": {},
   "source": [
    "Since $\\mathbf{A}$ is real symmetric, the eigenvectors are also real, and we have $\\mathbf{R}^{-1} = \\mathbf{R}^\\top$, i.e., $\\mathbf{R}$ is orthogonal. Therefore, we should have $\\mathbf{A} = \\mathbf{R} \\cdot \\mathbf{D} \\cdot \\mathbf{R}^\\top$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "232da084-b69c-4783-b96a-f84669efe47f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.33222557  0.36450267 -0.05915374]\n",
      " [ 0.36450267 -0.3405165  -0.15128737]\n",
      " [-0.05915374 -0.15128737 -0.83804111]]\n"
     ]
    }
   ],
   "source": [
    "D = np.diag(w)\n",
    "R = v\n",
    "A_decomp = np.dot(R, np.dot(D, R.T))\n",
    "print(A_decomp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e731e96f-80eb-4574-8d9f-f1c647ca10de",
   "metadata": {},
   "source": [
    "Notice that here we used the transpose `R.T` instead of the inverse `np.linalg.inv(R)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4317b0fa-fb87-4607-b0c9-8c340d639ac4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}