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", "Now let's consider some random matrices, check their ranks, and the ranks of their products" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The rank of A is 10\n", "The rank of B is 12\n", "The rank of AB is 10\n" ] } ], "source": [ "m = 10\n", "n = 15\n", "p = 12\n", "A = np.random.random((m,n))\n", "B = np.random.random((n,p))\n", "\n", "\n", "print('The rank of A is ', np.linalg.matrix_rank(A))\n", "print('The rank of B is ', np.linalg.matrix_rank(B))\n", "print('The rank of AB is ', np.linalg.matrix_rank(A@B))\n", "\n", "\n", "# you could extend this by inventing some matrices that don't have full rank" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework - A simple mixed determined problem [solved by hand and using the pseudo-inverse]\n", "\n", "Consider the problem from the lecture\n", "\n", "$$\n", "\\begin{pmatrix}\n", "1 & 0 & 0 \\\\\n", "1 & 0 & 0 \\\\\n", "0 & 2 & 2 \\\\\n", "0 & 3 & 3\n", "\\end{pmatrix}\n", "\\begin{pmatrix}\n", "m_1\\\\\n", "m_2\\\\\n", "m_3\n", "\\end{pmatrix}\n", "=\n", "\\begin{pmatrix}\n", "1\\\\\n", "2\\\\\n", "2\\\\\n", "3\n", "\\end{pmatrix}.\n", "$$\n", "\n", "Can you come up with a sensible looking \"solution\" to this problem by considering the under- and over-determined components separately?\n", "\n", "Hint: consider problems for $m_1$ and both $m_2$ and $m_3$ separately, and think about what the least squares and minimum norm solution would be for each.\n", "\n", "Once you have found the solution \"by-hand\", try using the pseudo-inverse $A^+$ which you can find using the function `np.linalg.pinv(A)`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution\n", "\n", "Since the two parts of this problem are uncoupled (this won't generally be the case), we can consider solving for $m_1$ separately to $m_2$ and $m_3$ (which are coupled to one another).\n", "\n", "The first two equations tell us that\n", "\n", "$$m_1 = 1 \\qquad \\textrm{and} \\qquad m_1=2$$\n", "\n", "without additional information we could simply take \n", "\n", "$$m_1 = \\frac{3}{2}$$\n", "\n", "Notice that this does correspond to the least squares solution to this problem.\n", "\n", "\n", "
\n", "\n", "The second part of the problem tells us that \n", "\n", "$$2m_2 + 2m_3 = 2 \\qquad \\textrm{and} \\qquad 3m_2 + 3m_3 = 3$$\n", "\n", "but these two just boil down to one equation\n", "\n", "$$m_2 + m_3 = 1$$\n", "\n", "So for any $\\alpha\\in \\mathbb{R}$, $m_2 = \\alpha$ and $m_3 = 1 - \\alpha$ is a solution.\n", "\n", "We have infinuitely many possible solutions, we've seen previously that seeking the minimum norm solution may be a sensible approach (if we know that $\\boldsymbol{m}$ is the perturbation from some a priori model guess for example).\n", "\n", "We say an example of how to find the mimimim norm solution from an infinite family of solution at the end of L3 - we form the norm as a function of $\\alpha$ and find the $\\alpha$ that minimises this - differentiate and set to zero:\n", "\n", "for simplicity form the square of the two norm of our solution:\n", "\n", "$$f(\\alpha) = \\alpha^2 + (1 - \\alpha)^2 =2\\alpha^2 - 2\\alpha + 1$$\n", "\n", "$$\\frac{df}{d\\alpha} = 4\\alpha - 2$$\n", "\n", "and so $\\alpha = 1/2$ and a sensible solution to this problem could be taken to be\n", "\n", "$$\\boldsymbol{m} = \\left(\\frac{3}{2} \\;\\; \\frac{1}{2} \\;\\; \\frac{1}{2} \\right)^T$$\n" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([1.5, 0.5, 0.5])\n" ] } ], "source": [ "G = np.array([\n", " [1, 0, 0],\n", " [1, 0, 0],\n", " [0, 2, 2],\n", " [0, 3, 3]])\n", "\n", "d = np.array([1, 2, 2, 3])\n", "\n", "m = np.linalg.pinv(G) @ d\n", "\n", "pprint(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework - A \"least squares\" and minimum norm solution to a simple $2\\times 2$ case\n", "\n", "In the L2 homework we considered the case \n", "\n", "$$\n", "\\left(\n", " \\begin{array}{rr}\n", " 2 & 3 \\\\\n", " 4 & 6 \n", " \\end{array}\n", "\\right)\\left(\n", " \\begin{array}{c}\n", " x \\\\\n", " y\n", " \\end{array}\n", "\\right) = \\left(\n", " \\begin{array}{c}\n", " 4 \\\\\n", " 7\n", " \\end{array}\n", "\\right),\n", "$$\n", "\n", "plotting the two constraints to make the point that this system has no solution:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU4AAAFSCAYAAACdRzaJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4uklEQVR4nO3deZyN9fvH8dc1m22sYYikRUqIBmXfsyZZQqQILbK0CPmmFe2lRYuUUqJUVLaIScqvb3ZKop2EkjW7z++Pa/hK92HOzDn3OWfmej4e52HmzOeccwnvPvf92cQ5hzHGmIyLi3QBxhgTayw4jTEmSBacxhgTJAtOY4wJkgWnMcYEyYLTGGOCZMFpchQR+UlE7gj0vTEZkRDpAkzOIiKnmjj8mnPuOj9qSVcd2OPj55lswILT+K3kcV+3Bsae8NxeP4txzm318/NM9mCX6sZXzrnfjz6A7UefA3YDvwLnH20rIhtEZM1x3zcVkT0ikpj+fRkReV9EdqU/3hOR0sHU43Hp7kSkj4i8k/5ZP4hItxNeU0pEJonIX+mP6SJS7rifnyEi00Rkm4j8LSLfikjnoP5DmahmwWmignNuN7AUaAiQHkQFgbIicrRH2gD4wjl3UEQEmAqkAI3SX3c6MDX9Z1kxHJgGXARMBl4RkTPT68oLzAf2AfWBmsAmYG76zwDGAHnTa7oQGEj6/yRM9mDBaaJJGunBiYbkQuC/6V8ffS4t/esmaLBd7Zz7yjm3GLgauBhonMU6Jjjn3nDOrQfuBg4BddN/1hkQoIdzbqVz7lvgBiAZvfUAcCaw0Dm3wjn3o3NulnNuVhZrMlHEgtNEkzSgdvqleAO0Z5cGNEjvzVXnf8F5AfCbc+6noy92zv0A/AZUyGIdK497z0PAVqB4+lOpwFnALhHZLSK7gR1AYeCc9Dajgf+IyCIReVBEUrNYj4kyNjhkoslnQC40IOsDT6E9uReB2sBBtAcK2usLNEKf1S2/Dnq839FORhywHO15nmgbgHNunIjMBlqiPeMvRGSUc+7eLNZlooT1OE3UOO4+Zx8gf/rXi4AyQFfS72+mN/8GKCUiZY++XkTORu9zfhPGMpcC5wJ/OOfWn/DYdtzvZYNz7iXn3FXoPdM+YazJ+MyC00SbNKAb8Jlz7rBzbh/wZfpzace1mwusAN4UkVQRqQa8iQbbvDDW9yawGZgmIvVF5CwRqScijx8dWReR0SLSXETOFpEqQHPCG+bGZxacJtrMB+L5Z0j+6zmnO3C3Re8/pqW3+R1o68K4O7dz7m+gHvAD8A7wLfAaeo/zr/RmccAzaFjOQYP22nDVZPwntgO8McYEx3qcxhgTJAtOY4wJkgWnMcYEyYLTGGOCZMFpjDFBivmVQ0WLFnVly5YN6jV79uwhX7584SkohKzO0IqVOiF2as3OdS5ZsuQP51wxzx8652L6kZqa6oI1f/78oF8TCVZnaMVKnc7FTq3ZuU5gsQuQO3apbowxQbLgNMaYIFlwGmNMkCw4jTEmSBacxhgTpJifjmRMVuzcuZMtW7Zw8OCJexdHVsGCBVmzZs2pG0ZYrNaZmJhI8eLFKVCgQKbez4LT5Fg7d+5k8+bNlCpVijx58pD1M95CZ9euXeTPnz/SZZxSLNbpnGPv3r1s3LgRIFPhGVWX6unHqs4XkTUi8rWIDIh0TSb72rJlC6VKlSJv3rxRFZomvESEvHnzUqpUKbZs2ZKp94i2Huch4Hbn3FIRyQ8sEZE5zjnbPduE3MGDB8mTJ0+kyzARkidPnkzfoomqHqdzbpNzbmn617uANUCp0L0/PPEE7N4dH6q3NDHOepo5V1b+7KOtx3lM+iFcVdHzZk78WR/SD79KSUkhLS0tQ++5fn0+7rwzlTPPrMzBg59TuHB0DQicaPfu3Rn+vUVSrNZZsGBBdu3aFbmCTuLw4cNRW9vxYr3Offv2Ze7vbqC1mJF8oEfCLgHanaptsGvVZ8xwLleuQ658eed++SWol/ouO68DjoQT6/zmm28iU0gG7Ny5M9IlZEis13myvwPE0lp1EUkE3gXedM69F+r3b9ECHnlkJZs2QZ068N13of4EY0x2F1XBKXrTYRywxjn3RLg+p3LlHaSlwd69ULcuLF8erk8yJjxGjRpF9erVKVCgAMWKFePyyy9n9erVvtfx0ksvUblyZQoUKECBAgWoWbMm06dP972O440cORIR4ZZbbgnbZ0RVcAK1gWuARiKyPP3RMhwfVLUqfPYZJCVBgwbw+efh+BRjwiMtLY2bb76ZL774gnnz5pGQkECTJk3Ytm1blt/7uuuu4957781Q21KlSvHwww+zdOlSFi9eTKNGjWjbti0rV67Mch2Z8X//93+MHTuWypUrh/Vzoio4nXMLnXPinKvsnKuS/pgRrs8rXx4WLoTixeGyy+Djj8P1ScaE1uzZs+nRowcVK1akUqVKTJgwga1bt/J5eg/gnXfeIVeuXPz888/HXjNgwADOOeccNm/eHLI6WrVqRYsWLTj33HM577zzGDFiBPnz52fRokXH2pQuXZonnvjnBeSqVavInTs333wTupmGO3bsoGvXrowbN47ChQuH7H29RFVwRsKZZ2rPs1w5aN0apkyJdEXGBG/Xrl0cOXLkWGB06NCBSpUq8eCDDwLw2GOP8dZbbzFr1ixSUlLCUsPhw4eZNGkSu3fvplatWseer1mzJl999dU/2g4cOJBevXpRoUKFfzw/cuRIkpOTT/r47LPPPD+/T58+dOjQgUaNGoX+N3eCqJ2O5KeUFEhL0+Ds1AnGjoWePSNdlfHbwIH+3++uUgWeeirr7zNgwACqVKlCzZo1AZ2jOHLkSFq1asU555zDiBEjmDdvHuXKlcv6h51g1apV1KxZk3379pGcnMz7779PpUqVjv28Zs2ajBkz5tj3U6dOZdmyZbz99tv/eq8bb7yRq6666qSfV6rUv6d2jx07lvXr1zNhwoQs/E4yzoIzXaFCMHs2tG8P118PO3bArbdGuipjTu22225j4cKFLFy4kPj4/y3uuOyyy6hevTr/+c9/+PDDD6levbrn60eOHMnIkSOPfb9//35EhMcee+zYczNnzqRu3bqery9fvjzLly9n+/btvPvuu1x77bWkpaVRsWJFAC699FJuv/12tm3bRr58+bjjjjsYPnw4p5122r/eq0iRIhQpUiSo3//atWu56667+Oyzz0hKSgrqtZkWaJ5SrDyCmsd54IBzAwe6z99+O2CTffuc69DBOXDu7rudO3Ik428farE6PzJaZcd5nAMHDnQlSpRwa9as+dfPPvnkE5c3b14XFxfnlixZEvA9/vzzT7du3bpjj3bt2rl+/fr947m///47w3U2btzY9ezZ89j3+/btc0lJSW727NnuoYcecuXLl3cHDhzwfL8RI0a4fPnynfSxYMGCf7zm1VdfdYCLj48/9gCciLj4+Hi3b9++kM/jzFk9zrVrYdw4qk6apMPq5577rya5csGkSXDDDfDAA7B9u15KxeX4u8Em2gwYMIBJkyaRlpbG+eef/4+frVixgnbt2vHMM88wffp0hg4dyuzZsz3f58ReXv78+SlSpAjnevz7yIgjR46wf//+Y9/nypWLqlWr8uGHH/Laa68xceJEEhMTPV+bmUv1tm3bUq1atX8816NHD8qVK8ddd91FUlISBw4cyNTvJZCcFZwVK8K8eSQ0bqyz3z/+GDymLcTH633OQoXg8cc1PF95BRJy1n8tE8X69u3LhAkTmDp1KoULF+b3338HIDk5mT///JOWLVty22230bNnT2rUqEHlypVJS0ujQYMGIa3jnnvu4corr+SMM85g165dTJw4kbS0tH/N5axZsyajR4+madOmtG7dOuD7ZeZSvVChQhQqVOgfz+XLl48iRYocu10QajkvCqpVY9no0dQYNgzq14cZMyD9hvrxRODRR6FwYfjPf2DnTu2J5s4dgZqNOcHRwZbGjRv/4/l+/foxZ84cWrduzfDhwwGoWLEiHTt2ZOjQof+YJhQKmzdvplu3bvz+++8ULFiQypUrM3PmTJo1a/aPdlWqVCEuLu5f05JiVc4LTuDvsmV1xnuTJvqYNk1/PYEIDBumPc9bboFWrWDqVIiBfVtNNqe34DJu8uTJGW47fvz4DLd94YUXMrSR8ZtvvskNN9zAhRdemOH3zopwbzqTc+/clS2rEzjPOUcT8f33Azbt2xcmTIBPP9V8DcHiDGOyvSNHjrB582YefvhhVq1adWxOaXaQc4MToGRJncB58cXQoQO89lrApt26wXvvwYoVUK8e/Pabf2UaE4sWLFhAyZIlGT9+PO+++27YV/P4KUdeqv9DkSIwZw5ceSVcd51O4Ozf37NpmzYwc6b+Wreuvuzss/0t15hY0aBBA44cORLpMsIiZ/c4j0pOho8+gnbtYMAAuO8+3S7eQ8OG8MknOtJepw58/bW/pRpjIs+C86hcuWDyZO113nuvLhsK8H/LGjX0fifoZfsJy3CNMdmcBefxEhJg3DjtdY4erWsvDx3ybFqxou6sVLAgNGoE8+f7XKsxJmIsOE8UFwdPPqm9zvHjddeP41ZBHO/sszU8zzxTd5b/4ANfKzXGRIgFpxcRuOceXWv53ntw+eWwZ49n09NP18v2iy7SW6RvvOFvqcYY/1lwnsyAAdrr/OQTaNoU/vrLs9lpp8HcuXq/85pr4Nln/S3TGOMvC85TufZa3d14yRI9YyN9TfCJ8ufX1Ztt2kC/fvDggwEH5o0xMc6CMyOuvBKmT4fvv9cJnD/95Nksd27N2G7d4O67YdAgC09jsiMLzoxq0kRnvP/xh07gXLPGs1lioi5AuuUW3Vmpd284fNjnWo0xYWXBGYyaNXUk6NAh7XkuWeLZLC4Onn5ad1UaNw46dw44MG+MiUEWnMGqXFnnICUn6zKiBQs8m4noRsiPPaaX71dcEXBg3pgs8eMc8ZPZtGkT1157LcWKFSN37txUqFCBT4+uEPFJ2bJlEZF/PVq1ahWWz7PgzIxzz9Vt6UqXhmbN9P5nALffDi+/rFf5zZrpUk1jQiUc54gHc6769u3bqV27Ns45pk+fzpo1a3jmmWcoXrx4yOrJiK+++opNmzYdeyxduhQROeVu8pllwZlZpUppb7NiRWjbFt56K2DT66/XTZD/+1/tpG7Z4l+ZJvs62Tnifp2rPnr0aEqWLMnrr79OjRo1OOuss2jcuDEXXHDBsTZ+nKterFgxSpQocewxY8YMChQoQMeOHUPy/iey4MyKokV1jmft2tC1K7zwQsCmHTvChx/qsUd168Ivv/hYp8mWTnaOuF/nqn/00UdccskldOrUieLFi1OlShWeffbZf2y07Ne56kc55xg3bhzdunUjb968Ifu9Hs+2lcuqAgV0r7mOHeGmm/RafMgQz6bNmukle6tWOjA/dy6cd56/5ZqTiKGD1U91jrhf56r/9NNPjBkzhltvvZUhQ4awfPly+vXrB3Dsnqsf56ofb86cOfz444/06tUr2N9OhllwhkKePLqD/LXXwtChGp6jRukI0Qlq19a9k5s10/CcPVsP3DQmozJ6jrgf56ofOXKEatWqMWrUKACqVq3KunXreO65544FZ7jPVT/R2LFjqV69OlWqVMnS+5yMBWeoJCbq+RoFC8LDD2t4PvecHpl5gipV9NSOpk11MdL06RqiJsIy0fOLhEWLFvHHH3/84wTHw4cPs2DBAl544QX27NlDrly5mDdvHitWrMA5d9LL8xN7eYMHD6ZUqVL0P25D70C9vBIlSvzrcvuCCy5g9OjRx75PTU0lKSmJxYsXs2zZMhISEujbt6/n+50Y4l4ChTjAli1bmDZtGs8999xJ3yOrLDhDKT4exozRozFHjdLd5F9/XUP1BOedp7OamjaFyy7TvUSaN49AzSbmZOQccb/OVb/kkktYu3btP5777rvvOPPMM499H+5z1Y83fvx4cuXKRefOnU9Ze1ZYcIaaCIwcqUdjDh6s5wq/8w543KQ+4wwdmG/eXNe4v/EGhGn2hMlGTnWO+M8//+zbuep9+/aladOmjBgxgk6dOrFs2TKefvrpf/Uaw3mu+lHOOV5++WU6d+6coZM3s8JG1cPlzjvhxRd14Kh5c+19eiheXDdBvuQSXWH08ss+12mylW3bttG8efOA56qHWmpqKlOnTuXtt9+mYsWKDBs2jAceeICbb775H+38OFc9LS2NdevW0bt377B9xlHW4wynPn30nme3brpN/KxZUKzYv5oVLKiDRO3b69r27dvhjjv8L9fEruPPEV/jsY9CuM5VB2jVqtUpV+j4ca56w4YNgz5vPrMsOMOtUyfdc659e53AOXeurjg6Qd68MG2a7uc5aJBu/dmkSQTqNSZEjhw5wtatWxk/fjyrVq0KKryjnV2q+6FlS/j4Y9i0SYfP163zbJaUBBMnQq9eept09Ohygc6LMybqZedz1S04/VK3rt7M3LNHv16xwrNZfDy89JJeqk+bVoru3eHgQZ9rNSYEjp6rvmbNGmrVqhXpckLKgtNPF1+sEzgTE3UC56JFns1E4JFH4Prrf+DNN/Uqf+9ef0s1xgRmwem388/XCZzFiv1vc2QPItCt2y8895yucW/RQmc2GWMiz4IzEs48U3ue5crpwvV33w3Y9OabdX7nwoXQuLFuQG+MiSwLzkhJSdF7ntWr66z3V18N2LRrV10Kv2oV1K8PGzf6WGc259f0FRN9svJnb8EZSYUL62h7kybQs+dJ10pffrlOA/3lFx1b+uEH/8rMrhITE9lrN49zrL179wZc+nkqFpyRli8ffPABdOgAt94K99wT8GjMBg20k7pzp85qWr3a31Kzm+LFi7Nx40b+/vtv63nmIM45/v77bzZu3JjpneptAnw0yJVLd5AvUADuv19nvwfofVarpuvbmzaFevV0Recll/hbbnZRoEABAH777TcORtmcr3379pE7d+5Il3FKsVpnYmIiKSkpx/4OBMuCM1okJOhC9YIF4cknYccOpHt3z6YVKuhgUZMmOmA0bZr+aoJXoECBTP/jCae0tDSqxsBGrTm1TrtUjyYiehj7/ffD669T4d57Yd8+z6ZnnaXhedZZujBp2jR/SzUmJ7PgjDYicPfd8PTTFFu4UKcr7drl2bRkST3mvWpVnSQf4BQFY0yIWXBGq379WDNkiCZj06awbZtnsyJFdN+QBg2ge3d45hl/yzQmJ7LgjGKbmzWDKVNg2TKdwLlpk2e75GT46CM9pbh/f3jggYAD88aYELDgjHZt28KMGfDjjzoH6ccfPZvlzq0bzV97LQwfDrffbuFpTLhYcMaCxo31/Pa//tLw/OYbz2YJCfDKK9rrfPJJuP56OHTI51qNyQEsOGPFJZfoBM4jR3QC5+LFns3i4nQK6D336CrOTp1g/35/SzUmu4u64BSRV0Rki4jYupgTVayoc5AKFNCjOI47LuF4InDvvdrrfO89Xa65Z4+vlRqTrUVdcALjATsoN5BzztGdlc44Qw+B+/DDgE0HDtRL908+0YH5v/7yr0xjsrOoC07n3ALAe+6NUaVK6WV75cpw5ZXw5psBm/boAW+/rVf2DRrA5s3+lWlMdhV1wWky6LTTtCtZr56e8DZmTMCm7dvrdKX163Vs6eeffazTmGxIonFXGBEpC3zknKsY4Od9gD4AKSkpqZMmTQrq/Xfv3k1ycnJWywy7jNQZd+AAFe67j6JffMEPvXrxy9VX601OD6tXF2Do0ErkyXOYxx5bSZkyf/tWZzSIlTohdmrNznU2bNhwiXOumucPnXNR9wDKAqsz0jY1NdUFa/78+UG/JhIyXOeBA8517eocODdokHNHjgRsumKFcykpzhUt6tySJT7XGWGxUqdzsVNrdq4TWOwC5I5dqmcHiYnw+uvQty88+ijccAMcPuzZtHJlHVvKlw8aNtRbpcaY4ERdcIrIW8AioLyIbBCR6yNdU0yIi9OF6sOGwdixcPXVcOCAZ9Ny5XRW0+mnQ7NmujDJGJNxUbcfp3OuS6RriFki8OCDUKgQDBqkuypNmQJ58/6raenS2tts3hyuuEIPhOvUyf+SjYlFUdfjNCFwxx3a65w1S7uUO3Z4NitWDObNg1q1oEsXfYkx5tQsOLOrXr1g8mT48ku9mblli2ezggU1X1u0gD594JFHfK7TmBhkwZmddeyoB8F9+60ejfnLL57N8uTR44c7d4bBg2HoUNtZyZiTseDM7po31yOIf/9dZ79/951ns6Qkvc95ww3w0ENw8826n4gx5t8sOHOCOnV0Q5B9+/TrZcs8m8XHw/PPw5Ah8MILuiApyg5/NCYqWHDmFFWr6gTO3Ln1nufnn3s2E4FRo7TXOXEitGsHe/f6XKsxUc6CMycpX14ncKak6HZJs2YFbDp4sPY+p0/XgaOdO32s05goZ8GZ05Qpoz3P8uWhTRs9byOAG2/UjZc+/1y3//zjDx/rNCaKWXDmRMWLw/z5UKOGDqWPGxewaZcuMHUqfP21bsS0YYN/ZRoTrSw4c6pChWD2bL1k79ULnngiYNNWrfSqfsMGHVtat86/Mo2JRhacOVm+fDrPs2NHPRbz7rsDTuCsX19XGe3erVNCV670uVZjoogFZ06XlARvvaVHYj74oB6RGWACZ7Vqur49IUGDdNEin2s1JkpE3SYfJgLi43WheuHC8Nhjurb9lVc0IU9QoYIOzDdtCk2a6P3PxET/SzYmkqzHaZSILlQfMQImTIAOHXTCvIeyZXVg/pxzoHVrWLCgqL+1GhNhFpzmf0Tgrrvguedg2jRo2VK3pvNQooQuRrr4YrjvvgsZP97XSo2JKAtO828336wL1xcs0OvxP//0bFakCMyZA1Wr/kWPHjB6tM91GhMhFpzGW9eu8N57sGKFjgT99ptns+RkGDlyFe3a6Tnu991nOyuZ7M+C0wTWpg3MnKnnCdepAz/84NksKckxeTJcdx3cey/ceqvtrGSyNwtOc3ING+r57Tt2aHiuXu3ZLCFBFyANHKiX7NdfD4cO+VuqMX6x4DSnVqPG/47DrF8f/vtfz2ZxcboA6b77YPx4uOoq2L/fvzKN8YsFp8mYCy/UCZyFCkHjxrqMyIMIDB+uvc7339fpSrt3+1uqMeFmwWky7uyzdQLnmWfqVKVp0wI27d9fe53z5ulk+W3b/CvTmHCz4DTBOf10+PRTqFwZ2rfXyfIBXHutnk68dCk0aKCndxiTHVhwmuCddpoOGNWrB927U+q99wI2vfJK3Qz5hx90bOmnn/wr05hwseA0mZM/P8yYAVdcQblnntENQgJM4GzSBObO1Xn0derAmjU+12pMiFlwmszLnRumTOH3pk11S7pBgwKG56WX6sD84cO6Ld3ixT7XakwIWXCarElI4NshQ+CWW+Dxx6F3b01HD5Uq6dhS/vx6FMenn/pcqzEhYsFpsi4uDp5+Wnud48bpcRwBJnCee67OaipdWo98nz7d51qNCQELThMaInD//drrnDIFrrgC9uzxbFqqlF62X3ghtG2r+ygbE0ssOE1o3Xab9jrnzIFmzWD7ds9mRYvqHM/atXU/kRde8LdMY7LCgtOEXs+eMHmyLs1s2BA2b/ZsVqCA7iHSsiXcdBM89JDPdRqTSRacJjw6dIAPP4TvvtP5nr/84tksTx5dmnn11TB0KAwZYtvSmehnwWnCp1kzvWTfvFkncK5d69ksMVEXIN10Ezz8sP4aYGDemKhgwWnCq1YtnXe0f79O4Fy2zLNZXJye2DF0KLz4InTrBgcP+lyrMRlkwWnC76KLdA5Snjy6aP2zzzybicDIkdrrnDRJR9z//tvXSo3JEAtO449y5TQ8S5bUS/iZMwM2vfNO7XXOnAktWugeysZEEwtO458zztDe5gUX6LEckycHbNqnj87v/OILXWW0dauPdRpzChacxl/FiukEzksvhS5dYOzYgE07ddItP7/5RgfmN2zwsU5jTsKC0/ivYEGYPVvXXPbpA48+GrBpy5bw8cd6yGadOrBunY91GhOABaeJjLx5YepUPZjozjvhrrsCTuCsWxfmz9cVnHXr6onFxkSSBaeJnKQkmDhRd1QaNQr69g14rvDFF+vt0cREPS/uiy98rtWY41hwmsiKj9ch9EGD4Pnn4ZprAk7gPP98HZgvXlzPMZozx+dajUlnwWkiT0Qnb44cqT3Q9u1h717PpmeeqT3PcuWgVSt4912fazUGC04TLUR02dCYMfDRRzoqtGuXZ9OUFEhLg+rV9Rbpq6/6W6oxFpwmutx0E7zxhnYrGzWCP/7wbFaokI62N2mimzE9+aS/ZZqczYLTRJ+rr9YR99WrdSRo40bPZvnywQcf6JX9bbfB8OG2s5LxhwWniU6tW8OsWfDrrzqB8/vvPZvlyqULkHr2hAcegAEDAg7MGxMyFpwmetWvr6uMdu3S8Fy1yrNZfDy8/DLceis88wz06AGHDvlcq8lRLDhNdKtWTQ8oiovTIP3yS89mInrc0f33w+uvQ8eOsG+fz7WaHMOC00S/ChV0AmeRItC4MXzyiWczET1o8+mn9RZpq1YBB+aNyZKoC04RaS4ia0VkvYgMiXQ9JkqcdZaOtJ99tk5Vmjo1YNN+/bTX+emnOlF+2zb/yjQ5Q1QFp4jEA88BLYAKQBcRqRDZqkzUKFlSJ3BefLGeafT66wGbXnONTo5ftkyv8Ddt8q9Mk/1FVXACNYD1zrkfnHMHgEnAFRGuyUSTIkV0rWWDBnDttXpdHsAVV8CMGfDjjzq2tGlTbv/qNNlatAVnKeDX477fkP6cMf+TnAzTp8OVV+r8o/vvDziB8+gt0b/+gv79q/LNNz7XarIlcVE0Y1hEOgLNnHO90r+/BqjhnOt3Qrs+QB+AlJSU1EmTJgX1Obt37yY5OTk0RYeR1Xlycvgw5R99lBKzZ/Nrhw58f9NNOvru4ccf83HHHZU4dCiehx9eyfnnR/eokf3Zh1Zm6mzYsOES51w1zx865076AM4/VZtQPYCawOzjvh8KDD3Za1JTU12w5s+fH/RrIsHqzIDDh53r3985cO6665w7eDBg0zffXOTOOsu55GTnov0/rf3Zh1Zm6gQWuwC5k5FL9WUiMlpECgcV15nzFVBORM4SkSSgM/CBD59rYlVcHDz1FNx7L4wfr+dt7N/v2fT00/excCGUKaObz3/4oZ+FmuwkI8FZA7gQWCci/dJHvsPCOXcIuAWYDawB3nbOfR2uzzPZhAjcc48G6HvvweWX63bxHk4/XefTV66st0jffNPfUk32cMrgdM6tcs41AXoB/YFVItIiXAU552Y4585zzp3jnBsRrs8x2dCAAbrH3Cef6ATOv/7ybHbaadqkXj2dtjRmjM91mpiX4VF159xUtOf5GjBJRGaIyPnhKsyYTLnuOnjnHViyRKcs/f67Z7P8+XWq0uWX64kdI0fazkom44KdjpQXWIKGZzNgpYg8LSIFQ16ZMZnVrp1uhrx+vZ7u9vPPns1y54YpU6BbNxg2DAYPtvA0GXPK4BSRgSLypoh8B/wJfAhUB0ajl+/lgW9E5JKwVmpMMI4eSvTHHzr7/dtvPZslJsJrr2mv89FH4YYb4PBhn2s1MScjPc7bgUTgeaAOUNA5V9M5d5tz7nXnXDPgaeCVMNZpTPBq1dIF6wcPQt26JH/3nWezuDjdju4//4GxY3Uf5QMHfK7VxJSEUzVwzp2Rgfd5FRiZ9XKMCbHKlXVzkKZNqXLbbXpUZr16/2omohshFyyoB27u3Klr3fPmjUDNJuqFasnlVqBRiN7LmNAqVw4WLmT/aadBs2Y6KhTAHXdor3P2bLjsMti+3b8yTewISXCmT7T/NBTvZUxYlC7N8qefhgsv1N0/TrJMt1cv/fF//wsNG8KWLT7WaWJCtG3yYUzYHCxYUI/iqFVLb2S++GLAtlddpQfBrV2rA/O//OJjoSbqWXCanKVAAT0ErmVLuPFGePjhgE2bN9cjiH//XQfmA4wtmRzIgtPkPHnywPvvQ5cuMGQIDB0acAJnnTq6d/K+fdrzXL7c10pNlLLgNDlTYiJMmKC9zocegptvDniucNWqOjCfK5cuRvr8c39LNdHHgtPkXPHxulB9yBB44QVdQnTwoGfT8uX1vLiUFJ1bP2uWz7WaqGLBaXI2ERg1Snudb72lWybt3evZtEwZ7XmWLw9t2uiSeJMzWXAaA7pQ/YUXdI5n8+Y6A95D8eIwfz7UqAGdO8O4cT7XaaKCBacxR91wg27Q+cUX0KiRrnP3UKiQjrZfdpnO+Xz8cX/LNJFnwWnM8bp00TPbv/5ah9E3bPBsljcvTJsGHTvqaqO777adlXISC05jTtSqla653LhR5yOtX+/ZLClJb4v26gUPPgj9+wccmDfZjAWnMV7q1dObmbt3a3iuXOnZLD4eXnpJe53PPqtHvQcYmDfZiAWnMYGkpuowekIC1K8PixZ5NhOBRx6BESPgjTegQwedMG+yLwtOY07mggt0AmfRotCkiW6O7EEE7roLnntO17i3bAm7ovvodpMFFpzGnErZstrzPPdcaN1al2sGcPPN2utcsEBz9s8//SvT+MeC05iMKFFCF62npuq1+PjxAZt27arZumKFXuH/9ptvVRqfWHAak1GFC+uleuPG0KMHjB4dsOnll8PMmXpOXJ068MMPPtZpws6C05hg5MsHH36oJ2kOHAj33htwAmfDhrr9544dGp6rV/taqQkjC05jgpUrF0yerGe433cf3HprwAmc1avr/U7QGU5ffulfmSZ8LDiNyYyEBF2oPnCgXrL37AmHDnk2vfBCHZgvXFiv8ufN87dUE3oWnMZkVlwcPPGE9jpfe03XXwaYwHn22TowX7asTlWaNs3fUk1oWXAakxUiMHy49jqnTtXpSrt3ezY9/XQ95v2ii6B9e91H2cQmC05jQqF/f+11pqXpTsfbtnk2O+00mDtXpyl1767LNE3sseA0JlS6d4cpU2DpUj1jY9Mmz2b588P06XpKcb9+ukGI7awUWyw4jQmltm11M+QfftBt6X76ybNZ7tyasddco1vSDRpk4RlLLDiNCbXGjfV6fNs2qF0bvvnGs1lCgi5AuuUW3Qy5d284fNjfUk3mWHAaEw6XXqojQUeO6ATOxYs9m8XFwdNPa69z3Dg9jmP/fp9rNUGz4DQmXCpV0jlI+fPrURyffurZTATuv19nNk2ZogfB7d1r/zSjmf3pGBNO556rs99Ll9ZD4D76KGDTW2/VXufcuXDnnRexfbt/ZZrgWHAaE26lSum6y4oV9fjhiRMDNu3ZE95+G779Nj8NG8LmzT7WaTLMgtMYPxQtCp98ooNF3brB888HbNq+PYwYsYrvvtOB+Z9/9rFOkyEWnMb4pUAB3WuuVSvd8XjUqIBNa9T4izlzYMsW3Vnp2299rNOckgWnMX7Kkwfeew+uvlrP2hg8OOAEzlq1dDzpwAHteS5d6nOtJiALTmP8lpioC9VvuklPebvxxoATOC+6SMeW8ubV/T0/+8znWo0nC05jIiEuTk92GzpUzxfu2lW7lh7KldPwLFkSmjXTq30TWRacxkSKCIwcqb3OyZN1xP3vvz2bnnGGDsyff77O85w82edazT9YcBoTaYMGaa9z5kyd67ljh2ez4sVh/nxdlNSlC4wd63Od5hgLTmOiQe/e8NZbsGgRNGxI4l9/eTYrWBBmz9ZL9j59tLNq/GfBaUy06NQJPvgAvv2WqgMGwK+/ejbLm1d3kL/qKh2Uv+su21nJbxacxkSTFi1g9myStm3TCZzr1nk2S0rSBUi9e+t00L59A54XZ8LAgtOYaFO3LsufeAL27tXwXLHCs1l8PLz4Itx5py5EuuYaOHjQ51pzKAtOY6LQ7vPO00mbuXLpbvJffOHZTgQeflgH5ydO1OWae/f6W2tOZMFpTLQqX14ncBYvrucYffxxwKZDh8KYMbr5UosWsHOnj3XmQBacxkSzMmW053neeXqC5pQpAZvedBO88YZmbePG8McfPtaZw0RNcIpIRxH5WkSOiEi1SNdjTNQ4OoGzenUdeX/llYBNr75aTylevVpP0ty40b8yc5KoCU5gNdAOWBDpQoyJOoUK6aV6kyZw/fW6XXwArVvDrFk6m6lOHfj+e//KzCmiJjidc2ucc2sjXYcxUStfPp3n2b493H47DB8ecAJn/fowbx7s2qXhuWqVz7Vmc1ETnMaYDMiVCyZN0q3iH3gA+vcPOIGzWjVd3x4fr0H6f//nc63ZmDgflxyIyFyghMePhjnnpqW3SQPucM55HwuobfoAfQBSUlJSJ02aFFQdu3fvJjk5OajXRILVGVqxUidkoFbnOOf55znjnXf4vWlT1g4ejIuP92z6+++5uf32i9i2LYkHH1xFaup2/+qMEpmps2HDhkucc97jLc65qHoAaUC1jLZPTU11wZo/f37Qr4kEqzO0YqVO5zJY65Ejzj3wgHPgXJs2zu3dG7Dpb785V6mSc0lJzr33ns91RoHM1AksdgFyxy7VjYlVIvCf/8Azz+i9z1at9Kamh5IlIS0NLr4YOnSA117zt9TsJmqCU0SuFJENQE1guojMjnRNxsSEW27RJPz0Ux1137bNs1mRIjBnjh7xft118PTT/paZnURNcDrn3nfOlXbO5XLOpTjnmkW6JmNiRvfu8O67sHw51KsHmzZ5NktO1tVF7drBgAFw//22s1JmRE1wGmOy6IordDPkn37SOUg//ujZLFcu3UH+uuvgnnvgtttsZ6VgWXAak500aqQTOLdv1zPcv/7as1lCAowbp7OZnnpK59QfOuRrpTHNgtOY7KZGDb3fCXrZ/tVXns3i4jQ0hw+H8eN1Nef+/b5VGdMsOI3JjipW1N0+ChbUXuj8+Z7NROC+++DJJ/W498svhz17fK41BllwGpNdnX22hmeZMrrX3AcfBGw6cKDuHfLJJ7qDXYAjj0w6C05jsrPTT9d1l5Ur61D6G28EbNqjB7zzDixZonsnb97sX5mxxoLTmOzutNO0K1mvnp6vMWZMwKbt2ul0pfXrdWD+5599rDOGWHAakxPkzw8zZkCbNnqy28iRASdwNm0Kc+fqRsi1a8OaNT7XGgMsOI3JKXLn1kny11wDw4bpKW8BwrNmTR2YP3RIO6pLlvhca5Sz4DQmJ0lI0LlHt9wCjz0GffrA4cOeTStX1lM78uWDhg31VqlRFpzG5DRxcbpQ/e674eWXoUsXOHDAs2m5cjowX6oUNGsG06f7XGuUsuA0JicS0YXqjz2mQ+lt2gScwFm6tPY8L7wQ2raFt97yt9RoZMFpTE52++0wdqyeZ9SsmS7V9FC0qK7krFULunaFF1/0t8xokxDpAowxEdarl64w6tpVb2bOnq0na56gQAE9BK5jR7jxRs3YSy7xv9xoYD1OY4ym4QcfwNq1ULcu/PKLZ7M8eeD99/W26JAh8NJLZ+fIbeksOI0xqnlzvWTfvFlnv6/1PnQ2MREmTNBe51tvleGmmwIOzGdbFpzGmP+pU0fP2Ni3T3uey5d7NouP1wVIV1/9My++CN26wcGDvlYaURacxph/qlJFh9Fz59ZF659/7tlMBHr3/pGHHtITi6+8Evbu9bXSiLHgNMb8W/nyOoEzJUXXYM4OfATY4MHwwgu6orN5c9i508c6I8SC0xjjrUwZ7XmWL68bdb7zTsCmN9wAEyfCF1/owPzWrT7WGQEWnMaYwIoX13uel1wCnTvrSqMAOneGqVPhm290ffuGDb5V6TsLTmPMyRUsqJfql10GvXvraqMAWrXSphs36jjTunU+1ukjC05jzKnlzQvTpsFVV8GgQbq7UoAJnPXq6Ukdu3frwPzKlT7X6gMLTmNMxiQl6Y3MXr10P89bbgl4rnBqqt4eTUiA+vVh0SKfaw0zC05jTMbFx8NLL2mvc8wYzn/ooYATOC+4QAfmixaFJk1gzhyfaw0jC05jTHBE4OGHYeRISsyZA+3b64R5D2XLas/z3HOhdWs9STM7sOA0xgRPBIYO5buBA/WQopYtYdcuz6YlSujAfGqqLokfP97PQsPDgtMYk2m/XXGFLlxfsAAaN4Y///RsV7iwXqo3bqynaT71lL91hpoFpzEma7p21S2TVq7UIfWNGz2b5csHH36oV/a33gr33htwYD7qWXAaY7Lu8st1s85fftE5SN9/79ksVy5d196jB9x3HwwcGHBgPqpZcBpjQqNBA90mfscOnf2+apVns4QEXYB066169FHPnnqaZiyx4DTGhE716nq/U0QncH75pWezuDh4/HHtdb72ms6r37/f51qzwILTGBNaF16oEzgLF9bRoE8+8WwmAsOHw+jReou0VStdbRQLLDiNMaF39tkanmedpVOVpk0L2LR/f+11pqXpDnbbtvlXZmZZcBpjwqNkSfj0U6haVYfSJ0wI2LR7d5gyBZYu1Sv8TZt8rDMTLDiNMeFTpAjMnatp2L07PPtswKZt2+pmyD/+qAPzP/3kW5VBs+A0xoRXcjJMnw5XXAH9+sGDDwacwHn0lui2bVC7tu7tGY0sOI0x4Zc7t16Ld+8Od98Nd9wRMDwvuUSv8I8c0fn0ixf7XGsGWHAaY/yRkACvvqq9ziee0O3pApwrXKmSji3lzw+NGunAUTSx4DTG+CcuTucfDR8Or7wCnToFnMB5zjkanqVL6yFwH33kc60nYcFpjPGXiM58f/JJePddaNMG9uzxbFqqlM6nr1RJjx+eONHnWgOw4DTGRMbAgdrrnDtXzzP66y/PZkWL6oBRnTrQrRs8/7y/ZXqx4DTGRE6PHvD22/DVV7rWffNmz2YFCuhUpdat4eab4aGH/C3zRBacxpjIat9eb2CuX68TOH/+2bNZnjx6Zd+1KwwdCoMHR25bOgtOY0zkXXaZ7nS8datek3/7rWezxER4/XW46SZ45BG48caAA/NhZcFpjIkOtWrpBM6DB7XnuXSpZ7O4OHjuOe11vvSS9kAPHPC3VAtOY0z0qFxZT3fLmxcaNtSvPYjoCcWPPAKTJ+tyzb//9q9MC05jTHQpV04ncJYsqZfwM2YEbDpokPY6Z82CZs10D2U/WHAaY6LPGWdob7NCBV3jPnlywKa9e+txHF9+qZ3UrVvDX54FpzEmOhUrpkdx1KwJXbpo1zKAq66CDz7QMaW6deHXX8NbmgWnMSZ6FSyo1+EtWsANN+hNzQCaN4ePP9a9POvUgXXrwldW1ASniDwqIt+KyEoReV9ECkW6JmNMFMibV8/W6NxZJ2/edVfACZx16sD8+bB3r369fHl4Soqa4ATmABWdc5WB74ChEa7HGBMtkpLgjTe01zlqFPTtG/Bc4Ysv1tujSUm6GOnzz0NfTtQEp3PuY+fc0UNC/w8oHcl6jDFRJj5eF6oPHqy/XnONzvn0UL68Bmbx4jow/9VXhUNaStQE5wl6AjMjXYQxJsqI6EL1UaN0q6R27fS63EOZMtrzLFcOhg2rxIYNISzD+bjYU0TmAiU8fjTMOTctvc0woBrQzgUoTkT6AH0AUlJSUidNmhRUHbt37yY5OTmo10SC1RlasVInxE6tkazz9A8+oNxTT7GjcmVWjRjB4Xz5PNvt3p3A55/noVmzXUG9f8OGDZc456p5/tA5FzUP4FpgEZA3o69JTU11wZo/f37Qr4kEqzO0YqVO52Kn1ojXOXGicwkJzlWr5tzWrQGbZaZOYLELkDtRc6kuIs2BwUAb55yPi6eMMTGrSxeYOhVWr9YDijZu9OVjoyY4gWeB/MAcEVkuIi9EuiBjTAxo1Urnem7YoHOQ1q8P+0dGTXA65851zp3hnKuS/rgx0jUZY2JE/fo6gXPXLg3PlSvD+nFRE5zGGJMlqak6jJ6QoEG6aFHYPsqC0xiTfVxwge6sVLQoNGmi5xmFgQWnMSZ7KVtWe57nnKP3P99/P+QfYcFpjMl+SpTQ3eRTU6FDB1JmzQrp21twGmOyp8KFdbukRo0o/8QThHLpkAWnMSb7Sk6Gjz5i+eOPQ+nQbX9hwWmMyd5y5WJnpUohfUsLTmOMCZIFpzHGBMmC0xhjgmTBaYwxQbLgNMaYIFlwGmNMkCw4jTEmSBacxhgTJAtOY4wJkgWnMcYEyddTLsNBRLYCPwf5sqLAH2EoJ9SsztCKlTohdmrNznWe6Zwr5vWDmA/OzBCRxS7QsZ9RxOoMrVipE2Kn1pxap12qG2NMkCw4jTEmSDk1OF+KdAEZZHWGVqzUCbFTa46sM0fe4zTGmKzIqT1OY4zJtBwZnCLSUUS+FpEjIhJ1I4Ii0lxE1orIehEZEul6AhGRV0Rki4isjnQtJyMiZ4jIfBFZk/7nPiDSNXkRkdwi8l8RWZFe532RrulkRCReRJaJyEeRruVkROQnEVklIstFZHEo3jNHBiewGmgHLIh0IScSkXjgOaAFUAHoIiIVIltVQOOB5pEuIgMOAbc75y4ALgX6Rul/0/1AI+fcRUAVoLmIXBrZkk5qALAm0kVkUEPnXJVQTUnKkcHpnFvjnFsb6ToCqAGsd8794Jw7AEwCrohwTZ6ccwuAbZGu41Scc5ucc0vTv96F/mMvFdmq/s2p3enfJqY/onIQQkRKA62AlyNdSyTkyOCMcqWAX4/7fgNR+I88VolIWaAq8GWES/GUfvm7HNgCzHHORWWdwFPAncCRCNeREQ74WESWiEifULxhQijeJBqJyFyghMePhjnnpvldTxDE47mo7HXEGhFJBt4FBjrndka6Hi/OucNAFREpBLwvIhWdc1F1D1lEWgNbnHNLRKRBhMvJiNrOud9EpDgwR0S+Tb9ayrRsG5zOuSaRriGTNgBnHPd9aeC3CNWSbYhIIhqabzrn3ot0PafinNsuImnoPeSoCk6gNtBGRFoCuYECIvKGc65bhOvy5Jz7Lf3XLSLyPno7LEvBaZfq0ecroJyInCUiSUBn4IMI1xTTRESAccAa59wTka4nEBEplt7TRETyAE2AbyNalAfn3FDnXGnnXFn07+e8aA1NEcknIvmPfg1cRgj+R5Qjg1NErhSRDUBNYLqIzI50TUc55w4BtwCz0UGMt51zX0e2Km8i8hawCCgvIhtE5PpI1xRAbeAaoFH6lJTl6b2laFMSmC8iK9H/gc5xzkX1VJ8YkAIsFJEVwH+B6c65WVl9U1s5ZIwxQcqRPU5jjMkKC05jjAmSBacxxgTJgtMYY4JkwWmMMUGy4DTGmCBZcBpjTJAsOI0xJkgWnCZbS9+0er+InHncc6NF5HsRSYlkbSZ22cohk62lr1P/CljmnOstIneg26HVds6ti2x1JlZl292RjAHdHFhE7kL3JPgeGIbusm6haTLNepwmRxCRL9DtxC53zs2MdD0mttk9TpPtiUgj4CJ0k+jNES7HZAPW4zTZmohcBHwK3IaekZPsnGsW2apMrLPgNNlW+kj6F8CLzrn7RaQisBK9x5kW0eJMTLPgNNmSiBQBPgcWOOduOO75yUAZ51zNiBVnYp4FpzHGBMkGh4wxJkgWnMYYEyQLTmOMCZIFpzHGBMmC0xhjgmTBaYwxQbLgNMaYIFlwGmNMkCw4jTEmSP8PHXSI67uIlGwAAAAASUVORK5CYII=\n", "text/plain": [ "
\n", "\n", "\n", "[The point here is that the minimiser of $ (4 - \\beta)^2 + (7-2\\beta)^2$ is not the same as the minimiser of $ (4 - \\beta)^2 + (3.5-\\beta)^2$! So our initial intuition in solution 1 was technically wrong. The former is effectively weighting the misfit in our satisfaction of the second equation more than the first.]\n", "\n", "\n", "
\n", "\n", "We can then repeat the solution above but starting from $2x+3y=3.6$ ...\n", "\n", "We can parametrise this line by introducing the arbitrary parameter $\\alpha\\in\\mathbb{R}$ and setting $x=\\alpha$, $y = (3.6 - 2\\alpha)/3$.\n", "\n", "To find the solution in this infinite family introduce the function of $\\alpha$\n", "\n", "\\begin{align}\n", "f(\\alpha) &:= \\alpha^2 + \\left((3.6 - 2\\alpha)/3\\right)^2\\\\\n", "&= \\alpha^2 + \\frac{4}{9}\\alpha^2 - \\frac{2\\times 2\\times 3.6}{9}\\alpha + 3.6^2 \\\\\n", "& = \\left(1+\\frac{4}{9}\\right) \\alpha^2 - \\frac{4\\times 3.6}{9}\\alpha + 3.6^2\n", "\\end{align}\n", "\n", "differentiate w.r.t. $\\alpha$:\n", "\n", "$$\n", "\\frac{df}{d\\alpha} = 2\\left(1+\\frac{4}{9}\\right) \\alpha - \\frac{4\\times 3.6}{9}\n", "$$\n", "\n", "and so $\\alpha = ...$ and our \"solution\" is ... [let's get some help with this:]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.55384615 0.88076923]\n", "1.0825\n", "3.75\n" ] }, { "data": { "text/plain": [ "[
\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": [ "### Solution\n", "\n", "
![](\"./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": [ "### Solution\n", "\n", "Recall that the null space of the sensitivity matrix $G$ is defined as the set of all column vectors $\\boldsymbol{v}$ such that $G\\boldsymbol{v}= \\boldsymbol{0}$. \n", "\n", "Let's start by writing our matrix $G$ as a numpy array and compute its rank." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[1, 1, 1, 0, 0, 0, 0, 0, 0],\n", " [0, 0, 0, 1, 1, 1, 0, 0, 0],\n", " [0, 0, 0, 0, 0, 0, 1, 1, 1],\n", " [1, 0, 0, 1, 0, 0, 1, 0, 0],\n", " [0, 1, 0, 0, 1, 0, 0, 1, 0],\n", " [0, 0, 1, 0, 0, 1, 0, 0, 1]])\n", "5\n" ] } ], "source": [ "# use numpy to work out the rank of G\n", "G = np.array([\n", " [1,1,1,0,0,0,0,0,0],\n", " [0,0,0,1,1,1,0,0,0],\n", " [0,0,0,0,0,0,1,1,1],\n", " [1,0,0,1,0,0,1,0,0],\n", " [0,1,0,0,1,0,0,1,0],\n", " [0,0,1,0,0,1,0,0,1]])\n", "\n", "pprint(G)\n", "\n", "print(np.linalg.matrix_rank(G))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By the rank-nullity theorem if the rank of $A$ is 5 we expect the dimention of the null space to be 4 (as $n=9$).\n", "\n", "Based on the two possible familes of solution given above (after rearrangement into $3\\times 3$) being\n", "\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", "\n", "you can hopefully see that we can spot that two possible vectors in the null space are\n", "\n", "$$\\boldsymbol{n}_1 = (1, \\; -1,\\; 0, \\;-1, \\;1, \\;0, \\;0, \\;0, \\;0)^T$$ \n", "\n", "$$\\boldsymbol{n}_1 = (0, \\;0, \\;0,\\; 0,\\; 1,\\; -1,\\; 0,\\; -1,\\; 1)^T$$ \n", "\n", "which we can check by multiplying by $G$:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([0, 0, 0, 0, 0, 0])\n", "array([0, 0, 0, 0, 0, 0])\n" ] } ], "source": [ "n_1 = np.array([1, -1, 0, -1, 1, 0, 0, 0, 0])\n", "pprint(G@n_1)\n", "\n", "n_2 = np.array([0, 0, 0, 0, 1, -1, 0, -1, 1])\n", "pprint(G@n_2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and with a bit of Sudoko like thinking we can come up with another two:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([0, 0, 0, 0, 0, 0])\n", "array([0, 0, 0, 0, 0, 0])\n" ] } ], "source": [ "n_3 = np.array([0, 1, -1, 0, -1, 1, 0, 0, 0])\n", "pprint(G@n_3)\n", "\n", "n_4 = np.array([0, 0, 0, 1, -1, 0, -1, 1, 0])\n", "pprint(G@n_4)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([[ 1, 0, 0, 0],\n", " [-1, 0, 1, 0],\n", " [ 0, 0, -1, 0],\n", " [-1, 0, 0, 1],\n", " [ 1, 1, -1, -1],\n", " [ 0, -1, 1, 0],\n", " [ 0, 0, 0, -1],\n", " [ 0, -1, 0, 1],\n", " [ 0, 1, 0, 0]])\n", "4\n" ] } ], "source": [ "# form a 9x4 matrix from these null vectors\n", "null_vecs = np.vstack((n_1,n_2,n_3,n_4)).T\n", "pprint(null_vecs)\n", "# and check its rank\n", "pprint(np.linalg.matrix_rank(null_vecs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The rank is 4 and so these 4 vectors are linearly independent, they are all in the null space of $G$, and hence they form a basis for the null space for this problem.\n", "\n", "i.e. given a particular solution such as \n", "\n", "$$\n", "\\begin{pmatrix}\n", "2 & 2 & 2 \\\\\n", "2 & 2 & 2 \\\\\n", "2 & 2 & 2\n", "\\end{pmatrix}\n", "$$\n", "\n", "we can add on any linear combination of the following matrices\n", "\n", "$$\n", "\\begin{pmatrix}\n", "1 & -1 & 0 \\\\\n", "-1 & 1 & 0 \\\\\n", "0 & 0 & 0 \n", "\\end{pmatrix},\\quad\n", "\\begin{pmatrix}\n", "0 & 1 & -1 \\\\\n", "0 & -1 & 1 \\\\\n", "0 & 0 & 0 \n", "\\end{pmatrix},\\quad\n", "\\begin{pmatrix}\n", "0 & 0 & 0 \\\\\n", "0 & 1 & -1 \\\\\n", "0 & -1 & 1 \n", "\\end{pmatrix},\\quad\n", "\\begin{pmatrix}\n", "0 & 0 & 0 \\\\\n", "1 & -1 & 0 \\\\\n", "-1 & 1 & 0 \n", "\\end{pmatrix}\n", "$$ \n", "\n", "and the result will always satisfy our problem.\n", "\n", "\n", "\n", "\n", "That is, any model in the null space of $G$ can be written as a linear combination of these four models. The significance of the null space is that, if $\\boldsymbol{v}$ is any vector in the null space, $\\alpha$ is any scalar, and $\\boldsymbol{m}$ is any model that satisfies $G\\boldsymbol{m} = \\boldsymbol{d}$, then the model parameters $\\boldsymbol{m} + \\alpha \\boldsymbol{v}$ also satisfies the observed data $\\boldsymbol{d}$.\n", "\n", "In other words, we can add any linear combination of vectors from the null space to a model that satisfies the data, and still satisfy the data. \n", "\n", "For all under-determined and mixed-determined problems with a solution (i.e. representing a consistent problem), the matrix $G$ will have a non-trivial null space, that is there will be solutions to $G\\boldsymbol{v}= 0$ for which $\\boldsymbol{v}\\ne \\boldsymbol{0}$. \n", "\n", "Note that the null space is purely a property of $G$, that is it is a property of the physics and the experimental geometry, and it does not depend upon either the data or the \"model\" (here the solution vector).\n" ] }, { "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": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0\n", "-7.6438855245442e-14\n" ] } ], "source": [ "G = np.array([\n", " [1,1,1,0,0,0,0,0,0],\n", " [0,0,0,1,1,1,0,0,0],\n", " [0,0,0,0,0,0,1,1,1],\n", " [1,0,0,1,0,0,1,0,0],\n", " [0,1,0,0,1,0,0,1,0],\n", " [0,0,1,0,0,1,0,0,1]])\n", "\n", "print(sl.det(G.T@G))\n", "print(sl.det(G@G.T))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are both zero so we can't use our least squares or minimum norm solutions.\n", "\n", "Let's use the pseudo inverse" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([2.01111111, 2.01111111, 2.04444444, 2.01111111, 2.01111111,\n", " 2.04444444, 1.91111111, 1.91111111, 1.94444444])\n" ] } ], "source": [ "G = np.array([\n", " [1,1,1,0,0,0,0,0,0],\n", " [0,0,0,1,1,1,0,0,0],\n", " [0,0,0,0,0,0,1,1,1],\n", " [1,0,0,1,0,0,1,0,0],\n", " [0,1,0,0,1,0,0,1,0],\n", " [0,0,1,0,0,1,0,0,1]])\n", "d = np.array([6.07,6.07,5.77,5.93,5.93,6.03])\n", "\n", "\n", "m = np.linalg.pinv(G) @ d\n", "\n", "pprint(m)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Some investigation will show that this solution both provides a least-squares fit to the data – \n", "that is, no other solution can provide a better fit.\n", "\n", "It also provides a minimal-norm solution - that is, no other least-squares solution has a smaller norm than this solution. \n", "\n", "Neither of these properties mean of course that the solution is correct, but in many problems it may be the most appropriate solution that we can discover without additional information. " ] }, { "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solution" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "array([2.01111111, 2.01111111, 2.04444444, 2.01111111, 2.01111111,\n", " 2.04444444, 1.91111111, 1.91111111, 1.94444444])\n", "\n", "mu= 1.0 \n", "\n", "array([1.72142857, 1.72142857, 1.74642857, 1.72142857, 1.72142857,\n", " 1.74642857, 1.64642857, 1.64642857, 1.67142857])\n", "\n", "mu= 0.1 \n", "\n", "array([1.97778953, 1.97778953, 2.01004759, 1.97778953, 1.97778953,\n", " 2.01004759, 1.88101534, 1.88101534, 1.9132734 ])\n", "\n", "mu= 0.010000000000000002 \n", "\n", "array([2.00772798, 2.00772798, 2.04095058, 2.00772798, 2.00772798,\n", " 2.04095058, 1.90806021, 1.90806021, 1.9412828 ])\n", "\n", "mu= 0.0010000000000000002 \n", "\n", "array([2.01077228, 2.01077228, 2.04409451, 2.01077228, 2.01077228,\n", " 2.04409451, 1.9108056 , 1.9108056 , 1.94412783])\n", "\n", "mu= 0.00010000000000000002 \n", "\n", "array([2.01107722, 2.01107722, 2.04440945, 2.01107722, 2.01107722,\n", " 2.04440945, 1.91108056, 1.91108056, 1.94441278])\n", "\n", "mu= 1.0000000000000003e-05 \n", "\n", "array([2.01110772, 2.01110772, 2.04444094, 2.01110772, 2.01110772,\n", " 2.04444094, 1.91110806, 1.91110806, 1.94444128])\n" ] } ], "source": [ "G = np.array([\n", " [1,1,1,0,0,0,0,0,0],\n", " [0,0,0,1,1,1,0,0,0],\n", " [0,0,0,0,0,0,1,1,1],\n", " [1,0,0,1,0,0,1,0,0],\n", " [0,1,0,0,1,0,0,1,0],\n", " [0,0,1,0,0,1,0,0,1]])\n", "d = np.array([6.07,6.07,5.77,5.93,5.93,6.03])\n", "\n", "# use the pinv to compute the solution\n", "m = np.linalg.pinv(G) @ d\n", "\n", "# print the solution obtained using the generalised inverse\n", "pprint(m)\n", "\n", "# now apply our \"least squqres\" analytical approach with \n", "# successively smaller regularisation parameters\n", "for k in range(6):\n", " mu = 0.1**k\n", " m = sl.inv( G.T@G + mu*np.eye(G.shape[1]) ) @ G.T @ d\n", " print('\\nmu=',mu,'\\n')\n", " pprint(m)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We do indeed find that as $\\mu$ gets smaller the solution to the damped least squares problem converges to that obtained using the generalised inverse." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }