{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "using LinearAlgebra" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Errors, Norms, and Condition Numbers\n", "\n", "(See also section 9.2 of Strang, *Introduction to Linear Algebra*. The book *Numerical Linear Algebra* by Trefethen and Bau is a good source for a much more in-depth discussion.)\n", "\n", "Throughout the semester, when we solve problems like $Ax=b$ on the computer, we notice that the answer is correct \"up to roundoff errors\" — little errors that appear in the solution, typically $\\lesssim 10^{-15}$:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "5-element Vector{Float64}:\n", " 0.0\n", " 6.938893903907228e-18\n", " 0.0\n", " 0.0\n", " 1.1102230246251565e-16" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = randn(5,5)\n", "b = randn(5)\n", "x = A \\ b # solve Ax = b\n", "b - A*x # this \"residual\" should be zero" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These arise because the computer [only stores about 15 decimal digits](https://en.wikipedia.org/wiki/Floating-point_arithmetic) for each number, so almost every arithmetic operation makes a tiny [round-off error](https://en.wikipedia.org/wiki/Round-off_error).\n", "\n", "In real applications, there are lots of other sources of error, too. If $A$ or $b$ comes from measured data, they probably include [measurement errors](https://en.wikipedia.org/wiki/Observational_error). Often, the equations you solve on the computer are only approximations for the *real* equations of nature (which may not even be fully known!), so there could also be **systematic errors** or **modelling errors** in your problem.\n", "\n", "Such (hopefully) tiny errors are almost inevitable in most real problems. But the key question is **do tiny errors in the *inputs* produce equally tiny errors in the *outputs*,** or **do errors get amplified?**?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Output error from input error\n", "\n", "Suppose that we we want to solve $Ax = b$, but we actually have a *little error Δb in b*. Does that mean **we have an equally little error in x**?\n", "\n", "Let's work it out: we solve $A(x+\\Delta x) = b + \\Delta b$, and we should obtain\n", "$$\n", "x + \\Delta x = A^{-1} (b + \\Delta b) = A^{-1} b + A^{-1} \\Delta b\n", "$$\n", "so our error should be $\\boxed{\\Delta x = A^{-1} \\Delta b}$.\n", "\n", "Let's check:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 0.3333333333333346\n", " 0.3333333333333327" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = [1 2\n", " 2.01 3.99]\n", "b = [1,2]\n", "x = A \\ b # \"exact\" answer" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 1.930000000000043\n", " -0.47000000000002146" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Δb = [-0.01, 0.004]\n", "x′ = A \\ (b + Δb)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 1.5966666666667084\n", " -0.8033333333333541" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x′ - x # the error Δx" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 1.5966666666667066\n", " -0.8033333333333534" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Δx = A \\ Δb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, great, our $\\Delta x$ formula worked. But look at what happened here — we put a \"tiny\" error $\\sim 0.01$ into b, and got a \"huge\" error $\\sim 1$ in x! Why?\n", "\n", "That is, if we compare $\\Vert \\Delta x \\Vert$ to $\\Vert \\Delta b \\Vert$, we get a big increase:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.010770329614269008" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(Δb) # the size of the error in the input" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.787369264838424" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(Δx) # the size of the error in the output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Error growth and a matrix norm\n", "\n", "In general, what is the relationship between $\\Vert \\Delta x \\Vert$ and $\\Vert \\Delta b \\Vert$? If $\\Delta b$ is a \"random\" error that could be in *any direction*, how big can $\\Vert \\Delta x \\Vert / \\Vert \\Delta b \\Vert$ be?\n", "\n", "That is, we would like to know the **maximum possible value** of\n", "$$\n", "\\frac{\\Vert \\Delta x \\Vert}{\\Vert \\Delta b \\Vert} = \\frac{\\Vert A^{-1} \\Delta b \\Vert}{\\Vert \\Delta b \\Vert}\n", "$$\n", "over **all possible Δb ≠ 0**.\n", "\n", "More generally, for *any* matrix $B$, we define the [induced or \"operator\" matrix norm](https://en.wikipedia.org/wiki/Matrix_norm)\n", "$$\n", "\\Vert B \\Vert = \\max_{y\\ne 0} \\frac{\\Vert B y \\Vert}{\\Vert y \\Vert}\n", "$$\n", "This is a measure of \"how big\" the matrix is, according to the *maximum* amount by which it can \"stretch\" a vector.\n", "\n", "By this definition, $\\boxed{\\Vert \\Delta x \\Vert \\le \\Vert A^{-1} \\Vert \\; \\Vert \\Delta b \\Vert}$: the **norm of A⁻¹ *bounds* how much the error can increase**.\n", "\n", "In Julia, the induced norm $\\Vert B \\Vert$ is computed by `opnorm(B)`, for example:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4.996014806077393" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(A) # ‖A‖" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "166.53382686925062" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(inv(A)) # ‖A⁻¹‖" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fact that our $A^{-1}$ from above has such a big norm explains why $\\Vert \\Delta x \\Vert$ could be 100× bigger than $\\Vert \\Delta b \\Vert$.\n", "\n", "But how can we get $\\Vert A^{-1} \\Vert$?\n", "\n", "One way is to simply compute $\\Vert A^{-1} y \\Vert / \\Vert y \\Vert$. Since $\\Vert A^{-1} y \\Vert / \\Vert y \\Vert = \\Vert A^{-1} \\alpha y \\Vert / \\Vert \\alpha y \\Vert$ for any scalar α, we can freely restrict ourselves to $\\Vert y \\Vert = 1$ (by choosing $\\alpha = 1 / \\Vert y \\Vert$), i.e. we can *equivalently* define\n", "\n", "$$\n", "\\Vert B \\Vert = \\max_{\\Vert y \\Vert = 1} \\Vert B y \\Vert\n", "$$\n", "\n", "So, we can just plot $\\Vert A^{-1} y \\Vert$ versus angle $\\theta$ for $y = (\\cos \\theta, \\sin \\theta)$ on the unit circle:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAHMCAYAAADWN6wLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABoC0lEQVR4nO3dd3gUVdsG8HvTC0noKSSEIEUglNCCobcAUpQiKIog6KcCIiJ2KVYUAbug76tgAUQFFASV0DskgUhoghAISAISII30zPfHeWeTJYXsZndnZuf+XddemexOdp/d7M4+c85zzjFIkiSBiIiISKeclA6AiIiISElMhoiIiEjXmAwRERGRrjEZIiIiIl1jMkRERES6xmSIiIiIdI3JEBEREekakyEiIiLSNSZDREREpGtMhoiIiEjXmAwRERGRrjEZIiJdW7x4Mdq3bw9XV1fMnTtX6XCISAFMhohI1wIDA/Haa6/h3nvvVToUIlKIi9IBEBEpSU6CfvnlF2UDISLFsGWISGHLli2DwWDAuXPnFL1PW8RhjkaNGmHu3Lk4d+4cDAYDtm/fXuG+x44dg6urKwwGAy5fvmy/IG/DnOdAVWer92ZGRgZmzpyJkJAQeHh4oHPnzti3b59VH4O0gckQkcIGDx6Mffv2ITAwUNH7tEUctjJt2jQUFhYCABISEpQNhjTp6tWr6Nq1K3bs2IEPPvgAa9euRVFREYYMGYLr168rHR7ZGZMhIoXVq1cPXbp0gbu7u6L3aYs4bOGnn37C1q1bMXjwYAAVJ0N9+/aFh4dHuZdZs2bZMWJtuHnzptIh2NWECRMgSRJ27NiBkSNHYtCgQViyZAmuXbvGLlMdYjJEujN37lwYDAYcOXIE9913H/z8/FC7dm3MmDEDhYWF+OuvvzBw4ED4+PigUaNGmD9/vsnf//3333jkkUfQtGlTeHl5oUGDBhg6dCgSExON++Tm5iIiIgJNmjRBenq68frU1FQEBASgV69eKCoqAlC2C6C68ZV3nwaDocKLvE95XRFyLMeOHcMDDzwAPz8/+Pv7Y+LEiSbPCxA1N23atIG7uzsaN26MDz/80Pj31pKTk2Ps1vj666/h7OxcYTK0ZcsW5Obmlnt54403qvyYI0aMQHBwcJnrCwsL0a5dO/Tv39/SpwMA2L17N/r27QsfHx94eXkhKioKGzZsMN7+888/w2AwYMuWLWX+dvHixcb3iuz06dMYO3Ys6tevD3d3d7Ro0QKffvqpyd/J/5dDhw5h1KhRqFWrFu64445K46zK+770fdvjPVOV51qerVu3YsOGDVi0aBG8vLyM1zdu3BgAcPbs2So9PjkOFlCTbo0ePRoPPfQQHn/8ccTExGD+/PkoKCjA5s2bMXnyZMycORMrVqzACy+8gCZNmmDEiBEAgEuXLqFOnTp45513UK9ePVy7dg1ff/01IiMjcfjwYTRv3hweHh744Ycf0KFDB0ycOBGrV69GcXExHnzwQUiShJUrV8LZ2dkm8ZXn1jqInJwcjBs3DkVFRahdu/ZtX6uRI0dizJgxmDRpEhITE/HSSy8BAL766isAwO+//44RI0agR48eWLVqFQoLC7FgwQKr1/O8++67OH/+PL7//nvUqVMHTZs2rXY3WWFhIQoLC1FUVITCwkLk5ubC1dXV+P/p0aMH1q5di/PnzyM0NNT4d4sWLcLJkyfx448/WvzYO3bsQP/+/dGmTRt8+eWXcHd3x2effYahQ4di5cqVGDNmDIYMGYL69etj6dKl6Nu3r8nfL1u2DO3bt0ebNm0AAMePH0dUVBQaNmyIhQsXIiAgAH/88QemTZuGq1evYs6cOSZ/P2LECNx///144oknkJ2dXWmsVXnfl2br94y5z7W0L774Ao0aNULv3r2N3a2AqCECAFdX1yrFQA5EItKZOXPmSACkhQsXmlzfrl07CYC0Zs0a43UFBQVSvXr1pBEjRlR4f4WFhVJ+fr7UtGlT6ZlnnjG5bdWqVRIA6YMPPpBmz54tOTk5SZs2bTLZZ+nSpRIAKSkpyWrx3Xqft8Z7zz33SDVq1JDi4+Mr/Rs5lvnz55vcx+TJkyUPDw+puLhYkiRJ6tSpkxQSEiLl5eUZ98nMzJTq1KkjVfUwExoaKs2ZM0dKSkqSAEjbtm0zuf38+fOSp6en1LNnT+N1o0ePlpycnKTs7OwqPUZ55OdY+rJ06VLj7fHx8RIAacWKFcbrzp49K3l5eUmvv/66Wc/hVl26dJHq168vZWZmGq8rLCyUwsPDpeDgYOPrO2PGDMnT01O6ceOGcb/jx49LAKSPP/7YeN2AAQOk4OBgKT093eRxpk6dKnl4eEjXrl0zec6zZ8+u2otUjore97Z4z5T33qzqc71VUVGRVLNmzTL/89KXb7/91qLXhLSL3WSkW0OGDDH5vUWLFjAYDBg0aJDxOhcXFzRp0gTnz583XldYWIi3334bLVu2hJubG1xcXODm5obTp0/jxIkTJvc5evRoPPnkk3juuefw5ptv4uWXX65yt4ql8d3O1KlTsWHDBvz4449o3759lf5m2LBhJr+3adMGubm5uHLlCrKzsxEXF4d7770Xbm5uxn1q1KiBoUOHVjmu25kxYwby8/Px0UcfmcRRXFxcpqvGHHPnzoUkSSaXCRMmGG9v164dfH19sWfPHuN1Tz75JEJCQvDCCy9Y/LjZ2dk4cOAARo0ahRo1ahivd3Z2xrhx43Dx4kX89ddfAICJEyciJycHq1atMu63dOlSuLu7Y+zYsQBE1+yWLVswfPhweHl5GVu8CgsLcffddyM3Nxf79+83iWHkyJFVjtec9z1g2/eMJc9V9tdff+HGjRt44403EBsba3IZP348AKBz585Vfl3IMTAZIt26tXvIzc0NXl5e8PDwKHN9bm6u8fcZM2Zg1qxZuPfee7F+/XocOHAAsbGxaNu2LXJycso8zsSJE1FQUAAXFxdMmzbN5vFV5s0338SSJUvw+eefY+DAgVWOpU6dOia/y0XWOTk5uH79OiRJgr+/f5m/K+86S2zbtg2rV6/GQw89hIYNG+LGjRu4ceOGscbDliPKnJycEBUVhb179wIAli9fjj/++ANLliwx+SI3l/y6lTd6LygoCACQlpYGAGjVqhU6deqEpUuXAgCKiorw3Xff4Z577jG+T9LS0lBYWIiPP/4Yrq6uJpe7774bgBhBVZo5IwfNfd/b8j1jyXOVyTVxkZGR6Nixo8nl6NGjuOOOO9CsWbPbxkCOhTVDRGb67rvv8PDDD+Ptt982uf7q1auoWbOmyXXZ2dkYN24cmjVrhsuXL+PRRx9VbKTKsmXLMGvWLMydOxcTJ0602v3WqlWrwvl+UlNTq33/RUVFxiTy66+/xtdff11mH1sPr+/RowdmzZqF5ORkzJgxA+PHj0evXr2qdZ+1atWCk5MTUlJSytx26dIlAEDdunWN1z3yyCOYPHkyTpw4gbNnzyIlJQWPPPKIyf3JrUpTpkwp9zHDwsJMfjenuN2c9/3tVPc9Y8lzlRUUFABAmZq9hIQExMfH4913373t45PjYTJEZCaDwVBm+PmGDRvwzz//oEmTJibXP/HEE0hOTsbBgwdx8uRJjBo1Cu+//z6eeeYZe4aM33//HY899hgmTpxYaWGpJby9vdGxY0f8/PPPWLBggbG1JCsrC7/++mu17//TTz/F0aNH8dprr6FHjx5lbh85cqRdkiF5DpqioiIsWLCg2vfp7e2NyMhIrFmzBgsWLICnpycAoLi4GN999x2Cg4NNWigeeOABzJgxA8uWLcPZs2fRoEEDREdHG2/38vJC7969cfjwYbRp06ZarVblMed9fzvVfc9U57nKrYmJiYno06cPANEFOG3aNISFhWHq1KlmPRdyDEyGiMw0ZMgQLFu2DHfeeSfatGmD+Ph4vPfee2WGX//3v//Fd999h6VLl6JVq1Zo1aoVpk6dihdeeAFdu3a1W11CUlIS7rvvPjRu3BiPPPJImVqKiIiIas8t9Prrr2Pw4MEYMGAAnn76aRQVFeG9995DjRo1cO3aNYvv9+rVq5g7dy6ioqIwa9asclsy2rZti4MHD6K4uBhOTrbp+e/UqRM8PT2RmJiIr776yqTFpjrmzZuH/v37o3fv3pg5cybc3Nzw2Wef4ejRo1i5cqXJ861ZsyaGDx+OZcuW4caNG5g5c2aZ5/vhhx+iW7du6N69O5588kk0atQImZmZ+Pvvv7F+/Xps3brV4lir+r6vquq+Zyx9ruHh4ejQoQPefPNN+Pv7w8/PDwsXLsTx48exZcsWk6H2pB9MhojM9OGHH8LV1RXz5s1DVlYW2rdvjzVr1uDVV1817pOYmIhp06Zh/PjxJsW4CxYswL59+zBmzBgcPnzY7O4FS5w/fx5ZWVk4deoUunfvXub2pKQkNGrUqFqPMXDgQKxevRqzZ8/GmDFjEBAQgMmTJ+PSpUv49ttvLb7fV155BZmZmfj8888r7NJp27Yttm3bhtOnT5cZ3m0tTk5OqFWrFjp27Gjy/6yunj17YuvWrZgzZw4mTJiA4uJitG3bFuvWrStTQA+IrrKVK1cCQLlxtGzZEocOHcIbb7yBV199FVeuXEHNmjXRtGlTYy2NparyvjdHdd8z1Xmuq1evxuOPP45JkybBw8MDgwYNQlxcXLU/B6RdBkmSJKWDICLHU1BQgHbt2qFBgwbYtGnTbfdv1KgRJkyYgAkTJiAsLAzbtm2rdl2OtSxYsACvvPIKEhIS0KJFiwr3U/Nz0AJz3zNE1sKWISKyikmTJqF///4IDAxEamoqlixZghMnTuDDDz9UOjSL3Lx5E3/++SdiY2Pxyiuv4K233qo0ESLzOdp7hrSLyRARWUVmZiZmzpyJf//9F66urmjfvj02btyIfv36KR2aRTZt2oThw4cjICAAL7/8MmbOnKl0SA7H0d4zpF1MhojIKn744QelQ7Cqe++9F6wisC1He8+QdrFmiIiIiHSNM1ATERGRrjEZIiIiIl1jzVAVFBcX49KlS/Dx8TFr+noiIiJSjiRJyMzMRFBQUKWTsjIZqoJLly4hJCRE6TCIiIjIAhcuXKh0tnQmQ1Xg4+MDQLyYvr6+CkdDREREVZGRkYGQkBDj93hFmAxVgdw15uvry2SIiIhIY25X4sICaiIiItI1JkNERESka0yGiIiISNeYDBEREZGuMRkiIiIiXWMyRERERLrGZIiIiIh0jckQERER6RqTISIiItI1JkNERESka0yGiIiISNeYDBEREZGucaFWhWVnZ1d4m7OzMzw8PKq0r5OTEzw9PS3a9+bNm5Akqdx9DQYDvLy8LNo3JycHxcXFFcbh7e1t0b65ubkoKiqyyr5eXl7GBfzy8vJQWFholX09PT3h5CTONfLz81FQUGCVfT08PODs7Gz2vgUFBcjPz69wX3d3d7i4uJi9b2FhIfLy8irc183NDa6urmX3LS4GbtwA/PyA/8VYet+ioiLk5uZWeL+urq5wc3Mze9/i4mLk5ORYZV8XFxe4u7sDACRJws2bN62yrzmfex4jyt+XxwgzjxH5+SgoKkJ+Ja+vPY4RipLottLT0yUAUnp6utXvG0CFl7vvvttkXy8vrwr37dmzp8m+devWrXDfjh07muwbGhpa4b4tW7Y02bdly5YV7hsaGmqyb8eOHSvct27duib79uzZs8J9vby8TPa9++67K33dShs1alSl+2ZlZRn3HT9+fKX7Xrlyxbjv5MmTK903KSnJuO/MmTMr3ffo0aPGfefMmVPpvgcPHjTuO3/+/Er33bZtm3HfTz75pNJ9f/31V+O+S5curXTfH374wbjvDz/8UOm+S5cuFTteuCD9+uKLle77SYsWkvTqq5J05oy0bdu2SvedP3++MYaDBw9Wuu+cOXOM+x49erTSfWfOnGncNykpqdJ9J0+ebNz3ypUrle47fvx4475ZWVmV7jtq1CiT93Bl+/IYIS48RpRczDpG+PlJkru7JAHSJ66ule5r82OEjVT1+5stQ0RkOykpwEMPAStXitagypw4Abz5pri0b2+f+Ij0LD29ZLuSFiQ9MEhSBe2ZZJSRkQE/Pz+kp6fD19fXqvfNJnDz92UTuAa6yf7+G3jnHbitWgXX/71fCtu1Q1779iLR6dABaNJEdJVduQJcvgy31FS4/vwzEBODIklCLgDUrg289howYQLwv9ceYDeZJfvyGGHZvpo9RhQVIf/991EwezYgv4+dnIC+fYFRo4CwMHjUqwfnWrUAX18UJCcjf/t2YO9eYNcucSIDAG5uwBNPwP3ll+Hi7w9AW91kVf3+ZjJUBbZMhogcSn4+8PzzwMcfl7QEDR0KzJ1b9daec+eAr74Sl3/+EdfdfTfw3/8CgYG2iJrIsZw4AUyaBOzbJ37v0AF45BFg9GigXr3b/70kAQcPAi++CGzfLq7z8wNefhl49lljnZ8WMBmyIiZDRFWQmgrcdx+we7f4fcgQkQR16GDZ/RUVAR98IA7A+fmilWjJEvEYRFRWQQEwfz7w+uviM+PjA7z3HvDYY6JVyFySBPzxB/DCC8CRI+K6wYOBFSsAjXwXVvX7m0Prif5n+PDhqFWrFkaNGqV0KNpz4IBIenbvFgfJdeuA9estT4QAcfb57LPAoUNARARw7Zo4sx0/XhzoiahEdrZIVF59VXw+7r4bOHYMePxxyxIhQHRNDxwIHD4sWmY9PIANG4CoKCApybrxK4zJENH/TJs2Dd98843SYWjPl18CPXoAly4BLVoAsbGia8xaWrUC9u8HZs0SCdI334iah0rqEIh05cYNIDoaiIkBvL2Bb78Ffv0VCAmxzv07OYlut507RVf1sWNA586itshBMBkiTdi+fTsaNWoEQBRhnjt3zuqP0bt3b/j4+Cj2+Jr03nvAo4+KM9Hhw0ULUbNm1n8cNzfR9L9xozg7Xb8eGDECqKRwmkgX/v0X6N1bFD7XrAls3ixGcJYacGA1nTqJk5327YGrV0Ux9ooV1n8cBTAZIodx+PBhuLq6onv37kqHog9Ll4piaQCYPRv46SdRo2BL0dHijNfTUyRG99xTMlKGSG8uXhStsgkJQP36oti5SxfbPmaDBqJFaNQoUaP08MPis6hxTIbIYUybNg0zZ87En3/+We7Q3g4dOiA8PLzM5dKlSwpEq3G//CJahADguefE8HdL6xLM1bcv8Ntvojtg0yZRqF3JUHUih5SaKhKhkydFd9iuXUDbtvZ5bC8vYNUqkQgVFYlBDXFx9nlsG2EyRA5hxYoVqFWrFqZMmYLMzEycPXu2zD7x8fE4evRomUtQUJACEWvYjh3AmDFi6PwjjwDvvmv/GHr2BH7/HahRA9i6FRg79vaTOhI5irw8YORIUcR8xx0iEbJF93RlnJyA//wH6N9fnIwMHgycOWPfGKyIyRBpXnZ2Nl5++WW8++67CA4Ohp+fHxISEpQOyzEdPgwMGyYOxsOGAV98YZvahKro1k0M+3V3Fy1V8+YpEweRPUkS8NRTokbIz090UYWGKhOLmxuwejXQrp2YPHXgQFHDpEFMhkjz3nrrLQwcOBAtWrQAALRs2dKiZGjAgAG47777sHHjRgQHByM2NtbKkWrc1auiSyojQzTPf/894KLwij5RUcCnn4rtWbNEckTkyJYsES0yTk7iM2jvFqFb+fiUJGR//63ZbmuuTUaadvbsWXzxxRc4evSo8brw8HCLkqE/+EVaMUkCJk4Uw+ebNxctMaWWa1DUpElittwvvgAeeACIjwfCwpSOisj6du4Epk0T2/PmiZYYNQgMFN3WXbuKz+LMmcBnnykdlVnYMkSa9swzzyAtLQ3BwcFwcXGBi4sL/vvf/7KbzNo++0wMZ3dzE2ejNWsqHZGpjz4S855cvy6G3GvwzJSoUsnJYgRXYaFI+p97TumITN15pzg2AMDixWLUp4YwGSLNiomJwZ49e3D48GEkJCQYL19++SUuXryItLQ0pUN0DEeOiJmgATGvULt2ioZTLnd3MbS/Xj0xzPjJJ0VrFpEjKCoC7r9f1ONERIjZoJWq1atM//7AM8+I7YkTgcuXlY3HDEyGSJMKCwvx9NNP47nnnkO7du1Mhsr37dsXANg6ZA03b4qDcF6eGC3y1FNKR1SxkBAx3NfJScxS/dNPSkdEZB0ffigWXfXxAdasEUPb1ertt4HWrUXiNnGiZk5KmAyRJq1fvx5paWmYOnVqmdtCQkLg5eXFZMgannlGrIAdGCgmWVTj2WhpvXuLtZkAYOpUsZ4ZkZadPg288orYXrgQ+N9M+Krl4SFmpXZ3F4XVixcrHVGVMBkiTRo+fDguX74Mb2/vMrcZDAZkZ2fjWblrhyzzyy8lQ+e//VZ0QWnByy+LNdKuXCnp3iPSouJiMblpbq6YbFSe6FTtwsNL5h979llxQqVyTIaIqKzs7JIusZkzxYFYK9zdS2oqli0TazURadHixWIEmbe3GE6v9pbZ0p56Siyfk5sLjBsn6p5UjMkQEZX11lvAhQti7pC5c5WOxnxRUcCUKWL7//5PJHdEWnLuHPDCC2L7nXe0N12Ek5M4GfHzE9Nd/Pe/SkdUKSZDpAmNGjXC9OnTAQBz5sxBTTsP7Vb68e3qr7+ABQvE9ocfqrtYszJvvy2KqpOSgDlzlI6GqOokCXjsMZHEd+8OTJ6sdESWCQwEXn9dbL/yiqpr+AxSeStakomMjAz4+fkhPT0dvr6+SodDZDuSBAwYAMTEAHffLeYK0VLT/K02bhSj4JycgP37gU6dlI6I6PZWrAAefFAUIx85AjRtqnRElissFNMBHD0qkjp5xng7qer3N1uGiKjETz+JRMjdXUxkqOVECBAJnbyI6+TJXMyV1C83F3jpJbH96qvaToQAsWTPxx+L7SVLxDxgKsRkiIiErKySCdNefFGshu0IFi0Sq9vHxQE//qh0NESV+/BDMdt0cDAwY4bS0VhHr17A6NHiZOSpp1Q59xCTISIS3ngD+OcfUagpF246An9/4PnnxfbLLwP5+crGQ1SRf/8VtW6A+KmW9f+sYcECUX+4e7foBlQZJkNEJEauvP++2P7oI8c6CAPiDDsgADh7VjTVE6nR668DGRlA+/aiZsiRhISIkxFArKuWmalsPLdgMkREwGuvAQUFQL9+wJAhSkdjfd7e4jkC4gsnPV3ZeIhu9ddfJYn6ggWi6N/RPPss0LgxkJJSMmJVJRzw1SYis/z1l1jLCwDefFPZWGxp4kSxsnZaGjB/vtLREJl64QUx8mroULGsjCPy8CiZmfr991U11J7JEJHezZ0rChuHDgUiI5WOxnZcXEoOxIsWARcvKhsPkWzHDrH8jbNzyXvUUY0YAbRpI7rJFi5UOhojVSdDO3fuxNChQxEUFASDwYCff/7Z5PYJEybAYDCYXLp06WKyT15eHp566inUrVsX3t7eGDZsGC7yIEgkHDkCfP+92JYnR3NkQ4cC3bqJ4cuciJHUQJJKCvz/7//EunqOzMmppMv6ww+Bq1eVjed/VJ0MZWdno23btvjkk08q3GfgwIFISUkxXjZu3Ghy+/Tp07F27Vp8//332L17N7KysjBkyBAUqXydFCK7kBOC0aOBdu0UDcUuDAbgvffE9rJlwPHjioZDhJgY4OBBMWhBLwn6PfeIiRizs0s+jwpTdTI0aNAgvPnmmxgxYkSF+7i7uyMgIMB4qV27tvG29PR0fPnll1i4cCH69euHiIgIfPfdd0hMTMRmLt5IehcbC/z8szhT0+L6Y5bq0gUYPlx0Db7zjtLRkN7JdXqPPy6mgdADg6GkJfqTT4DLl5WNBypPhqpi+/btqF+/Ppo1a4bHHnsMV65cMd4WHx+PgoICREdHG68LCgpCeHg49u7dW+F95uXlISMjw+RC5HBmzRI/H3rI8ZvmbyUP8V2xQqxdRqSEnTuBXbsANzdg5kylo7GvwYPF8jg3b6piQIOmk6FBgwZh+fLl2Lp1KxYuXIjY2Fj06dMHeXl5AIDU1FS4ubmhVq1aJn/n7++P1NTUCu933rx58PPzM15CQkJs+jyI7G7XLuCPP0RRsV6a5kvr2BGIjgaKilQ3xJd05K23xM+JE4EGDZSNxd5Ktw599pkYbq8gTSdDY8aMweDBgxEeHo6hQ4fit99+w6lTp7Bhw4ZK/06SJBgqWXPppZdeQnp6uvFy4cIFa4cumujXrwdmz7b+fRPdTumDcOPGysaiFHn9py+/BCo5OSKyidhYYNMmMYJMLqDWmwEDgLvuEgMaFO6y1nQydKvAwECEhobi9OnTAICAgADk5+fj+vXrJvtduXIF/pX0zbq7u8PX19fkYnVnz4oisjfeEPO8ENnLkSOiVcjJybGW3TBXz57iQJyXVzL7NpG9yCckDz0klsDRo9KtQ59/ruh0Fw6VDKWlpeHChQsIDAwEAHTo0AGurq6IiYkx7pOSkoKjR48iKipKqTCFJk3EMF8A+OADRUMhnVm0SPwcOVK/rUKAOBDLtUOLFwO3nDQR2UxiophXyGAoaaHUq759gWHDxGfRFg0PVaTqZCgrKwsJCQlISEgAACQlJSEhIQHJycnIysrCzJkzsW/fPpw7dw7bt2/H0KFDUbduXQwfPhwA4Ofnh0mTJuHZZ5/Fli1bcPjwYTz00ENo3bo1+vXrp+Az+x95ReJly1Qz1wI5uH/+KVkkUW8Fm+UZPBho3VpMAPfpp0pHQ3ohL8Z6331A8+bKxqI0g0EkhrNnK5oMQVKxbdu2SQDKXMaPHy/dvHlTio6OlurVqye5urpKDRs2lMaPHy8lJyeb3EdOTo40depUqXbt2pKnp6c0ZMiQMvvcTnp6ugRASk9Pt+bTk6TiYknq0EGSAEl64w3r3jdReV54QbzfevRQOhL1WLFCvCZ16khSVpbS0ZCj++svSTIYxHvuzz+VjsbhVfX72yBJkqRcKqYNGRkZ8PPzQ3p6uvXrh1asEKsT+/uLlcM9PKx7/0SyzEyxcnR6OrBuXUk3rd4VFoo1y86cEbVD06crHRE5ssmTRbfskCFiEA3ZVFW/v1XdTaYL990HBAeLSadWrlQ6GnJk//2vSISaNxfdQyS4uJSM5nn/fZEcEdnCjRsliyLLZRKkCkyGlObqCkybJrYXLRLr1BBZW0FBSaH+s8+KkWRU4uGHgTp1gORk4NdflY6GHNXSpWIJivBwoFcvpaOhUnhEVIPHHgNq1ACOHhXr1BBZ208/iS/6+vWBceOUjkZ9PDzE5xAQywMQWVtRUcl7a9o0UThMqsFkSA1q1gQmTRLbCxcqGgo5IEkqmWX5qadYl1aRJ54QLWZbtnABV7K+334T88vVqiXqRElVmAypxdNPiwPxpk2ihYjIWvbsAQ4dEqtiP/mk0tGoV2iomAgV4DB7sr6PPhI/H30U8PJSNhYqg8mQWoSFASNGiG1OwkjWtGSJ+Dl2rKiLoYpNnSp+fv21KDYnsobjx0UJhJOTGE1GqsNkSE2eflr8XLmSB2KyjqtXgR9/FNtPPKFsLFrQuzfQooUocpVH/RBVl1wrNGwY0KiRoqFQ+ZgMqUnXruJAfPNmySzBRNWxbBmQnw906CBWaqfKGQwlrUOffCIWVCaqjtLD6eWRw6Q6TIbUxGAAHn9cbH/+OYfZU/UUF4v3EcBWIXOMGwf4+ACnToliaqLq4HB6TWAypDbjxonRPn/+CRw8qHQ0pGVbtwJ//y3W+7n/fqWj0Q4fH+CRR8T2xx8rGwtpW3FxSTH+U09xOL2KMRlSm9q1xazUAPDFF8rGQtomF06PGyfmsaKqk4tcf/0VSEpSNhbSru3bxTIvvr4cTq9yTIbUSO4q+/57FlKTZS5dAn7+WWzL7yequubNgf79RVf1l18qHQ1plfzeGTsW8PZWNhaqFJMhNYqKAlq1EoXU332ndDSkRV99JWa87doVaN1a6Wi0SZ6Retky8VoSmeP6dWD1arEtT6pLqsVkSI0MBuD//k9ss5CazFVUVNLFysJpyw0bJrqt//mHy+SQ+VasAPLygDZtxGhOUjUmQ2olF1InJgIHDigdDWnJb78BFy6IL/JRo5SORrvc3UvqPJYuVTYW0h65i2ziRBZOawCTIbWqVQsYPVpsy8OjiapCbhWaMIHrkFXXxIni588/A2lpioZCGnL4sLi4uQEPPaR0NFQFTIbUTC58XbVKTNxFdDtXroiWIYB1CtbQrh0QESEmruREqFRVcqvQvfdyCRyNYDKkZnfdJQqpc3JKllQgqszKlUBhoZhtumVLpaNxDHLr0FdfKRsHaUNODrB8udjmCYlmMBlSM4MBePhhsf3tt8rGQtogT/s/fryycTiSsWNFd0dCguj6IKrM2rWiJb9hQ6BfP6WjoSpiMqR2Dz4okqJduzj5G1Xu6FHg0CHA1ZUzTltT7dqiuwNgITXdntyC+MgjYpV60gT+p9SuQQOgb1+xzTmHqDJyq9DgwUDdusrG4mjkrrLvvgNyc5WNhdQrKUmsZ2cwlCzpQprAZEgL5K6yb77hnENUvqKikmRZfr+Q9fTrBwQHi4n01q1TOhpSq6+/Fj/79QNCQ5WNhczCZEgLhg8HvLzEopucc4jKs3kzkJIiunQGD1Y6Gsfj7FxSh8VCaiqPJJUUTrNmT3OYDGlBjRrAiBFim4XUVB65i+yBB0SxL1nfhAniZ0wMkJqqaCikQnFx4oTVywu45x6loyEzMRnSCrnr4/vvxZwnRLKMDDGCBeAZqS01aQJERgLFxcAPPygdDamN3Cp0zz3iBJY0hcmQVvTpAwQFAdeuARs3Kh0NqclPP4m5Te68U8wvRLYzdqz4uXKlsnGQuhQViRNVoGQJF9IUJkNa4exc8iGTu0SIANO5hbgGkm2NHi2GS+/fD5w9q3Q0pBZbtwKXL4vZpqOjlY6GLMBkSEvGjRM/f/1VtBARnTsH7NghkiCugWR7AQGilRYoaQkgkpdque8+Mc8XaQ6TIS1p3VqslVRQwJoFElatEj979xZDv8n2HnhA/GRXGQGii3r1arHNLjLNYjKkNXLrECdgJKAkGRozRtk49GTECDFi7+hRIDFR6WhIaRs2AJmZYvmNqCiloyELMRnSmjFjRJfInj3AxYtKR0NKOn1arJXl7Fwy9QLZXs2awN13i222DpHcRfbAA1x+Q8P4n9OaBg2Abt3E9k8/KRsLKUvuKu3Xj8tv2FvpUWWcFV6/btwQLUMAu8g0jsmQFo0eLX6ybkjf2EWmnCFDxFwy586JkWWkT6tXi3nfwsNFTSdpFpMhLRo5UnSV7dsHJCcrHQ0p4cQJUa/i6lqyojrZj6enWCYHKOkmIf2R//dySyFpFpMhLQoMBHr0ENvsKtMnuVUoOhqoVUvZWPRKHlX2ww9AYaGysZD9paQA27aJbfm9QJrFZEir2FWmX5LELjI1kGu1rlwRk+6RvqxdKz6LkZFAo0ZKR0PVxGRIq0aMEF1lBw6IugXSj8RE4ORJwN2dC0IqydVVTLIH8KREj+S5hUaNUjYOsgomQ1oVEAD07Cm22VWmL3Kr0KBBgK+vsrHonfxF+Msv7CrTk3//BbZvF9sjRyoaClkHkyEtY1eZ/rCLTF169BDrUV29CuzapXQ0ZC+//AIUFwPt2wNhYUpHQ1bAZEjLRowQk3zFxrKrTC8OHwbOnBGjmYYMUToacnEp6aqUu03I8cmt8WwVchiqToZ27tyJoUOHIigoCAaDAT///LPxtoKCArzwwgto3bo1vL29ERQUhIcffhiXLl0yuY9evXrBYDCYXO6//347PxMb8fcHevUS2z/+qGgoZCdyq9DgwWKeG1Ke/IW4Zo1oLSDHdv06sGWL2Ga9kMNQdTKUnZ2Ntm3b4pNPPilz282bN3Ho0CHMmjULhw4dwpo1a3Dq1CkMGzaszL6PPfYYUlJSjJfPP//cHuHbB7vK9EOSSpJe+f9OyuvbV9RupaRwAkY9WLdO1IeFhwPNmikdDVmJi9IBVGbQoEEYNGhQubf5+fkhJibG5LqPP/4YnTt3RnJyMho2bGi83svLCwEBATaNVTEjRgCTJwNxccDZs0DjxkpHRLaSmAgkJQEeHiVrY5Hy3N1Fl+WKFaKrjIt1Oja5i4ytQg5F1S1D5kpPT4fBYEDNmjVNrl++fDnq1q2LVq1aYebMmcjMzKz0fvLy8pCRkWFyUa169YDevcX22rXKxkK2Jf9/o6MBb29lYyFTpbvKuFaZ48rIADZtEtusF3IoDpMM5ebm4sUXX8TYsWPhW2q48YMPPoiVK1di+/btmDVrFlavXo0Rt1nhe968efDz8zNeQkJCbB1+9cjLApSqqSIHJP9/5f83qcfAgYCXlxjIcPiw0tGQrfz6q1iLrHlzoFUrpaMhK3KIZKigoAD3338/iouL8dlnn5nc9thjj6Ffv34IDw/H/fffj59++gmbN2/GoUOHKry/l156Cenp6cbLhQsXbP0Uqkeuk9qzR8yGS44nKQlISBCjBzmKTH28vMS8TwBHlTky+X8rrw9JDkPzyVBBQQFGjx6NpKQkxMTEmLQKlad9+/ZwdXXF6dOnK9zH3d0dvr6+JhdVCwkBOnYUzfPr1ysdDdmC3CrUo4dYAoLUR25xXr2aXWWOKDsb+O03sc16IYej6WRIToROnz6NzZs3o06dOrf9m2PHjqGgoACBgYF2iNCO5JXL2VXmmOT/K1eoV68hQwA3N+Cvv4Djx5WOhqztt9+AnBwxyWK7dkpHQ1am6mQoKysLCQkJSEhIAAAkJSUhISEBycnJKCwsxKhRoxAXF4fly5ejqKgIqampSE1NRX5+PgDgzJkzeP311xEXF4dz585h48aNuO+++xAREYGuXbsq+MxsQP6SjIkBsrIUDYWs7N9/gd27xTaTIfXy9QX69xfba9YoGwtZn/w/ZReZQ1J1MhQXF4eIiAhEREQAAGbMmIGIiAjMnj0bFy9exLp163Dx4kW0a9cOgYGBxsvevXsBAG5ubtiyZQsGDBiA5s2bY9q0aYiOjsbmzZvh7Oys5FOzvpYtgSZNgLw84I8/lI6GrGn9+pKp/0NDlY6GKiOPMGLdkGMpKAA2bhTbHMDgkFQ9z1CvXr0gVdL3XtltABASEoIdO3ZYOyx1MhhEq8GCBaJLhcM+HYc8pJ6tQuo3bBjg7Az8+acoeue6VY5h1y4gPV1MZRIZqXQ0ZAOqbhkiM8lrJP36qziTIe3LyhJdnwDPSLWgTh2gWzexzcEMjmPdOvFzyBCR7JLDYTLkSO66S5y53LgB7NypdDRkDb//Lro+77iD85poxdCh4ieTIccgSSXJUDnLPZFjYDLkSJydSz6sHFXmGOQusuHDWbSpFfJncPt20bVC2nbsmOjydHcvKZAnh8NkyNGUHmLPuU60LT8f2LBBbLNeSDuaNhUzFBcWcjCDI5Bbhfr14zI4DozJkKPp21d8YC9eBCqZZZs0YMcO0bLg7y+6QEk72FXmONhFpgtMhhyNp6dYJwlgV5nWyV+kQ4eKZThIO+Qvzg0bRAsRaVNqKnDggNjmMjgOjUdYR8TZqLVPksSoQIAHYS266y6gdm3g+nXgf/OekQbJn8FOnYCgIGVjIZtiMuSI7r5btCQcPQokJysdDVni5ElRtOnmJro+SVtcXMTnEGBXmZaxi0w3mAw5otq1S2pM5AJc0hb5/9a7N1CjhrKxkGXkL1D5C5W05ebNkjm+mAw5PCZDjmrwYPGTyZA2sYtM+wYMAFxdgVOnxOKtpC1btgC5uWIJnNatlY6GbIzJkKOSk6GtW8VKy6QdN26ULMwq/x9Je3x9gV69xDa7yrSndBcZ5/hyeEyGHFXr1kBwsEiEtm1TOhoyxx9/AEVFYvFdrm2lbRxir03FxSX/M3aR6QKTIUdlMLCrTKvk/xdbhbRPTob27AHS0pSNhaouLg64fFm07vXooXQ0ZAdMhhxZ6WSIs1FrQ1ERsHGj2Ga9kPY1aiRaaYuKgN9+Uzoaqir5MxgdLUZ0ksNjMuTI+vQR6+mcPw8cP650NFQVBw+KFoSaNYGoKKWjIWuQW4fkonhSPzlxHTRI2TjIbpgMOTJvbzE0G2BXmVbIX5gDB4q5akj75PmGNm0SLUSkbv/+C8TGim15Nn9yeEyGHB3rhrSFQ+odT2SkaOm7fl20/JG6/fGHKCto146zTusIkyFHJydDe/aIgzGp14ULwJEjYvZwnpE6DhcXUXsCsG5IC9hFpktMhhxdWBjQooVont+0SeloqDJy691ddwF16igbC1mXnNwyGVK3oiLg99/FNpMhXWEypAfsKtMGDql3XHIyFBcHXLmibCxUsdhY4No1wM+vZEkj0gUmQ3ogf7n+9hsLONUqN1dM/w8wGXJEgYGiBgUQNSmkTnLLXXQ0BzDoDJMhPejaVZzpXL1aMkqC1GX3bjFbeFAQ10FyVHK3i9wNQ+ojzy/ELjLdYTKkB66uLOBUO/kLcsAAroPkqOQvWHm5FVKXy5dFNybAAQw6xGRIL+QPN5vo1Un+vwwYoGwcZDtduojlHdLSSr50ST3kz2BEhOjWJF1hMqQXcsuQXCBI6vHPP8DRo2JIfb9+SkdDtuLqCvTvL7bZQqs+HFKva0yG9CI4WKyCXlwMbN6sdDRUmnxG2qkTh9Q7OvmLlsmQuhQVlXwO5RnDSVeYDOmJ3AXD+YbUhV1k+iF3V8fGigENpA4HDohJaWvWFDOGk+4wGdIT+ctWnm6elFdUBMTEiG0WbTq+Bg2ANm3E548nJerBIfW6x2RIT3r0ADw8gIsXgRMnlI6GANFCIJ+RduqkdDRkD+wqUx/OOq17TIb0xNNTJEQAR5Wphfx/6NePZ6R6UXpkZ3GxsrGQGN0XHy+25YEmpDtMhvSmdFcZKU/+P7CLTD+6dgV8fIB//wUOHVI6GtqyRXRbhodzlXodYzKkN/KZz44dYsZjUs7166JwE2DxtJ64ugK9e4ttuV6MlCPXbsnTHpAuMRnSm1atRBFnbi6wa5fS0ejb5s2im6RlSzH1AemH/MXLZEhZklTyP2AXma4xGdIbg6HkQ8+uMmWxi0y/5GRozx7g5k1lY9Gz06eB5GTAza2knpJ0icmQHnG+IeVJkul6ZKQvzZoBISFAfj6wc6fS0eiXfAzs1g3w8lI2FlIUkyE96tdPtBAdPSqWgiD7O35cvPYeHkD37kpHQ/ZmMLCrTA3k1571QrrHZEiP6tQpmdOGrUPKkLvIevYUUx6Q/sjd1UyGlFFQAGzbJrZZL6R7TIb0ikPslSWvD8eDsH717StaiBITgdRUpaPRnwMHgMxMoG5doF07paMhhTEZ0qvSZ6Wc+M2+SteJcJV6/apbF4iIENtcPNn+5Fbxfv0AJ34V6h3fAXrVpYuY+O3aNSAhQelo9OXAASA7G6hfX0z0RvrFuiHlsF6ISlF1MrRz504MHToUQUFBMBgM+Pnnn01ulyQJc+fORVBQEDw9PdGrVy8cO3bMZJ+8vDw89dRTqFu3Lry9vTFs2DBcvHjRjs9CpVxcRL0KwLNSe5Nf7759eUaqd6WTIS6ebD/XrwMHD4ptJkMElSdD2dnZaNu2LT755JNyb58/fz4WLVqETz75BLGxsQgICED//v2RmZlp3Gf69OlYu3Ytvv/+e+zevRtZWVkYMmQIioqK7PU01EvuotmyRdk49EZOhthFRl27ihGFKSlihCHZx7ZtojzgzjvFFAdEkkYAkNauXWv8vbi4WAoICJDeeecd43W5ubmSn5+ftGTJEkmSJOnGjRuSq6ur9P333xv3+eeffyQnJyfp999/r/Jjp6enSwCk9PT06j8RNUlMlCRAkjw9JSk3V+lo9CE9XZKcncXrfu6c0tGQGkRHi/fD++8rHYl+PPGEeM2nTVM6ErKxqn5/q7plqDJJSUlITU1FdKnROO7u7ujZsyf27t0LAIiPj0dBQYHJPkFBQQgPDzfuU568vDxkZGSYXBxSq1aAv79Yo2zfPqWj0YcdO4CiIqBJEyA0VOloSA1YN2R/XI+MbqHZZCj1f0NR/f39Ta739/c33paamgo3NzfUqlWrwn3KM2/ePPj5+RkvIY7ajGowiLoVgF1l9sIuMrqV/IW8Y4cYaUi2dfasuLi6Ar16KR0NqYRmkyGZwWAw+V2SpDLX3ep2+7z00ktIT083Xi5cuGCVWFVJ/lJmEbV9yEknkyGStW4tRhZmZ7OF1h7kY12XLkCNGsrGQqqh2WQoICAAAMq08Fy5csXYWhQQEID8/Hxcv369wn3K4+7uDl9fX5OLw5JbhmJjAUftDlSLlBTg2DHRIte7t9LRkFo4OZUkx+wqs72tW8VP+dhHBA0nQ2FhYQgICEBMqYNHfn4+duzYgaioKABAhw4d4OrqarJPSkoKjh49atxH9xo2FPUrRUWimZ5sR24Vat8eqF1b2VhIXeSuMrbQ2lZxMZMhKpeL0gFUJisrC3///bfx96SkJCQkJKB27dpo2LAhpk+fjrfffhtNmzZF06ZN8fbbb8PLywtjx44FAPj5+WHSpEl49tlnUadOHdSuXRszZ85E69at0Y/dFCX69QP+/lsciIcOVToax8V6IapInz7iZ1ycaKF15NZoJR07Bvz7r1ihvnNnpaMhFVF1MhQXF4fepboTZsyYAQAYP348li1bhueffx45OTmYPHkyrl+/jsjISGzatAk+Pj7Gv3n//ffh4uKC0aNHIycnB3379sWyZcvg7Oxs9+ejWn37AkuWsIjaliSJyRBVTG6h/ftvYNcuYPBgpSNyTHKrUPfugJubsrGQqhgkidOe3k5GRgb8/PyQnp7umPVDaWlAvXriCzslBfhfPRZZ0cmTQIsWgLu7mP2WK9XTrR5/HPjiC2DGDGDhQqWjcUzDhgHr1wPz5wPPPad0NGQHVf3+1mzNEFlRnTolC0aydcg25Ne1WzcmQlQ+uatMbr0g6yosLKmLlF9rov9hMkQCl+awLXaR0e3Ic94kJIjWWrKuQ4dEPVbNmkC7dkpHQyrDZIgEeWTF5s1cMNLaCgvFWkgAR7BQxfz9gfBwsb19u6KhOCS5xa1XL4A1o3QLJkMkdOsmCgovXBBFnGQ9hw8D6emAn58YVk9UEXnACLvKrE9u9eYJCZXDotFkycnJFj1YzZo1HbMA2RF4eQFRUeKMdPNmoGlTpSNyHHKrUM+ePCOlyvXpA3z8MZMha8vLA3bvFtusF6JyWJQMNWrUyOy/MRgMmDNnDmbPnm3JQ5I99O0rkqFt24Ann1Q6GschJ0OcdZpup2dPMUP5yZPApUtAUJDSETmG/fuB3FwxUrZFC6WjIRWyKBkqLi62dhykBnIB5/btom7oNmu8URUUFIh5YwAmQ3R7tWqJrtT4eJFEP/ig0hE5BrmLrE8fHteoXBYlQ2FhYbddDLU806dPx7Rp0yx5SLKHzp3FsO9//wWOHwdatVI6Iu2LjxcLcNauLRbkJLqdPn3E+2brViZD1iJ3O7KLjCpgUTK0bNkyix7Mku41siM3N6BrV1EztG0bkyFrKF0v5MTxClQFffoA773HuiFrycoCDhwQ20yGqAIWJUM9e/Y0bl+4cAEhISFWC4gU1ru3SIa2bwemTlU6Gu1jvRCZq1s3wMUFOHcOSEoCwsKUjkjbdu0S01uEhfG1pApV+1T1zjvvxKxZs5CdnW2NeEhp8pf29u1ihWeyXH4+sGeP2GYyRFVVowYQGSm22TpUfewioyqodjIUExODTZs2oWnTpli6dKk1YiIldewIeHuLGXCPHlU6Gm07eBC4eVOs+8YuRzIHl+awHiZDVAXVToaioqJw4MABvPPOO5g9ezYiIiKwnbOnaperq2imBzgLbnXJXWS9enEEC5mndDLEGeEtd+OGmPQUKBktS1QOq1V0Pvzwwzh16hSGDh2KwYMHY/jw4fibMxlrk3zQkL/MyTKsFyJLdekCeHgAqaliziGyzK5dIpls1oxzNlGlrDq8RZIkREdH4//+7/+wbt06hIeH49lnn0VmZqY1H4ZsTf7y3rGDdUOWys0F9u4V20yGyFweHmJGeIAttNUhr1JfatAPUXmqnQwtWbIEkyZNQps2beDn54d+/fphz549mDJlCj777DMkJCSgZcuWiIuLs0a8ZA/t24sizuvXgSNHlI5Gm/bvF0sABAQAzZsrHQ1pkfwFLn+hk/nkRJJdZHQbFg2tL+2tt95Cly5dMH78eHTp0gUdO3aEu7u78faJEyfi7bffxoQJE3CUBbna4OoKdO8O/PabOJi0a6d0RNrDeiGqLjkZ4ozwlklPL6kXYssQ3Ua1k6ELFy7cdp9JkyZh1qxZ1X0osqdevUQytG0bMH260tFoD+uFqLoiIwF3d+DyZeDUKbYwmmv3btHN36QJ0KCB0tGQytllStz69etjK4eIakvpuqGiImVj0ZqbN0U3GcBkiCzn4SEKqQF2lVmCXWRkBqslQ/Hx8RXeZjAYTGatJg2IiAB8fUVT859/Kh2NtuzbJxZobdBAnJUSWap0VxmZh8kQmcFqydDw4cOtdVekBi4uom4I4BB7c5XuImOdB1WH/EW+YwfnGzJHRgZw6JDY5ok4VYFZNUOjR48u93pJknDt2jWrBEQq0rs3sGGDOMN69lmlo9EOuUuDZ6RUXV26iAWUL10CzpxhS2NVyfVCd9wBBAcrHQ1pgFnJ0ObNm/Htt9+iRo0aJtdLkoSdO3daNTBSAfnLfOdOsdChS7Xr7R1fTo5YhgPgGSlVn6cn0Lmz+HLfvp3JUFWxi4zMZNa3W69evVCjRo1y638iIiKsFhSpRLt2gJ9fSd1Qhw5KR6R+Bw6IBVoDA8VZKVF19eolkqEdO4BHH1U6Gm3gZItkJrNqhtasWVNhIfTvv/9ulYBIRZydS9YpY8tf1civU8+erBci6yg9+SLrhm4vMxOQB/QwGaIqqlYBdWpqqrXiILXq0UP8ZDJUNfIZqfy6EVXXXXeJLuoLF4CkJKWjUb89e8R0II0bAw0bKh0NaUS1kqHo6GhrxUFqVToZ4jpllcvPF8PqAZ6RkvV4ewOdOoltzjd0e3K9ED+DZIZqJUMSm2wdX4cOgJcXcO0acPy40tGoW1ycKKCuWxdo0ULpaMiRlB5iT5Vj8TRZoFrJkIE1EY7P1bVk9Wx2lVVOfn169GC9EFkXJ1+smqwscVICsGWIzGKX5ThI41g3VDWlkyEia4qKEgMazp8XFyqfXC/UqBEQGqp0NKQhTIbo9uQvd45mqVhhoRj+DPCMlKzPxwfo2FFss6usYqVHcxKZoVrJkJubm7XiIDXr3FnMgpuaCvz9t9LRqNOff4ohvX5+QOvWSkdDjohdZbfH1lmyULWSoTi5b5Ycm6cnEBkpttlVVj75bL17d9GdQWRtpecborJyc0tmf5fXVSSqInaTUdWwbqhyPCMlW+vaVRTmnz0r1iojUwcPiuktAgK4bAmZjckQVU3puiEyVVwM7NoltpkMka34+QFt24pt+f1GJeTXpHt3juYks9kkGbp+/Tq2bduG999/3xZ3T0q46y6OZqnIsWNiHiZvb6B9e6WjIUcmd/8wGSqrdDJEZKZqL0OelJSEhIQEk8vFixchSRK8vb3xzDPPWCNOUpqPj5iA8eBBcdDhsNUScmtZVJSYl4nIVrp3Bz7+mMnQrQoLxbB6gK2zZBGLW4Z69uyJmjVr4o477sDEiROxefNmBAQE4J9//sGXX36J8+fPIzMz05qxktJYN1Q+Ducle5FbPRITgRs3FA1FVf78U0y46OcHhIcrHQ1pkMXJ0L59+zBlyhRcuHAB169fx549e/D555/DYDCgc+fOCAkJsWacpAZMhsqSJBZPk/0EBABNm4r3ndwSQiUtZd26cTQnWcTiZOjAgQPYtWsXpkyZglOnTlkzpipr1KgRDAZDmcuUKVMAABMmTChzW5cuXRSJ1SF06yYKE//6S8w5RMDp08Dly4C7u5iPicjWWDdUlnxCwnohspDFyVBERAR27tyJ0aNHY8CAAZgyZQquXLlizdhuKzY2FikpKcZLTEwMAOC+++4z7jNw4ECTfTZu3GjXGB1KrVpAmzZimwdiQX4dIiNFQkRka/IXPltoBUli8TRVW7VHk40dOxbHjh1DzZo10apVKxQXF6OoqMgasd1WvXr1EBAQYLz8+uuvuOOOO9CzVO2Gu7u7yT61a9e2S2wOi0PsTclLcHTrpmwcpB/yF35cHJCTo2wsanDyJHD1KuDhUbJkCZGZrDK03svLC2+99RYOHDiAIUOGoG/fvliwYAFy7PhBzc/Px3fffYeJEyfCUGqOie3bt6N+/fpo1qwZHnvssSq1XuXl5SEjI8PkQv8jH4jlJEDveEZK9ta4MRAYCBQUAAcOKB2N8uTPYJcuYtkgIgtYdZ6hxo0b45dffsHy5cuxdOlSNG7c2Jp3X6mff/4ZN27cwIQJE4zXDRo0CMuXL8fWrVuxcOFCxMbGok+fPsjLy6v0vubNmwc/Pz/jhcXgpcgtIEeOAOnpysaitJQU4MwZUUd1111KR0N6YTCUtNCyu5onJGQVBkkyfxny5OTk2+5TUFCAX3/9FcOHDzdeV7NmTfj6+pr7cFUyYMAAuLm5Yf369RXuk5KSgtDQUHz//fcYMWJEhfvl5eWZJEwZGRkICQlBenq6zeLXlCZNRBLw22/AwIFKR6OcH38ERo8G2rUDDh9WOhrSk08/BaZOBfr3BzZtUjoaZYWGAsnJQEwM0K+f0tGQymRkZMDPz++2398WTbrYqFGjKu8rT7poMBgwZ84czJ4925KHrNT58+exefNmrFmzptL9AgMDERoaitOnT1e6n7u7O9xZDFuxbt1EMrR7t76TIdYLkVLkVpB9+8SEgy7Vnj9Xm86fF4mQs7PoJiOykEWfoOLiYmvHUS1Lly5F/fr1MXjw4Er3S0tLw4ULFxAYGGinyBxUt27A11+zbqj03CZE9hQeDtSsKSZeTEjQb+Gw/Bls3x6oUUPZWEjTLEqGwsLCTIqUq2r69OmYNm2aJQ9ZoeLiYixduhTjx4+HS6mzo6ysLMydOxcjR45EYGAgzp07h5dffhl169Y16bojC8hf/gcOiFWi9Vi0mJEhZr0FmAyR/Tk5iVXsN2wQCYHekyFOeErVZFEytGzZMosezJzutaravHkzkpOTMXHiRJPrnZ2dkZiYiG+++QY3btxAYGAgevfujVWrVsHHx8fqcehK8+ZAnTpAWhoQH6/P4uF9+8Rq9WFhQIMGSkdDetS9u0iGdu4E9LoGJIunyUosSoZ6qmgNpujoaJRXA+7p6Yk//vhDgYh0wGAQrSG//CK6yvSYDMldhDwIk1Lk1pDdu8XEgxa01mva1avAiRNiu2tXZWMhzbPq0HrSEb3PN8TiaVJahw6Ap6dICk6eVDoa+9u7V/xs0QKoW1fZWEjzmAyRZeQkYM8e0V2kJ/n5wP79YpvJECnFzU0sAwPoc74heaFafgbJCpgMkWUiIsRZaVqaWLhVTw4dAnJzRd3UnXcqHQ3pWemTEr2RW2fZRUZWwGSILKPns9LSQ+r1VqdB6iInAnpLhnJzxdpsAFuGyCqYDJHl5IOQ3uqGWDxNanHXXSIhP3MGSE1VOhr7iYsT3dX+/mKtNqJqYjJEltNjEXVxMYunST38/IDWrcW2nlqHSn8G2TpLVsBkiCzXpYuY/C0pCfjnH6WjsY+TJ4Fr10S9VESE0tEQ6bOrTH6urBciK2EyRJbz9QXathXbejkQy/VCXbroc+ZtUh+9dVcXF3MkGVkdkyGqHr0diNlFRmojt44cPgzcvKlsLPZw8iRw/Trg5QW0a6d0NOQgmAxR9chJgV5GlPGMlNSmYUOxJExhIXDwoNLR2J58QhIZCbi6KhsLOQwmQ1Q9clJw5AiQnq5sLLaWkiLqowwG0U1GpAby8jiAPlpoeUJCNsBkiKonKEgsVlpcXDIrs6OSD8KtW4t6KSK10FMRNSdbJBtgMkTVJx+U5LWCHJX8/HgQJrWRW0n27gWKipSNxZZSUoCzZ8UoVj0uEE02w2SIqk8vZ6Uczktq1bo1UKMGkJEBHDumdDS2w9ZZshEmQ1R9UVHi54EDoojTEd28KdYkA5gMkfq4uJTUsTnySQnrhchGmAxR9bVqJc7SsrKAxESlo7GN2FiR6AUGAqGhSkdDVJYeiqg5tQXZCJMhqj5n55KzUketGypdL8Tp/0mNHL27OjtbzKUEsHWWrI7JEFmHox+IWS9EahcZKQqLz58HLl5UOhrrO3BAFIc3bAiEhCgdDTkYJkNkHXLdkCO2DBUXcyQZqZ+PT8mMzI54UsITErIhJkNkHaXPSh1t0VZ5+n9PT07/T+rmyC20PCEhG2IyRNbh41OyaKujtQ7Jz6dzZ07/T+omJwqOVkRdXAzs2ye25VZoIitiMkTWIx+kHO2slM3zpBXyKKs//xSjOx3FiRNiuR9vbzHHEJGVMRki63HUJnomQ6QVDRqIAuPiYsdatFX+DEZGijmViKyMyRBZj9wydPiwGAbrCK5cAU6fFtuc/p+0wBEHM8jPhV1kZCNMhsh6GjYUZ6ZFRWKSQkcg1ym0bAnUqqVsLERVISftTIaIqozJEFmPweB4i7ayi4y0Rk4Y9u8X3WVa9++/Ja2z8uSuRFbGZIisy9GKqJkMkda0bSumgbh+HfjrL6WjqT62zpIdMBki65KThn37tH9WmpsLxMWJbSZDpBWurmIaCMAxWmjZRUZ2wGSIrKttW8DLS5yVnjypdDTVc+gQkJ8P1KsH3HGH0tEQVZ0jFVEzGSI7YDJE1lX6rFTrXWVcnJW0ylGSofz8ksEYTIbIhpgMkfU5St2QXKvAIfWkNXKh8cmTQFqasrFUR0KC6K6uXRto1kzpaMiBMRki6ytdN6RVklRyVs1kiLSmbt2S5GH/fmVjqY7SXWRsnSUbYjJE1ieflZ46BVy9qmwsljp/HkhNFbPdduyodDRE5nOErjLWC5GdMBki66tdG2jeXGxr9axUbtWKiBDDlIm0Rk4gtNpCK0klXe1MhsjGmAyRbWj9QMwzUtI6+b174ABQWKhsLJZITgYuXQKcnYFOnZSOhhwckyGyDbnORqvJEIunSetatAD8/ICbN4EjR5SOxnzyCUlEhJiug8iGmAyRbchJxMGD2jsrzc4Wo1gAJkOkXU5O2l6njK2zZEdMhsg2WrYEfH1FYpGYqHQ05omLE4vNNmgAhIQoHQ2R5ZgMEVUJkyGyDSenklFlWusqK91FxuG8pGVaHVGWnQ38+afYZjJEdsBkiGxHq3VDcrw8CJPWde4sTkzOnxfFyFoRG8vWWbIrTSdDc+fOhcFgMLkEBAQYb5ckCXPnzkVQUBA8PT3Rq1cvHDt2TMGIdUaLTfScbJEcia8v0Lq12NbSSQkHMJCdaToZAoBWrVohJSXFeEksVZ8yf/58LFq0CJ988gliY2MREBCA/v37IzMzU8GIdSQyUnQznT0LXLmidDRVc+aMmCjSzU2MYiHSOi2elDAZIjvTfDLk4uKCgIAA46VevXoARKvQBx98gFdeeQUjRoxAeHg4vv76a9y8eRMrVqxQOGqdqFlTFFID2jkrlb8wOnYE3N2VjYXIGrTWXS1JJZO1MhkiO9F8MnT69GkEBQUhLCwM999/P86ePQsASEpKQmpqKqKjo437uru7o2fPnth7mzOkvLw8ZGRkmFzIQlo7EPOMlByN/F4+dAjIy1M2lqo4exb491/ROtu+vdLRkE5oOhmKjIzEN998gz/++AP/+c9/kJqaiqioKKSlpSE1NRUA4O/vb/I3/v7+xtsqMm/ePPj5+RkvISzgs5zWmuhZPE2OpkkToE4dkQjJ82epmfwZbN+erbNkN5pOhgYNGoSRI0eidevW6NevHzZs2AAA+Prrr437GG4ZGi1JUpnrbvXSSy8hPT3deLlw4YL1g9cLOamIiwMKCpSN5XYyM0vmRGLLEDkKg0Fb01zIMcoxE9mBppOhW3l7e6N169Y4ffq0cVTZra1AV65cKdNadCt3d3f4+vqaXMhCzZoBtWoBOTkl84ao1cGDQHEx0KgREBiodDRE1qOl7mp2VZMCHCoZysvLw4kTJxAYGIiwsDAEBAQgJibGeHt+fj527NiBKHaB2E/pyRfV3lXGIfXkqOT3tFyYrFbZ2SXrqPFzSHak6WRo5syZ2LFjB5KSknDgwAGMGjUKGRkZGD9+PAwGA6ZPn463334ba9euxdGjRzFhwgR4eXlh7NixSoeuL1o5K+UZKTmqTp3EiYm8ErxacSkcUoiL0gFUx8WLF/HAAw/g6tWrqFevHrp06YL9+/cjNDQUAPD8888jJycHkydPxvXr1xEZGYlNmzbBx8dH4ch1Rm6JU3MyxOG85Mh8fIDwcNHqsm8fMHKk0hGVjyckpBBNJ0Pff/99pbcbDAbMnTsXc+fOtU9AVL5blwQIClI6orJOnQKuXwc8PIC2bZWOhsj67rpLJEP79zMZIrqFprvJSCPks1JAva1DcqtQx46Aq6uysRDZgtq7qyWJI8lIMUyGyD7kA/GBA8rGURE5GeJBmByV/N6OiwPy85WNpTzyZIuurpxskeyOyRDZh9rnOWHzPDm6Zs2A2rXF5ItqnOZCPiFp3150VxPZEZMhso/SZ6Vqm3wxK6tkskW2DJGjUvvkizwhIQUxGSL7kCdfzM0tmUdELeLixGSLISHqLO4mshY11w0xGSIFMRki+3ByAiIjxbbaJn7jQZj0Qm4ZUttnMDu7pOuOn0NSAJMhsh+1HohZPE160bmz6C47dw64zYLVdiVPthgUBAQHKx0N6RCTIbIfNdYrlJ5skckQOTpfX3VOc1G6dfY2C2kT2QKTIbKfzp3FzzNnxBBaNUhKAq5cEcN5IyKUjobI9tTYQsvZ30lhTIbIfmrVAlq0ENtqmW+Iw3lJb9RWRC1JJccDts6SQpgMkX2p7ayUXWSkN3IypJZpLpKTRf2SiwsnWyTFMBki+2IyRKSsZs2AmjWBnJyS+bWUJH8G27UDPD0VDYX0i8kQ2ZecdBw4IEaPKCknBzh8WGyzVoH0Qm3TXPCEhFSAyRDZV6tWQI0aYtbn48eVjeXQIaCwEAgIABo2VDYWIntSUwstkyFSASZDZF/OziWjypQ+EJc+CHM4L+mJWpKhvDxxUgIwGSJFMRki+1PLgVgeTcODMOmNfEJy+jSQlqZcHAkJQH4+ULcu0LixcnGQ7jEZIvtTSzLEuU1Ir2rXFoXUAHDwoHJxsHWWVILJENmfXLx5/Dhw44YyMVy8CPzzj+i269BBmRiIlKSGkxLOL0QqwWSI7K9+feCOO8S2Umel8hdAmzaAt7cyMRApSQ3JkPzY8gkSkUKYDJEylD4QcwQL6V3paS6Ki+3/+Jcvi+VwDAagUyf7Pz5RKUyGSBlMhoiU1bq1mOQwPR04dcr+jy93kbVsCfj52f/xiUphMkTKKJ0MSZJ9H7ugAIiPF9tsnie9cnEBOnYU20qclPCEhFSEyRApo00bwN0duH5dDO+1pyNHgNxcsXBs06b2fWwiNVGyhZbJEKkIkyFShptbySgue69gLz9e585iaQIivVIqGSoqKhk8wWSIVIDfBKQcpdZH4hkpkSB/BhITxRI59nLsGJCdDfj4AC1a2O9xiSrAZIiUU3o0iz1xOC+REBQEBAeL0WRyHZ09lG6ddXa23+MSVYDJEClHTkb+/FOsIG8PaWklNUrykgREeqZEVxlbZ0llmAyRcho2BPz9xcrx8mKNtibXKTRtCtSpY5/HJFIzJZMhts6SSjAZIuUYDPbvKuP0/0Sm7D3NRXo6cOKE2GYyRCrBZIiUZe8iap6REplq317MOZSaCiQn2/7xDh4USVdYmFiah0gFmAyRsuzZMlRczOG8RLfy9ATathXb9vgcsnWWVIjJECmrY0fRXZacDKSk2PaxTp8Wkzx6eIhJH4lIkBOTffts/1hyMsTWWVIRJkOkLB8foFUrsW3rs1L5/jt0AFxdbftYRFoiJya2/gxKEruqSZWYDJHy7NVVxoMwUfnkz8ShQ0B+vu0eJykJuHpVzEAfEWG7xyEyE5MhUp69iqhZq0BUvqZNxVp9eXli7T5bkT+D7dqJtQmJVILJEClPToZiY8WaRbZw86aY3LH04xGRYDDYp6uM9UKkUkyGSHktWwI1aoi1io4ds81jxMeLRCswEAgJsc1jEGkZkyHSMSZDpDxnZ6BTJ7FtqwNx6YOwwWCbxyDSMlsnQ3l5JTPNMxkilWEyROpg6yJqroVEVDl5rb5Tp8QUFNb255+iOLtOHeCOO6x//0TVwGSI1MHWRdRsnieqXJ06QJMmYluenNSa2DpLKqbpZGjevHno1KkTfHx8UL9+fdx7773466+/TPaZMGECDAaDyaULWwfUR05Sjh8HMjKse9+XLgEXLwJOTmKSRyIqny27ynhCQiqm6WRox44dmDJlCvbv34+YmBgUFhYiOjoa2dnZJvsNHDgQKSkpxsvGjRsVipgqFBAAhIaKSdliY6173/JBuFUrUahNROWzZQstkyFSMRelA6iO33//3eT3pUuXon79+oiPj0ePHj2M17u7uyMgIMDe4ZG5IiOB8+fFQbNvX+vdLw/CRFUjt5rLi6laqzsrLQ34+2+xLdcmEamIpluGbpWeng4AqF27tsn127dvR/369dGsWTM89thjuHLlSqX3k5eXh4yMDJML2YGtmuiZDBFVTdu2YjLEtDTgzBnr3a/8GWzeXEzuSKQyDpMMSZKEGTNmoFu3bggPDzdeP2jQICxfvhxbt27FwoULERsbiz59+iAvL6/C+5o3bx78/PyMlxDOS2MfpZMhSbLOfRYVAXFxpvdPROUrvUyGNU9KeEJCKucwydDUqVNx5MgRrFy50uT6MWPGYPDgwQgPD8fQoUPx22+/4dSpU9iwYUOF9/XSSy8hPT3deLlw4YKtwycAaN8ecHEBLl8Wq9hbw/HjQFaWqBVq2dI690nkyGzRQstkiFTOIZKhp556CuvWrcO2bdsQHBxc6b6BgYEIDQ3F6dOnK9zH3d0dvr6+JheyA09PoE0bsW2tA7F8Px07iskdiahy1k6GJKlkqD6TIVIpTSdDkiRh6tSpWLNmDbZu3YqwsLDb/k1aWhouXLiAwMBAO0RIZpOLK601zwnPSInMI39WEhLErNHVdfq0mMTRw6PkZIdIZTSdDE2ZMgXfffcdVqxYAR8fH6SmpiI1NRU5OTkAgKysLMycORP79u3DuXPnsH37dgwdOhR169bF8OHDFY6eymXts1ImQ0TmCQsD6tYVs0UnJFT//uRh+h06AK6u1b8/IhvQdDK0ePFipKeno1evXggMDDReVq1aBQBwdnZGYmIi7rnnHjRr1gzjx49Hs2bNsG/fPvj4+CgcPZVLTlri44GCgurdV1ZWycKvTIaIqsbaK9jzhIQ0QNPzDEm3GXHk6emJP/74w07RkFU0bw74+QHp6cDRoyUjWywRFwcUFwPBwUBQkPViJHJ0XboAGzaIVp1p06p3X0yGSAM03TJEDsjJyXor2PMgTGQZa7UM5eSIBVpL3yeRCjEZIvWx1oGYyRCRZeQTkrNngX//tfx+EhKAwkLA3x9o2NAqoRHZApMhUh8mQ0TKqlkTuPNOsV2dtQLlz2DnzlypnlSNyRCpj5y8nDwpaocs8c8/YrV6Z2cxioWIzGONkxKekJBGMBki9alfH2jUSEzWJi+lYS75IBweDnh7Wy00It2Q5/yqTjLEyRZJI5gMkTpV90BcunmeiMwnJzDyCvbm+vdfUXMElNQgEakUkyFSp+o20bN5nqh62rQRK9hfvw78/bf5fy+3Ct15p5gug0jFmAyROlVnBXuuVE9Ufa6uYvFkwLKTEp6QkIYwGSJ1qs4K9seOAdnZYqX6Fi1sEx+RHlSnhZb1QqQhTIZInaqzgr28f6dOXKmeqDpK1w2ZgyvVk8YwGSL1svSslM3zRNZh6Qr2pVeqb93aJqERWROTIVIvS89KeUZKZB2NGlm2gr18QtK+PVeqJ01gMkTqZckK9qVXqueweqLqKb2CvTknJTwhIY1hMkTq1ayZGJKbkyNWsK+K+HiuVE9kTZZ0V3OeL9IYJkOkXqVXsK/qWSnrhYisy9wJUHNzS7rU+DkkjWAyROpm7lmpnDTxjJTIOuTP0t9/A9eu3X7/P/8U3dr16omaIyINYDJE6iYfiNkyRKSMWrVElzVQtc9h6c8gV6onjWAyROomJ0PHjwMZGZXve+kScPGi6F7jSvVE1mNOCy3rhUiDmAyRugUEAA0biknc4uMr31c+a23VSsw+TUTWYU4LLVtnSYOYDJH6VfWslAdhItuo6lqBaWnAmTNimy1DpCFMhkj9qpoMsXiayDbatgXc3ESyc/ZsxfvJn8FmzYCaNe0SGpE1MBki9atKE31RERAbK7bZMkRkXW5uQESE2K7sc8jWWdIoJkOkfu3biwVX5QLp8pw8CWRmAl5eQMuW9o2PSA+q0kLLmadJo5gMkfp5ewPh4WK7orNS+fqOHQEXF/vERaQnXbqIn3v2lH87V6onDWMyRNpwu7NSNs8T2Vb37uLnoUOiFfZWZ86ImiJ3d6BNG/vGRlRNTIZIG25XN8TiaSLbCg4GGjcWa//t3Vv2dvmEJCJC1BgRaQiTIdIGucUnLk4US5d28yZw5IjpfkRkfT16iJ87dpS9jV1kpGFMhkgbWrQQEylmZYnZqEs7fFgkSAEB4uyViGyjZ0/xc+fOsrdx5mnSMCZDpA3OzqI4GijbVca1kIjsQ24ZOngQyMkpuT4vT5yUAGwZIk1iMkTaUVERNYuniewjLAxo0ECsSr9/f8n1R44A+flAnTqirohIY5gMkXaUV0R96RKwe7fp7URkGwZD+V1lpbvI2DpLGsRkiLRDbvlJTASuXQPmzRPT/l+6BNSqBXTqpGx8RHogd5WVlwyxdZY0irPTkXY0aCAu//wjkqC0NHF9ZCTw6aeAr6+y8RHpgZwM7dsnusbc3Di1BWkeW4ZIW+SDbVoaEBgIfPONmPOkQwdl4yLSizvvBOrVEwXUcXHA9evAqVPiNiZDpFFMhkhbJk0CQkOBl14SB+Bx4wAnvo2J7MZgMO0qk1uFmjQRBdREGsRuMtKWwYPFhYiU06MHsHq1mHyxoEBcx1Yh0jAmQ0REZB65Zaj0oq0sniYNY/8CERGZp3VroGZNsWDrpk3iOiZDpGFMhoiIyDzOzkC3bmK7uBhwdQXatlU2JqJq0E0y9NlnnyEsLAweHh7o0KEDdu3apXRIRETaJXeVAUC7doCHh2KhEFWXLpKhVatWYfr06XjllVdw+PBhdO/eHYMGDUJycrLSoRERaVPpZIhdZKRxukiGFi1ahEmTJuHRRx9FixYt8MEHHyAkJASLFy9WOjQiIm1q3x7w9hbbHElGGufwyVB+fj7i4+MRHR1tcn10dDT27t1b7t/k5eUhIyPD5EJERKW4ugLPPgt07MjpLkjzHD4Zunr1KoqKiuDv729yvb+/P1JTU8v9m3nz5sHPz894CQkJsUeoRETa8tprQGwsULu20pEQVYvDJ0Mywy0rKUuSVOY62UsvvYT09HTj5cKFC/YIkYiIiBTg8JMu1q1bF87OzmVaga5cuVKmtUjm7u4Od3d3e4RHRERECnP4liE3Nzd06NABMTExJtfHxMQgKipKoaiIiIhILRy+ZQgAZsyYgXHjxqFjx46466678MUXXyA5ORlPPPGE0qERERGRwnSRDI0ZMwZpaWl4/fXXkZKSgvDwcGzcuBGhoaFKh0ZEREQKM0iSJCkdhNplZGTAz88P6enp8PX1VTocIiIiqoKqfn87fM0QERERUWWYDBEREZGuMRkiIiIiXWMyRERERLrGZIiIiIh0jckQERER6RqTISIiItI1JkNERESka7qYgbq65HkpMzIyFI6EiIiIqkr+3r7d/NJMhqogMzMTABASEqJwJERERGSuzMxM+Pn5VXg7l+OoguLiYly6dAk+Pj4wGAxWu9+MjAyEhITgwoULXObDhvg62wdfZ/vha20ffJ3tw5avsyRJyMzMRFBQEJycKq4MYstQFTg5OSE4ONhm9+/r68sPmh3wdbYPvs72w9faPvg624etXufKWoRkLKAmIiIiXWMyRERERLrGZEhB7u7umDNnDtzd3ZUOxaHxdbYPvs72w9faPvg624caXmcWUBMREZGusWWIiIiIdI3JEBEREekakyEiIiLSNSZDREREpGtMhmzss88+Q1hYGDw8PNChQwfs2rWr0v137NiBDh06wMPDA40bN8aSJUvsFKm2mfM6r1mzBv3790e9evXg6+uLu+66C3/88Ycdo9Uuc9/Psj179sDFxQXt2rWzbYAOwtzXOS8vD6+88gpCQ0Ph7u6OO+64A1999ZWdotU2c1/r5cuXo23btvDy8kJgYCAeeeQRpKWl2Sla7dm5cyeGDh2KoKAgGAwG/Pzzz7f9G0W+ByWyme+//15ydXWV/vOf/0jHjx+Xnn76acnb21s6f/58ufufPXtW8vLykp5++mnp+PHj0n/+8x/J1dVV+umnn+wcubaY+zo//fTT0rvvvisdPHhQOnXqlPTSSy9Jrq6u0qFDh+wcubaY+zrLbty4ITVu3FiKjo6W2rZta59gNcyS13nYsGFSZGSkFBMTIyUlJUkHDhyQ9uzZY8eotcnc13rXrl2Sk5OT9OGHH0pnz56Vdu3aJbVq1Uq699577Ry5dmzcuFF65ZVXpNWrV0sApLVr11a6v1Lfg0yGbKhz587SE088YXLdnXfeKb344ovl7v/8889Ld955p8l1jz/+uNSlSxebxegIzH2dy9OyZUvptddes3ZoDsXS13nMmDHSq6++Ks2ZM4fJUBWY+zr/9ttvkp+fn5SWlmaP8ByKua/1e++9JzVu3Njkuo8++kgKDg62WYyOpCrJkFLfg+wms5H8/HzEx8cjOjra5Pro6Gjs3bu33L/Zt29fmf0HDBiAuLg4FBQU2CxWLbPkdb5VcXExMjMzUbt2bVuE6BAsfZ2XLl2KM2fOYM6cObYO0SFY8jqvW7cOHTt2xPz589GgQQM0a9YMM2fORE5Ojj1C1ixLXuuoqChcvHgRGzduhCRJuHz5Mn766ScMHjzYHiHrglLfg1yo1UauXr2KoqIi+Pv7m1zv7++P1NTUcv8mNTW13P0LCwtx9epVBAYG2ixerbLkdb7VwoULkZ2djdGjR9siRIdgyet8+vRpvPjii9i1axdcXHioqQpLXuezZ89i9+7d8PDwwNq1a3H16lVMnjwZ165dY91QJSx5raOiorB8+XKMGTMGubm5KCwsxLBhw/Dxxx/bI2RdUOp7kC1DNmYwGEx+lySpzHW327+868mUua+zbOXKlZg7dy5WrVqF+vXr2yo8h1HV17moqAhjx47Fa6+9hmbNmtkrPIdhzvu5uLgYBoMBy5cvR+fOnXH33Xdj0aJFWLZsGVuHqsCc1/r48eOYNm0aZs+ejfj4ePz+++9ISkrCE088YY9QdUOJ70GertlI3bp14ezsXOYM48qVK2WyXllAQEC5+7u4uKBOnTo2i1XLLHmdZatWrcKkSZPw448/ol+/frYMU/PMfZ0zMzMRFxeHw4cPY+rUqQDEl7YkSXBxccGmTZvQp08fu8SuJZa8nwMDA9GgQQP4+fkZr2vRogUkScLFixfRtGlTm8asVZa81vPmzUPXrl3x3HPPAQDatGkDb29vdO/eHW+++SZb761Aqe9BtgzZiJubGzp06ICYmBiT62NiYhAVFVXu39x1111l9t+0aRM6duwIV1dXm8WqZZa8zoBoEZowYQJWrFjB/v4qMPd19vX1RWJiIhISEoyXJ554As2bN0dCQgIiIyPtFbqmWPJ+7tq1Ky5duoSsrCzjdadOnYKTkxOCg4NtGq+WWfJa37x5E05Opl+bzs7OAEpaL6h6FPsetGl5ts7Jwza//PJL6fjx49L06dMlb29v6dy5c5IkSdKLL74ojRs3zri/PKTwmWeekY4fPy59+eWXHFpfBea+zitWrJBcXFykTz/9VEpJSTFebty4odRT0ARzX+dbcTRZ1Zj7OmdmZkrBwcHSqFGjpGPHjkk7duyQmjZtKj366KNKPQXNMPe1Xrp0qeTi4iJ99tln0pkzZ6Tdu3dLHTt2lDp37qzUU1C9zMxM6fDhw9Lhw4clANKiRYukw4cPG6cvUMv3IJMhG/v000+l0NBQyc3NTWrfvr20Y8cO423jx4+XevbsabL/9u3bpYiICMnNzU1q1KiRtHjxYjtHrE3mvM49e/aUAJS5jB8/3v6Ba4y57+fSmAxVnbmv84kTJ6R+/fpJnp6eUnBwsDRjxgzp5s2bdo5am8x9rT/66COpZcuWkqenpxQYGCg9+OCD0sWLF+0ctXZs27at0uOtWr4HDZLEtj0iIiLSL9YMERERka4xGSIiIiJdYzJEREREusZkiIiIiHSNyRARERHpGpMhIiIi0jUmQ0RERKRrTIaIiIhI15gMERERka4xGSIiIiJdYzJERA7lww8/RFhYGLy8vHDvvfciPT1d6ZCISOW4NhkROYyXX34ZP/74I7788kvUqFEDw4cPx8iRI7Fo0SLjPhMmTEBAQAD8/PywZs0anDx5Ep6enoiKisK7776L5s2bK/gMiEgJTIaIyCHExsaiS5cuiI2NRfv27QEAb7/9NpYtW4ZTp04BAIqLi+Hv749169bhtddew/33349OnTqhsLAQr7zyChITE3H8+HF4e3sr+VSIyM7YTUZEDmHBggXo06ePMRECgHr16uHq1avG3/fs2QMnJydERkbi999/x4QJE9CqVSu0bdsWS5cuRXJyMuLj443733nnnTAYDOVePvroI7s+PyKyHSZDRKR5eXl5WL9+PYYPH25yfU5ODvz8/Iy/r1u3DkOHDoWTU9lDn1xbVLt2beN1a9euBQBs2bIFKSkpSE5OhouLC3788Uc8/vjjtngqRKQAJkNEpHmHDh1CTk4Onn32WdSoUcN4ee6550xqgNatW4d77rmnzN9LkoQZM2agW7duCA8PN16fmpoKFxcXdO3aFQEBAUhLS0NhYSG6d+8Od3d3uzw3IrI9F6UDICKqrlOnTsHDwwOJiYkm1w8bNgxdu3YFAJw4cQIXL15Ev379yvz91KlTceTIEezevdvk+sTERDRr1syY+CQkJKBevXrw9/e30TMhIiUwGSIizcvIyED9+vXRpEkT43XJyck4efIkRo4cCUC0CvXv3x+enp4mf/vUU09h3bp12LlzJ4KDg01uO3LkCFq3bm38PSEhAW3atLHhMyEiJbCbjIg0r27dusjIyEDpwbFvvfUW7r77brRs2RIA8Msvv2DYsGHG2yVJwtSpU7FmzRps3boVYWFhZe73yJEjJskPkyEix8RkiIg0r0+fPsjNzcU777yDc+fO4e2338a6deuwePFiAMCVK1cQGxuLIUOGGP9mypQp+O6777BixQr4+PggNTUVqampyMnJASCG4R87dswk+Tl79ixCQ0Pt++SIyOaYDBGR5vn7+2PZsmVYvHgxWrZsib1792L37t0ICQkBAKxfvx6RkZGoX7++8W8WL16M9PR09OrVC4GBgcbLqlWrAABnzpzBzZs3TZKhtm3bYu7cudi5c6d9nyAR2RQnXSQihzds2DB069YNzz//vNKhEJEKsWWIiBxet27d8MADDygdBhGpFFuGiIiISNfYMkRERES6xmSIiIiIdI3JEBEREekakyEiIiLSNSZDREREpGtMhoiIiEjXmAwRERGRrjEZIiIiIl1jMkRERES69v9EvKan26kAswAAAABJRU5ErkJggg==", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "PyObject Text(0.5, 1.0, 'maximizing $\\\\Vert A^{-1} y \\\\Vert$ over angle $\\\\theta$')" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using PyPlot\n", "θ = range(0,2π,length=100)\n", "plot(θ/(2π), [norm(A \\ [cos(θ), sin(θ)]) for θ in θ], \"r-\")\n", "plot(θ/(2π), ones(length(θ))*norm(inv(A)), \"k--\")\n", "xlabel(L\"\\theta / 2\\pi\")\n", "ylabel(L\"\\Vert A^{-1} y \\Vert\")\n", "text(0.16,opnorm(inv(A))*0.95, L\"\\Vert A^{-1} \\Vert\")\n", "title(L\"maximizing $\\Vert A^{-1} y \\Vert$ over angle $\\theta$\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yup, the maximum certainly seems to be matching what Julia's `norm` function computes for $\\Vert A^{-1} \\Vert$.\n", "\n", "But there must be a better way to compute it. Maybe $\\Vert A^{-1} \\Vert$ is related to $|\\lambda|$?" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " -166.5334932692097\n", " 0.20015993587216604" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "eigvals(inv(A))" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 166.5334932692097\n", " 0.20015993587216604" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs.(eigvals(inv(A)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hey, it certainly *looks* like $\\Vert A^{-1} \\Vert$ is the magnitude of the biggest $\\lambda$ of $A^{-1}$. Let's check it more carefully:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-0.0003336000409035478" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs.(eigvals(inv(A)))[1] - opnorm(inv(A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hmm, this difference of 0.0003 is *much* bigger than the $\\sim 10^{-15}$ we normally get from mere roundoff errors. Let's try another matrix $B$, chosen at random:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.8720563752626838" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B = randn(2,2)\n", "opnorm(B)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 1.7338448966661324\n", " 1.7338448966661324" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "abs.(eigvals(B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Whoops, no, in this case the eigenvalues have *much* smaller magnitudes than $\\Vert B \\Vert$. So the near-match in the case of $A^{-1}$ was just a coincidence?\n", "\n", "In fact, the right answer involves the **singular values** of the matrix:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 2.8720563752626838\n", " 1.0467127844662298" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "svdvals(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The **norm is the biggest singular value** of the matrix!" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "svdvals(B)[1] - opnorm(B)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "svdvals(inv(A))[1] - opnorm(inv(A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(In fact, computing the singular values is precisely how Julia's `norm` function works.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Rayleigh quotient and the min–max theorem\n", "\n", "Why is $\\Vert B \\Vert$ equal to the largest singular value of $B$? We can derive it!\n", "\n", "To avoid square roots, it is convenient to look at $\\Vert B \\Vert^2$:\n", "\n", "$$\n", "\\Vert B \\Vert^2 = \\max_{y\\ne 0} \\frac{\\Vert By \\Vert^2}{\\Vert y \\Vert^2} \\\\ = \\max_{y\\ne 0} \\frac{(By)^H (By)}{y^H y} = \\max_{y\\ne 0} \\frac{y^H (B^H B) y}{y^H y}\n", "$$\n", "\n", "This final ratio is called the [Rayleigh quotient](https://en.wikipedia.org/wiki/Rayleigh_quotient) of the *Hermitian* matrix $S = B^H B$:\n", "$$\n", "R_S(y) = \\frac{y^H S y}{y^H y}\n", "$$\n", "\n", "It is a remarkable and important fact, known as the [min–max theorem](https://en.wikipedia.org/wiki/Min-max_theorem) that the *maximum* of the Rayleigh quotient is given by the *biggest λ of S*!\n", "\n", "To prove this, since $S$ is an $m\\times m$ Hermitian matrix, we use that it's eigenvalues $\\lambda_1,\\ldots,\\lambda_m$ are real and the corresponding eigenvectors $q_1,\\ldots,q_m$ can be chosen orthonormality. Any vector $y$ can be written in this basis as $y = c_1 q_1 + \\cdots + c_m q_m$ for some coefficients $c_1,\\ldots,c_m$, and then orthonormality means:\n", "$$\n", "R_S(y) = \\frac{(c_1 q_1 + \\cdots + c_m q_m)^H (\\lambda_1 c_1 q_1 + \\cdots \\lambda_m c_m q_m)}{(c_1 q_1 + \\cdots + c_m q_m)^H (c_1 q_1 + \\cdots + c_m q_m)} = \\frac{\\lambda_1 |c_1|^2 + \\cdots + \\lambda_m |c_m|^2}{ |c_1|^2 + \\cdots + |c_m|^2}\n", "$$\n", "which is just a [weighted average](https://en.wikipedia.org/wiki/Weighted_arithmetic_mean) of the eigenvalues λ. If we number the eigenvalues in order $\\lambda_1 \\ge \\lambda_2 \\ge \\cdots \\lambda_m$, we immediately get:\n", "$$\n", "R_S(y) = \\frac{\\lambda_1 |c_1|^2 + \\cdots + \\lambda_m |c_m|^2}{ |c_1|^2 + \\cdots + |c_m|^2} \\le \\frac{\\lambda_1 |c_1|^2 + \\cdots + \\lambda_1 |c_m|^2}{ |c_1|^2 + \\cdots + |c_m|^2} = \\lambda_1\n", "$$\n", "so $\\lambda_1$ is an upper bound, achieved by $R_S(q_1) = \\lambda_1$.\n", "\n", "So,\n", "$$\n", "\\Vert B \\Vert = \\sqrt{\\max\\;\\lambda\\;\\mathrm{of}\\;B^T B}\n", "$$\n", "Let's check:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2.8720563752626838" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(B)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 1.0467127844662298\n", " 2.8720563752626838" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sqrt.(eigvals(B'B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But we already learned, in a previous lecture, that the **eigenvalues of BᴴB** are precisely the **the squares σ² of the singular values σ of B**." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 2.8720563752626838\n", " 1.0467127844662298" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "svdvals(B)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key properties of the matrix norm\n", "\n", "Note that $\\Vert B\\Vert$ is defined above for *any* matrix, including non-square matrices. Some key properties of this norm are:\n", "\n", "* Positivity: $\\Vert B \\Vert \\ge 0$, and $= 0$ only if $B = 0$. (Obvious from the definition and the SVD.)\n", "\n", "* Scaling: $\\Vert \\alpha B \\Vert = |\\alpha|\\;\\Vert B \\Vert$. (Obvious from the definition.)\n", "\n", "* [Triangle inequality](https://en.wikipedia.org/wiki/Triangle_inequality): $\\Vert B_1 + B_2 \\Vert \\le \\Vert B_1 \\Vert + \\Vert B_2 \\Vert$. (Not so obvious, but easy to show from the definition and the triangle inequality for vector norms.)\n", "\n", "* Products: $\\Vert B_1 B_2 \\Vert \\le \\Vert B_1 \\Vert \\; \\Vert B_2 \\Vert$. (Not so obvious, but easy to show from the definition above.)\n", "\n", "It's easy to try some of these:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "10.0" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(10*B) / opnorm(B)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.6003396253535439" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "B₂ = randn(size(B))\n", "opnorm(B + B₂) / (opnorm(B) + opnorm(B₂))" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.8142010561693229" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(B * B₂) / (opnorm(B) * opnorm(B₂))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Norm of the inverse\n", "\n", "One thing to be careful of is the inverse. Even if $A$ is invertible, $\\Vert A \\Vert^{-1} \\ne \\Vert A^{-1} \\Vert$ in general! However, there is a cute trick to get $A^{-1}$. Let's look at the definition:\n", "\n", "$$\n", "\\Vert A^{-1} \\Vert = \\max_{y\\ne 0} \\frac{\\Vert A^{-1} y \\Vert}{\\Vert y \\Vert} = \\max_{z\\ne 0} \\frac{\\Vert z \\Vert}{\\Vert Az \\Vert}\n", "$$\n", "\n", "where we have substituted $z = A^{-1} y$, i.e. $y = Az$. But then\n", "\n", "$$\n", "\\frac{1}{\\Vert A^{-1} \\Vert} = \\min_{z\\ne 0} \\frac{\\Vert Az \\Vert}{\\Vert z \\Vert}\n", "$$\n", "\n", "since the maximum is 1 over the minimum. But, exactly analogous to the discussion above, $\\Vert Az \\Vert / \\Vert z \\Vert$ is has a minimum for the *minimum* singular value of $A$! So\n", "\n", "* $\\Vert A^{-1} \\Vert$ is the *inverse* of the *minimum* singular value of $A$ (assuming $A$ is invertible) \n", "\n", "Let's check:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "166.53382686925062" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(inv(A))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2-element Vector{Float64}:\n", " 0.20015953491241695\n", " 166.53382686924056" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "1 ./ svdvals(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yup, they match!" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.0061285138363019e-11" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(inv(A)) - 1/svdvals(A)[2]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(At least up to roundoff errors.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compared to what? Relative errors\n", "\n", "There is still a problem with our presentation above. How do we decide whether an error $\\Vert \\Delta x \\Vert$ is big or small? Compared to what?\n", "\n", "Above, we compared to $\\Vert \\Delta b \\Vert$. But this is often *not* the right thing to do. We can see that in a couple of ways:\n", "\n", "* In a physical system, $x$ and $b$ may have *different units*. e.g. $b$ may be a current and $x$ may be a voltage. $\\Vert \\Delta x \\Vert / \\Vert \\Delta b \\Vert$ is therefore a *dimensionful* quantity, and we can't say whether it is big or small without comparing to some other dimensionful quantity with the same units.\n", "\n", "* If we multiply $A$ by 1000, it is easy to see from above that $\\Vert A \\Vert$ multiples by 1000 and $\\Vert A^{-1} \\Vert$ is *divided* by 1000. But simply scaling the problem shouldn't change whether the errors are big or small! (You can't reduce the errors in your experiment by changing units from meters to millimeters!)\n", "\n", "The right thing to do is generally to compare $\\Vert \\Delta x \\Vert$ to $\\Vert x \\Vert$ and $\\Vert \\Delta b \\Vert$ to $\\Vert b \\Vert$. The ratio\n", "$$\n", "\\frac{\\Vert \\Delta x \\Vert}{\\Vert x \\Vert}\n", "$$\n", "is called the [relative error](https://en.wikipedia.org/wiki/Approximation_error), as opposed to the *absolute* error $\\Vert \\Delta x \\Vert$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Condition number of a matrix\n", "\n", "How much can the relative error grow when we solve $Ax = b$? Well, from above:\n", "\n", "$$\n", "\\frac{\\Vert \\Delta x \\Vert}{\\Vert x \\Vert} \\le \\frac{\\Vert A^{-1} \\Vert \\; \\Vert \\Delta b \\Vert}{\\Vert x \\Vert} = \\frac{\\Vert A^{-1} \\Vert \\; \\Vert \\Delta b \\Vert}{\\Vert b \\Vert} \\frac{\\Vert Ax \\Vert}{\\Vert x \\Vert} \\le \\boxed{\\Vert A \\Vert \\; \\Vert A^{-1} \\Vert \\frac{\\Vert \\Delta b \\Vert}{\\Vert b \\Vert}}\n", "$$\n", "\n", "The quantity $\\kappa(A) = \\Vert A \\Vert \\; \\Vert A^{-1} \\Vert$ is called the [condition number](https://en.wikipedia.org/wiki/Condition_number) of the matrix A. It gives an **upper bound** on how much the **relative error can grow** when comparing input to output. Since $\\kappa(A) = \\kappa(A^{-1})$, it is the same for computing $x = A^{-1} b$ or for computing $b = Ax$ (i.e. x from b or b from x)! In Julia, it is computed by `cond(A)`." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "832.005464751455" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cond(A)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "832.0054647515053" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "opnorm(A) * opnorm(inv(A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the definition above and the relationship of the matrix norm to the singular values, the **condition number is the ratio of the largest to the smallest singular value** of A:\n", "\n", "$$\n", "\\boxed{\\kappa(A) = \\Vert A \\Vert \\; \\Vert A^{-1} \\Vert = \\frac{\\sigma_\\max}{\\sigma_\\min}}\n", "$$" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "832.005464751455" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "σ = svdvals(A)\n", "maximum(σ) / minimum(σ)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this, we can see that $\\kappa(A) \\ge 1$. **The smallest possible condition number is 1**. This happens for $\\kappa(I) = 1$, or for any multiple of $I$, and in fact **any unitary matrix** has condition number 1." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ill-conditioned matrices\n", "\n", "If the condition number is $\\gg 1$, we say that the matrix is **ill-conditioned** (or \"badly conditioned\"). When you solve an ill-conditioned problem, **any error can be greatly magnified**. This includes both measurement errors and things like roundoff errors during the calculations.\n", "\n", "Our matrix $A$ above pretty badly conditioned — it can magnify erros by a factor of $\\sim 1000$:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "832.005464751455" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cond(A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An **ill-conditioned matrix is a matrix that is *almost* singular**. A singular matrix has a condition number of \"∞\": the condition number blows up as one of the singular values goes to zero.\n", "\n", "It is easy to see that our matrix $A$ from above is almost singular: the second row is *almost* twice the first row:" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 1.0 2.0\n", " 2.01 3.99" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's why our error increased by so much:" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "787.1845290391896" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "norm(Δx)/norm(x) / (norm(Δb)/norm(b))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(In fact, I didn't pick `Δb` at random: I picked it to lie almost parallel to the singular vector for the smallest singular value, to magnify the error as much as possible.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Just because a matrix is ill-conditioned does not mean that you can't get accurate answers, but it indicates that you need to be **much more careful** about errors." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Nearly defective matrices\n", "\n", "In class, we discussed defective matrices, in which we are \"missing\" some eigenvectors (due to repeated roots) and hence don't have a basis. In practice, *exactly* defective matrices are extremely uncommon.\n", "\n", "However, it is much more common to have a *nearly* defective matrix $A$, in which two or more eigenvectors are *nearly* parallel. Such matrices are still diagonalizable, $A = X \\Lambda X^{-1}$, so why worry?\n", "\n", "The problem with a *nearly* defective matrix is that the eigenvector matrix $X$ is *nearly* singular — it is **ill-conditioned**! That means that *working with eigenvectors of a nearly-defective matrix is very sensitive to errors (roundoff, measurement, …)*!\n", "\n", "For example, let's construct a nearly defective matrix $M$ with eigenvectors $(1,0)$ and $(1,10^{-14})$, and eigenvalues $\\lambda = 1$ and $\\lambda = 1 + 10^{-14}$." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 1.0 0.999201\n", " 0.0 1.0" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X = [1 1\n", " 0 1e-14]\n", "M = X * Diagonal([1,1+1e-14]) / X" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, let's compute the matrix exponential $e^M$, both by the built-in Julia function `expm` and by the 18.06 formula involving the diagonalization $X e^\\Lambda X^{-1}$:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 2.71828 2.71611\n", " 0.0 2.71828" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "exp(M)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×2 Matrix{Float64}:\n", " 2.71828 2.70894\n", " 0.0 2.71828" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "X * exp(Diagonal([1,1+1e-14])) / X" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The two don't match very well, and in fact the diagonalization-based formula gives horribly wrong results (off by about 10%). Computing matrix exponentials (and many other matrix functions) robustly for arbitrary matrices (including nearly defective matrices) is a tricky problem. There is a famous paper called [Nineteen Dubious Ways to Compute the Exponential of a Matrix](http://epubs.siam.org/doi/abs/10.1137/S00361445024180) on the many possible approaches." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The \"opposite\" of a nearly defective matrix is a matrix with an *orthonormal* basis $X=Q$ of eigenvectors (e.g. a Hermitian, unitary, or other [\"normal\"](https://en.wikipedia.org/wiki/Normal_matrix) matrix). In this case, the eigenvector basis has condition number $= 1$, the best possible conditioning, and working with eigenvectors is great!" ] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.8.0", "language": "julia", "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.3" } }, "nbformat": 4, "nbformat_minor": 2 }