chapter16.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第十六讲：投影矩阵和最小二乘\n",
    "\n",
    "上一讲中，我们知道了投影矩阵$P=A(A^TA)^{-1}A^T$，$Pb$将会把向量投影在$A$的列空间中。\n",
    "\n",
    "举两个极端的例子： \n",
    "* 如果$b\\in C(A)$，则$Pb=b$；\n",
    "* 如果$b\\bot C(A)$，则$Pb=0$。\n",
    "\n",
    "一般情况下，$b$将会有一个垂直于$A$的分量，有一个在$A$列空间中的分量，投影的作用就是去掉垂直分量而保留列空间中的分量。\n",
    "\n",
    "在第一个极端情况中，如果$b\\in C(A)$则有$b=Ax$。带入投影矩阵$p=Pb=A(A^TA)^{-1}A^TAx=Ax$，得证。\n",
    "\n",
    "在第二个极端情况中，如果$b\\bot C(A)$则有$b\\in N(A^T)$，即$A^Tb=0$。则$p=Pb=A(A^TA)^{-1}A^Tb=0$，得证。\n",
    "\n",
    "向量$b$投影后，有$b=e+p, p=Pb, e=(I-P)b$，这里的$p$是$b$在$C(A)$中的分量，而$e$是$b$在$N(A^T)$中的分量。\n",
    "\n",
    "回到上一讲最后提到的例题：\n",
    "\n",
    "我们需要找到距离图中三个点 $(1, 1), (2, 2), (3, 2)$ 偏差最小的直线：$y=C+Dt$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAFRCAYAAACot654AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VFW+9vHnVFIJAZIQAZkFB0AUHJgaRZyDBFtFJkXm\nRmn10g003aDSF0REUBGHV1AvtFdtB+4S00RAbERAJBihGSIqIKIIBklnIAkhCalK1ftHOsEQhqRS\nVbuG72ctl8mpU+f8dirkqb1rn7Mtt9vtFgAAMMJmugAAAMIZQQwAgEEEMQAABhHEAAAYRBADAGAQ\nQQwAgEGR/jqR01mmo0eL/HW6gJKQUD9s296tW2fZbJa2bt1luhRjwvn1l2h/OLc/nNsuSU2bxtZo\nP7/1iCMjI/x1qoATzm0Hrz/tD9/2h3Pba4OhaQAADCKIAQAwiCAGAMAgghgAAIMIYgAADCKIAQAw\niCAGAMAgghgAAIMIYgAADCKIAQAwiCAGAMAgghgAAIMIYgAADCKIAQAwiCAGAMAgghgAAIMIYgAA\nDCKIAQAwiCAGAMAgghgAAIMIYgAADCKIAQAwiCAGAMAgghgAAIMIYgAADCKIAQAwiCAGAMAgghgA\nAIMIYgAADCKIAQAwiCAGAMAgghgAAIMIYgAADCKIAQAwiCAGAMAgghgAAIMiPX2iy+XSX//6V/34\n44+y2WyaNWuWLrnkEm/WBgBAyPO4R7xu3TpZlqX33ntPEydO1IIFC7xZFwAAYcHjHvGtt96qm2++\nWZKUkZGh+Ph4rxUFAEC48DiIJclms+mRRx7R2rVr9dJLL3mrJgAAwobldrvddT1ITk6OhgwZoo8+\n+kj16tXzRl0IEe3atZMkHThwwGgdABCoPO4Rp6SkKDMzU+PHj1d0dLRsNptstrN/5JyVdczT0wW1\npk1jw7btLpdbNpsVtu2Xwvv1l2h/OLc/nNsulbe/JjwO4r59++rRRx/ViBEj5HQ6NX36dEVFRXl6\nOAAAwpLHQRwTE6MXXnjBm7UAABB2uKEHAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMA\nYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGebz6EgDgpKIi6Y037Nq6NUIjRjiUk2Np\n164IDR0qXXihtHJlpOx26YsvIvTssydkWaYrRqCgRwwAXrB6daRGj3YoK8tSSYmloUOdGj26VJMn\nl4fv7t0RGjTIqfT0CO3ezZ9enESPGAC8oG9fp+x26cABm5KSnJKkjAybcnKkxMQyXXttmU6cKN/3\n4otdBitFoOFtGQB4QWystH17hLp2LZPtP39Z16+PVGJi+df5+ZZeeCFKM2acUHS0uTr95fuj+7Ri\n7wrTZQQFghgAvCQ1NUKdOpX3drOzLa1ZE6FZs8ofa9nSrWnTSvXCC1HKzTVYpI99f3SfHvrkft25\n/DYdKTxiupygQBADgJds3hwhSUpOjtSLL0ZpyZISXXBB1X2aNnUrLS30PhX8dQB3PO9SfTl8px7o\n9oDpsoJC6P02AIABDoe0e7dNy5YVy7KkgQOdlY898USU2rZ1a/RohzIyLLVpEzqfEX9/dJ+e+9fT\n+uzndRp/xcN65oYFio2KM11WUCGIAcALtm0rH5Y+3WVJ99zj1Lff2vTWW3bdcYdTXboEfxATwN5D\nEANAHe3ebdP8+VHKzbW0YUOEbryxrMrjHTu61LFj8IevRAD7AkEMAHXUqZNLy5YVmy7Dpwhg3yGI\nAQBnRAD7HkEMAKiGAPYfghgAUIkA9j+CGABAABtEEAOAj5SVlZ17J8MIYPO4sxYA+MjatWtMl3BG\np7sT1qRufyaEDaBHDABhhB5w4CGIASAMEMCBiyAGgBBGAAc+ghgAQhABHDwIYgAIIQRw8CGIASAE\nEMDBiyAGgCBGAAc/ghgAghABHDoIYgAIIgRw6CGIASAIEMChiyAGgABGAIc+ghgAvKi4uFgxMTHV\nthcVFal+/fo1Pg4BHD5Y9AEAvOjgwZ+0e/e3p9n2TY2ez2IM4YcgBgAv6tjxUn333Z4q23bs2KZu\n3Xqc9XkEcPhiaBoAvOzSSy/Tt9+W94B/+ukntWrV+oz7MgQNesQA4GUdO16qffv2SpK+/PJLde/e\ns9o+9IBRgSAGAB+47LLO2r59m9q0aVNlOwGMUzE0DQA+0L59B8XE1NM111yjrKxjDEHjjAhiAPCB\n3Nw8ffPNxepy0yvKu+LvKmrxnR686r8IYFRDEAOADzw8412ti9gl9VgjpU1U/0OjNWncfabLQgAi\niAHAiyqGoDde/LG0aZp2/WO7urin6/BVyyVJEbu+Uswbf5M7IUEqLZXtaK6OzXtOatDAcOUwhSAG\nAC849TPgvntu0kef/16dNV2SW23bFih66TuK+Z9XlP/eB3I3ayZJinn5RUVt2qjS25LMNgDGEMQA\nUAdnmoR1dF6e7K6/SynSXXf9XS8Ni1PsyPuVt/yjyhCO3LFN9rRUlQwbYbgVMIkgBgAPnGsWdEJC\nIy1efLeUMlqLF9+t+KED5G6UoOjVqxS9aoXkcqmsQ0cVvPq61LChJMn+RapsR36RVVQk+6aNKhk2\nQo7rbzTUQvgLQQwAteDRZUhOZ3mwjhit4zOeOONucb8bocLH56hk+Ci54uIVP+peZX+zn8+PQxw3\n9AAMyc/P05gx98nhcEiS/vWvLXrggVFKTLxew4cP1sqVKTU+VnFxsRYseFp3391f/fvfounT/6Ls\n7CxJksPh0Nix9ykvL88n7QgXdbkRh5WbK5WVqazthWfdL2/5ap248+7yb9wuyen0RukIcB4FsdPp\n1NSpUzV8+HANHTpU69at83ZdQMh79dWXNWjQUNntdh06dFDTpk3WDTfcrDfeeFejR9+vBQue0ebN\nm2p0rBdfnK/09B168slntHDhYp04cUKPPjpFkmS32zVkyDAtWvSiL5sTsrxxJyx348Zyx8ZJZdWD\n1crMVNTKDyVJZR0vlf6zhGL0qg9VNGUaveEw4FEQf/jhh0pISNA777yjxYsXa/bs2d6uCwhpmZlH\ntGHDOvXrd7skad26T9ShQ0eNGDFGrVq1Vt++/dSvX3+tWbP6nMdyOp1au/afmjBhsi6/vLMuvPAi\nPfLIDO3Zs1sHD/4kSerbN0mbNm1UZuYRr7YjL++oSkpKvHrMQOHVW1FGRKhk1FhFrVtbZXPkzu1q\nOHuGHNffcHLbjm2q//yzctdvoKKH/lDXZiAIePQZcVJSkvr16ydJcrlciozko2agNlasWK7u3XvK\nbrdLkm65pa+uuaZ3lX0sy1Jh4bEaHW/u3OfUpcuV1bZXPD8yMlI9evRUSkqyxo9/uI7Vn1SvXow+\n+2y9HA6HbrzxJjVsGOu1Y5viq1tRHp8+Uw1mz1Tsg7+Tq2VrqfSEyi5ur2MvvSLZTvaJnFd3k/Pq\nbqr35utqdMdtyktZLdWvX+fzI3B51COOiYlR/fr1VVhYqIkTJ2ry5MnergsIaWlpm9Wjx28qv2/d\nuo06dLi08vvc3Bx9+umaKvucSXnI/kb16tWr3Pb+++8pPr6R2rfvWLmtR4/fKC1ts5daUK5evXq6\n7bYk9evXX1u2pOnDD/+h7OxsSVJ+bq5WPjBGK3r21MoHRiv/aK5Xz+1t3l6MoaL9ksrbX5Cv448/\nqWOvvq7jM57Q8SefVsnY+ytDOHLbVjW+/BLZDh2UJDl691HkVzur9aIRejzuyv7yyy+aMGGCRowY\nof79+3uzJiCkuVwu7du3V23PMHGnpKRE06f/RU2bnq8BAwbX+vgbNnyqpUvf1iOP/Hdlj1uS2rW7\nSPv375PT6awyinXsWIHy8vLkcJSqtNQhh8Pxq69L5XCUyvmfSUNut2RZ5f+vYFmqfCwqKlovvDBf\nHTp01Hmff6YxKcmyJLm1VW/I0m8Xv1Hr9viar3rAn0/7k8akJEuSxqT849ztj4iQ89LL5Greovzb\nAz9IUVFyXt65zrUgsHkUxNnZ2Ro3bpxmzJihXr161fh5TZsG/7CVp8K17TZb+V/pcG1/hV+3Pzc3\nVy6XSxde2LLaz6VilCkz84jee+89tW7dpFbnWbVqlWbN+qvGjRunUaOGVXmsXbuWcrvdiox0qmnT\nhMrtpaUFiowsU/36MYqKildUVJSioqJkt9srv46MjJRVkbhn8O9//1vr16/X2LEj1atXL634v7dV\n8QxLUpPDhwLq92Bv9l7N3jhba/av0aRek/T6oMWKi/beYgxNDh+qXfsTb5CyHlDT/3uz/N1Naqq0\napUa96z+kUMwCaTXPFB5FMSvvfaaCgoKtGjRIi1cuFCWZWnJkiWKioo66/Oysmr2eVeoado0Nmzb\n7nK5ZbNZYdt+qfrrn59/XJKUnX1McXG/3p6nyZP/S3l5eXrppdcUHR1fq5/bihXLNX/+XN1zz3CN\nGjW+2nNzcgolSUePFkk6+VhUVJxatKgeQE6n5HS6VFR09slYBQX52rBhnc47r7FuuilJllX+eme3\nbC23tv6nRyxlt2wTEL8Hp/aA0+7bodioOJ0okLLkvfqyW7bWt9qqy1SL9ifecfLr4ePK/x8APzNP\nhfPfPqnmb0I8CuLp06dr+vTpnjwVCHvx8Y0UERGhgoL8ym1Op1N/+cskHTt2TAsXLlaLFi1rdczP\nPluvZ599SqNG/U733//gaffJy8uTZVlq1CjhtI97orDwmLZu/VK3336nIiIiqjzW55nn9YYsjU1J\n1ht3DVSfZxZ47bye8Pd6wH2eeV5jtm/T4NiGim9/qfH2I3Ax3RkwoH37Dtq/f5+uuqqrJGnp0rf1\n3Xd7tGDBy4qOjlZubo4kKTLSrri4ODmdThUU5KtRowTZbFXnWBYXF+vZZ+eod+8+GjhwSOVzJSku\nLr7y8+D9+/epQ4eO5xxiro2GDWN1yy19T/tYfMJ55Z+JpiQb/WzY3wFc4cDBnzRw0l/UvGUT3Xxz\nf6/+3BFaCGLAgF69eis9facGDbpHkrRhwzq5XC5NmlT10qIuXa7UwoWLtWtXuv74xwf1/vsr1Lx5\n8yr77NixTQUFBUpN/VypqeUr+LjdblmWpeefX6hu3XpIkr76aqeuueY6P7QuMJgK4AoHD/6k888/\nX4mJt+jTT9fppptu8du5EVwIYsCA22+/U2PGDFNxcbFiYmK0ZMlbZ93/6qu7KSnpt4qOrj4P49pr\nr9PGjVvO+vzi4mKlpW3WQwFwgwhfr8drOoAl6fDhDDVv3kJ5eUfVqFEjFRcXq6SkpMolZkAF7jUN\nGNC8eQtdf/1NWr16ZY32z8j4WXl5R5WQcJ5H5/v441Xq0+cGNWvW/Nw7+1D00ncUO/FhHZ/6mI7/\n9XEdf+IpOTt2UtSmjXU+trevA66Lbdu2qnv3nmrQoKEKCwt16619tXbtGr/XgeBAEAOGPPTQH5SS\nkly56MPZtGrVWvPmeTbZx+FwaPnyD/Tww3/06PneErn1S8VO+aMKn36u2nq8ju49PT5uIAVwhZYt\nW8myLDVs2FDHjh1TVFSUunS5wlg9CGyW2/3rS/N9K1ynsYfzFP5u3TrLZrO0desu06UYE86vvyQ1\nPT9OWf8uUPzQAYr85muV3HOf5HJVrsdbMmBQ5Xq8kdv/JfuXabKOFci+9UsV/WmqHKfc+rPCqUPQ\n47qMNxq+p/PDD98rLq6emjRpbboUI8L+d9+Xly8BQK3UZD3e4mJFr16l49NnSpKiVixX/LBByv1y\np1y/GlIPhM+Aa6pBg1gdO5arJrW7LwvCDEPTAHyuJuvxRvz4g2L+3/OyHfhRkuS46RapuFiRW9Ik\nBeYQ9LlUDE0DZ0OPGIDPnWs9XvvWL1X62zuVt3KNXO3Kw9qWkSFZlg40sWv2J/cHRQ/4VPXr11dR\nUZHpMhDg6BED8L0arsfr/NWkrbL5s5SS1F637vlD0PSAT8VNPFAT9IgB+EVN1+P9/ug+7Zj/eznz\nv1HmtL/oywCchAV4E0EMwD8iInT88SfP+HDFJKyoTz7WoIQbdfXbyxXrjpYtM1OuNsEbxJdddpnp\nEhDgCGIARv16FvRT6q8BnWfI6jdAOlos+7aNcjVrJlebC0yX6bGLL744rC/hwbkRxACMOPUypAUX\nTlLrvrfJKnpbmjpVcrsly1L29z+bLhXwKYIY8JKiIumNN+zaujVCI0Y4lJNjadeuCA0dKnXpUp4r\nM2dG64knTpgu1aizXQec80OG4eoA/2PWNOAlq1dHavRoh7KyLJWUWBo61KnRo0s1ebJUWCi99ppd\nX3wRce4DhahgvA4Y8AeCGPCSvn2dstulAwdsSkoqv142I8OmnJzyOzg++KBDsbF+u6Oscbm5eXrg\ngX9Ikro+eZt+m9w36AM4Pz9PY8bcV3l/8NTUzzV69L265Zbe+t3vhistbXONj1VcXKwFC57W3Xf3\nV//+t2j69L8oOztLUvn9wceOvU95eXk+aQcCC0EMeElsrLR9e4S6di2rvBpn/fpIJSaarcuUh2e8\nq5SI8hWHft5xu3ptmRO0AVzh1Vdf1qBBQ2W323XgwI/67/9+RAMGDNbbb7+vxMQkPfbYn/XLL4dr\ndKwXX5yv9PQdevLJZ7Rw4WKdOHFCjz46RZJkt9s1ZMgwLVr0oi+bgwBBEANelJoaoU6dXJKk7GxL\na9ZEaNYsw0X5WcUQ9MaLn5KyLtPXVifp8+k6/GMz06XVSWbmEW3YsE79+t0uSfr3vzM1ePA9uvvu\nwWrRoqWGDRuhevVi9M03517gxOl0au3af2rChMm6/PLOuvDCi/TIIzO0Z89uHTz4kySpb98kbdq0\nUZmZR3zaLphHEANetHlz+WfAycmRevHFKC1ZUqILgvfKm1o59TPgvnuelj6fpi7ubyW51bZtgekS\n62TFiuXq3r2n7Ha7JKlnz16VS0s6nU6tXLlcDodDnTvXbLnDuXOfU5cuV1bbXlhYfqlTZGSkevTo\nqZSUZC+1AIGKWdOAlzgc0u7dNi1bVizLkgYOrH5fZf8tOuo/Z5oFfXRenr7dkajY2KG66KL6euaZ\nm0yXWidpaZt15513V9t+8OABjRx5j9xutx58cIKaN29xzmOVh+xvqmx7//33FB/fSO3bd6zc1qPH\nb5ScvEzjxz9c9wYgYBHEgJds21Y+LH262wuXlEhvvWXXvn02vfaaXWPGOBQd7f8avelcyxEeP35M\no0bdrMsvP0833ZRksNK6c7lc2rdvr9qeZvWoxo2baMmSv2vXrnS9/PLzatWqjW64oXZvOjZs+FRL\nl76tRx7578oetyS1a3eR9u/fJ6fTqchI/lyHKoamAS/Yvdum+fOjlJtracOG6pco1asnjR/v0Ndf\nH9fvfx/cIVzTy5C2bduq9u076JJLLtEPP+w3VK135Ofny+VyqVGjRtUea9Cgodq376CBA4eof/87\n9MEH/1erY69d+0/NmvVXDRs2UklJv63yWFxcvNxut/LzmT0dyghiwAs6dXJp2bJirV9fpBtvLDNd\njk/U5jrgQ4cOqlWr1pKkrl27ateudH+X61U2W/kwR1nZydd2//7v9fXXX1XZr127C2sVmitWLNfs\n2TM0ZMgwPfjghGqPu//zWYbNxp/qUMarC+CsPLkRx7ZtW9W9e8/KYfp27S7Ujz/+4KeKvS8+vpEi\nIiJUUJBfuW3duk80f/68Kvvt3bvntMPXp/PZZ+v17LNPaeTIsZWTvk6Vl5cny7LUqFGC58Uj4BHE\nAE7L0zthFRTk64IL2ko6OTntyiuvrrwsJ1i1b99B+/fvq/z+9tvvVEbGz1q8+BX9/PMhvf/+Un36\n6ScaOXKspPKZ1Lm5OXK5XNWOVVxcrGefnaPevfto4MAhys3NqfzP6Tw5yW///n3q0KEj6xqHOIIY\nQBV1vRVlXFy8unbtXm17bScwBZpevXorPX1n5fctW7bSc8+9pC1bvtCYMcO0YsU/NGfO02rfvoMk\nadeudN11Vz8dOVL9OuAdO7apoKBAqamfa8CAJA0YkKS77uqnAQOSlJ6+o3K/r77aqWuuuc73jYNR\nTMMDIOncs6A9EUodudtvv1NjxgxTcXGxYmJiJElXXHGVFi9+67T7X311NyUl/VbR0dE6tVN87bXX\naePGLWc9X3FxsdLSNuuhh/7glfoRuOgRAz50pqHJQMJiDDXTvHkLXX/9TVq9emWN9s/I+Fl5eUfV\nuHFjj8738cer1KfPDWrWrLlHz0fwIIgBH/rqq/TKBQICjT8CONRuYPLQQ39QSkpyjV7TVq1aa968\nBR6dx+FwaPnyD844iQuhhaFpwIcsywq4HrEvhqDDRULCeXrzzfdqvL+nlx3Z7fZanQfBjSAGfMhm\ns1VeC2oaAQwEJoIY8CGbzWa8R0wAA4GNIAZ86OTQtP+nDxPAQHAgiAEfOjk07b8gDqQAbtKkiZHz\nAsGEIAZ8yLIqhqZ9f4FCIAVwhW7dehg9PxAMCGLAh2w238+aDsQABlBzBDHgQ768fIkABkIDQQz4\nUMWsaW+uYkcAA6GFIAZ8yJvXERPAQGgiiAEfqhiajojw/BgEMBDauNc04AOHD2dIqtojPnToYK2O\nwWIMQHigRwz4wJYtaRowYJCk8h6x2+3W9u3/Ups2F5zzucHaAy4qkt54w66tWyM0YoRDOTmWdu2K\n0NChUufO0uuv21VcXL7vhAmBuRAGYAI9YsAHrrjiKu3Ysa1ystbmzZvOucB7sPeAV6+O1OjRDmVl\nWSopsTR0qFOjR5dq8mRp7doI9e/v1IQJDm3fHqFdu/jTA1TgXwPgAxdddLF++umAbDabysrKlJ2d\npfPPP/+0+wZ7AFfo29cpu106cMCmpCSnJCkjw6acnPJtycnlA3Dt2rmUkeH/W34CgYqhacBHunbt\nrn/9a4t++aWBevW6ttrjwToEfSaxsVJaWoS6di2rvFxr/fpIJSZKY8Y4VFpavm337gj9/vcMTQMV\nCGLARy64oK1WrfpQJ0401M03N6/cHmoB/GupqRHq1Kn8BibZ2ZbWrInQp59Kdnv5f1u22HTttWVq\n1iwwloYEAgFBDPjQzTcnqkmTWEmhHcAVNm+OUPfuZUpOjtSOHRFasqREF1zQQFlZUmGhtHlzpCZN\nKjVdJhBQCGLAR3Jz8zR//m7tzz+uvCv+rqIW3+nBq/4rJANYkhwOafdum5YtK5ZlSQMHOqs8npxs\n14QJpXI4pC++iND115cZqhQILAQx4CMPz3hX6yJ2ST3WSGkT1f/QaE0ad58idn2lmDf+JndCglRa\nKtvRXB2b95zUoIHpkutk27byYWnrNPOwUlIi9cQT0Zo7N0oul6UPPyzyf4FAgCKIAS+rGILeePHH\n0qZp2vWP7erinq7DVy1X9NJ3FPM/ryj/vQ/kbtZMkhTz8ouK2rRRpbclGa7cc7t32zR/fpRycy1t\n2BChG2+s2tu96y6n7rqr0FB1QGAjiAEvOfUz4L57btJHn/9enTVdklt9Y7crdso85S3/qDKEI3ds\nkz0tVSXDRpgtvo46dXJp2bJi02UAQYkgBuroTJOwjs7Lk931dylFuuuuv2tGzmdyN0pQ9OpVil61\nQnK5VNahowpefV1q2PDkAQsLFTfxYRXOnitXy1bmGgbALwhiwEPnmgWdkNBIixffLaWM1uJX7lD9\nC8apZMRoHZ/xxBmPWe+dt2Q7nKGoVR9Ks+b4oxkADCOIgVry5DIkKzdXKitTWdsLz7pfyfBRkqT6\n8+d5rV4AgY1bXAI1VJdbUbobN5Y7Nk4qc1Z7zMrMVNTKD31RMoAgUKcgTk9P18iRI71VCxCQvHIv\n6IgIlYwaq6h1a6tsjty5XQ1nz5Dj+hu8XDWAYOHx0PSSJUuUkpKiBkF+7SNwJt6+E9bx6TPVYPZM\nxT74O7latpZKT6js4vY69tIrqrw5M4Cw43EQt23bVgsXLtTUqVO9WQ9gnLcCOD83V59P+5PGSlr5\nwGj1eeZ56fEnvV8wgKDm8dvwxMRERUREeLMWwChvL0f4+bQ/aUxKsiRpTMo/9PnUP3mzXAAhwq+z\npps2jfXn6QJKuLbdZiu/32Egt39v9l7N3jhba/av0aRek/T6oMWKi677vaCbHD6kirs9Wv/5/pw/\nh3fflTZtkixLjZ99UrruOunhh+tci2mB/Pr7Qzi3P5zbXlN1DmK3u+bLmWVlHavr6YJS06axYdt2\nl8stm80KyPafOgSddt8OxUbF6USBlKW615vdsrXc2ipLkltSdss25/45JN5R/t+sp09uC8CfXW2E\n8++/FN7tD+e2SzV/E1LnILZOd4d3IID5aznCPs88r0UOpzK/260LL79CfZ5Z4PVzAAh+dQriVq1a\naenSpd6qBfApf68HHJ9wnhoOHab96VvVf9pM2ZgZDeA0+MuAkOftSVg15XK5VFZWph49emjDhk99\nei4AwYsgRsgyFcAVNm/epN69r1PDhg114kSpiopYgxdAdQQxQo7pAK5w9OhRnXdeY0nSrbf21aef\nrvHr+QEEBxZ9QMjw92fAZ7N797fq2PHSyu/tdrsaN26iI0d+UfPmLYzUBCAwEcQIeoEUwBWio6N1\n0UUXSzp5id811/TW8eOFJssCEIAIYgStQAzgChUhLJ28xM+yLDVsyM0NAFRFECPoBHIAA0BtEcQI\nGgQwgFBEECPgEcAAQhlBjIBFAAMIBwQxAk6oBXBtFkYBEH4IYgSMUAvgCtddd53KykxXASBQEcQw\nLlQDuMJ5550X1kvBATg7ghjGhHoAA0BNEMTwOwIYAE4iiOE3BDAAVMfqS/A5R5wjIFZD8ob8/DyN\nGXOfHA5Hle0//3xIt9zSWy6Xy6PjPv30k1qy5NXK7x0Oh8aOvU95eXl1qhdA4COI4TPfH92nnBuy\ndaT/kaAP4AqvvvqyBg0aKrvdXrktM/OIpk6dVC2ca+qdd97UypUpVbbZ7XYNGTJMixa9WKd6AQQ+\nghhe9+v1gO15drX6oFXQB7BUHrgbNqxTv363V27buHGD7r9/lKKjo2t9vKKi4/rrX6fq3XffUrNm\nzas93rdvkjZt2qjMzCN1qhtAYCOI4TW/DuCKHnBcerxsjtD4NVuxYrm6d+9ZpTf8xRepGj/+Yf3x\nj1NqfbyHDVVFAAAKSElEQVTDhw/L4XDob397Ry1atKz2eGRkpHr06KmUlOQ61Q0gsIXGX0gYdboA\nDoUe8KnS0jarR4/fVNk2bdp03XHHAI+Od8kl7fX008+refPqveEKPXr8Rmlpmz06PoDgwKxpeCyc\nZkG7XC7t27dXbdte6Nfztmt3kfbv3yen06nISP65AqGIf9motXAK4Ar5+flyuVxq1KiRX88bFxcv\nt9ut/Pw8NW7cxK/nBuAfBDFqLBwDuILNZkmSyvx80+iKBSNsNj5FAkIVQYxzCucArhAf30gREREq\nKMj363nz8vJkWZYaNUrw63kB+A9vs3FG4TIJq6bat++g/fv31Xh/p9Op3Nwcj2/yIUn79+9Thw4d\nZVmWx8cAENgIYlRDAJ9er169lZ6+s8b779qVrrvu6qcjR859HfCZgvarr3bqmmuuq/E5AQQfhqZR\niSHos7v99js1ZswwFRcXKyYmpspjV1/dTRs3bqm2LSnpt4qOjta5OsUvvfRqtW3FxcVKS9ushx76\nQ51rBxC46BGDHnANNW/eQtdff5NWr15Zo/0zMn5WXt5RNW7c2KPzffzxKvXpc8Np77oFIHQQxGGM\nAK69hx76g1JSkmt0X+lWrVpr3rwFHp3H4XBo+fIP9PDDf/To+QCCB0PTYYghaM8lJJynN998r8b7\ne3rZkd1ur9V5AAQvgjiMEMAAEHgI4jBAAANA4CKIQxgBDACBjyAOQQQwAAQPgjiEEMAAEHwI4hBA\nAANA8CKIgxgBDADBjyAOQgQwAIQOgjiIEMAAEHoI4iBAAANA6CKIAxgBDAChjyAOQAQwAIQPgjiA\nEMAAEH4I4gBAAANA+CKIDSKAAQAEsQEEMACgAkHsRwQwAOBUBLEf7M3eq+mfzCCAAQDVEMQ+5na7\nNW7FON3Q8hYCGABQDUHsY5ZladPvNikr65jpUgAAAchmugAAAMIZQQwAgEEEMQAABhHEAAAYRBAD\nAGCQR7Om3W63Hn/8ce3du1dRUVGaM2eO2rRp4+3aAAAIeR71iNeuXavS0lItXbpUU6ZM0dy5c71d\nFwAAYcGjIN62bZv69OkjSbryyiv19ddfe7UoAADChUdD04WFhYqNjT15kMhIuVwu2WxnzvV27drJ\n5XJ7crqgZ7NZYdv2w4czJEndunU2XIk54fz6S7Q/nNsfzm2XpIMHf6rRfh4FccOGDXX8+PHK788V\nwhVsNsuT04WEcG67RPtpP+0PV+Hc9pryKIi7du2q9evXq1+/ftq5c6c6dOhwzuccOHAgbG/z2LRp\nbNi2vVu3zrLZLG3dust0KcaE8+sv0f5wbn84t702PArixMREpaam6t5775UkJmsBAOAhj4LYsizN\nmjXL27UAABB2uKEHAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAG\nEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACA\nQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMA\nYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQA\nABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQxAAAGEcQAABhEEAMAYBBBDACAQQQx\nAAAGEcQAABhEEAMAYBBBDACAQXUK4k8++URTpkzxVi0AAISdSE+fOGfOHKWmpqpTp07erAcAgLDi\ncY+4a9euevzxx71YCgAA4eecPeJly5bpzTffrLJt7ty5SkpK0pYtW3xWGAAA4eCcQTx48GANHjzY\nKydr2jTWK8cJRuHadpvNkhS+7a9A+2l/uArntteUx58ReyIr65g/TxcwmjaNDdu2u1xu2WxW2LZf\nCu/XX6L94dz+cG67VPM3IVy+BACAQXXqEffs2VM9e/b0Vi0AAIQdesQAABhEEAMAYBBBDACAQQQx\nAAAGEcQAABhEEAMAYJDldrvdposAACBc0SMGAMAgghgAAIMIYgAADCKIAQAwiCAGAMAgghgAAIP8\nEsSFhYV68MEHNXLkSN17773auXOnP04bcD755BNNmTLFdBl+43a7NXPmTN17770aNWqUDh06ZLok\nv0tPT9fIkSNNl+F3TqdTU6dO1fDhwzV06FCtW7fOdEl+5XK59Nhjj2nYsGEaPny4vv/+e9MlGZGT\nk6Mbb7xRP/74o+lS/G7gwIEaNWqURo0apccee+ys+9ZpGcSa+t///V9de+21GjVqlH788UdNmTJF\nycnJ/jh1wJgzZ45SU1PVqVMn06X4zdq1a1VaWqqlS5cqPT1dc+fO1aJFi0yX5TdLlixRSkqKGjRo\nYLoUv/vwww+VkJCgZ555Rvn5+RowYIBuvvlm02X5zbp162RZlt577z1t2bJFCxYsCKvffan8zdjM\nmTNVr14906X4XWlpqSTprbfeqtH+fukRjx07Vvfee6+k8hcnOjraH6cNKF27dtXjjz9uugy/2rZt\nm/r06SNJuvLKK/X1118brsi/2rZtq4ULF5ouw4ikpCRNnDhRUnnvMDLSL+/5A8att96q2bNnS5Iy\nMjIUHx9vuCL/e/rppzVs2DCdf/75pkvxuz179qioqEjjxo3TmDFjlJ6eftb9vf6vY9myZXrzzTer\nbJs7d646d+6srKwsTZ06VdOnT/f2aQPGmdqflJSkLVu2GKrKjMLCQsXGxlZ+HxkZKZfLJZstPKYm\nJCYmKiMjw3QZRsTExEgq/x2YOHGiJk+ebLgi/7PZbHrkkUe0du1avfTSS6bL8avk5GQ1btxYvXv3\n1quvvmq6HL+rV6+exo0bpyFDhujAgQN64IEH9M9//vOMf/u8HsSDBw/W4MGDq23fu3ev/vznP2va\ntGnq3r27t08bMM7U/nDUsGFDHT9+vPL7cAphSL/88osmTJigESNGqH///qbLMWLevHnKycnRkCFD\n9NFHH4XNMG1ycrIsy1Jqaqr27NmjadOm6ZVXXlHjxo1Nl+YX7dq1U9u2bSu/btSokbKystSsWbPT\n7u+X8aLvv/9ekyZN0gsvvKCOHTv645QIAF27dtX69evVr18/7dy5Ux06dDBdkhHheDv37OxsjRs3\nTjNmzFCvXr1Ml+N3KSkpyszM1Pjx4xUdHS2bzRZWb0Lffvvtyq9HjhypJ554ImxCWJI++OADfffd\nd5o5c6YyMzN1/PhxNW3a9Iz7+yWIFyxYoNLSUs2ZM0dut1txcXFh+9lZOElMTFRqamrl/IC5c+ca\nrsgMy7JMl+B3r732mgoKCrRo0SItXLhQlmVpyZIlioqKMl2aX/Tt21ePPvqoRowYIafTqenTp4dN\n208Vjr//gwcP1qOPPqr77rtPNptNTz311FnfiLH6EgAABoXPWAkAAAGIIAYAwCCCGAAAgwhiAAAM\nIogBADCIIAYAwCCCGAAAgwhiAAAM+v89Qjx0rPOhJgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1173d5d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import linear_model\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "x = np.array([1, 2, 3]).reshape((-1,1))\n",
    "y = np.array([1, 2, 2]).reshape((-1,1))\n",
    "predict_line = np.array([-1, 4]).reshape((-1,1))\n",
    "\n",
    "regr = linear_model.LinearRegression()\n",
    "regr.fit(x, y)\n",
    "ey = regr.predict(x)\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.axis('equal')\n",
    "plt.axhline(y=0, c='black')\n",
    "plt.axvline(x=0, c='black')\n",
    "\n",
    "plt.scatter(x, y, c='r')\n",
    "plt.scatter(x, regr.predict(x), s=20, c='b')\n",
    "plt.plot(predict_line, regr.predict(predict_line), c='g', lw='1')\n",
    "[ plt.plot([x[i], x[i]], [y[i], ey[i]], 'r', lw='1') for i in range(len(x))]\n",
    "\n",
    "plt.annotate('(1, 1)', xy=(1, 1), xytext=(-15, -30), textcoords='offset points', size=14, arrowprops=dict(arrowstyle=\"->\"))\n",
    "plt.annotate('(2, 2)', xy=(2, 2), xytext=(-60, -5), textcoords='offset points', size=14, arrowprops=dict(arrowstyle=\"->\"))\n",
    "plt.annotate('(3, 2)', xy=(3, 2), xytext=(-15, -30), textcoords='offset points', size=14, arrowprops=dict(arrowstyle=\"->\"))\n",
    "\n",
    "plt.annotate('$e_1$', color='r', xy=(1, 1), xytext=(0, 2), textcoords='offset points', size=20)\n",
    "plt.annotate('$e_2$', color='r', xy=(2, 2), xytext=(0, -15), textcoords='offset points', size=20)\n",
    "plt.annotate('$e_3$', color='r', xy=(3, 2), xytext=(0, 1), textcoords='offset points', size=20)\n",
    "\n",
    "plt.annotate('$p_1$', xy=(1, 7/6), color='b', xytext=(-7, 30), textcoords='offset points', size=14, arrowprops=dict(arrowstyle=\"->\"))\n",
    "plt.annotate('$p_2$', xy=(2, 5/3), color='b', xytext=(-7, -30), textcoords='offset points', size=14, arrowprops=dict(arrowstyle=\"->\"))\n",
    "plt.annotate('$p_3$', xy=(3, 13/6), color='b', xytext=(-7, 30), textcoords='offset points', size=14, arrowprops=dict(arrowstyle=\"->\"))\n",
    "plt.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据条件可以得到方程组 \n",
    "$\n",
    "\\begin{cases}\n",
    "C+D&=1 \\\\\n",
    "C+2D&=2 \\\\\n",
    "C+3D&=2 \\\\\n",
    "\\end{cases}\n",
    "$，写作矩阵形式\n",
    "$\\begin{bmatrix}1&1 \\\\1&2 \\\\1&3\\\\\\end{bmatrix}\\begin{bmatrix}C\\\\D\\\\\\end{bmatrix}=\\begin{bmatrix}1\\\\2\\\\2\\\\\\end{bmatrix}$，也就是我们的$Ax=b$，很明显方程组无解。\n",
    "\n",
    "我们需要在$b$的三个分量上都增加某个误差$e$，使得三点能够共线，同时使得$e_1^2+e_2^2+e_3^2$最小，找到拥有最小平方和的解（即最小二乘），即$\\left\\|Ax-b\\right\\|^2=\\left\\|e\\right\\|^2$最小。此时向量$b$变为向量$p=\\begin{bmatrix}p_1\\\\p_2\\\\p_3\\end{bmatrix}\\\\$（在方程组有解的情况下，$Ax-b=0$，即$b$在$A$的列空间中，误差$e$为零。）我们现在做的运算也称作线性回归（linear regression），使用误差的平方和作为测量总误差的标准。\n",
    "\n",
    "注：如果有另一个点，如$(0, 100)$，在本例中该点明显距离别的点很远，最小二乘将很容易被离群的点影响，通常使用最小二乘时会去掉明显离群的点。\n",
    "\n",
    "现在我们尝试解出$\\hat x=\\begin{bmatrix}\\hat C\\\\ \\hat D\\end{bmatrix}$与$p=\\begin{bmatrix}p_1\\\\p_2\\\\p_3\\end{bmatrix}$。\n",
    "\n",
    "$$\n",
    "A^TA\\hat x=A^Tb\\\\\n",
    "A^TA=\n",
    "\\begin{bmatrix}3&6\\\\6&14\\end{bmatrix}\\qquad\n",
    "A^Tb=\n",
    "\\begin{bmatrix}5\\\\11\\end{bmatrix}\\\\\n",
    "\\begin{bmatrix}3&6\\\\6&14\\end{bmatrix}\n",
    "\\begin{bmatrix}\\hat C\\\\\\hat D\\end{bmatrix}=\n",
    "\\begin{bmatrix}5\\\\11\\end{bmatrix}\\\\\n",
    "$$\n",
    "\n",
    "写作方程形式为$\\begin{cases}3\\hat C+16\\hat D&=5\\\\6\\hat C+14\\hat D&=11\\\\\\end{cases}$，也称作正规方程组（normal equations）。\n",
    "\n",
    "回顾前面提到的“使得误差最小”的条件，$e_1^2+e_2^2+e_3^2=(C+D-1)^2+(C+2D-2)^2+(C+3D-2)^2$，使该式取最小值，如果使用微积分方法，则需要对该式的两个变量$C, D$分别求偏导数，再令求得的偏导式为零即可，正是我们刚才求得的正规方程组。（正规方程组中的第一个方程是对$C$求偏导的结果，第二个方程式对$D$求偏导的结果，无论使用哪一种方法都会得到这个方程组。）\n",
    "\n",
    "解方程得$\\hat C=\\frac{2}{3}, \\hat D=\\frac{1}{2}$，则“最佳直线”为$y=\\frac{2}{3}+\\frac{1}{2}t$，带回原方程组解得$p_1=\\frac{7}{6}, p_2=\\frac{5}{3}, p_3=\\frac{13}{6}$，即$e_1=-\\frac{1}{6}, e_2=\\frac{1}{3}, e_3=-\\frac{1}{6}$\n",
    "\n",
    "于是我们得到$p=\\begin{bmatrix}\\frac{7}{6}\\\\\\frac{5}{3}\\\\\\frac{13}{6}\\end{bmatrix}, e=\\begin{bmatrix}-\\frac{1}{6}\\\\\\frac{1}{3}\\\\-\\frac{1}{6}\\end{bmatrix}$，易看出$b=p+e$，同时我们发现$p\\cdot e=0$即$p\\bot e$。\n",
    "\n",
    "误差向量$e$不仅垂直于投影向量$p$，它同时垂直于列空间，如 $\\begin{bmatrix}1\\\\1\\\\1\\end{bmatrix}, \\begin{bmatrix}1\\\\2\\\\3\\end{bmatrix}$。\n",
    "\n",
    "接下来我们观察$A^TA$，如果$A$的各列线性无关，求证$A^TA$是可逆矩阵。\n",
    "\n",
    "先假设$A^TAx=0$，两边同时乘以$x^T$有$x^TA^TAx=0$，即$(Ax)^T(Ax)=0$。一个矩阵乘其转置结果为零，则这个矩阵也必须为零（$(Ax)^T(Ax)$相当于$Ax$长度的平方）。则$Ax=0$，结合题设中的“$A$的各列线性无关”，可知$x=0$，也就是$A^TA$的零空间中有且只有零向量，得证。\n",
    "\n",
    "我们再来看一种线性无关的特殊情况：互相垂直的单位向量一定是线性无关的。\n",
    "\n",
    "* 比如$\\begin{bmatrix}1\\\\0\\\\0\\end{bmatrix}\\begin{bmatrix}0\\\\1\\\\0\\end{bmatrix}\\begin{bmatrix}0\\\\0\\\\1\\end{bmatrix}$，这三个正交单位向量也称作标准正交向量组（orthonormal vectors）。\n",
    "* 另一个例子$\\begin{bmatrix}\\cos\\theta\\\\\\sin\\theta\\end{bmatrix}\\begin{bmatrix}-\\sin\\theta\\\\\\cos\\theta\\end{bmatrix}$\n",
    "\n",
    "下一讲研究标准正交向量组。"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}