Eigen analysis and reconstruction in Python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performing eigenvalue analysis/reconstruction in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are new to Python, this jumps in a bit quickly. I'd suggest looking at Numpy for Matlab Users before reading this. The following is simply a brief demonstration of using Python 3.5 for eigenanalysis. A few things to note:\n",
    "\n",
    "1. If you use an earlier version of Python (than 3.5), the @ operator isn't defined yet. To perform matrix multiplication, you neet to type use numpy.dot(A,B) in place of A@B. That will require you to start with import numpy before doing any of this.\n",
    "\n",
    "2. I demonstrate for a non-symmetric matrix. For a symmetric matrix, you should use eigh. The H in eigh means Hermetian, and is a more general definition than symmetric. If you don't deal with complex valued matrices, it's irrelevant to you. Using eigh leverages the properties of Hermetian matrices in the solution process resulting in potentially faster and more accurate results than the more general eig code. For a small matrix, this is irrelevant, but it becomes important for more substantial calculations.\n",
    "3. I show a couple of tips later that matter for later matrices (avoiding the inverse). Be aware of them. Why may or may not matter to you, but when you get to big or sensitive problems, they do.\n",
    "\n",
    "The first thing I need to do is import a couple of tools (scipy, and its linear algebra package)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to attempt to solve for\n",
    "\n",
    "$$A r = r v$$\n",
    "\n",
    "where A is the matrix, and $r$ represents the right eigenvectors, while\n",
    "$v$ represents the eigenvalues. We are also going to obtain the left\n",
    "eigenvalues as well, for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.46454727 -0.88290596  0.40824829]\n",
      " [-0.57079553 -0.23952042 -0.81649658]\n",
      " [-0.67704379  0.40386512  0.40824829]]\n",
      "[[  1.61168440e+01+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j]\n",
      " [  0.00000000e+00+0.j  -1.11684397e+00+0.j   0.00000000e+00+0.j]\n",
      " [  0.00000000e+00+0.j   0.00000000e+00+0.j  -9.75918483e-16+0.j]]\n",
      "[[-0.23197069 -0.78583024  0.40824829]\n",
      " [-0.52532209 -0.08675134 -0.81649658]\n",
      " [-0.8186735   0.61232756  0.40824829]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[1,2,3],[4,5,6],[7,8,9]])  # Defining the array (I'm avoiding using the matrix class)\n",
    "# It's a personal preference, and I'm still not locked into it myself. \n",
    "(v, l, r) = la.eig(A, left = True) # You can read the help, buy the left eigenvectors don't get created without this. \n",
    "v = np.diag(v)  # by default, eig puts the eigenvalues in a 1-D array. We will need a diagonal matrix in a moment.\n",
    "print(l)\n",
    "print(v)\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.61168440e+01+0.j,   0.00000000e+00+0.j,   0.00000000e+00+0.j],\n",
       "       [  0.00000000e+00+0.j,  -1.11684397e+00+0.j,   0.00000000e+00+0.j],\n",
       "       [  0.00000000e+00+0.j,   0.00000000e+00+0.j,  -9.75918483e-16+0.j]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These should be identical based on the eigen equation. They are to the\n",
    "default precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  1.61168440e+01+0.j   0.00000000e+00+0.j   0.00000000e+00+0.j]\n",
      " [  0.00000000e+00+0.j  -1.11684397e+00+0.j   0.00000000e+00+0.j]\n",
      " [  0.00000000e+00+0.j   0.00000000e+00+0.j  -9.75918483e-16+0.j]]\n",
      "[[  1.61168440e+01   3.55271368e-15   1.77635684e-15]\n",
      " [ -2.49800181e-15  -1.11684397e+00  -2.77555756e-17]\n",
      " [  2.79947848e-15   4.62304004e-16   9.86076132e-32]]\n"
     ]
    }
   ],
   "source": [
    "print(v)\n",
    "print(la.inv(r)@A@r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In reality, one should never ever use the inverse function, unless the\n",
    "actual answer you want is the inverse itself (which I've never seen for\n",
    "a real problem). What you typically want is the inverse of a matrix\n",
    "times another matrix or vector, which is the solution to a linear\n",
    "algebra problem. We can use the `solve` function to obtain this. This is\n",
    "arguably no better, but illustrative. For a larger problem, the benefit\n",
    "is easier to demonstrate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.61168440e+01,   3.05037570e-15,  -4.48960353e-17],\n",
       "       [ -2.74605686e-15,  -1.11684397e+00,   3.95339943e-16],\n",
       "       [  3.86524889e-15,  -7.16944125e-16,  -3.70074342e-17]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "la.solve(r,A)@r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is \"rebuilding\" the original matrix from the eigensolution. Looks\n",
    "pretty good."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.+0.j,  2.+0.j,  3.+0.j],\n",
       "       [ 4.+0.j,  5.+0.j,  6.+0.j],\n",
       "       [ 7.+0.j,  8.+0.j,  9.+0.j]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r@v@la.inv(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Avoiding the inverse is a bit uncomfortable in this case, but a bit of\n",
    "doodling yields that\n",
    "\n",
    "$$A B^{-1}  = \\left(\\left(A B^{-1}\\right)^T\\right)^T = \\left(\\left(B^{-1}\\right)^T A^T  \\right)^T$$\n",
    "\n",
    "Noting that\n",
    "\n",
    "$$\\left(B^{-1}\\right)^T  = \\left(B^{T}\\right)^{-1}$$\n",
    "\n",
    "we can use"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.+0.j,  2.+0.j,  3.+0.j],\n",
       "       [ 4.+0.j,  5.+0.j,  6.+0.j],\n",
       "       [ 7.+0.j,  8.+0.j,  9.+0.j]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r@la.solve(r.T,v.T).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the left eigenvectors, they are actually simply the right\n",
    "eigenvectors of the transpose of the matrix, so\n",
    "\n",
    "$$A^T l = l  v$$\n",
    "\n",
    "Below I lazily obtain the eigenvalues using the left eigenvectors, with\n",
    "an inverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.61168440e+01,   1.77635684e-15,   8.88178420e-16],\n",
       "       [  1.94289029e-15,  -1.11684397e+00,  -4.99600361e-16],\n",
       "       [  1.66684734e-15,   1.64791705e-16,   9.06493304e-17]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "la.inv(l)@A.T@l"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}