TB Lectures 25-26

TB-Lectures-25-26-julia.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 you call `eig`, you get the complete EVD. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "λ = [34.0,8.94427,-8.94427,2.05764e-15]\n",
      "X = [-0.5 -0.823607 0.376393 -0.223607; -0.5 0.423607 0.0236068 -0.67082; -0.5 0.0236068 0.423607 0.67082; -0.5 0.376393 -0.823607 0.223607]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " -0.5  -0.823607    0.376393   -0.223607\n",
       " -0.5   0.423607    0.0236068  -0.67082 \n",
       " -0.5   0.0236068   0.423607    0.67082 \n",
       " -0.5   0.376393   -0.823607    0.223607"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [\n",
    "    16     2     3    13\n",
    "     5    11    10     8\n",
    "     9     7     6    12\n",
    "     4    14    15     1\n",
    "\n",
    "];\n",
    "(λ,X) = eig(A);\n",
    "@show λ;\n",
    "@show X;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want the full decomposition, call with two output arguments."
   ]
  },
  {
   "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": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank_X = rank(X) = 1\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = eye(4) + diagm([1;1;1],1)   # Jordan block\n",
    "(λ,X) = eig(A);\n",
    "@show rank_X = rank(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.8459253727671276e-16"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = diagm(λ);\n",
    "residual = norm(A*X-X*D)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can request the eigenvalues alone using `eigvals`. This can be significantly faster than the full EVD."
   ]
  },
  {
   "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": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 65.0    5.3192e-15   -2.099e-15   2.71171e-15   3.63509e-15\n",
       "  0.0  -21.2768       -2.58883     2.18707      -3.48932    \n",
       "  0.0    0.0         -13.1263     -3.38447      -2.82394    \n",
       "  0.0    0.0           0.0        21.2768        2.62871    \n",
       "  0.0    0.0           0.0         0.0          13.1263     "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [\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",
    "(T,Q) = schur(A);\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": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Complex{Float64},1}:\n",
       "        65.0+0.0im    \n",
       " 1.55431e-15+21.2768im\n",
       " 1.55431e-15-21.2768im\n",
       "         0.0+13.1263im\n",
       "         0.0-13.1263im"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = A[:,end:-1:1]';\n",
    "λ = eigvals(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 65.0   -4.9109e-15    2.68933e-15    1.39058e-15   4.36516e-15\n",
       "  0.0    1.55431e-15  20.9795         6.8026        2.53686e-14\n",
       "  0.0  -21.5782        1.55431e-15    5.24025e-14  -1.70351    \n",
       "  0.0    0.0           0.0            0.0          13.4714     \n",
       "  0.0    0.0           0.0          -12.79          0.0        "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(T,Q) = schur(A);\n",
    "T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Complex{Float64},1}:\n",
       " 1.55431e-15+21.2768im\n",
       " 1.55431e-15-21.2768im"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(T[2:3,2:3])"
   ]
  },
  {
   "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": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       "  15.0     -12.2967  -25.744    3.85678      1.07758\n",
       " -30.4959   34.543    16.0572  -7.03914     -1.13121\n",
       "   0.0      33.2719   17.3227   3.97494     -2.24645\n",
       "   0.0       0.0     -10.6355  -0.0166871  -17.4078 \n",
       "   0.0       0.0       0.0     13.8061      -1.84898"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = hessfact(A)[:H]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Complex{Float64},1}:\n",
       "        65.0+0.0im    \n",
       " 1.55431e-15+21.2768im\n",
       " 1.55431e-15-21.2768im\n",
       "         0.0+13.1263im\n",
       "         0.0-13.1263im"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(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": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       "  30.0     -50.0999    0.0        0.0           0.0        \n",
       " -50.0999  104.382     8.18347   -7.07535e-15   7.89787e-16\n",
       "   0.0       8.18347   0.578209  -1.64647       4.82947e-15\n",
       "   0.0       0.0      -1.64647   -5.10605      -3.04241    \n",
       "   0.0       0.0       0.0       -3.04241       0.145366   "
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = A+A';\n",
    "hessfact(A)[:H]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

TB-Lectures-25-26-matlab.ipynb