TB-Lecture-08-julia.ipynb

TB-Lecture-08-matlab.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 8: Gram-Schmidt orthogonalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modified Gram-Schmidt\n",
    "\n",
    "The classical Gram-Schmidt process is numerically unstable, as will be seen in the next lecture. A mathematically equivalent formulation, which we call MGS, turns out to be a lot better. \n",
    "\n",
    "In classical GS the key step is to project a new column of $A$ onto the orthogonal complement of the span of all the previous $q_j$s. We express this projection as $P_j=I - Q_jQ_j^*$, where $Q_j=\\bigl[ q_1 \\; q_2\\; \\cdots \\; q_j\\bigr]$. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "    1.0000\n",
      "    1.0000\n",
      "    1.0000\n",
      "    0.0000\n",
      "    0.0000\n",
      "    0.0000\n"
     ]
    }
   ],
   "source": [
    "[Qj,~] = qr(rand(6,3),0);   % get a random Q with 3 ON columns\n",
    "Pj = eye(6) - Qj*Qj';\n",
    "svd(Pj)       % all 1 or 0, as required for an orthogonal projector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's intuitively clear that we get the same result by subtracting off projections onto each column, one at a time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   1.2505e-16\n"
     ]
    }
   ],
   "source": [
    "Mj = eye(6) - Qj(:,1)*Qj(:,1)' - Qj(:,2)*Qj(:,2)' - Qj(:,3)*Qj(:,3)';\n",
    "norm(Pj-Mj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also factor this into a product of orthogonal projectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ans =\n",
      "\n",
      "   2.3849e-16\n"
     ]
    }
   ],
   "source": [
    "Gj = eye(6);\n",
    "for k = 1:3\n",
    "    Gj = Gj*(eye(6)-Qj(:,k)*Qj(:,k)');\n",
    "end\n",
    "norm(Gj-Pj)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works because the \"cross-terms\" you get in the products have inner products between different $q_k$ and thus are zero. In MGS, we apply each of these factors, as soon at it is found, to \\emph{all} the columns of $A$ at once. By the time we are ready to choose a new column from $A$, it is already orthgonalized against the entire ON basis set found so far."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operation counts\n",
    "\n",
    "It's standard in numerical analysis to compare the performance of algorithms by counting the number of floating point operations, or flops, that are performed as a function of the input size. By either form of QR factorization, the number of flops is $\\sim 2mn^2$ in asymptotic notation. Let's see if we observe the implicit proportionality with respect to $n^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAABcSAAAXEgFnn9JSAAAA\nB3RJTUUH4AkTDBk4LzTDwAAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAxOS1TZXAtMjAxNiAwODoyNTo1Ni38EfEAACAA\nSURBVHic7d1/UFT3vf/xDxdJVsEYSM11wR9QLAjGzoQJkwlxvoI64CQ2HS0avr2molTbYL5pKLa3\nE+53gNzSXHVMNTeaTkuKMBnDELzftNe1Q0oE52oUzODUDjEqyDqgIG1XIRF23EG+f5z25OTsgqi7\n53zO2efjr93PObv7wZw5r3zO+ZzPO2J8fFwAAGC2fzK7AwAACEEgAQAkQSABAKRAIAEApEAgAQCk\nQCABAKRAIAEApEAgAQCkQCABAKRAIAEApEAgAQCkQCABAKRAIAEApEAgAQCkQCABAKRAIAEApEAg\nAQCkQCABAKRAIAEApEAgAQCkMM3sDtyF06dPV1dXd3d3x8XFLV269Ac/+MGDDz5odqcAAMFhmRHS\n//zP/7zwwgsjIyOFhYWZmZnV1dXFxcXj4+Nm9wsAEBwRVjmnf+tb34qKimpsbPynf/onIcR///d/\nb9++/be//e3TTz9tdtcAAEFgjRHSX//61wsXLnzrW99S0kgIsWrVqgceeOD48ePmdgwAECzWCKTe\n3l4hRGJiotoSFRWVkJBw5coV0/oEAAgqawTSyMiIECI6OlrbOGPGjJs3b5rUIwBAkFkjkJQbXbrb\nXePj41a5AQYAuCNDp30PDQ2tWbNm8+bNGzZs8N/a3Nx8+PDhrq6umJiYlJSUoqKiBQsWKJumT58u\nhBgdHdXuPzo66nQ6Deg2AMAAho6Qqqurr1y5EvA6W1VV1bZt244ePRoTEzM8PNzQ0PDcc8+pcxaU\n4HG73er+Y2NjfX19BBIA2EbIA2lsbMztdh87dqykpOTXv/51wH3a2trq6urmz5//hz/8ob6+/siR\nI3v37vX5fGVlZV6vVwgRHx/vdDr/+Mc/qh9pbW31+XxPPPFEqPsPIMy5Pd6kqo91LWZ1xt5Cfsnu\n0qVLq1evnnyfQ4cOCSG2b9+ekJCgtOTl5a1atcrlcp04cWLFihVCiO9973s7duyoqqpat27dhQsX\ndu/ePW/evJycnFD3H0A4q2zqqfiwRwgRUXpUbUyMcxQ+4SzPSzKvX/YU8kByOp1vvvmm8vrs2bPV\n1dX++7S1tUVGRi5btkzbuGLFCpfLderUKSWQNm7cODQ0VFNTU1dXJ4RYsmTJ66+/7nA4Jvnp1NTU\noP0ZAMLShdVv+ze6Pd6KD3t2vtc09+Qvje/S/Tt//rzZXQgs5IEUExOTl5f39x+bFuDnfD7f4ODg\n3LlzdemSnJwshOjr61PeRkZGlpSUvPzyy1euXImNjZ05c+ZUfl3af/cpSk1NNf5PCPqP3ucX3sPH\n7+ojU9n5/vcx5T9lcBn/J8hwKGoHRjojj6RM5dum/qNT3POOu93xUJxKZ0xh/rTv4eHh27dvz5o1\nS9eutAwNDWkbIyMj58+fP8U0AoB7s6n+XETpUW0aVeQmVeQmPXLBVZHLlbpQMX+1b5/PJwINnpQW\nZWvYMuX/qYP+o/f5hffw8bv6yFR2DtY+lmb8H2jiodjafV19nRjnaHkxIzHOIYQoz9sthDjwSb8y\nryGi9GhinKOnLCsoPzrFPe+4m3UPRfNHSJGRkSJQ8CgtylYAMJJ2Hl1PWZaSRqqWFzO0e7Z23zCu\nZ7ZmfiAFfOhVbVG23rPU1FSZL5gifFj3f1rD0IHT/errwkynEMLt8WojKjHOkZ0cq7491nVdWIH8\n50PzAykmJmbmzJm9vb1jY2PaduUx2Pt89PX8+fOcCADcldrTA+rrjZlOIcSm+nM5b3doM0k75/vA\nJ/3CCuQ/H5ofSEKItLS0W7dudXZ2ahs7OjqEEOnp6SZ1CkDY+VvD7gv58dnJD6stlU09OfvPtHZf\nd3u8OW93VDb1KO21mlFUOdMcgkSKQMrNzRVC7Nq1S23p7+8/ePBgZGTk8uXLzesXgDAy3NLwt4bd\nQoj//Zsvy362dl9XJzhobxdpL+tpL9/hfpg/y04IUVBQUF9f397evnbt2mefffbatWsul2tkZGTL\nli3x8fFm9w6A/Q23NAzse0V9e/7yhtQF7+r2SYxz1BSkuT3eTfXndO1GdDEMSBFIUVFRtbW1r732\nWnNzs3LhLjo6urS0tKioyOyuAbC/0c6T2jQSQszZtqdlfkbO/g5tY2Ls9NrT/cpKQipl1gOCIkKq\nkkI+n8/tdkdHRzudzoiIiPv8NnU+ieT38QCYaLTzZG/5d7Qtc7bt6Zifu6n+08kXUR3fbbEbCvKf\nEqUYIamioqK+8Y1vBPELpf13ByAD/zR6a9ba/zz8NSE6JvqIdSnnQ5lnfksxqQEAjOefRv8v5n/9\n58Nrp/JZFhAKBblGSABgjIBp9LNHtup2064bpEiq+rimIF07NRzBwggJQNjxDfbq0qjdkeb936/7\nT+CuKUjXTaLrKcsijUKEQAIQXnyDvT3FT2pbBp2PL33TtTHTqV1TVZGzv0NXLhahQyABCCP+adTu\nSFv6n67EOIfu6SIVy6caxuaBJP9iggCMdG1fifZtuyPto9VvCSGU9YHU9uzkWCsunzo5+c+HNg8k\n+RcTBGCkuZWN6ut2R9oL/1y2MdNZ2dSjL4BU/LgVl0+dnPznQ2bZAQgvKY1XL+THX5k2+4V/LhNC\n+D8AW1OQLlg+1Qw2HyEBgL8Hfn1pecIvlde6NGopzlAm0bF8qvEIJABhR1dhT6Wkkdvjzdl/Rre/\nUV0LawQSgHDUUvy4LpMS4xzHuq5XNvUkVX2svaXE8qmG4R4SAJtTauv1lGVpWzbVn9M9deT2eHUr\neVtu+VSrs3kgKXMcJZ9YAiB0Kpt6lJiJKD1qdl9MJvmcb2H7QCKKgDCnG/TcxQdtN7NO/tW+bR5I\nADC5luKMxFgHy6fKgEACYEOb6s9p520LISpyk9zXvbpGdZK3lvZuE4zELDsANqRbeaGnLKs8L8m/\nAixrp0qFQAJgE6OdJ9XX2uzpKctKjHPoVqvT7snaqZIgkADYgVJw70J+vPjqIgvKU0T+q9XZb+1U\nGyCQAFietvzrhfz42tMD6qaNmc7W7hvauXZKEVj7rZ1qAwQSAGvzL0b+w4t71deb6j/N2d+h3aoU\ngWXtVAnZPJDkr/8B4H74p9FD2euf/L+/Ut+ydqpK/vOhzad982AsYGO+wV5dGs1YnDXnpT1CiOzk\nWP8pDBW5ScraqbrisGGydqr8D8bafIQEwK78i5HPWJyl1t/zXztVwdqpMrP5CAmALU2eRsJvWp2C\ntVMlxwgJgMX4Bnuv7SvRtujSSDetDlbBCAmAlShpNNL55fIK2jRye7zu617dtLqA7Ld2qg0QSACs\nRJdGV6bNXv5FsfhHaYnEOIduWl1NQVp2cixrp1oCgQTAMvrK8/VplPBL7Q66NCrMdPrPWWDtVGlx\nDwmAZWhvFPmnkU52cmxNQVroO4WgIZAAWMnr+R8JIa5Mm/2zR7YqLRW5SRW5SbpniRLjHC3Fj5vQ\nP9wHAgmAlbR2X09d8O7PHtna7khT60osWxiru1jn9nipK2E5Nr+HpDyTzHoNgG0owdPuSBP/uBvU\n2n0j4LQ6pa4EkxdUMq/RoLB5IBFFgJ3415Vwe7yb6j+daP9jXdcJJBVLBwFA0OjqSgghNtWf016s\ny06ObSnOUN9SV8JaCCQAlqEd7lQ29eiKwGYnx7YUP05dCesikABYxkbNQ0Wt3dd1RWCVaXVhWFfC\nNggkAJahKz2uVVOQ7vZ4c/af0e1vSL8QHDaf1ADAZlqKH/efVpedHHus63rOVxdUpa6E5RBIAKwk\n4CRv3eU76kpYFJfsAEhhtPPkHfdxe6a0kjcsikACYL7RzpO95d+5kB8/+W660uMBUVfCurhkB8Bk\nShopry/kx6c0Xg24m/8k75qCNOpK2AkjJABm8g32qmmkOP5/ntUtQ6dMn/N/5Eg3ia6nLIs0sjRG\nSABM4xvs7Sl+UtvS7kh74YFS4fFG/KPmnj9W8rYrAgmAOXyDvX3l+dqWGYuzXvii+I4frClID1mn\nYCYu2QEwgW+w99q+Et9fetWWGYuztPX3JtJSnMF1ObuyeSClpqbKvLQtEJ6UNNIWI293pM37olh7\nma4iN8l/UQbS6H7Ifz6MGB8fN7sPoZKamkr5CUBCfeX52jSKmj1v+dxfqot2J8Y5Wl7MEEIErLCn\nFOUzpp+2JPOJ0eYjJACy8U+jpLfbtCUkesqyEuMcEz1ypJTdC3kvYQYCCYBxAqaRf9k9/0ne2st3\nx7q+3AQ7IZAAGMQ/jea8tEf4ld3TpZEyybs878v1Fyi7Z1cEEgCD+KfR9MVPia+W3cvZ36FPoxcz\nhBCU3QsHBBIAg6Q0XlWXBVLTSHy17J5WYpyjpiBdWY6BsnvhgAdjARgqpfHqaOdJNY3EP8ruaQdG\nCmVhOrfHq5vgQNk9uyKQABhNm0YTUWruHeu6XkHZvbBBIAEw2ab6c/7DI13NPUHZvTDAPSQAZtpU\nf057fwjhjEACYBr/NJroihxl98IBl+wAmKOyqUeXRspSdeW5SZTdC08EEgATHDjdr5utoC6c6l92\nz9CewTxcsgMQQm6P17/864HT/bqZ3CzjDcEICUDoVDb1KMMgbV2JxDiHdilVIURNQRppBMEICUBQ\n+AZ7/Rt1F+UU/mnEo0VQEEgA7pdvsLen+MkL+fF3+0HSCFoEEoD7oqSR8lrJpE315yJKj96x/Ctp\nBB0CCcC9U4qRa1v6yvN1y3X3lGUtWxhgqbpN9ecC1oRF2LJ5IMlfQx6wLiWNtEUlZizOmlvZqCv/\n6r7uzdnfEfAbKP9qJPnPhzYPpPPnz0tbPR6wOl0aRc2eN7eyUVf+1e3Rp1FhppPyr6aQ/3xo80AC\nECIBi5ELv/KvuueNCjOdNQVplH9FQAQSgLs2URqJScu/ZifH1hSkCcq/YgIEEoC7459Gc17ao76d\nqPxrdnJsS/HjymvKvyIgAgnAXQiYRv7lX3WfSoxzKGMjt8ebs/+MblMo+wsrYekgAFOlSyMhhC6N\nFDUFabr53Imx05XLdJR/xSQIJABTMvDWK7o0mld5KGAxct1EBkH5V0wNl+wATMlwa4P27URplLP/\njP8zsMBUEEgApiSl8ar6eqI0qmzq0S3TEPCrKP+KgLhkB2CqUhqvXsiPnyiNWrtvaG8RJcY5Wl7M\nEF+NJcq/YhIEEoC7oB0nafmvyFBTkO4/QqL8KybBJTsAQaCbyEDNPdwDAgnA/dJNZMhOjmU+N+4B\ngQTgvugmMmhXZADuCoEE4N4FmMhAGuFeEUgA7lFr9w3/iQxmdQY2QCABuEeb6j/Vvm0pzmAiA+4H\ngQTgXuTsP6OtDJudHEsa4T4RSACEb7D3rvZnIgNCgUACwp1vsPfavpIL+fH+m9wer27dbrfHy0QG\nhAgrNQDh7tq+EmUZ7wv58dqFGCqbepTgiSg9qjYmxjm0V+oEExkQPIyQgLCmK3HU8+KT6mtd7SKF\nLo2YyIAgsl4gnT179qmnnhoZGTG7I4Dl+Zd/TXq7beofr8hNIo0QRBa7ZHft2rWf//znHo9nfHzc\n7L4A1jZRGm2qP3fgdL92T6VahG7AlJ0cW55HFQkEk2UC6f3336+rq+vq6rp9+7bZfQEszz+N5ry0\nR3mtK2jU8mJGYpyjtfuG+GogtXZfT6r6mNW7EUSWCaQlS5Zs3bpVCHHq1KnGxkazuwNYWMA0Uksc\nae8SKXnjvyKDumdr9w2u2iFYLBNIixYtWrRokRDC6/USSMA9G3jrFW0aCSG0aaS9WKes2O2fRtqJ\ndse6rhNICBbrTWoAcM8G3npluLVB26Ir/1p7ekB9vTHT6Z9GFblJ2nneBz75yt0m4H5IN0L6/PPP\nDx8+rL6dPXv2ypUrTewPYBt/a9g9eRoJIbKTH1bvIW2q/1Q3yVuZyKCtxVeey7wGBE2oAmloaGjN\nmjWbN2/esGGD/9bm5ubDhw93dXXFxMSkpKQUFRUtWLBA2eTxeHbu3KnumZGRQSAB92+4peFvDbu1\nLf5pJITYmOlUZ9P5p5GyIoP2sl52cmxIuouwFKpAqq6uvnLlys2bN/03VVVV1dXVPfjgg+np6cPD\nww0NDb/73e/27du3dOlSIcSCBQvOnDkTol4BYWtg3yvatwHTSAiRGOfITo7VTrRTFGY6awrS3B6v\nrlR5Ypwj6F1F2ApmII2NjfX29l6+fPmDDz44cuRIwH3a2trq6urmz59/4MCBhIQEIURTU1NJSUlZ\nWVlTU5PDwcENhERK41V1tbo52/YETCPhNypSJcY61JWEVNQpR3AFM5AuXbq0evXqyfc5dOiQEGL7\n9u1KGgkh8vLyVq1a5XK5Tpw4sWLFiiD2B4CWkklztu15KGd9wB3cHm/O2x0BM0kXReO7l4ekiwhv\nwQwkp9P55ptvKq/Pnj1bXV3tv09bW1tkZOSyZcu0jStWrHC5XKdOnQp6IKWmpqqvz58/H9wvByxH\nu3aqzkQPG8EGtGdCmQUzkGJiYvLy8v7+vdMCfLPP5xscHJw7d67u0lxycrIQoq+vbyq/sm7dunXr\n1k2xS4QQMBUBHzbKTo7VrSGkqGBmndVoz4Qyh5Oh076Hh4dv3749a9YsXbvSMjQ0ZGRnACj8bw6p\nE+rKc5O00xaSqj6uKUjnSViEiKGB5PP5RKDBk9KibAVgJP+lVLXlX3WT6Fi5DiFlaCBFRkaKQMGj\ntChbARgmZ/8Z3QxvZXq3Wf1BmDN06aDp06cLIUZHR3XtSouyFYAB3B6vfxpV5CaRRjCRoSOkmJiY\nmTNn9vb2jo2NacdDbrdbCOF0Bv+ZBuX2HVMbAK2A07up/Wp7Mk9nUBi9uGpaWtqtW7c6Ozu1jR0d\nHUKI9PT0CT50786fP08aAVqt3TeSqj7WplFinIM0Cgfynw+NDqTc3FwhxK5du9SW/v7+gwcPRkZG\nLl/Oo3ZAaAWc3s3EOUjC6NW+CwoK6uvr29vb165d++yzz167ds3lco2MjGzZsiU+Pt7gzgD20Fee\nP7fyzkXCJpneDcjA6BFSVFRUbW1tXl7eZ599tnPnztra2tHR0dLS0pKSEoN7AtiDUv5VWafO7fEm\nVX2l+J56aW5T/TnSCJKLGB8fN+WHfT6f2+2Ojo52Op0RERGh+An1Dp7kl02Be6YrRp664F3dDolx\njsInnK3dN5jeDflPiaYFkgFSU1Ol/XcH7p8uja5Mm/2zR7a2O+4cMxW5SeV5LP8TpmQ+MUpXMRbA\nVOjSSAgxxTRiQh2kRSAB1jPw1iu6NHrhn8vUNFIWP/W/TMeEOkiOQAIsZuCtV4ZbG7Qtaholxjla\nXsxIjHO0dt/QTWFQkEaQmdGz7ADcj+GWhonSSAjRU5aVGOc4cLo/YGUjt8fb2n3DiF4C94RAAixj\nuKVhYN8r2paPVr+lppFSUPzA6f5N9ecm+oZjXdcn2gSYzuaBlJqaKv/yTcBUjHae1KXRvMpDDSNf\nTpbbmOn0T6PCTGdLcYb69sAnAQruIUzIfz60+T0kaWc3AndltPNkb/l3tC1ztu2Zvvip7L4edebC\npvpPdeulKtO7tRFVTrHXMKacD2XOJJuPkAB78E+jh3LWCyE2Zn65Rr4ujWoK0pSHjbT197KTY0Pb\nUeA+EEiABaQ0XlVfP7K+VEkjIURinCNgxtQUpBVmOpWiR9p2XQVYQCoEEmANSiY9lL3+kfWl2nb/\nNRcS4xyXPd7Kpp6kqo+1jyIVZga/5BgQRDa/hwTYiXacpPAvJyGEcHu8uoeQxndT2wUWwAgJsKqA\naQRYF4EEWJJ/Gk00YaGCmXWwCJtfslMmODL5Gzbjn0ZKOQm3x6udtpBU9TGL10El84RvBeUnAIuZ\nKI3M6g+sReYTI5fsACsJeKWONII9EEiAZQRMI8qQwzYIJMAaSCPYHoEEWIDb4yWNYHsEEmCyvvL8\nyXdwe7xJVV+pD0sawZYIJMBMfeX5I50fX8iPn2gH0gjhg0ACTDPw1isjnX8Pm4CZRBohrPBgLGCO\ngbde0RUjH+08OX3xU+pb0gjBJf+DsTYPJKIIcvpbw25dGs2rPKRLI13tV9II94kCfQD0hlsa/taw\nW9sSMI20lSNII4QDAgkw1HBLw8C+V7Qtd0yjxDgHaYRwQCABxhntPKlLoznb9twxjXrKsozrImAe\nAgkwyGjnyd7y72hb5mzboxYjF6QRwh6BBBjBP40eWV+qTSMhBGmEMEcgASHnn0YPZa9/ZH2ptiVn\n/xnSCGGOQAJCK2AazXlpj7bFP41aXswwqH+ANAgkILSmL34qpfGq+nbG4qw7plFNQbq28CsQJggk\nwAhKJs1YnDW3slHbHjCNKDqO8GTzEubKC9ZrgJxIIxhJ/lMiSwcB5iCNYDCWDgIQgC6NhBCkEUAg\nAaHiv1y32+MVgdKopTiDNAJsfskOMEtlU0/Fhz1CiIjSo2qjMndOiSUVaQQoCCQgJJQ00tFFkSCN\nAA0u2QGmIY0ALUZIQDBtqj934HS/tqUiN0kEGjCRRoAOIyQgmPzXoyvPS3Jf11+pE0Jsqv/UwH4B\nFkAgAfdl4K2v1DfS3iXqKctKjHP4j5nUPVu7b4S8f4B1EEjAvRt465Xh1oYL+fHKW23wFGY6hRCV\nTT26NNIuUnes6yuTv4EwRyAB92i4pWG4tUF5rWRS7ekBdevGTOeB0/26W0ctxRk1Benq2wOfBBg5\nAWGLQALuxXBLg64Y+WjnSe0khU31n26qP6fdQZnFUKsZMJXnJoW6n4CFEEjAXRvtPKlLo3mVh6Yv\nfmpjplNt0T1yVFOQpsSV9gpednJsiHsKWInNAyk1NVXmlQRhRf4F9+Zs2zN98VNCiMQ4R8CMqSlI\nK8x0uj3enP1ntO0UPYKR5D8f2vw5JFb7RnAFTKOHctarb8vzklr3f2WqQmKc47LHq64kpCrUDKcA\nA8i/2rfNAwkIIv80emR9qTaNWrtv5Ozv0H3K7fHqomh89/LQdRKwLptfsgOCxT+NHspe/8j6UvVt\nwDQCMHUEEnBnAdNozkt71Lf+aTTRhIUKZtYBE+CSHXAHvsFeXRrNWJw1eRoVZjprCtLcHq922kJS\n1cdU4QMmQSABk/EN9vYUP6ltmbE4a25lo/p2ojQSfpPoesqyQtlTwPK4ZAdMJurRedq3d0yj7ORY\nJY0A3C0CCbiDlMaryouppFFL8eOGdg6wEQIJuLOUxqukERBqBBIwJdo0cnu8pBEQdAQScHfcHm9S\n1cfaFtIICAoCCbgLpBEQOgQSMFWkERBSBBIwJW6PV1ffiDQCgotAAu5MSaPW7i+X8SaNgKAjkIA7\n8E+jxDgHaQQEHYEETCZgGrEIEBAKBBLC2nBLw+Q7kEaAYQgkhK/hloaBfa9cyI+faIec/WdII8Aw\nEePj42b3IVTUSr0UMoc/XYkjdcE6FWkEm5H/lGjzEdL58+el/aeHifwL7umu3fmnUU1BukGdA0JD\n/vOhzQMJ8OefRo+sL30oZ736NmAaUVgPCDUCCeElYDHyR9aXqm9JI8AsBBLCSMA00hYjJ40AExFI\nCBe+wV5dGs1YnDVJGgkhSCPASAQSwoJvsLen+Elti67gnn8atRRnkEaAkQgk2B9pBFgCgQSb8w32\nXttXom3RpZFuLQZBGgEmIZBgZ0oajXR+WcTIP40OnO7XfoQ0AsxCIMHOoh6dp02jqNnzSCNAWgQS\nbE5dEyhq9rykt9vU9sqmHtIIkAqBBPtLabyqS6MDp/srPuzR7kMaAaYjkBAWdGmkK0ZOGgEyIJAQ\nXvzTqKYgjTQCZEAgIYwETKPCTKdZ/QGgRSAhXLR23yCNAJkRSAgLrd03cvZ3aFtII0A2BBLszz+N\nKnKTSCNANgQSbM4/jQozneV5SWb1B8BECCTYWcA0qilIM6s/ACZBIMG2SCPAWggkWNVo58kL+fET\nbSWNAMshkGBJajHygJnkn0bZybGkESA5KwXSsWPHfvjDH65cufL555+vrq6+deuW2T2COdQ0UlzI\njx946xX1bcA0ail+3Lj+Abgnlgmk9957b+vWrdOnTy8qKlqyZMnevXtffvnl8fFxs/sFo+nSSAgx\nY3HWnJf2KK/dHi9pBFjUNLM7MCUXL16sqKioqanJyspSWtatW/fcc8+5XK7Vq1eb2zcYZvJK5G6P\nN6nqY+1WogiwFmsEUkJCQmNj48KFC9WWP//5z0KIpCSeJgkXvsHevvJ8bQtpBNiMNQJpxowZS5Ys\nUd9evHjxjTfeePLJJxcvXmxir2AYpRK57y+9aosujXLe5jIdYHnSBdLnn39++PBh9e3s2bNXrlyp\nvh0fH3///ff/4z/+4+tf//revXvN6CBMcG1fibYSuS6NNtWfc3u86lbSCLCoUAXS0NDQmjVrNm/e\nvGHDBv+tzc3Nhw8f7urqiomJSUlJKSoqWrBggbLJ4/Hs3LlT3TMjI0MNpKtXr7766quffPLJli1b\niouLo6KiQtR5SKWvPF+bRlGz5+nSqLX7urqVNAKsK1SBVF1dfeXKlZs3b/pvqqqqqqure/DBB9PT\n04eHhxsaGn73u9/t27dv6dKlQogFCxacOXPG/1MXL14sLCyMj4///e9///Wvfz1E3YZs/NNIrf3q\nn0aJcQ7SCLCuYAbS2NhYb2/v5cuXP/jggyNHjgTcp62tra6ubv78+QcOHEhISBBCNDU1lZSUlJWV\nNTU1ORyOgJ+6ffv29u3bFy5c+Jvf/OaBBx4IYp8huYnGRkII/zTqKcsytHMAgiqYgXTp0qU7TsI+\ndOiQEGL79u1KGgkh8vLyVq1a5XK5Tpw4sWLFioCfOnny5GefffbTn/60vb1d2/7YY489/DDFp+0s\npfGqshZD1Ox5c17aE/XoPKU9Z/8Z0giwmWAGktPpfPPNN5XXZ8+era6u9t+nra0tMjJy2bJl2sYV\nK1a4XK5Tp05NFEh/+tOfhBDae0uKd955R7nQN5HU1FT19fnz56fwR0A6KY1Xe158cs5Le6Yvfkpp\nIY2Au6I9E8osmIEUExOTl5f39++dFuCbfT7f4ODg3LlzdZfmkpOThRB9fX0TeksKJwAADndJREFU\nfXNxcXFxcfE9dIkQsi5lMreSNMp9I7fHmxjn8E+jmoJ003oJWIH2TChzOBk67Xt4ePj27duzZs3S\ntSstQ0NDRnYGMqts6qn4sEcIEVF6VG1MjHMIIbQzvJU0yk7msi1gB4YGks/nE4EGT0qLshUQQihp\npKONIkEaAbZj6OKqkZGRIlDwKC3KVmAqSCPAfgwdIU2fPl0IMTo6qmtXWpStwaVcLeVOklVsqj93\n4HS/tqUiN0kEGjCRRsDdkvnukcLQEVJMTMzMmTN7e3vHxsa07W63WwjhdDqD/ovnz58njSzEf+5c\neV5Sa/cN/z031X9qYL8AO5D/fGh0PaS0tLRbt251dnZqGzs6OoQQ6enMlQp32rtEPWVZ/nPqtHsG\nDCoA1mV0IOXm5gohdu3apbb09/cfPHgwMjJy+fLlBncGUtFerCvMdAohWrtv6NJImWinONYVIKgA\nWJfRgVRQULBw4cL29va1a9e+8847v/jFL/Lz80dGRjZv3hwfH29wZyCV2tMD6uuNmU4hRHbywzUF\naWpjS3GG9pGjA5985W4TAKszOpCioqJqa2vz8vI+++yznTt31tbWjo6OlpaWlpSUGNwTyEY7SaGy\n6e+zGAoznUomtRRnZCc/XKsZRZXnUp4RsJWI8fFxU37Y5/O53e7o6Gin0xkRERGKn1CnlEh+Hw8K\nbdXXiVYD0j4nq9xkMqhzgPXJf0o0eoSkioqK+sY3vhEfHx+iNFLIP6sEqsQ4R3ZyrPLavyS52+PN\n2X9Gt79xnQOsT/7zoXQVYxHOlGpGyjDI7fGqF+6E36NIyqwHAHZCIEFeARcQGt/NbEzAnky7ZAdM\npILZCkBYIpAgnfK8pIkyiawCbMy0WXYGkH9KCSahVD9S3yZVfcz6dcD9kP+UaPNAkvbfHQBMIfOJ\nkUt2AAApEEgAACkQSAAAKRBIAAApEEgAACkQSAAAKdh86SBl3r20cxwBwDDqc0jSsnkgEUUAoFDO\nhzLHEpfsAABSIJAAAFIgkAAAUiCQAABSIJAAAFIgkAAAUiCQAABSsPlzSDwYCwAKmZ9AUtg8kIgi\nAFDwYCwAAFNCIAEApEAgAQCkQCABAKRAIAEApEAgAQCkQCABAKRAIAEApEAgAQCkYPOVGlg6CAAU\nMq/RoLB5IBFFAKBg6SAAAKaEQAIASIFAAgBIgUACAEiBQAIASIFAAgBIgUACAEiBQAIASIFAAgBI\ngUACAEiBQAIASIFAAgBIweaLq7LaNwAoZF5WVWHzQCKKAEDBat8AAEwJgQQAkAKBBACQAoEEAJAC\ngQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEEAJACgQQAkAKBBACQAoEE\nAJACgQQAkAKBBACQAoEEAJDCNLM7EFpK9XilkjwAhDPlfCgzmwcSUQQACuV8KHMscckOACAFAgkA\nIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAF\nAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAghWlmd+AuuFyu999///Lly48+\n+ugzzzzz3e9+NyoqyuxOAQCCwzIjpNra2tLS0oSEhO9///sJCQk7duyorKw0u1MAgKCxxghpbGxs\nz549RUVFP/nJT4QQ//Iv/zI2NvZf//Vf//Zv/+ZwOMzuHQAgCKwxQurv74+Pj1+5cqXa8s1vflMI\ncevWLfM6BQAIJmuMkObOnetyudS33d3d77777qpVqx566CETewUACCJrjJBUp06devrpp5955pnY\n2NidO3ea3R0AQNBEjI+Pm92Hr/j8888PHz6svp09e7b2St1f//rXs2fPut3u3/72t4mJifv27Zs1\na9ZEX5Wamnr+/PnQdhcALEXmE2OoLtkNDQ2tWbNm8+bNGzZs8N/a3Nx8+PDhrq6umJiYlJSUoqKi\nBQsWKJs8Ho926JORkaENpK997WvLly8XQmRlZX3729/+6KOP1q5dG6I/AQgimc8CgCRCdcmuurr6\nypUrN2/e9N9UVVW1bdu2o0ePxsTEDA8PNzQ0PPfcc8ePH1e2Lliw4IzGO++8I4Robm4uKioaHh5W\nv2TRokVxcXEdHR0h6j8AwGDBDKSxsTG3233s2LGSkpJf//rXAfdpa2urq6ubP3/+H/7wh/r6+iNH\njuzdu9fn85WVlXm93om+edq0acePH//www/VloGBgevXr6vjKrtKTU21wY/e5xfew8fv6iNT2TlY\n+1ia8X+gbIfivX3D1D8yxT3vuJt1D8VgXrK7dOnS6tWrJ9/n0KFDQojt27cnJCQoLXl5eatWrXK5\nXCdOnFixYkXAT2VlZS1cuHDPnj0ejycvL6+vr2/Hjh2zZs26488BAKwimJMavvjiixMnTiivz549\nW11d/eMf//gHP/iBdp9ly5b95S9/6ejo0D7Q6nK5fvzjH3/ve98rKyub6Mt7e3srKyuPHz+udPiJ\nJ5549dVXFy9ePEl/rPu/CQAQOtLezgzmCCkmJiYvL+/v3zstwDf7fL7BwcG5c+fqlldITk4WQvT1\n9U3y5fPmzauurvZ6vQMDA48++uiMGTPu2B9p/9EBAP4MfTB2eHj49u3b/hO1lZahoaE7foPD4UhM\nTAxF3wAA5jL0wVifzycCDZ6UFmUrACA8GRpIkZGRIlDwKC3KVgBAeDI0kKZPny6EGB0d1bUrLcpW\nAEB4MjSQYmJiZs6c2dvbOzY2pm13u91CCKfTaWRnAABSMXpx1bS0tFu3bnV2dmoblQUX0tPTDe4M\nAEAeRgdSbm6uEGLXrl1qS39//8GDByMjI5VF6gAA4cnoekgFBQX19fXt7e1r16599tlnr1275nK5\nRkZGtmzZEh8fb3BnAADyMDqQoqKiamtrX3vttebmZuXCXXR0dGlpaVFRkcE9AQBIxbR6SD6fz+12\nR0dHO53OiIgIU/oAAJCHdAX6AADhyWIlzIPr7NmzTz311MjIiNkdQfg6duzYD3/4w5UrVz7//PPV\n1dW3bt0yu0cIXy6Xq7CwMCcn5/nnn6+trTV+9ZzwDaRr1679/Oc/93g8jBFhlvfee2/r1q3Tp08v\nKipasmTJ3r17X375ZQ5ImKK2tra0tDQhIeH73/9+QkLCjh07KisrDe5DOF6ye//99+vq6rq6um7f\nvi2E6OjoiI6ONrtTCDujo6NLly7Ny8v7xS9+obS89957FRUVhw4deuyxx8ztG8LN2NjYE0888d3v\nfvcnP/mJ0vKjH/3oj3/8o65UUKiF4whpyZIlW7du3blzZ35+vtl9Qfjq6+v74osvtEUmMzIyxD8W\nLgGM1N/fHx8fv3LlSrXlm9/8phDC4GvIRk/7lsGiRYsWLVokhPB6vY2NjWZ3B2EqISGhsbFx4cKF\nasuf//xnIURSUpJ5nUKYmjt3rsvlUt92d3e/++67q1ateuihh4zsRjiOkAAZzJgxY8mSJeqawhcv\nXnzjjTeefPLJyesgAyF16tSpp59++plnnomNjd25c6fBv04gASYbHx9vaGh4/vnn4+Pj9+7da3Z3\nENYWLlz47//+7//6r/86ODhYWFg4lbqpQWTVS3ZDQ0Nr1qzZvHnzhg0b/Lc2NzcfPny4q6srJiYm\nJSWlqKhowYIFxncS4eA+D8WrV6+++uqrn3zyyZYtW4qLi6OioozqOGzo/k+MX/va15RlRbOysr79\n7W9/9NFHa9euNaLrQgjrjpCqq6uvXLly8+ZN/01VVVXbtm07evRoTEzM8PBwQ0PDc889d/z4ceM7\niXBwP4fixYsX161bd/Pmzd///vc/+tGPSCPcp3s+Gpubm4uKioaHh9X9Fy1aFBcXp5RiMIyVRkhj\nY2O9vb2XL1/+4IMPjhw5EnCftra2urq6+fPnHzhwICEhQQjR1NRUUlJSVlbW1NRk5PxF2FhQDsXb\nt29v37594cKFv/nNbx544AFj/wLYR1COxmnTph0/fvzDDz9U5x4PDAxcv37d4GtLVgqkS5cuaefI\nBnTo0CEhxPbt25V/dCFEXl7eqlWrXC7XiRMnVqxYEfJeIgwE5VA8efLkZ5999tOf/rS9vV37wcce\ne+zhhx8OUc9hP0E5GrOyshYuXLhnzx6Px5OXl9fX17djx45Zs2bd8ZuDy0qB5HQ633zzTeX12bNn\nq6ur/fdpa2uLjIxctmyZtnHFihUul+vUqVMEEoIiKIfin/70JyGE/0Smd955Z+nSpaHpOGwoKEfj\nAw888Ktf/aqysvKNN97YvXu3EOKJJ554/fXXDS7kbaVAiomJycvLU15Pmxag5z6fb3BwcO7cubpL\nc8nJyUKIvr4+3f7r1q1bt25daDoLOwvKoVhcXFxcXBz6zsLmgnVinDdvXnV1tdfrHRgYePTRR2fM\nmBHijgdgpUC6o+Hh4du3b8+aNUvXrrQYPH8R4YxDEfK4q6PR4XAkJiYa1jcdq86yC0hZm9b//xGU\nFuNXrkXY4lCEPCx0NNoqkCIjI0Wgf1+lRdkKGIBDEfKw0NFoq0BSVmEZHR3VtSst6hotQKhxKEIe\nFjoabRVIMTExM2fO7O3tHRsb07YryycbPF0E4YxDEfKw0NFoq0ASQqSlpd26dauzs1PbqDxsnJ6e\nblKnEI44FCEPqxyNdguk3NxcIcSuXbvUlv7+/oMHD0ZGRioLNAHG4FCEPKxyNNpq2rcQoqCgoL6+\nvr29fe3atc8+++y1a9dcLtfIyMiWLVvi4+PN7h3CCIci5GGVo9FugRQVFVVbW/vaa681Nzcr49Po\n6OjS0tKioiKzu4bwwqEIeVjlaIwYHx83uw8h4fP53G53dHS00+mMiIgwuzsIXxyKkIfkR6NtAwkA\nYC12m9QAALAoAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAgBQIJACAFAgkAIAUCCQAghf8PrVRvvgRu\nDAoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_ = 50:50:500;\n",
    "time_ = [];\n",
    "for k = 1:length(n_)\n",
    "    n = n_(k);\n",
    "    A = rand(1200,n);\n",
    "    Q = zeros(1200,n);  R = zeros(600,600); \n",
    "    \n",
    "    tic\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\n",
    "    time_(k) = toc;\n",
    "end\n",
    "loglog(n_,time_,'-o',n_,(n_/500).^2,'--')\n",
    "axis tight, xlabel('n'), ylabel('elapsed time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flops aren't everything. On massively parallel computers, for example, the time needed for memory access and communication often competes with or dwarfs the floating point arithmetic. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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
}