Table of Contents
\n", "- 1 Homework
- 1.1 Homework - Rank of the product of two matrices (and the outer-product of two vectors)
- 1.2 Homework - A simple mixed determined problem [solved by hand and using the pseudo-inverse]
- 1.3 Homework - A \"least squares\" and minimum norm solution to a simple $2\\times 2$ case
- 1.4 Homework - Lagrange multiplier derivation of minimal-norm solution
- 1.5 Homework - Differentiation of inner products
- 1.6 Homework - Curve-fitting - response to outliers [read-through]
- 1.7 Homework - Tomography example - rank and null space
- 1.8 Homework - Tomography example - solution via the pseudo-inverse
- 1.9 Homework - Damped least squares applied to the tomography example
\n", "\n", "The answer (if you work it out component by component) turns out to be equivalent to\n", "\n", "$$\\frac{\\partial}{\\partial \\boldsymbol{x}} \\left(\\boldsymbol{a}^T\\boldsymbol{b}\\right) \n", "=\\left(\\frac{\\partial \\boldsymbol{a}}{\\partial \\boldsymbol{x}}\\right)^T\\boldsymbol{b} +\n", "\\left(\\frac{\\partial \\boldsymbol{b}}{\\partial \\boldsymbol{x}}\\right)^T\\boldsymbol{a}$$\n", "\n", "The differentials $\\partial \\boldsymbol{a}/\\partial \\boldsymbol{x}$ and $\\partial \\boldsymbol{a}/\\partial \\boldsymbol{x}$ are both $m\\times n$ matrices, so that their\n", "transposes are $n\\times m$. \n", "\n", "Thus the products \n", "$(\\partial \\boldsymbol{a}^T/\\partial \\boldsymbol{x}) \\boldsymbol{b}$ and $(\\partial \\boldsymbol{a}^T/\\partial \\boldsymbol{x}) \\boldsymbol{a}$\n", "are both column vectors of length $n$ as \n", "required. Note that it does not matter if we differentiate a vector and then transpose the \n", "result, or if we transpose the vector before differentiation - both generate the same outcome. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework - Curve-fitting - response to outliers [read-through]\n", "\n", "Here we are going to fit a *linear* line to some invented data, and see what happens if we create an outlier - how much is the slope of the best-fit line impacted?\n", "\n", "As numpy's polyfit function only has the option to minimise the 2 norm, we have to do some work ourselves to create an approach that minimses other norms - so just read through the following solution." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGKCAYAAAARyQg4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmPklEQVR4nO3de5xddXnv8c8TSMRJUDlAok3MjEdBwBqrREARmaBW4qV4O1XcSrFyRltB21fxUqcHUTpeqiJ6KtIRkUMZjVXBUsQblwEUUUKLIFcDZIYEhOI9DBgSnvPH3omTYU+yJpnstfeaz/v12q/JWuu393rmmcl897rstSIzkSRJ1TGr7AIkSdL0MtwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdAiLiTyPiWxHxi4h4KCJui4iPRcQeZde2SUScHRGrS1hvb0RkRPS2et1liIiexvd77HY8928i4jU7oSxpSgx3zXgR8X7gO8BDwHHAS4EzgGOBayLiyeVVt4VTgFeXXYS26m8Aw12l27XsAqQyRcQy4B+B0zLzb8ctujwizgeuBc4Blk3T+h6Tmb/fnudm5u3TUUM72pG+SHo0t9w1070H+CXw9xMXZOadwEeB3og4GCbfZdts13VEDEfE9yPilRHxXxHxe+CvJyskIt7YGLcuIn4TETdExNvGLd9it/y4Wt4eER+JiJ9HxO8i4tyI6IqIp0XEdxqvtyoi/mLC+k5uPP+ZEXFZRIxFxD0R8aGI2Obfhoh4TURc3XjeryPiqxGxuMDzJu1LRBwfET+MiF82XvPqiHj5hOf/NCLOHDf9+IjYGBFrJoz7QUT82zZq6YqI0xuHY9ZFxAXAoibjnhsRX4uINRHxYETcGhEfjojHjhuzGugGao2+ZkSc3Vj2tIj414i4s/H8OyLic+102EfVYrhrxoqIXYHDge9l5kOTDLug8fWI7VzNvsBngP9LfXf/JZPU8gLgXOBy4FXA/wI+DzyhwDr+Hvgj4C+Ak4DXUz+scD7wTeq78q8HvhgRz2jy/G8AFzfW+yXg/zReZ1IR8Xbg68BNwOuAtwF/TH2Px+4Fap6sLz3AmdS//9cDK4ELI2L5uOdeypY/j17g98DCiNi3Ud9c4LnAZduo41+oH4o5lfru9Fup92CixcB1wNuBI4FPA38JfHHcmFcDP6d+iOd5jccpjWV/BKyhvtv+pcCHgBcBF22jPmn7ZKYPHzPyASwAEvjIVsbs1hhzemO6pzF97IRxvY35vePmDQOPAH9SoJYTgV9uY8zZwOpx05tquXTCuPMa8980bt4ewAbgA+PmndwY974Jz/888DvgCc2+N2Ae8BvgrAnP6wHWA3+zje+jUF+ob3zsCnwX+Pdx81/dqKe7MX0a9TdhPwPe1ph3ZGPMflt5/acDG5t8/59r9jMetzwadb2p8X3sOW7ZauDcAj/vXYEXNNbz7LL/L/io3sMtd81k0YJ1rM7M6wqMuwbYo7FL/RUR8YQprONbE6ZvaXz9zqYZmfkr4D6g2cmBE3ddr6Ae4H88yfqeBzwOGIqIXTc9qG+Z3gK8sEDNTfsSEQdGxIURcS/1NyMPAy+hHsSbXE49VDdtvR9BfWv+0gnz7snMW5jcwdTfQDT7/ifW9bjGpydup76X4GHgX6n/Du2zlXVsev6ciHh/RNwSEQ82nn9lY/HTt/JUabsY7prJ7gcepL7FOZlNy+7aznXcU2RQZl5OfVf0k6nvTv/viLg4IpYUePqvJkyv38r83Zo8/95JphdOsr75ja8XUw+p8Y9nAntuo15o0pfGpxIuAf4HcALwfOq71r89vu7M/CXwE2BZROxF/U3IZY1Hb2PYMra9S/5Jja+Tff/jfZH6LvnPUH+z8VzgHY1lzXo60Ueo7yk5F3g5cBB/OKu+yPOlKfFsec1YmbkhIq4AXhIRu2Xz4+5/1vh6aePrpjFzJoybLNAK31M5M78GfC0i5lEPqY8B346IRZn5SNHX2Q4LgDsmTAOsnWT8LxpfjwVubLL8dwXW2awvRwKPB/48MzefHBcRXU3GXkb9mPyyRj3XU3/DMD8iDgWeTf14+tZseoMx2fe/af27AUcBJ2fmp8fNf+Y2Xn+8NwDnZOY/jnv+vCk8X5oSt9w1032cejB/eOKCiHgK8F7gisz8UWP2vdR3y07cZf1ypklmrsvMC6mH05MotiW8I/58wvQbgHXATycZfxX1AH9aZq5s8rh1O+vYFOIPb5rROEHu0CZjL6O+Z+FtwHDW3Uf9zcYHgV34wxuyyfyI+u79Zt//eI9pvN7DE+Yf2+Q1fw88tsn8ribPf8s26pO2m1vumtEy85KIOAn4UET0UP9M+6+A5wDvo37i2JvHjc+I+Arw1oi4jfrZ1S/nD7uDt0tEfIj6FuNlwN3UP471TuC6zPzvHXntAv5346Nv11A/k/s46lupv242ODN/GxHvBj4bEXtTP+b/G+phezj1sG12xvm2XEz9OPs5EfFJ6m9sPgiM8ugNkSuonwz3Iv6wexzq/TseGM3MO9iKzLw1Ir5E/We/6ft/CfCyCeN+ExFXA38XEfdQP5zzlzQ/bHETcFhEvIL6mfP3Z+Zq6ocW/iIibgBWUd8l//yt1SftCLfcNeNl5inAcmAu9WOr36X+uetzgKWZOTrhKe+ifkb6ycBXqB8zPWEHy/gR9eP7nwK+R32X/OVM4x6BrTiKeqhdQP0M8H/kDx/haioz/4X6IYunUz+x7FvUg3hX6h8Zm7LMvBGoUf+s+AXUr0HwPupBPnHsb6lfYAi23ELf9O9tHW/f5G3AF6h/WuF8YD/gjU3GHd1Y32epf2rh59R/Dyb6e+pv+P6N+puFkxvzT6D+PQ1Q/53ZvfGa0k4RmYUPCUqqkIg4GfgAMDszN5RcjqRp5Ja7JEkVY7hLklQx7paXJKli3HKXJKliDHdJkiqmMp9z32uvvbKnp6fsMrbwwAMPMHfu3LLL6Aj2qhj7VJy9KsY+FdOufbr22mvvz8y9J86vTLj39PSwcuXKssvYwvDwML29vWWX0RHsVTH2qTh7VYx9KqZd+xQRI83mu1tekqSKMdwlSaqYlod7RJwVEfdFRNObUkRELSKubzyuiohntbpGSZI6WRlb7mdTv7XjZO4EDs/MJdSvbz3YiqIkSaqKlp9Ql5lXNO6+Ndnyq8ZNXk397liSJKmgUq5Q1wj3CzNz4j2xJ447EdgvM4+bZHkf0AewYMGCA1esWDHpa82ZM4fddtuNiNjuuqcqM1u6vk42Wa8yk4ceeoj169eXUFX7WbduHfPmzSu7jI5gr4qxT8W0a5+WLVt2bWYunTi/bT8KFxHLgLcCL5hsTGYO0thtv3Tp0pzsYwqjo6NEBAsWLGD27NktC9zf/e537L777i1ZV6dr1qvM5OGHH+bee+8lM1m8eHFJ1bWPdv04TjuyV8XYp2I6rU9tebZ8RCwBzgSOysxf7OjrPfDAAyxcuJA5c+a4Jd1BIoI5c+awcOFCHnjggbLLkaSO0XbhHhGLgfOAN2fmbdP1urNmtd23qoL82UnS1LR8t3xEfBnoBfaKiDXAB4DZAJl5BnASsCdwemMre0Oz4wmSJKm5Ms6WP3oby48Dmp5AJ0mSts39nSrka1/7mucrSFKHMNy10/T09PCJT3yi7DIkacYx3KdgaGiInp4eZs2aRU9PD0NDQy1dv5/1liQVYbgXNDQ0RF9fHyMjI2QmIyMj9PX17dSA7+3t5a/+6q848cQT2XvvvTn00EMBOPXUU1myZAlz585l4cKFHHfccfz617/e/LwnPvGJfOUrX9k8feihh7L77ruzYcMGAH72s58REaxdu3bSdZ9zzjl0d3fT1dXFK17xCu69994tlt9+++0cddRRPPGJT2Tu3Lk85znP4cILL9yi9pGREd797ncTEZt36f/iF7/g6KOPZtGiRTz2sY/lGc94Bl/84hd3uFeSpD8w3Avq7+9nbGxsi3ljY2P09/fv1PWee+65ZCZXXnkl55xzDlD/aNhpp53GjTfeyJe+9CV+/OMfc8IJJ2x+zuGHH85ll122ucaVK1fymMc8ZvP97oeHh3na057GwoULm67zRz/6Ecceeyx9fX1cd911vPKVr+Skk07aYsy6detYvnw53/ve9/jJT37Ca1/7Wl7zmtdwyy23AHDeeeexaNEiTjrpJO655x7uueceAB566KHNbwRuvPFG3vWud/G2t72N4eHhae2bJM1omVmJx4EHHpiTuemmmyZdVlREJPCoR0RM+pzf/va3O7TOww8/PJ/5zGduc9y3vvWtnDNnTm7cuDEzM08//fTcd999MzPzu9/9bu633355zDHH5Ic//OHMzHzjG9+Yxx133KSvd/TRR+eLX/ziLea99a1vzfqvy+QOPvjgPOWUUzZPd3d358c//vFt1v/6178+jznmmK2OmY6fYRVcdtllZZfQMexVMfapmHbtE7Aym2SiW+4FTXbp0519SdQDDzzwUfMuvfRSXvKSl7Bo0SJ23313XvOa17B+/Xp+/vOfA/Vd4rfddht33303w8PDLFu2jN7e3s1bx5dffvlWL6N4880387znPW+LeROnH3jgAd7znvdwwAEHsMceezBv3jxWrlzJ6OjoVr+fjRs3MjAwwJIlS9hzzz2ZN28e5513HnfddVeBbkiSijDcCxoYGKCrq2uLeV1dXQwMDOzU9c6dO3eL6ZGREV7+8pez//7789WvfpVrr72Ws846C/jDCXf7778/CxYsYHh4eHO4L1u2jB/84AfcdNNNrF27dqvhngVuJnTiiSfy1a9+lVNOOYXLL7+c6667joMOOmibJ/194hOf4JOf/CTvfve7ueSSS7juuut41atexcMPP7zNdUqSimnbG8e0m1qtBtSPvY+OjrJ48WIGBgY2z2+VlStXsn79ej71qU+xyy67AGxxItsmhx9+ON/85jdZuXIlhx9+OPPnz2evvfbin/7pn7Z6vB3ggAMO4Oqrr95i3sTp73//+xxzzDG89rWvBerH0m+//Xb23XffzWPmzJnDxo0bH/W8V77ylbz5zW8G6m8kbrvtNm+wI0nTyC33KajVaqxevZpHHnmE1atXtzzYAfbZZx8eeeQRTjvtNO68806+/OUvc9pppz1qXG9vL1/5ylfYZ599mD9/PlAP/HPPPXebdzZ65zvfycUXX8xHPvIRfvazn/H5z3+e888/f4sx++67L+effz7/+Z//yQ033MCb3vQmHnrooS3G9PT0cOWVV7J27Vruv//+zc+75JJL+P73v88tt9zC8ccfz5133rn9DZEkPYrh3mGWLFnCpz/9aU499VQOOOAAzjzzzKYXilm2bBkbN27cIsibzWvmkEMO4Qtf+AKf+9znWLJkCeeddx4nn3zyFmNOPfVU5s+fz2GHHcby5cs55JBDOOyww7YY86EPfYi77rqLpz71qey9994A/MM//AMHHXQQy5cv54UvfCFz584t5U2SJFVZFDm+2gmWLl2amz7qNdHNN9/M/vvv3+KKvJ/7VGyrV2X9DNtNp91Tukz2qhj7tHVDQ0OlH47dmoi4NpvcXM1j7pIkNbHp4mWbrnGy6eJlQFsFfDPulpckqYmyLl42HQx3SZKamOy6Hdu6nkc7MNwlSWqirIuXTYcZE+5VOXFwJvJnJ6kMZV28bDrMiHCfPXs2Dz74YNllaDs9+OCDzJ49u+wyJM0wtVqNwcFBuru7iQi6u7sZHBxs+5PpYIaE+/z581m7di1jY2NuBXaQzGRsbIy1a9duvhCPJLXSpouXXXrppaVdvGx7zIiPwj3ucY8D4O67727pNcwfeughdtttt5atr5NN1qvZs2ezYMGCzT9DSdK2zYhwh3rAtzoghoeHefazn93SdXYqeyVJ02dG7JaXJGkmMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdkqSKMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdkqSKMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdkqSKMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdklSKoaEhenp6mDVrFj09PQwNDZVdUmXsWnYBkqSZZ2hoiL6+PsbGxgAYGRmhr68PgFqtVmZpleCWuySp5fr7+zcH+yZjY2P09/eXVFG1tDzcI+KsiLgvIn46yfKIiM9ExKqIuD4intPqGiVJO9fo6OiU5mtqythyPxs4civLlwP7NB59wOdaUJMkqYUWL148pfmampaHe2ZeAfxyK0OOAs7JuquBJ0TEk1pTnSSpFQYGBujq6tpiXldXFwMDAyVVVC2Rma1faUQPcGFm/nGTZRcCH83M7zemLwHem5krm4zto751z4IFCw5csWLFTq17qtatW8e8efPKLqMj2Kti7FNx9qqYMvt08cUXc+aZZ3Lfffcxf/58jjvuOF784heXUsu2tOvv07Jly67NzKUT57djuH8T+MiEcH9PZl67tddcunRprlz5qPwv1fDwML29vWWX0RHsVTH2qTh7VYx9KqZd+xQRTcO9Hc+WXwM8edz0IuDukmqRJKnjtGO4XwAc0zhr/hDgN5l5T9lFSZLUKVp+EZuI+DLQC+wVEWuADwCzATLzDOAi4GXAKmAMeEura5QkqZO1PNwz8+htLE/gHS0qR5KkymnH3fKSJGkHGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JFTQ0NERPTw+zZs2ip6eHoaGhsktSC7X8CnWSpJ1raGiIvr4+xsbGABgZGaGvrw+AWq1WZmlqEbfcJali+vv7Nwf7JmNjY/T395dUkVrNcJekihkdHZ3SfFWP4S5JFbN48eIpzVf1GO6SVDEDAwN0dXVtMa+rq4uBgYGSKlKrGe6SVDG1Wo3BwUG6u7uJCLq7uxkcHPRkuhnEs+UlqYJqtZphPoO55S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JUsWUEu4RcWRE3BoRqyLifU2WPz4i/iMifhIRN0bEW8qoU5KkTtTycI+IXYDPAsuBA4CjI+KACcPeAdyUmc8CeoFPRsSclhYqSVKHKmPL/SBgVWbekZnrgRXAURPGJLB7RAQwD/glsKG1ZUqS1JnKCPeFwF3jptc05o33z8D+wN3ADcC7MvOR1pQnSVJn27WEdUaTeTlh+qXAdcARwFOB70XElZn52y1eKKIP6ANYsGABw8PD017sjli3bl3b1dSu7FUx9qk4e1WMfSqm0/pURrivAZ48bnoR9S308d4CfDQzE1gVEXcC+wE/Hj8oMweBQYClS5dmb2/vzqp5uwwPD9NuNbUre1WMfSrOXhVjn4rptD6VsVv+GmCfiHhK4yS5NwAXTBgzCrwIICIWAE8H7mhplZIkdaiWb7ln5oaIOB74DrALcFZm3hgRb28sPwM4BTg7Im6gvhv/vZl5f6trlSSpE5WxW57MvAi4aMK8M8b9+27gT1tdlyRJVeAV6iRJqhjDXZKkijHcJUmqGMNdkqSKMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdkqSKMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdkqSKMdwltY2hoSF6enqYNWsWPT09DA0NlV2S1JF2LbsASYJ6sPf19TE2NgbAyMgIfX19ANRqtTJLkzqOW+6S2kJ/f//mYN9kbGyM/v7+kiqSOpfhLqktjI6OTmm+pMkZ7pLawuLFi6c0X9LkDHdJbWFgYICurq4t5nV1dTEwMFBSRVLnMtwltYVarcbg4CDd3d1EBN3d3QwODnoynbQdPFteUtuo1WqGuTQN3HKXJKliDHdJkirGcJckqWIMd0maAi+Rq07gCXWSVJCXyFWncMtdkgryErnqFIa7JBXkJXLVKQx3SSrIS+SqUxjuklSQl8hVp9hmuEfEG1pRiCS1Oy+Rq05R5Gz5cyKiDzg+M2/a2QVJUjvzErnqBEV2yx8IzAb+KyI+ERHzdnJNkiRpB2wz3DPzhsw8DOgD3gTcGhFH7/TKJEnSdil8Ql1m/j/g6cA3gH+NiMsi4hk7qzBJkrR9pnS2fGb+JjPfATwX2Iv6rvpPRsTuO6U6SZI0ZYXCPSJmR8RBEfHOiPgS8HXgGdRPyHsHcEtE/NlOrFOSJBVU5KNwVwG/BX4IfBLYF/gP4PXAImA+sAL4WkS8feeVKkmSiijyUbh1wEeAHwBXZ+YDTcb8XUTcC7wfOGMa65MkSVO0zXDPzD8t+FpXAB/dsXIkSdKOms7Lz/4EOGoaX0+SJG2Habufe2Y+SP1YvCRJKpE3jpEkqWIMd0mSKsZwlySpYkoJ94g4MiJujYhVEfG+Scb0RsR1EXFjRFze6holSepU03ZCXVERsQvwWeAlwBrgmoi4YPztZCPiCcDpwJGZORoR81tdpyRJnaqMLfeDgFWZeUdmrqd+dbuJH6F7I3BeZo4CZOZ9La5RkqSOVUa4LwTuGje9pjFvvH2BPSJiOCKujYhjWladJEkdruW75YFoMi8nTO8KHAi8CHgs8MOIuDozb9vihSL6qN9nngULFjA8PDz91e6AdevWtV1N7cpeFWOfirNXxdinYjqtT2WE+xrgyeOmFwF3Nxlzf+M69g9ExBXAs4Atwj0zB4FBgKVLl2Zvb+/Oqnm7DA8P0241tSt7VYx9Ks5eFWOfium0PpWxW/4aYJ+IeEpEzAHeAFwwYcy/A4dFxK4R0QUcDNzc4jolSepILd9yz8wNEXE88B1gF+CszLxx0+1iM/OMzLw5Ir4NXA88ApyZmT9tda2SJHWiMnbLk5kXARdNmHfGhOmPAx9vZV2SJFWBV6iTJKliDHdJkirGcJckqWIMd0mSKsZwlySpYgx3aYYYGhqip6eHWbNm0dPTw9DQUNklSdpJSvkonKTWGhoaoq+vj7GxMQBGRkbo6+sDoFarlVmapJ3ALXdpBujv798c7JuMjY3R399fUkWSdibDXZoBRkdHpzRfUmcz3KUZYPHixVOaL6mzGe7SDDAwMEBXV9cW87q6uhgYGCipIkk7k+EuzQC1Wo3BwUG6u7uJCLq7uxkcHPRkOqmiPFtemiFqtZphLs0QbrlLklQxhrskSRVjuEuSVDGGuyRJFWO4S5JUMYa7JEkVY7hLklQxhrskSRVjuEs7yPukS2o3XqFO2gHeJ11SO3LLXdoB3iddUjsy3KUd4H3SJbUjw13aAd4nXVI7MtylHeB90iW1I8Nd2gHeJ11SO/JseWkHeZ90Se3GLXdJkirGcJckqWIMd0mSKsZwlySpYgx3SZIqxnCXJKliDHdJkirGcJckqWIMd0mSKsZwlySpYgx3SZIqxnCXJKliDHdJkirGcJckqWIMd0mSKsZwlySpYgx3SZIqxnCXJKliDHdJkiqmlHCPiCMj4taIWBUR79vKuOdGxMaIeF0r65MkqZO1PNwjYhfgs8By4ADg6Ig4YJJxHwO+09oKJUnqbGVsuR8ErMrMOzJzPbACOKrJuBOArwP3tbI4SZI63a4lrHMhcNe46TXAweMHRMRC4NXAEcBzJ3uhiOgD+gAWLFjA8PDwdNe6Q9atW9d2NbUre1WMfSrOXhVjn4rptD6VEe7RZF5OmD4NeG9mboxoNrzxpMxBYBBg6dKl2dvbO00lTo/h4WHaraZ2Za+KsU/F2ati7FMxndanMsJ9DfDkcdOLgLsnjFkKrGgE+17AyyJiQ2Z+oyUVSpLUwcoI92uAfSLiKcBa4A3AG8cPyMynbPp3RJwNXGiwS5JUTMvDPTM3RMTx1M+C3wU4KzNvjIi3N5af0eqaJEmqkjK23MnMi4CLJsxrGuqZeWwrapIkqSq8Qp0kSRVjuEuSVDGGuyRJFWO4S5JUMYa7JEkVY7hLklQxhrskSRVjuEuSVDGGuyRJFWO4S5JUMYa7JEkVY7hLklQxhrskSRVjuEuSVDGGuyRJFWO4S5JUMYa7JEkVY7hLklQxhrskSRVjuEuSVDGGuyRJFWO4S5JUMYa7JEkVY7hLklQxhrskSRVjuM9gQ0ND9PT0MGvWLHp6ehgaGiq7JEnSNNi17AJUjqGhIfr6+hgbGwNgZGSEvr4+AGq1WpmlSZJ2kFvuM1R/f//mYN9kbGyM/v7+kiqSJE0Xw32GGh0dndJ8SVLnMNxnqMWLF09pviSpcxjuM9TAwABdXV1bzOvq6mJgYKCkiprzpD9JmjrDfYaq1WoMDg7S3d1NRNDd3c3g4GBbnUy36aS/kZERMnPzSX8GvCRtneE+g9VqNVavXs0jjzzC6tWr2yrYwZP+JGl7Ge5qW570J0nbx3BX2/KkP0naPoa72lannPQnSe3GcFfb6oST/iSpHXn5WbW1Wq1mmEvSFLnlLklSxRjukiRVjOEuSVLFGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFVMKeEeEUdGxK0RsSoi3tdkeS0irm88roqIZ5VRpyRJnajl4R4RuwCfBZYDBwBHR8QBE4bdCRyemUuAU4DB1lYpSVLnKmPL/SBgVWbekZnrgRXAUeMHZOZVmfmrxuTVwKIW1yhJUseKzGztCiNeBxyZmcc1pt8MHJyZx08y/kRgv03jJyzrA/oAFixYcOCKFSt2XuHbYd26dcybN6/sMjqCvSrGPhVnr4qxT8W0a5+WLVt2bWYunTi/jFu+RpN5Td9hRMQy4K3AC5otz8xBGrvsly5dmr29vdNU4vQYHh6m3WpqV/aqGPtUnL0qxj4V02l9KiPc1wBPHje9CLh74qCIWAKcCSzPzF+0qDZJkjpeGcfcrwH2iYinRMQc4A3ABeMHRMRi4DzgzZl5Wwk1SpLUsVq+5Z6ZGyLieOA7wC7AWZl5Y0S8vbH8DOAkYE/g9IgA2NDsmIIkSXq0Uj7nnpkXZea+mfnUzBxozDujEexk5nGZuUdm/knj0VHBPjQ0RE9PD0cccQQ9PT0MDQ2VXZIkaQYp45h7pQ0NDdHX18fY2BgAIyMj9PX1AVCr1cosTZI0Q3j52WnW39+/Odg3GRsbo7+/v6SKJEkzjeE+zUZHR6c0X5Kk6Wa4T7PFixdPab4kSdPNcJ9mAwMDdHV1bTGvq6uLgYGBkiqSJM00hvs0q9VqDA4O0t3dTUTQ3d3N4OCgJ9NJklrGs+V3glqtRq1W67jLFUqSqsEtd0mSKsZwlySpYgx3SZIqxnCXJKliDHdJkirGcJckqWIMd0mSKsZwlySpYgx3SZIqxnCXJKliDHdJkirGcJckqWIMd0mSKsZwlySpYgx3SZIqxnCXJKliDHdJkirGcJckqWIMd0mSKsZwlySpYgx3SZIqxnCXJKliIjPLrmFaRMR/AyNl1zHBXsD9ZRfRIexVMfapOHtVjH0qpl371J2Ze0+cWZlwb0cRsTIzl5ZdRyewV8XYp+LsVTH2qZhO65O75SVJqhjDXZKkijHcd67BsgvoIPaqGPtUnL0qxj4V01F98pi7JEkV45a7JEkVY7hPg4g4MiJujYhVEfG+JssjIj7TWH59RDynjDrLVqBPtUZ/ro+IqyLiWWXU2Q621atx454bERsj4nWtrK9dFOlTRPRGxHURcWNEXN7qGttFgf9/j4+I/4iInzR69ZYy6ixbRJwVEfdFxE8nWd4Zf88z08cOPIBdgNuB/wnMAX4CHDBhzMuAbwEBHAL8qOy627RPzwf2aPx7+UzsU9FejRt3KXAR8Lqy627HPgFPAG4CFjem55dddxv36v3Axxr/3hv4JTCn7NpL6NULgecAP51keUf8PXfLfccdBKzKzDsycz2wAjhqwpijgHOy7mrgCRHxpFYXWrJt9ikzr8rMXzUmrwYWtbjGdlHkdwrgBODrwH2tLK6NFOnTG4HzMnMUIDPt1eS9SmD3iAhgHvVw39DaMsuXmVdQ/94n0xF/zw33HbcQuGvc9JrGvKmOqbqp9uCt1N8dz0Tb7FVELAReDZzRwrraTZHfqX2BPSJiOCKujYhjWlZdeynSq38G9gfuBm4A3pWZj7SmvI7SEX/Pdy27gAqIJvMmfgShyJiqK9yDiFhGPdxfsFMral9FenUa8N7M3Fjf0JqRivRpV+BA4EXAY4EfRsTVmXnbzi6uzRTp1UuB64AjgKcC34uIKzPztzu5tk7TEX/PDfcdtwZ48rjpRdTf+U51TNUV6kFELAHOBJZn5i9aVFu7KdKrpcCKRrDvBbwsIjZk5jdaUmF7KPp/7/7MfAB4ICKuAJ4FzLRwL9KrtwAfzfqB5VURcSewH/Dj1pTYMTri77m75XfcNcA+EfGUiJgDvAG4YMKYC4BjGmdZHgL8JjPvaXWhJdtmnyJiMXAe8OYZuGU13jZ7lZlPycyezOwBvgb89QwLdij2f+/fgcMiYteI6AIOBm5ucZ3toEivRqnv4SAiFgBPB+5oaZWdoSP+nrvlvoMyc0NEHA98h/oZqWdl5o0R8fbG8jOon838MmAVMEb9HfKMUrBPJwF7Aqc3tkg3ZAfdqGG6FOzVjFekT5l5c0R8G7geeAQ4MzObfsSpygr+Tp0CnB0RN1Df9fzezGzHu6DtVBHxZaAX2Csi1gAfAGZDZ/099wp1kiRVjLvlJUmqGMNdkqSKMdwlSaoYw12SpIox3CVJqhjDXZKkijHcJUmqGMNdkqSKMdwlTUlEPC0iHo6ID06Y/7mI+F1EzLirCkrtxnCXNCWZuYr6zX3+NiL2AoiIk4C/BF6dmSvLrE+Sl5+VtB0i4onA7cDpwC3AIHB0Zv5bqYVJArxxjKTtkJk/j4jTgL+j/nfknQa71D7cLS9pe/0MeAzww8z8bNnFSPoDw13SlEXEEcC/AD8EDo2IZ5VckqRxDHdJUxIRzwG+Qf2kul5gFPhwiSVJmsBwl1RYRDwN+BbwXeCEzFwPfBB4WUS8sNTiJG3m2fKSCmmcIX8V9S31l2bm7xvzdwF+CvwqM59fYomSGgx3SZIqxt3ykiRVjOEuSVLFGO6SJFWM4S5JUsUY7pIkVYzhLklSxRjukiRVjOEuSVLFGO6SJFXM/weGpySmwtD6TwAAAABJRU5ErkJggg==\n", "text/plain": [ "
![](\"./figures/straight_ray_tomog_schematic.png\")
\n",
"\n",
"\n",
"Now, let the slowness of each block be $x_j$, let's simplify by assuming that the size $h=1$, and the total travel time across the model be $T_i$, then the following equations relate the travel times to the slownesses \n",
"\n",
"$$\\begin{align*}\n",
"T_1 &= x_1 + x_2 + x_3\\\\\n",
"T_2 &= x_4 + x_5 + x_6\\\\\n",
"& \\vdots \\\\\n",
"T_6 &= x_3 + x_6 + x_9\n",
"\\end{align*}\n",
"$$\n",
"\n",
"Given the six measurements, the inverse problem is to determine information \n",
"about the nine slownesses $x_1, \\ldots x_9$.\n",
"\n",
"The equation of condition for this system is \n",
"\n",
"\n",
"$$\n",
"\\begin{pmatrix}\n",
"1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\\n",
"0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\\\\n",
"\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots \\\\\n",
"0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 1\n",
"\\end{pmatrix}\n",
"\\begin{pmatrix}\n",
"x_1\\\\\n",
"x_2\\\\\n",
"\\vdots \\\\\n",
"x_9\n",
"\\end{pmatrix}\n",
"=\n",
"\\begin{pmatrix}\n",
"T_1\\\\\n",
"T_2\\\\\n",
"\\vdots \\\\\n",
"T_6\n",
"\\end{pmatrix}\n",
"$$\n",
"\n",
"\n", "\n", "We can see here why the matrix $G$ is sometimes called the *sensitivity matrix*. \n", "\n", "The element in the $i$-th row and $j$-th column gives the sensitivity $\\partial T_i/\\partial x_j$ of the $i$-th measurement to a change in the $j$-th model variable. \n", "\n", "So, for example, the sixth measurement is only sensitive to $x_3$, $x_6$ and $x_9$, and it is equally sensitive to each of these variables. \n", "\n", "Note that when $\\partial T_i/\\partial x_j=0$ a change in the slowness $x_j$ will not affect the value of the travel time $T_i$, thus we can find no information about the value of $x_j$ from the measured value of $T_i$.\n", "\n", "
\n", "\n", "Let's suppose that there are no errors in the observations, and suppose further than the travel time for all six observations, i.e. along every row and every column, is equal to 6 units ($T_i=6\\; \\forall i$). \n", "\n", "Then clearly a homogeneous \"model\" (or solution vector), for which the slowness is each of the nine blocks is 2 will satisfy the data exactly. We can write this as \n", "\n", "$$\n", "\\begin{pmatrix}\n", "2 & 2 & 2 \\\\\n", "2 & 2 & 2 \\\\\n", "2 & 2 & 2\n", "\\end{pmatrix}\n", "$$\n", "\n", "\n", "However, note that the following solutions also satisfy the data exactly \n", "\n", "$$\n", "\\begin{pmatrix}\n", "1 & 2 & 3 \\\\\n", "2 & 2 & 2 \\\\\n", "3 & 2 & 1 \n", "\\end{pmatrix},\\quad\n", "\\begin{pmatrix}\n", "-2 & 0 & 8 \\\\\n", "-2 & 6 & 2 \\\\\n", "10 & 0 & -4 \n", "\\end{pmatrix},\\quad\n", "\\begin{pmatrix}\n", "1 & 2 & 3 \\\\\n", "2 & 2 & 2 \\\\\n", "3 & 2 & 1 \n", "\\end{pmatrix},\n", "$$\n", "\n", "we've reordered the 9 solution values which we write as a column vector when specifying our linear system into $3\\times 3$ matrices here so you can more easily map them onto the 2D physical solution domain.\n", "\n", "The following are also solutions\n", "$$\n", "\\begin{pmatrix}\n", "2+\\alpha & 2-\\alpha & 2 \\\\\n", "2-\\alpha & 2+\\alpha & 2 \\\\\n", "2 & 2 & 2 \n", "\\end{pmatrix}, \\quad\n", "\\begin{pmatrix}\n", "2 & 2 & 2 \\\\\n", "2 & 2+\\beta & 2-\\beta \\\\\n", "2 & 2-\\beta & 2+\\beta\n", "\\end{pmatrix},\\quad\n", "$$\n", "where $\\alpha$ and $\\beta$ are arbitrary constants.\n", "\n", "\n", "In this problem therefore, there are infinitely many model parameters that satisfy the data. \n", "\n", "Some of these may not satisfy other constraints on the model parameters - for example, the slowness cannot be negative. \n", "\n", "But even with external constraints, in a problem such as this there will still be infinitely many models that satisfy the data exactly. \n", "\n", "
\n", "\n", "**Questions**\n", "\n", "Use `numpy` to compute the rank of this matrix. What dimension do we expect the null space to be?\n", "\n", "Starting from the last two of these example solutions, can you figure out an appropriate number of linearly independent vectors that span the the null space of $G$?\n", "\n", "Can you find a basis for the null space?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework - Tomography example - solution via the pseudo-inverse\n", "\n", "In class we introduced the concept of the data potentially containing errors, i.e. if instead of all data values being 6, consider the case with\n", "\n", "$$ \\boldsymbol{T} = (6.07, 6.07, 5.77, 5.93, 5.93, 6.03)^T $$\n", "\n", "We stated tht now, even though there are fewer data than model parameters, there is no model that fits the data exactly. We can see this because it should be the case from $G$ that\n", "\n", "$$ T_1 + T_2 + T_3 = T_4 + T_5 + T_6$$ \n", "\n", "whereas for this data\n", "\n", "$$ T_1 + T_2 + T_3 = 17.91, \\qquad T_4 + T_5 + T_6 = 17.89 $$\n", "\n", "so that, with these data, there can be no solution. \n", "\n", "
\n", "\n", "Use the pseudo inverse computed using `np.linalg.pinv` to find a solution. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework - Damped least squares applied to the tomography example\n", "\n", "Consider again the tomography example from the previous exercise with noisy data.\n", "\n", "Now try solving this problem using the damped least squares approach and show that as $\\mu$ tends to zero we recover the same solution as found in the previous question using the generalised inverse." ] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "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.7.13" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": true }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 1 }