Trefethen & Bau lecture notes in Julia

TB-Lecture-01-julia.ipynb

TB-Lecture-02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 2: Orthogonal vectors and matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "With real vectors and matrices, the transpose operation is simple and familiar. It also happens to correspond to what we call the **adjoint** mathematically. In the complex case, one also has to conjugate the entries to keep the mathematical structure intact. We call this operator the **hermitian** of a matrix and use a star superscript for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2x4 Array{Complex{Float64},2}:\n",
       " 0.725337+0.164674im  0.288889+0.794957im  …  0.825082+0.589749im\n",
       " 0.188618+0.51475im   0.878958+0.840362im      0.25438+0.839209im"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(2,4) + 1im*rand(2,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4x2 Array{Complex{Float64},2}:\n",
       " 0.725337-0.164674im   0.188618-0.51475im \n",
       " 0.288889-0.794957im   0.878958-0.840362im\n",
       " 0.384954-0.0683305im  0.777636-0.230011im\n",
       " 0.825082-0.589749im    0.25438-0.839209im"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Aadjoint = A'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get plain transpose, use a `.^` operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4x2 Array{Complex{Float64},2}:\n",
       " 0.725337+0.164674im   0.188618+0.51475im \n",
       " 0.288889+0.794957im   0.878958+0.840362im\n",
       " 0.384954+0.0683305im  0.777636+0.230011im\n",
       " 0.825082+0.589749im    0.25438+0.839209im"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Atrans = A.'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inner products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If **u** and **v** are column vectors of the same length, then their **inner product** is $\\mathbf{u}^*\\mathbf{v}$. The result is a scalar. (In Julia, though, the scalar inner product is `dot(u,v)`, whereas the matrix multiplication result is a 1-by-1 matrix that is not equivalent for all purposes.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dot(u,v) gives -2 - 3im\n",
      "u'*v gives Complex{Int64}[-2 - 3im]\n"
     ]
    }
   ],
   "source": [
    "u = [ 4; -1; 2+2im ]\n",
    "v = [ -1; 1im; 1 ]\n",
    "println(\"dot(u,v) gives \", dot(u,v))\n",
    "println(\"u'*v gives \",u'*v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inner product has geometric significance. It is used to define length through the 2-norm,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Complex{Int64},1}:\n",
       " 25+0im"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length_u_squared = u'*u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum( abs(u).^2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm_u = norm(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It also defines the angle between two vectors as a generalization of the familiar dot product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Complex{Float64},1}:\n",
       " -0.23094-0.34641im"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cos_theta = (u'*v) / ( norm(u)*norm(v) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The angle may be complex when the vectors are complex! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1-element Array{Complex{Float64},1}:\n",
       " 1.79019+0.34786im"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta = acos(cos_theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The operations of inverse and hermitian commute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4x4 Array{Complex{Float64},2}:\n",
       "   1.82916+2.59705im   -0.816702-1.96154im   …  -0.530187+0.312639im\n",
       " -0.303269-0.356645im    1.02432+0.929259im     -0.350343-1.09856im \n",
       "  -1.69821-0.467787im   0.649194+0.784488im       1.01278-1.67234im \n",
       "  0.684181-1.39733im   -0.779362+0.402974im     -0.378864+2.93732im "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(4,4)+1im*rand(4,4);  (inv(A))'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4x4 Array{Complex{Float64},2}:\n",
       "   1.82916+2.59705im   -0.816702-1.96154im   …  -0.530187+0.312639im\n",
       " -0.303269-0.356645im    1.02432+0.929259im     -0.350343-1.09856im \n",
       "  -1.69821-0.467787im   0.649194+0.784488im       1.01278-1.67234im \n",
       "  0.684181-1.39733im   -0.779362+0.402974im     -0.378864+2.93732im "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv(A')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So we just write $\\mathbf{A}^{-*}$ for either case. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Orthogonality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Orthogonality, which is the multidimensional extension of perpendicularity, means that $\\cos \\theta=0$, i.e., that the inner product between vectors is zero. A collection of vectors is orthogonal if they are all pairwise orthogonal. \n",
    "\n",
    "Don't worry about how we are creating the vectors here for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x3 Array{Float64,2}:\n",
       " -0.520352    0.113586    0.187403\n",
       " -0.169743   -0.61493     0.140415\n",
       " -0.0291219  -0.763508   -0.268636\n",
       " -0.746073    0.135121   -0.53895 \n",
       " -0.378085   -0.0880756   0.763238"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q = qr(rand(5,3))[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $\\mathbf{Q}^*\\mathbf{Q}$ is a matrix of all inner products between columns of $\\mathbf{Q}$, those columns are orthogonal if and only if that matrix is diagonal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3x3 Array{Float64,2}:\n",
       "  1.0          5.09384e-17  -1.42308e-17\n",
       "  5.09384e-17  1.0           6.8765e-17 \n",
       " -1.42308e-17  6.8765e-17    1.0        "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QhQ = Q'*Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact we have a stronger condition here: the columns are **orthonormal**, meaning that they are orthogonal and each has 2-norm equal to 1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given any other vector of length 5, we can compute its inner product with each of the columns of $\\mathbf{Q}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3x1 Array{Float64,2}:\n",
       " -0.862349\n",
       " -0.715461\n",
       "  0.514809"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = rand(5,1);  c = Q'*u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then use these coefficients to find a vector in the column space of $\\mathbf{Q}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x1 Array{Float64,2}:\n",
       " 0.463935\n",
       " 0.658623\n",
       " 0.433078\n",
       " 0.269245\n",
       " 0.781978"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = Q*c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As explained in the text, $\\mathbf{r} = \\mathbf{u}-\\mathbf{v}$ is orthogonal to all of the columns of $\\mathbf{Q}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3x1 Array{Float64,2}:\n",
       " -3.85016e-18\n",
       " -8.15e-18   \n",
       " -1.01795e-16"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = u-v;  Q'*r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consequently, we have decomposed $\\mathbf{u}=\\mathbf{v}+\\mathbf{r}$ into the sum of two orthogonal parts, one lying in the range of $\\mathbf{Q}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1x1 Array{Float64,2}:\n",
       " -1.55597e-17"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v'*r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unitary matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We just saw that a matrix whose columns are orthonormal is pretty special. It becomes even more special if the matrix is also square, in which case we call it **unitary**. (In the real case, such matrices are confusingly called _orthogonal_. Ugh.) Say $\\mathbf{Q}$ is unitary and $m\\times m$. Then $\\mathbf{Q}^*\\mathbf{Q}$ is an $m\\times m$ identity matrix---that is, $\\mathbf{Q}^*=\\mathbf{Q}^{-1}$! It can't get much easier in terms of finding the inverse of a matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x5 Array{Float64,2}:\n",
       " 2.28878e-16  1.46376e-16  1.11886e-16  2.48253e-16  2.83052e-16\n",
       " 1.66533e-16  1.24127e-16  2.77556e-17  1.31656e-16  1.57009e-16\n",
       " 5.55112e-17  3.51083e-16  8.32667e-17  1.14439e-16  1.14439e-16\n",
       " 5.55112e-17  2.08167e-16  2.00148e-16  1.24127e-16  1.66533e-16\n",
       " 2.16778e-16  1.8619e-16   1.11022e-16  3.92523e-16  1.57009e-16"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q = qr(rand(5,5)+1im*rand(5,5))[1]\n",
    "abs( inv(Q) - Q' )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rank of $\\mathbf{Q}$ is $m$, so continuing the discussion above, the original vector $\\mathbf{u}$ lies in its column space. Hence the remainder $\\mathbf{r}=\\boldsymbol{0}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x1 Array{Complex{Float64},2}:\n",
       " 3.33067e-16-1.66533e-16im \n",
       "                  0.0-0.0im\n",
       "   2.77556e-17-1.249e-16im \n",
       " 5.55112e-16-1.11022e-16im \n",
       " 1.11022e-16+2.22045e-16im "
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c = Q'*u; \n",
    "v = Q*c;\n",
    "r = u - v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is another way to arrive at a fact we already knew: Multiplication by $\\mathbf{Q}^*=\\mathbf{Q}^{-1}$ changes the basis to the columns of $\\mathbf{Q}$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.3",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}

TB-Lecture-03.ipynb

TB-Lecture-04.ipynb

TB-Lecture-05.ipynb

TB-Lecture-06.ipynb

TB-Lecture-07.ipynb

TB-Lecture-08.ipynb

TB-Lecture-09.ipynb

TB-Lecture-10.ipynb

TB-Lecture-11.ipynb

TB-Lecture-12.ipynb