TB Lecture 7: QR factorization

TB-Lecture-07-julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 7: The QR factorization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gram-Schmidt orthogonalization\n",
    "\n",
    "The classical way to prove the existence of a QR factorization is to construct it via the Gram-Schmidt process. Let's look at it for a small example. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6x4 Array{Float64,2}:\n",
       " 0.0   2.0  9.0  1.0\n",
       " 7.0   8.0  1.0  4.0\n",
       " 8.0   9.0  9.0  1.0\n",
       " 1.0   5.0  5.0  1.0\n",
       " 7.0  10.0  4.0  2.0\n",
       " 0.0   5.0  4.0  4.0"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = round(10*rand(6,4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we want an orthonormal basis for the span of the first column. That's pretty easy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(Q,R) = (zeros(6,4),zeros(4,4));\n",
    "R[1,1] = norm(A[:,1]);\n",
    "Q[:,1] = A[:,1]/R[1,1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we add another orthonormal vector to become a basis for the span of the first _two_ columns. We can do that via orthogonal projection..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "v' * Q[:,1] = [-4.6629367034256575e-15]\n"
     ]
    }
   ],
   "source": [
    "P = Q[:,1]*Q[:,1]';  \n",
    "v = A[:,2] - P*A[:,2];\n",
    "@show v'*Q[:,1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we'll organize the process in a particular way. \n",
    "\n",
    "$$\n",
    "r_{22}q_2 = v = a_2 - (QQ^*)a_2 = a_2 - Q(Q^*a_2) = a_2 - q_1(r_{12})\n",
    "$$\n",
    "\n",
    "Doing so allows us to rearrange to\n",
    "\n",
    "$$\n",
    "a_2 = q_1r_{12} + q_2 r_{22}, \n",
    "$$\n",
    "\n",
    "which is the second column of $A=QR$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R = [12.767145334803704 15.900187134755534 0.0 0.0\n",
      " 0.0 6.7958847164850145 0.0 0.0\n",
      " 0.0 0.0 0.0 0.0\n",
      " 0.0 0.0 0.0 0.0]\n",
      "Q = [0.0 0.29429575153747534 0.0 0.0\n",
      " 0.5482823149915702 -0.10562148138001427 0.0 0.0\n",
      " 0.6266083599903659 -0.14173138954412173 0.0 0.0\n",
      " 0.07832604499879574 0.5524815949108433 0.0 0.0\n",
      " 0.5482823149915702 0.18867427015746108 0.0 0.0\n",
      " 0.0 0.7357393788436885 0.0 0.0]\n"
     ]
    }
   ],
   "source": [
    "R[1,2] = dot(Q[:,1],A[:,2]);\n",
    "v = A[:,2] - Q[:,1]*R[1,2];\n",
    "\n",
    "R[2,2] = norm(v);\n",
    "@show R;\n",
    "Q[:,2] = v/R[2,2];\n",
    "@show Q;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the relevant projector is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "v' * Q[:,1:2] = [-2.852267617604043e-15 3.2753572043548736e-15]\n"
     ]
    }
   ],
   "source": [
    "P = Q[:,1:2]*Q[:,1:2]';\n",
    "v = A[:,2] - P*A[:,2];\n",
    "@show v'*Q[:,1:2];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we break it into steps as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R = [12.767145334803704 15.900187134755534 8.772517039865122 0.0\n",
      " 0.0 6.7958847164850145 7.727520347118983 0.0\n",
      " 0.0 0.0 9.128437657679335 0.0\n",
      " 0.0 0.0 0.0 0.0]\n"
     ]
    }
   ],
   "source": [
    "(R[1,3],R[2,3]) = dot(Q[:,1],A[:,3]),dot(Q[:,2],A[:,3]);\n",
    "v = A[:,3] - Q[:,1:2]*R[1:2,3];\n",
    "R[3,3] = norm(v);\n",
    "@show R;\n",
    "Q[:,3] = v/R[3,3];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Q$ is unitary and $R$ is upper triangular by construction, and we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm(Q * R[:,1:3] - A[:,1:3]) = 1.275549143317629e-15\n"
     ]
    }
   ],
   "source": [
    "@show norm( Q*R[:,1:3] - A[:,1:3] );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rather than finish things manually, let's start over and automate everything in a loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "(m,n) = size(A);\n",
    "(Q,R) = (zeros(m,n),zeros(n,n));\n",
    "R[1,1] = norm(A[:,1]);\n",
    "Q[:,1] = A[:,1]/R[1,1];\n",
    "for j = 2:n\n",
    "    R[1:j-1,j] = Q[:,1:j-1]'*A[:,j];\n",
    "    v = A[:,j] - Q[:,1:j-1]*R[1:j-1,j];\n",
    "    R[j,j] = norm(v);\n",
    "    Q[:,j] = v/R[j,j];\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q = [0.0 0.29429575153747534 0.7367989840260739 0.3107599849234479\n",
      " 0.5482823149915702 -0.10562148138001454 -0.32794481561090993 0.6096183641478815\n",
      " 0.6266083599903659 -0.14173138954412198 0.5037334814158353 -0.11637297026426416\n",
      " 0.07832604499879574 0.5524815949108433 0.004773069733776502 -0.44844297969083946\n",
      " 0.5482823149915702 0.18867427015746083 -0.24843245882629786 -0.41255740103288735\n",
      " 0.0 0.7357393788436885 -0.184636307262656 0.3833357516263034]\n"
     ]
    }
   ],
   "source": [
    " \n",
    "@show Q;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5207720726002076e-15"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(Q'*Q-eye(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5036864803660879e-15"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(A-Q*R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This only produces the \"thin\" version of the QR factorization. Conceptually at least, the full version just tacks on an orthonormal basis for whatever remains of $\\mathbb{C}^m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear system solution\n",
    "\n",
    "We can apply QR to solve a linear system $Ax=b$ for $m\\times m$ matrix $A$. We have $QRx=b$, so $Rx=Q^*b$ is an equivalent triangular system, solvable by back substitution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = [0.40145622419481397,0.6999528823342192,1.4627482274682255,-0.4412557347859269,1.155606318206312,1.1186980122017582,-2.074662732712543,-0.17133966158239913,-0.9745719879262823]\n"
     ]
    }
   ],
   "source": [
    "A = round(10*rand(9,9));   b = collect(1:9);\n",
    "m = 9;\n",
    "(Q,R) = qr(A);\n",
    "z = Q'*b;\n",
    "x = zeros(m);\n",
    "x[m] = z[m]/R[m,m];\n",
    "for i = m-1:-1:1\n",
    "    y = z[i] - R[i,i+1:m]*x[i+1:m];\n",
    "    x[i] = y[1] / R[i,i];\n",
    "end\n",
    "@show x;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0880185641326534e-14"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(A*x-b)"
   ]
  }
 ],
 "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": 1
}

TB-Lecture-07-matlab.ipynb