Trefethen & Bau lecture notes in MATLAB

TB-Lecture-01.ipynb

TB-Lecture-02.ipynb

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

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 10: Householder triangularization\n",
    "\n",
    "The stable way to get a QR factorization is to operate orthogonally (unitarily) in order to produce a triangular $R$. It's a lot like Gaussian elimination, proceeding one column at a time. Instead of elementary triangular row operations, though, we can choose from reflections or rotations. This lecture describes the reflection approach.\n",
    "\n",
    "The key step is to find, given $x$, a unitary $F$ such that $Fx=\\alpha e_1$ for a scalar $\\alpha$. Since $F$ preserves the 2-norm, $\\alpha= \\pm \\|x\\|_2$. One can simply exhibit the solution. Define $v=\\alpha e_1 - x$ and set $$ F = I - 2 \\frac{vv^*}{v^*v}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   4.9032e-16\n"
     ]
    }
   ],
   "source": [
    "x = randn(5,1);  alpha = norm(x);\n",
    "v = alpha*eye(5,1) - x;\n",
    "F = eye(5) - 2*(v*v')/(v'*v);\n",
    "norm(F'*F-eye(5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "    3.0983\n",
      "    0.0000\n",
      "   -0.0000\n",
      "    0.0000\n",
      "    0.0000\n"
     ]
    }
   ],
   "source": [
    "F*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In context, $x$ is drawn from rows $j$ to $m$ of column $j$, so that $F$ (applied to the lower rows) puts zeros below the diagonal as desired. For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "    4.9308   -0.4355   -0.8077\n",
      "   -0.0000    0.4842    1.2494\n",
      "    0.0000    0.1275    1.5434\n",
      "    0.0000    2.7054    1.9893\n",
      "         0    1.3357   -0.1876\n",
      "   -0.0000   -0.8751    0.2201\n"
     ]
    }
   ],
   "source": [
    "A = randn(6,3);  \n",
    "x = A(:,1);  \n",
    "v = [x(1)+sign(x(1))*norm(x);x(2:end)]; \n",
    "v = v/norm(v);\n",
    "F = eye(6) - 2*(v*v');\n",
    "F*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once $j$ sweeps from 1 to $n$, the matrix will be transformed into the $R$ we seek. If we accumulate the actions of these reflectors, we end up with the $Q$ as well (the full one or the thin one, as we choose).  \n",
    "\n",
    "Writing codes for the algorithm is part of the exercises, so I won't reproduce them here. There is an some important shortcut to know. In applications we rarely want the actual $Q$ or even the thin $\\hat{Q}$; instead we want the capability of applying $Q$ or $Q^*$ to a given vector. It's more efficient to store the $v$ vectors of the reflectors and apply them on demand, using \n",
    "\n",
    "$$Fz = z - 2\\frac{v^*z}{v^*v}v.$$\n",
    "\n",
    "One form of the `qr` command will actually return those $v$-vectors for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "QR =\n",
      "\n",
      "    4.9308   -0.4355   -0.8077\n",
      "    0.0695   -3.1811   -1.8046\n",
      "   -0.0549    0.0348   -2.1743\n",
      "   -0.5736    0.7381   -0.0733\n",
      "   -0.4439    0.3644   -0.3601\n",
      "    0.2164   -0.2387    0.2628\n"
     ]
    }
   ],
   "source": [
    "QR = qr(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The upper triangle of the result is $R$. Each reflector is built from the $v$ constructed having $v_1=1$ and the remaining elements from the lower triangle in the appropriate column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "HH = @(v) eye(length(v)) - 2*(v*v')/(v'*v);\n",
    "v1 = [1;QR(2:6,1)];  F1 = HH(v1);\n",
    "v2 = [1;QR(3:6,2)];  F2 = HH(v2);\n",
    "v3 = [1;QR(4:6,3)];  F3 = HH(v3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   1.5456e-15\n"
     ]
    }
   ],
   "source": [
    "Q = F1*blkdiag(1,F2)*blkdiag(eye(2),F3);\n",
    "R = triu(QR);\n",
    "norm(Q*R-A)"
   ]
  }
 ],
 "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
}

TB-Lecture-11.ipynb

TB-Lecture-12.ipynb