TB Lecture 7: QR factorization

TB-Lecture-07-julia.ipynb

TB-Lecture-07-matlab.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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A =\n",
      "\n",
      "    35     1     6    26\n",
      "     3    32     7    21\n",
      "    31     9     2    22\n",
      "     8    28    33    17\n",
      "    30     5    34    12\n",
      "     4    36    29    13\n"
     ]
    }
   ],
   "source": [
    "A = magic(6);  A = A(:,1: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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R =\n",
      "\n",
      "   56.3471\n",
      "\n",
      "\n",
      "Q =\n",
      "\n",
      "    0.6211\n",
      "    0.0532\n",
      "    0.5502\n",
      "    0.1420\n",
      "    0.5324\n",
      "    0.0710\n"
     ]
    }
   ],
   "source": [
    "R = norm(A(:,1)),  Q = A(:,1)/R"
   ]
  },
  {
   "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": [
      "innerprod =\n",
      "\n",
      "  -2.2204e-15\n"
     ]
    }
   ],
   "source": [
    "P = Q*Q';  \n",
    "v = A(:,2) - P*A(:,2);\n",
    "innerprod = v'*Q"
   ]
  },
  {
   "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 =\n",
      "\n",
      "   56.3471   16.4693\n",
      "         0   54.2196\n",
      "\n",
      "\n",
      "Q =\n",
      "\n",
      "    0.6211   -0.1702\n",
      "    0.0532    0.5740\n",
      "    0.5502   -0.0011\n",
      "    0.1420    0.4733\n",
      "    0.5324   -0.0695\n",
      "    0.0710    0.6424\n"
     ]
    }
   ],
   "source": [
    "R(1,2) = Q(:,1)'*A(:,2);\n",
    "v = A(:,2) - Q(:,1)*R(1,2);\n",
    "R(2,2) = norm(v)\n",
    "Q(:,2) = v/R(2,2)"
   ]
  },
  {
   "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": [
      "innerprod =\n",
      "\n",
      "   1.0e-14 *\n",
      "\n",
      "    0.2228   -0.2821\n"
     ]
    }
   ],
   "source": [
    "P = Q(:,1:2)*Q(:,1:2)';\n",
    "v = A(:,2) - P*A(:,2);\n",
    "innerprod = v'*Q"
   ]
  },
  {
   "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 =\n",
      "\n",
      "   56.3471   16.4693   30.0459\n",
      "         0   54.2196   34.8797\n",
      "         0         0   32.4907\n",
      "\n",
      "\n",
      "Q =\n",
      "\n",
      "    0.6211   -0.1702   -0.2070\n",
      "    0.0532    0.5740   -0.4500\n",
      "    0.5502   -0.0011   -0.4460\n",
      "    0.1420    0.4733    0.3763\n",
      "    0.5324   -0.0695    0.6287\n",
      "    0.0710    0.6424    0.1373\n"
     ]
    }
   ],
   "source": [
    "R(1:2,3) = Q(:,1:2)'*A(:,3);\n",
    "v = A(:,3) - Q(:,1:2)*R(1:2,3);\n",
    "R(3,3) = norm(v)\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": [
      "ans =\n",
      "\n",
      "   3.5527e-15\n"
     ]
    }
   ],
   "source": [
    "norm( Q*R - 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": true
   },
   "outputs": [],
   "source": [
    "R = norm(A(:,1));\n",
    "Q = A(:,1)/R(1,1);\n",
    "for j = 2:size(A,2)\n",
    "    R(1:j-1,j) = Q'*A(:,j);\n",
    "    v = A(:,j) - Q*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 =\n",
      "\n",
      "    0.6211   -0.1702   -0.2070    0.4998\n",
      "    0.0532    0.5740   -0.4500    0.2106\n",
      "    0.5502   -0.0011   -0.4460   -0.4537\n",
      "    0.1420    0.4733    0.3763    0.5034\n",
      "    0.5324   -0.0695    0.6287   -0.2096\n",
      "    0.0710    0.6424    0.1373   -0.4501\n",
      "\n",
      "\n",
      "R =\n",
      "\n",
      "   56.3471   16.4693   30.0459   39.0969\n",
      "         0   54.2196   34.8797   23.1669\n",
      "         0         0   32.4907   -8.9182\n",
      "         0         0         0    7.6283\n"
     ]
    }
   ],
   "source": [
    "Q,R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   1.2510e-15\n"
     ]
    }
   ],
   "source": [
    "norm(Q'*Q-eye(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   7.9441e-15\n"
     ]
    }
   ],
   "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 =\n",
      "\n",
      "    0.0024\n",
      "    0.0024\n",
      "    0.0024\n",
      "    0.0024\n",
      "    0.1024\n",
      "    0.0024\n",
      "    0.0024\n",
      "    0.0024\n",
      "    0.0024\n"
     ]
    }
   ],
   "source": [
    "A = magic(9);   b = (1:9)';\n",
    "[Q,R] = qr(A);\n",
    "z = Q'*b;\n",
    "x(9,1) = z(9)/R(9,9);\n",
    "for i = 8:-1:1\n",
    "    x(i) = (z(i) - R(i,i+1:9)*x(i+1:9)) / R(i,i);\n",
    "end\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   3.1792e-15\n"
     ]
    }
   ],
   "source": [
    "norm(A*x-b)"
   ]
  }
 ],
 "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
}