TB Lectures 25-26

TB-Lectures-25-26-julia.ipynb

TB-Lectures-25-26-matlab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lectures 25-26: Eigenvalues and algorithmic considerations\n",
    "\n",
    "## Eigenvalue decomposition\n",
    "\n",
    "When it exists, the eigenvalue decomposition (EVD) expresses a square matrix $A$ as a change of basis away from a diagonal matrix: $A=XDX^{-1}$. There are two major syntaxes. If only the eigenvalues are desired, call `eig` with one or zero output arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      "\n",
      "    16     2     3    13\n",
      "     5    11    10     8\n",
      "     9     7     6    12\n",
      "     4    14    15     1\n",
      "\n",
      "\n",
      "lam =\n",
      "\n",
      "   34.0000\n",
      "    8.9443\n",
      "   -8.9443\n",
      "   -0.0000\n"
     ]
    }
   ],
   "source": [
    "A = magic(4)\n",
    "lam = eig(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want the full decomposition, call with two output arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X =\n",
      "\n",
      "   -0.5000   -0.8236    0.3764   -0.2236\n",
      "   -0.5000    0.4236    0.0236   -0.6708\n",
      "   -0.5000    0.0236    0.4236    0.6708\n",
      "   -0.5000    0.3764   -0.8236    0.2236\n",
      "\n",
      "\n",
      "D =\n",
      "\n",
      "   34.0000         0         0         0\n",
      "         0    8.9443         0         0\n",
      "         0         0   -8.9443         0\n",
      "         0         0         0   -0.0000\n"
     ]
    }
   ],
   "source": [
    "[X,D] = eig(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Technically, the output is guaranteed only to satisfy $AX=XD$, which holds for defective as well as diagonalizable matrices. (But defectiveness is a \"probability zero\" event and is not robust to floating point perturbations.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      "\n",
      "     1     1     0     0\n",
      "     0     1     1     0\n",
      "     0     0     1     1\n",
      "     0     0     0     1\n",
      "\n",
      "\n",
      "rank_X =\n",
      "\n",
      "     1\n",
      "\n",
      "\n",
      "residual =\n",
      "\n",
      "   3.8459e-16\n"
     ]
    }
   ],
   "source": [
    "A = eye(4) + diag([1 1 1],1)   % Jordan block\n",
    "[X,D] = eig(A);\n",
    "rank_X = rank(X)\n",
    "residual = norm(A*X-X*D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use `poly` to get the characteristic polynomial of a matrix. However, as pointed out in Lecture 15, using it to compute eigenvalues is an unstable algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "     0\n"
     ]
    }
   ],
   "source": [
    "A = eye(2) + 1e-10*[0 0;0 1];\n",
    "norm( eig(A) - diag(A) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   2.1003e-08\n"
     ]
    }
   ],
   "source": [
    "norm( roots(poly(A)) - diag(A) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(As an interesting note, the connection is actually in the other direction: eigenvalue computations are used to find polymomial roots!) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Schur factorization\n",
    "\n",
    "The EVD has some serious drawbacks as a target for numerical computation. First, it does not exist for every square matrix. And second, if it exists but the condition number of the eigenvalue matrix is large, it will be very sensitive to roundoff. A better alternative is the **Schur factorization** $A=QTQ^*$, in which $Q$ is unitary and $T$ is upper triangular."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      "\n",
      "    17    24     1     8    15\n",
      "    23     5     7    14    16\n",
      "     4     6    13    20    22\n",
      "    10    12    19    21     3\n",
      "    11    18    25     2     9\n",
      "\n",
      "\n",
      "ans =\n",
      "\n",
      "   1.2786e-15\n",
      "\n",
      "\n",
      "T =\n",
      "\n",
      "   65.0000    0.0000   -0.0000   -0.0000    0.0000\n",
      "         0  -21.2768   -2.5888    2.1871   -3.4893\n",
      "         0         0  -13.1263   -3.3845   -2.8239\n",
      "         0         0         0   21.2768    2.6287\n",
      "         0         0         0         0   13.1263\n"
     ]
    }
   ],
   "source": [
    "A = magic(5)\n",
    "[Q,T] = schur(A);\n",
    "norm( Q'*Q-eye(5) )\n",
    "T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Schur form always exists, it still reveals the eigenvalues of $A$, and it involves only stable unitary transformations. There is one quirk with it. If $A$ is real but has complex eigenvalues, there is an alternate form in which $T$ is \"quasitriangular\". Each real eigenvalue appears along the diagonal of $T$, but complex conjugate pairs are represented as $2\\times 2$ real blocks along the diagonal. This allows the entire factorizaton to be done in real arithmetic, which can improve speed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lam =\n",
      "\n",
      "  65.0000 + 0.0000i\n",
      "   0.0000 +21.2768i\n",
      "   0.0000 -21.2768i\n",
      "  -0.0000 +13.1263i\n",
      "  -0.0000 -13.1263i\n"
     ]
    }
   ],
   "source": [
    "A = rot90(A);\n",
    "lam = eig(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T =\n",
      "\n",
      "   65.0000    0.0000   -0.0000   -0.0000   -0.0000\n",
      "         0    0.0000   20.9795    6.8026   -0.0000\n",
      "         0  -21.5782    0.0000   -0.0000   -1.7035\n",
      "         0         0         0   -0.0000   13.4714\n",
      "         0         0         0  -12.7900   -0.0000\n"
     ]
    }
   ],
   "source": [
    "[Q,T] = schur(A);\n",
    "T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   0.0000 +21.2768i\n",
      "   0.0000 -21.2768i\n",
      "\n",
      "\n",
      "ans =\n",
      "\n",
      "  -0.0000 +13.1263i\n",
      "  -0.0000 -13.1263i\n"
     ]
    }
   ],
   "source": [
    "eig(T(2:3,2:3))\n",
    "eig(T(4:5,4:5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hessenbergification\n",
    "\n",
    "Eigenvalue algorithms are iterative. The iterative steps can be performed a lot faster if the matrix has structural zeros. Hence before iteration begins, an orthogonal similarity transformation is used to get an **upper Hessenberg** matrix with the same eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H =\n",
      "\n",
      "   15.0000  -12.2967  -25.7440    3.8568    1.0776\n",
      "  -30.4959   34.5430   16.0572   -7.0391   -1.1312\n",
      "         0   33.2719   17.3227    3.9749   -2.2465\n",
      "         0         0  -10.6355   -0.0167  -17.4078\n",
      "         0         0         0   13.8061   -1.8490\n"
     ]
    }
   ],
   "source": [
    "H = hess(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "  65.0000 + 0.0000i\n",
      "   0.0000 +21.2768i\n",
      "   0.0000 -21.2768i\n",
      "  -0.0000 +13.1263i\n",
      "  -0.0000 -13.1263i\n"
     ]
    }
   ],
   "source": [
    "eig(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $A=QHQ^*$ is \"nearly\" a Schur factorization. All that remains to be done is to make the subdiagonal zero. \n",
    "\n",
    "If $A$ is hermitian, then since the transformation preserves symmetry the resulting $H$ is actually hermitian and tridiagonal. This leads to even greater speedup in eigenvalue computations; for this and other reasons, one might even consider the hermitian and nonhermitian problems to be of distinct species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "    0.1631    3.2854         0         0         0\n",
      "    3.2854    4.8747    1.5240         0         0\n",
      "         0    1.5240    0.4987   -7.5820         0\n",
      "         0         0   -7.5820  102.4634  -54.0925\n",
      "         0         0         0  -54.0925   22.0000\n"
     ]
    }
   ],
   "source": [
    "A = A+A';\n",
    "hess(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "matlab",
   "version": "0.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}