lectures/lecture-1/lecture-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 1: Floating-point arithmetic, vector norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Syllabus\n",
    "\n",
    "**Today:** \n",
    "- Part 1: floating point, vector norms\n",
    "- Part 2: matrix norms and stability concepts\n",
    "\n",
    "**Tomorrow:** Matrix norms and unitary matrices\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Representation of numbers\n",
    "\n",
    "- Real numbers represent quantities: probabilities, velocities, masses, ...\n",
    "\n",
    "<img src=\"./Babylonian_numerals.png\" width=500>\n",
    "\n",
    "- It is important to know, how they are represented in the computer, which only knows about bits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fixed point\n",
    "\n",
    "- The most straightforward format for the representation of real numbers is **fixed point** representation, also known as **Qm.n** format.\n",
    "\n",
    "- A **Qm.n** number is in the range $[-(2^m), 2^m - 2^{-n}]$, with resolution $2^{-n}$.\n",
    "\n",
    "- Total storage is $m + n + 1$ bits.\n",
    "\n",
    "- The range of numbers represented is fixed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Floating point\n",
    "The numbers in computer memory are typically represented as **floating point numbers** \n",
    "\n",
    "A floating point number is represented as  \n",
    "\n",
    "$$\\textrm{number} = \\textrm{significand} \\times \\textrm{base}^{\\textrm{exponent}},$$\n",
    "\n",
    "where *significand* is integer, *base* is positive integer  and *exponent* is integer (can be negative), i.e.\n",
    "\n",
    "$$ 1.2 = 12 \\cdot 10^{-1}.$$\n",
    "\n",
    "- This format has a long history. It was already used in the world's first working programmable, fully automatic digital computer [Z3](https://en.wikipedia.org/wiki/Z3_(computer)) designed in 1935 and completed in 1941 in Germany by [Konrad  Zuse](https://en.wikipedia.org/wiki/Konrad_Zuse)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Floating point: formula\n",
    "\n",
    "$$f = (-1)^s 2^{(p-b)} \\left( 1 + \\frac{d_1}{2} + \\frac{d_2}{2^2}  + \\ldots + \\frac{d_m}{2^m}\\right),$$\n",
    "\n",
    "where $s \\in \\{0, 1\\}$ is the sign bit, $d_i \\in \\{0, 1\\}$ is the $m$-bit mantissa, $p \\in \\mathbb{Z}; 0 \\leq p \\leq 2^e$, $e$ is the $e$-bit exponent, commonly defined as $2^e - 1$\n",
    "\n",
    "Can be thought as a uniform $m$-bit grid between two sequential powers of $2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fixed vs Floating\n",
    "\n",
    "**Q**: What are the advantages/disadvantages of the fixed and floating points?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**A**:  In most cases, they work just fine.\n",
    "\n",
    "- However, fixed point represents numbers within specified range and controls **absolute** accuracy.\n",
    "\n",
    "- Floating point represent numbers with **relative** accuracy, and is suitable for the case when numbers in the computations have varying scale (i.e., $10^{-1}$ and $10^{5}$).\n",
    "\n",
    "- In practice, if speed is of no concern, use float32 or float64."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## IEEE 754\n",
    "In modern computers, the floating point representation is controlled by [IEEE 754 standard](https://en.wikipedia.org/wiki/IEEE_floating_point) which was published in **1985** and before that point different computers behaved differently with floating point numbers. \n",
    "\n",
    "IEEE 754 has:\n",
    "- Floating point representation (as described above), $(-1)^s \\times c \\times b^q$.\n",
    "- Two infinities, $+\\infty$ and $-\\infty$\n",
    "- Two zeros: +0 and -0\n",
    "- Two kinds of **NaN**: a quiet NaN (**qNaN**) and signalling NaN (**sNaN**) \n",
    "    - qNaN does not throw exception in the level of floating point unit (FPU), until you check the result of computations\n",
    "    - sNaN value throws exception from FPU if you use corresponding variable. This type of NaN can be useful for initialization purposes\n",
    "    - C++11 proposes [standard interface](https://en.cppreference.com/w/cpp/numeric/math/nan) for creating different NaNs \n",
    "- Rules for **rounding**\n",
    "- Rules for $\\frac{0}{0}, \\frac{1}{-0}, \\ldots$\n",
    "\n",
    "Possible values are defined with\n",
    "- base $b$\n",
    "- accuracy $p$ - number of digits\n",
    "- maximum possible value $e_{\\max}$\n",
    "\n",
    "and have the following restrictions \n",
    "- $ 0 \\leq c \\leq b^p - 1$\n",
    "- $1 - e_{\\max} \\leq q + p - 1 \\leq e_{\\max}$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The two most common format, single & double\n",
    "\n",
    "The two most common formats, called **binary32** and **binary64** (called also **single** and **double** formats). Recently, the format **binary16** plays important role in learning deep neural networks.\n",
    "\n",
    "| Name | Common Name | Base | Digits | Emin | Emax |\n",
    "|------|----------|----------|-------|------|------|\n",
    "|binary16| half precision | 2 | 11 | -14 | + 15 |\n",
    "|binary32| single precision | 2 | 24 | -126 | + 127 |  \n",
    "|binary64| double precision | 2 | 53 | -1022 | +1023 |  \n",
    "\n",
    "<img src=\"./double64.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples\n",
    "\n",
    "- For a number +0\n",
    "    - *sign* is 0\n",
    "    - *exponent* is 00000000000\n",
    "    - *fraction* is all zeros\n",
    "- For a number -0\n",
    "    - *sign* is 1\n",
    "    - *exponent* is 00000000000\n",
    "    - *fraction* is all zeros\n",
    "- For +infinity\n",
    "    - *sign* is 0\n",
    "    - *exponent* is 11111111111\n",
    "    - *fraction* is all zeros\n",
    "\n",
    "**Q**: what about -infinity and NaN ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Accuracy and memory\n",
    "\n",
    "The **relative accuracy** of single precision is $10^{-7}-10^{-8}$, while for double precision is $10^{-14}-10^{-16}$.\n",
    "\n",
    "<font color='red'> Crucial note 1: </font> A **float16** takes **2 bytes**, **float32** takes **4 bytes**, **float64**, or double precision, takes **8 bytes.**\n",
    "\n",
    "<font color='red'> Crucial note 2: </font> These are the only two floating point-types supported in hardware (float32 and float64) + GPU/TPU different float types are supported.\n",
    "\n",
    "<font color='red'> Crucial note 3: </font> You should use **double precision** in computational science and engineering and **float** on GPU/Data Science.\n",
    "\n",
    "\n",
    "Also, half precision can be useful in training deep neural network, see this [paper](https://arxiv.org/pdf/1905.12334.pdf).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How does number representation format affect training of neural networks (NN)?\n",
    "\n",
    "- Weights in layers (fully-connected, convolutional, activation functions) can be stored with different accuracies\n",
    "- It is important to improve energy efficiency of the devices that are used to train NNs\n",
    "- Project [DeepFloat](https://github.com/facebookresearch/deepfloat) from Facebook demonstrates how re-develop floating point operations in a way to ensure efficiency in training NNs, more details see in this [paper](https://arxiv.org/pdf/1811.01721.pdf)\n",
    "- Affect of the real numbers representation on the gradients of activation functions\n",
    "- Typically, the first digit is one.\n",
    "- Subnormal numbers have first digit 0 to represent zeros and numbers close to zero.\n",
    "- Subnormal numbers fill the gap between positive and negative\n",
    "- They have performance issues, often flushed to zero by default.\n",
    "\n",
    "<img width=500, src=\"./grad_norm_fp16.png\">\n",
    "\n",
    "- And on the learning curves\n",
    "\n",
    "<img width=500, src=\"./train_val_curves.png\">\n",
    "\n",
    "Plots are taken from [this paper](https://arxiv.org/pdf/1710.03740.pdf%EF%BC%89%E3%80%82)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## bfloat16 (Brain Floating Point)\n",
    "\n",
    "- This format occupies 16 bits\n",
    "    - 1 bit for sign\n",
    "    - 8 bits for exponent\n",
    "    - 7 bits for fraction\n",
    "    <img src=\"./bfloat16.png\">\n",
    "- Truncated single precision format from IEEE standard\n",
    "- What is the difference between float32 and float16?\n",
    "- This format is utilized in Intel FPGA, Google TPU, Xeon CPUs and other platforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor Float from Nvidia ([blog post about this format](https://blogs.nvidia.com/blog/2020/05/14/tensorfloat-32-precision-format/))\n",
    "\n",
    "- Comparison with other formats\n",
    "\n",
    "<img src=\"./tensor_float_cf.png\">\n",
    "\n",
    "- Results\n",
    "\n",
    "<img src=\"./TF32-BERT.png\">\n",
    "\n",
    "- PyTorch and Tensorflow supported this format are available in [Nvidia NCG](https://ngc.nvidia.com/catalog/all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Mixed precision ([docs from Nvidia](https://docs.nvidia.com/deeplearning/performance/mixed-precision-training/index.html))\n",
    "\n",
    "- Main idea:\n",
    "    - Maintain copy of weights in single precision\n",
    "    - Then in every iteration\n",
    "        - Make a copy of weights in half-precision\n",
    "        - Forward pass with weights in half-precision\n",
    "        - Multiply the loss by the scaling factor $S$\n",
    "        - Backward pass again in half precision\n",
    "        - Multiply the weight gradient with $1/S$\n",
    "        - Complete the weight update (including gradient clipping, etc.)\n",
    "    - Scaling factor $S$ is a hyper-parameter\n",
    "    - Constant: a value so that its product with the maximum absolute gradient value is below 65504 (the maximum value representable in half precision).\n",
    "    - Dynamic update based on the current gradient statistics\n",
    "- Performance comparison\n",
    "<img src=\"./mixed_precision_res.png\" width=500>\n",
    "\n",
    "- Automatic mixed-precision extensions exist to simplify turning this option on, more details [here](https://developer.nvidia.com/automatic-mixed-precision)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Alternative to the IEEE 754 standard\n",
    "\n",
    "Issues in IEEE 754:\n",
    "- overflow to infinity or zero\n",
    "- many different NaNs\n",
    "- invisible rounding errors\n",
    "- accuracy is very high or very poor\n",
    "- subnormal numbers – numbers between 0 and minimal possible represented number, i.e. significand starts from zero\n",
    "\n",
    "Concept of **posits** can replace floating point numbers, see [this paper](http://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf)\n",
    "\n",
    "<img width=600 src=\"./posit.png\">\n",
    "\n",
    "- represent numbers with some accuracy, but provide limits of changing\n",
    "- no overflows!\n",
    "- example of a number representation \n",
    "\n",
    "<img width=600 src=\"./posit_example.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Division accuracy demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9259246\n",
      "       inf\n",
      "nan\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "from jax.config import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "#c = random.random()\n",
    "#print(c)\n",
    "c = jnp.float32(0.925924589693)\n",
    "print(c)\n",
    "a = jnp.float32(1e-50)\n",
    "b = jnp.float32(c / a)\n",
    "print('{0:10.16f}'.format(b))\n",
    "print(abs(a * b - c)/abs(c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Square root accuracy demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "float32\n",
      "       nan\n"
     ]
    }
   ],
   "source": [
    "a = jnp.float32(1e-40)\n",
    "b = jnp.sqrt(a)\n",
    "print(b.dtype)\n",
    "print('{0:10.64f}'.format(abs(b * b - a)/abs(a)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exponent accuracy demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "float32\n",
      "nan\n"
     ]
    }
   ],
   "source": [
    "a = jnp.float32(1e99)\n",
    "b = jnp.exp(a)\n",
    "print(b.dtype)\n",
    "print(jnp.log(b) - a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of demos\n",
    "\n",
    "- For some values the inverse functions give exact answers\n",
    "- The relative accuracy should be preserved due to the IEEE standard\n",
    "- Does not hold for many modern GPU\n",
    "- More details about adoptation of IEEE 754 standard for GPU you can find [here](https://docs.nvidia.com/cuda/floating-point/index.html#considerations-for-heterogeneous-world) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Loss of significance\n",
    "\n",
    "- Many operations lead to the loss of digits [loss of significance](https://en.wikipedia.org/wiki/Loss_of_significance)\n",
    "- For example, it is a bad idea to subtract two big numbers that are close, the difference will have fewer correct digits\n",
    "- This is related to algorithms and their properties (forward/backward stability), which we will discuss later"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summation algorithm\n",
    "\n",
    "However, the rounding errors can depend on the algorithm.\n",
    "\n",
    "- Consider the simplest problem: given $n$ floating point numbers $x_1, \\ldots, x_n$  \n",
    "\n",
    "- Compute their sum\n",
    "\n",
    "$$S = \\sum_{i=1}^n x_i = x_1 + \\ldots + x_n.$$\n",
    "\n",
    "- The simplest algorithm is to add one-by-one \n",
    "\n",
    "- What is the actual error for such algorithm? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Na&iuml;ve algorithm\n",
    "\n",
    "Na&iuml;ve algorithm adds numbers one-by-one: \n",
    "\n",
    "$$y_1 = x_1, \\quad y_2 = y_1 + x_2, \\quad y_3 = y_2 + x_3, \\ldots.$$\n",
    "\n",
    "- The **worst-case** error is then proportional to $\\mathcal{O}(n)$, while **mean-squared** error is $\\mathcal{O}(\\sqrt{n})$.\n",
    "\n",
    "- The **Kahan algorithm** gives the worst-case error bound $\\mathcal{O}(1)$ (i.e., independent of $n$).  \n",
    "\n",
    "- <font color='red'> Can you find the better algorithm?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Kahan summation\n",
    "The following algorithm gives $2 \\varepsilon + \\mathcal{O}(n \\varepsilon^2)$ error, where $\\varepsilon$ is the machine precision.\n",
    "\n",
    "- The reason of the loss of significance in summation is operating with numbers of different magnitude\n",
    "- The main idea of Kahan summation is to keep track of small errors and aggregate them in separate variable\n",
    "- This approach is called *compensated summation*\n",
    "\n",
    "```python\n",
    "s = 0\n",
    "c = 0\n",
    "for i in range(len(x)):\n",
    "    y = x[i] - c\n",
    "    t = s + y\n",
    "    c = (t - s) - y\n",
    "    s = t\n",
    "```\n",
    "\n",
    "- There exists more advanced tricks to process this simple operation that are used for example in ```fsum``` function from ```math``` package, implementation check out [here](https://github.com/python/cpython/blob/d267006f18592165ed97e0a9c2494d3bce25fc2b/Modules/mathmodule.c#L1087)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error in np sum: -8.3e-07\n",
      "Error in Kahan sum Numba: 4.7e-08\n",
      "Error in Kahan sum JAX: 0.0e+00\n",
      "Error in dumb sum: -1.0e-03\n",
      "Error in math fsum: 1.3e-11\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "import jax.numpy as jnp\n",
    "import numpy as np\n",
    "import jax\n",
    "from numba import jit as numba_jit\n",
    "\n",
    "n = 10 ** 7\n",
    "sm = 1e-10\n",
    "x = jnp.ones(n, dtype=jnp.float32) * sm\n",
    "x = x.at[0].set(1)\n",
    "#x = jax.ops.index_update(x, [0], 1.)\n",
    "true_sum = 1.0 + (n - 1)*sm\n",
    "approx_sum = jnp.sum(x)\n",
    "math_fsum = math.fsum(x)\n",
    "\n",
    "\n",
    "@jax.jit\n",
    "def dumb_sum(x):\n",
    "    s = jnp.float32(0.0)\n",
    "    def b_fun(i, val):\n",
    "        return val + x[i] \n",
    "    s = jax.lax.fori_loop(0, len(x), b_fun, s)\n",
    "    return s\n",
    "\n",
    "\n",
    "@numba_jit(nopython=True)\n",
    "def kahan_sum_numba(x):\n",
    "    s = np.float32(0.0)\n",
    "    c = np.float32(0.0)\n",
    "    for i in range(len(x)):\n",
    "        y = x[i] - c\n",
    "        t = s + y\n",
    "        c = (t - s) - y\n",
    "        s = t\n",
    "    return s\n",
    "\n",
    "@jax.jit\n",
    "def kahan_sum_jax(x):\n",
    "    s = jnp.float32(0.0)\n",
    "    c = jnp.float32(0.0)\n",
    "    def b_fun2(i, val):\n",
    "        s, c = val\n",
    "        y = x[i] - c\n",
    "        t = s + y\n",
    "        c = (t - s) - y\n",
    "        s = t\n",
    "        return s, c\n",
    "    s, c = jax.lax.fori_loop(0, len(x), b_fun2, (s, c))\n",
    "    return s\n",
    "\n",
    "k_sum_numba = kahan_sum_numba(np.array(x))\n",
    "k_sum_jax = kahan_sum_jax(x)\n",
    "d_sum = dumb_sum(x)\n",
    "print('Error in np sum: {0:3.1e}'.format(approx_sum - true_sum))\n",
    "print('Error in Kahan sum Numba: {0:3.1e}'.format(k_sum_numba - true_sum))\n",
    "print('Error in Kahan sum JAX: {0:3.1e}'.format(k_sum_jax - true_sum))\n",
    "print('Error in dumb sum: {0:3.1e}'.format(d_sum - true_sum))\n",
    "print('Error in math fsum: {0:3.1e}'.format(math_fsum - true_sum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More complicated example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0\n",
      "0.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "test_list = [1, 1e20, 1, -1e20]\n",
    "print(math.fsum(test_list))\n",
    "print(jnp.sum(jnp.array(test_list, dtype=jnp.float32)))\n",
    "print(1 + 1e20 + 1 - 1e20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of floating-point \n",
    "- You should be really careful with floating point numbers, since it may give you incorrect answers due to rounding-off errors.\n",
    "\n",
    "- For many standard algorithms, the stability is well-understood and problems can be easily detected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Vectors\n",
    "- In NLA we typically work not with **numbers**, but with **vectors**\n",
    "- Recall that a vector in a fixed basis of size $n$ can be represented as a 1D array with $n$ numbers \n",
    "- Typically, it is considered as an $n \\times 1$ matrix (**column vector**)\n",
    "\n",
    "**Example:** \n",
    "Polynomials with degree $\\leq n$ form a linear space. \n",
    "Polynomial $ x^3 - 2x^2 + 1$ can be considered as a vector $\\begin{bmatrix}1 \\\\ -2 \\\\ 0 \\\\ 1\\end{bmatrix}$ in the basis $\\{x^3, x^2, x, 1\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Vector norm\n",
    "\n",
    "- Vectors typically provide an (approximate) description of a physical (or some other) object \n",
    "\n",
    "- One of the main question is **how accurate** the approximation is (1%, 10%)\n",
    "\n",
    "- What is an acceptable representation, of course, depends on the particular applications. For example:\n",
    "    - In partial differential equations accuracies $10^{-5} - 10^{-10}$ are the typical case\n",
    "    - In data-based applications sometimes an error of $80\\%$ is ok, since the interesting signal is corrupted by a huge noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Distances and norms\n",
    "\n",
    "- Norm is a **qualitative measure of smallness of a vector** and is typically denoted as $\\Vert x \\Vert$.\n",
    "\n",
    "The norm should satisfy certain properties:\n",
    "\n",
    "- $\\Vert \\alpha x \\Vert = |\\alpha| \\Vert x \\Vert$\n",
    "- $\\Vert x + y \\Vert \\leq \\Vert x \\Vert + \\Vert y \\Vert$ (triangle inequality)\n",
    "- If $\\Vert x \\Vert = 0$ then $x = 0$\n",
    "\n",
    "The distance between two vectors is then defined as\n",
    "\n",
    "$$ d(x, y) = \\Vert x - y \\Vert. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Standard norms\n",
    "The most well-known and widely used norm is **euclidean norm**:\n",
    "\n",
    "$$\\Vert x \\Vert_2 = \\sqrt{\\sum_{i=1}^n |x_i|^2},$$\n",
    "\n",
    "which corresponds to the distance in our real life. If the vectors have complex elements, we use their modulus."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $p$-norm\n",
    "Euclidean norm, or $2$-norm, is a subclass of an important class of $p$-norms:\n",
    "\n",
    "$$ \\Vert x \\Vert_p = \\Big(\\sum_{i=1}^n |x_i|^p\\Big)^{1/p}. $$\n",
    "\n",
    "There are two very important special cases:\n",
    "- Infinity norm, or Chebyshev norm is defined as the element of the maximal absolute value: \n",
    "\n",
    "$$ \\Vert x \\Vert_{\\infty} = \\max_i | x_i| $$\n",
    "\n",
    "<img src=\"chebyshev.jpeg\">\n",
    "\n",
    "- $L_1$ norm (or **Manhattan distance**) which is defined as the sum of modules of the elements of $x$: \n",
    "\n",
    "$$ \\Vert x \\Vert_1 = \\sum_i |x_i| $$\n",
    "  \n",
    "<img src=\"manhattan.jpeg\">  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We will give examples where $L_1$ norm is very important: it all relates to the **compressed sensing** methods \n",
    "that emerged in the mid-00s as one of the most popular research topics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Equivalence of the norms\n",
    "All norms are equivalent in the sense that\n",
    "\n",
    "$$ C_1 \\Vert x \\Vert_* \\leq  \\Vert x \\Vert_{**} \\leq C_2 \\Vert x \\Vert_* $$  \n",
    "\n",
    "for some positive constants $C_1(n), C_2(n)$, $x \\in \\mathbb{R}^n$ for any pairs of norms $\\Vert \\cdot \\Vert_*$ and $\\Vert \\cdot \\Vert_{**}$. The equivalence of the norms basically means that if the vector is small in one norm, it is small in another norm. However, the constants can be large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Computing norms in Python\n",
    "\n",
    "The NumPy package has all you need for computing norms: ```np.linalg.norm``` function.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative error in L1 norm: 0.0008608277121789923\n",
      "Relative error in L2 norm: 0.0010668749128008221\n",
      "Relative error in Chebyshev norm: 0.0025285461541625647\n"
     ]
    }
   ],
   "source": [
    "n = 100\n",
    "a = jnp.ones(n)\n",
    "b = a + 1e-3 * jax.random.normal(jax.random.PRNGKey(0), (n,))\n",
    "print('Relative error in L1 norm:', jnp.linalg.norm(a - b, 1) / jnp.linalg.norm(b, 1))\n",
    "print('Relative error in L2 norm:', jnp.linalg.norm(a - b) / jnp.linalg.norm(b))\n",
    "print('Relative error in Chebyshev norm:', jnp.linalg.norm(a - b, jnp.inf) / jnp.linalg.norm(b, jnp.inf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Unit disks in different norms\n",
    "\n",
    "- A unit disk is a set of point such that $\\Vert x \\Vert \\leq 1$\n",
    "- For the euclidean norm a unit disk is a usual disk\n",
    "- For other norms unit disks look very different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Unit disk in the p-th norm, $p=0.5$')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6cklEQVR4nO29e5gc1Xnn/3m7Z0ZGtpAGcdFlNBIyIIOE15YGIRbHxkb8FmUxSiQwt/xinhjL3ti7683PWROTaFkl2cXZxxuSX7RxhOL1JUZgkGxkFmIjLGzjZYQ0WjtoIBJi0EijO2IkC0swmu53/6iq1umaqurq6Z5LT7+f55lnuqtOVZ3qmf6eU+95L6KqGIZhGGOfzEh3wDAMwxgeTPANwzDqBBN8wzCMOsEE3zAMo04wwTcMw6gTTPANwzDqBBN8wzCMOsEE3zAMo04wwa8TRKRTRK6r4PhviMiflXMuEdkjIouHo3+Dve5wMNr6Y9QvJvg1goioiFwS2na/iPxDmuNVda6qPucfV5EAueeqFoM952gT09HWn1pBRM4Tke+JyK9FpFtE7kxo+5yIvC0ib/k/O4ezr7WMCb5hjEFEpGGk+1Amq4E+4CLgLuBvRWRuQvvPq+p7/J85w9LDMYAJ/hjBn1l+UUT+SUROiMijIvKu0P7FIvJtoBX4gT87+o8x5/ugiGwXkZMi8igw4FzO+y+JyH6/7U4RuT7ifJeLyOsickdC/xenuRfnmKR7+UDCZzFNRNaLyFG/T/+uxOf6RyLysoj0isj/jOpLJf2JuF7S3/Fyf4Z73DeD3Rw69ksi8k/Ar0Wkwd/2h/75fi0ify8iF4nI0/7fa5OINMfdf6hvd4nIC36fDorIPhFZkubYEud9N7Ac+BNVfUtVnwc2Av9vpec2Qqiq/dTAD6DAJaFt9wP/4L/eA7wITAPOA14BPuu03QMsDr+OuVYT0A38B6ARuAU4A/xZxLnmAPuAaf77WcB73XbAfGAvcFPCNcP9i72XuONC2yKPx5vkdAAr/fucDXQB/yrh/DuAGf65fh58DtXoT8zxcX1vBHYDX/b7/jHgJDDHOfYXfl/Pcba1482cpwNHgO3AB/EG8R8D/ynl/+ADwGlgmf85fhHojmj3JHA85ufJiPYfBE6Ftn0R+EFMP54DjgJv+H+P60b6+1krPzbDH1v8taoeUNU3gR8AHxjkeRbhicuDqnpGVR8Htsa0zQHjgCtEpFFV96jqa87+38Cbrf2uqj5ZRh8qvZe4468CLlDVVarap6pdwEPA7Qnn+htV3eef68+ByKeUQfannLaLgPcAD/h9/zGeuN4ROnafqp52tv3/qnpYVfcDPwO2qOr/UdW3ge/hCW4argT+UlU3qGoe+BbQGn5aUdWbVHVSzM9NEed9D/Cr0LYTwISYfnwJb6CeDqzBe6J6b8p7qGtM8GuHHJ4IuzTizbwDDjmvT+F9kQbDNGC/+tMpn+6ohqq6G/gC3tPGERF5RESmOU0+C/xvLX9BttJ7iTt+JjDNN4kcF5HjeDPmixLOtc953Y33+QQmjmDh8OlB9qecttOAfb7Yuv2ZHtPXgMPO69MR79N+tlcCjzvvLwTe8geOSngLODe07Vy8p5cBqOoWVT2pqu+o6jfxZvm/WWEf6gIT/NphL565xOViYoS4BKWKIBwEpouIONtaY0+m+rCqfghPTBX4irP7s3izwL8cRD/TUG5Bh33A66FZ5wRVTRKMGc7rVuAAgKp+R88uHAa27KEsMHEAmCEi7ve2FdjvvB+S64vIJLzP4aiz+RZgwEDnrw+8FfMTNTDuAhpE5FJn278AOlN2TwEp2cowwa8hHgX+WERaRCTjL3B+nOIZV1oO4z0Sx/EC0A/8OxFpFJFlwMKohiIyR0Q+JiLjgLfxZozuDPQkcCPwYRF5YBB9LUWpewnzInDSX9w8R0SyIjJPRK5KOOZz/ud+HnAf3t+iWv0phy14M/7/6P9drsP7H3ikWhcQL97iGxG7rsR7yrzTXwz+18Dv4z3ZFaGqS5yBMPwzYJFXVX8NbABWici7ReRaYCnw7Yj+TRKRfyUi7/L7cRfwYeAfK7jtusEEv3ZYBfxv4HmgF/gL4C5V3TGIc/1XvMHjuIh8MbxTVfvwFubuBt4EbsP7QkYxDm8x7w08U8SFwB+FznccuAFYIiJ/Ooj+JpF4L2FUNQfchGcXfx2v32uBiQmHPQz8CG9x9zXgz6rVn3Lw/y4fB5bg9ft/4K2N/HMVLzMDz0QS5krgO8A1eP9//xn4LVV9uUrX/X3gHLxF5XXAv1HVTig8MXzZb9eI9/kHi7b/1u/Hrir1Y0wjxWZawzBcRGQPcI+qbhrpvgw1ItIE/BJ4v6qeCe37W2CXqg6Vac4YBmyGbxgG4D1BqOrlYbH3uRLPRdSoYUzwDcNIwzygmqYjYwQwk45hGEadYDN8wzCMOmHUJlg6//zzddasWSPdDcMwjJqio6PjDVW9IGrfqBX8WbNmsW3btpHuhmEYRk0hIrHBmGbSMQzDqBNM8A3DMOoEE3zDMIw6wQTfMAyjTjDBNwzDqBNM8A3DMOoEE3zDKJOO7l5Wb95NR3fvSHfFMMqiKn74InIj8FdAFlirqg+E9rcC3wQm+W3uVdWnqnFtwxhOOrp7uWttO339eZoaMnznnkUsmJmqBrhhjDgVz/BFJAusxsvRfQVwh4hcEWr2x8B3VfWDeLVD/0el1zWMkaC96xh9/XnyCmf687R3HRvpLhlGaqph0lkI7FbVLr9AwyN41WpclLM1Kyfil4gzjFpj0ezJNDVkyAo0NmRYNHvySHfJMFJTDZPOdIoLJ/cAV4fa3A/8SET+LfBuYHHUiURkBbACoLU1toSqYYwYC2Y28517FtHedYxFsyebOceoKYYrl84dwDdU9asicg3wbRGZp6pu7VNUdQ2wBqCtrc3yNhujkgUzm03ojZqkGiad/Xh1MANa/G0unwK+C6CqLwDvAs6vwrUNwzCMlFRD8LcCl4rIxX5NzNuBjaE2e4HrAUTkcjzBP1qFaxvGsGIumUYtU7FJR1X7ReTzwA/xXC6/rqqdIrIK2KaqG4H/D3hIRP4D3gLu3Wqltowaw1wyjVqnKjZ836f+qdC2lc7rl4Frq3EtwxgpNmzv4Z0zeZSzLpkm+EYtYZG2hpGCju5eHtu2j+CxNJs1l0yj9jDBN4wUtHcdoz/vyb0Atyxosdm9UXOY4BtGCtyAq3GNGZbPbxnpLhlG2YzamraGMZqwgCtjLGCCbxgpsYAro9Yxk45hDALzxzdqEZvhG0aZmD++UavYDN8wysRSJBu1igm+YZRJOEVy8/gmM+8YNYGZdAyjTFyPnebxTax6stPMO0ZNYDN8wxgEC2Y287mPXkLvqT4z7xg1gwm+YVSAVcAyagkz6RhGBVhAllFLmOAbRoVYQJZRK5hJxzAMo04wwTcMw6gTqiL4InKjiOwUkd0icm9Mm0+IyMsi0ikiD1fjuoZhGEZ6Krbhi0gWWA3cAPQAW0Vko1/lKmhzKfBHwLWq2isiF1Z6XcMYbXR099rirTGqqcai7UJgt6p2AYjII8BS4GWnzaeB1araC6CqR6pwXcMYNVh+HaMWqIZJZzqwz3nf429zuQy4TER+LiLtInJj1IlEZIWIbBORbUePHq1C1wxjeLD8OkYtMFyLtg3ApcB1wB3AQyIyKdxIVdeoapuqtl1wwQXD1DXDSEdSSmQLwDJqgWqYdPYDM5z3Lf42lx5gi6qeAV4XkV14A8DWKlzfMIaMwC5fKmdOEIC1YXtPodC5YYw2qiH4W4FLReRiPKG/Hbgz1Ob7eDP7/yki5+OZeLqqcG3DGDJcu3xGhFxeUc6abKJs9Ou399DXn2fD9h6z4xujjopNOqraD3we+CHwCvBdVe0UkVUicrPf7IfAMRF5GdgM/KGqmpHTGNW4dvm8KtmMJJpszI5vjHaqklpBVZ8CngptW+m8VuAP/B/DqAkCu/yZ/jyNDRlW3jSX3lN9sW6X4fZmxzdGG+Jp8eijra1Nt23bNtLdMOqccn3rzRffGGlEpENV26L2WfI0w0jATYyWRswtkZoxmjHBN4wUhAOrSpl3DGM0YoJvGClwF2T7+vOsfGIHeVWLqjVqCsuWaRgpWDR7Mg0ZQQABcnklr/DOGc8F0zBqARN8w0jBzkMnySkokMkI2awA3vvHtu2LjL41jNGGCb5hlKCju5eVT+wgl/c82nJ5Zdbkdxf25/JqPvdGTWCCbxglaO86VhB7ABFhz7FfF95nM2I+90ZNYIJvGCVYNHsy4xozZICGjHD9+y4k7w8AAtzaNsMWbY2awLx0DKMEQWK0wAcf4KevHi1E1C6b3zLCPTSMdJjgG0aIqACrcEDVypvm8vSOgyyZNxWA1Zt3m0++MeoxwTcMhzSVqzq6ewupkrd0HQMR+nNW6coY/ZgN3zAc0mS8LGqTU85YhkyjRjDBNwyHUpWrOrp7+eW+4wjel6cxKzSWqHSVVCnLMIYTM+kYhkN4gdY1z3R093LHmhfoy3keOtmM8HvXXsyv3ulHgGXzW1gws7loDQDgjofaCwu893/ccvAYI0dVBN8vSv5XQBZYq6oPxLRbDjwOXKWqlvvYGJXEZbxs7zrGmdxZf/x8Xln7/OuFnDrL5rcU1gDeOZMnmxHmt06irz8PeDl47vv+S6AwrtHs/cbwU7FJR0SywGpgCXAFcIeIXBHRbgLw74EtlV7TMKpNGrNL8/imonq1Imdz6gT2+/auY7xzJo8C/Xll657i86mfnqHvjNn7jeGnGjP8hcBuVe0CEJFHgKXAy6F2fwp8BfjDKlzTMKpGkmeOa57pPdWH4Am2AIsvv6jIHz8w4WQzQr8TmZsVyIXqDGUsOtcYAaoh+NOBfc77HuBqt4GIzAdmqOr/EhETfGNUEeWZE9jiwznwxzWeLWH4mY+8l8985L0D7P2rls4rSp+88qa57Dhwgsc7eujPeQXRVy2dZ+YcY9gZ8kVbEckA/x24O0XbFcAKgNbW1qHtmGH4xNWiDQ8Evaf6Ihd0w8J959WtzJkyYUC75fNbrPyhMaJUQ/D3AzOc9y3+toAJwDzgOREBmAJsFJGbwwu3qroGWANeTdsq9M0wBhCOpI3zzKmkKHnUwq+VPzRGmoqLmItIA7ALuB5P6LcCd6pqZ0z754AvlvLSsSLmxlCQJpI23N5cLI1aYkiLmKtqv4h8Hvghnlvm11W1U0RWAdtUdWOl1zCMahFnr4/DnZXf972Xilws/+SJHaiVOTRqiKrY8FX1KeCp0LaVMW2vq8Y1DWMwuGaabDbD/uOn6ejuHRBgFWVrDz8LBzny0wwchjEasNQKRl0R2OtvX9gKqjzy4l7uWtte8L8PTD5f/dHOou3gLbo2Zb26tg0ZCikY3IHDxVIqGKMNS61g1B0LZjbT3nWM/lAh8mB7nMlnwcxm1q24psimv2F7D49t28cjL+7l8W37uLVtBsvmt7Dz0Mki18wkk0/cE4VhVBsTfKMuWTR7Mg3ZDH39+UIh8mXzW0p65oQ9bdyBoy+nfGfLXh55cS8KBLFXfU4WTVfYO7p7CwNGf97WAoyhxwTfqEsWzGzmlgUtrNviiXN/Tnlw0y6+sPiySBfNuFn4otmTachIIaEaRETVitA8vmlAENeqJzsLaRjA1gKMocds+Ebdsnx+S6FWbR74+e43uOOhdtZv7xkg9nF2/QUzm7m1bUb0BfDSLHz0fRfSeeBEkano6R0HC08X4KVqKNfX3zDKxQTfGPPELZ4GC7jXXno+GfFMMH39edZtKV7ILVUUZdn8Ft7VmEH89wI0ZIVLLnwPIvDsK4d5bNs+GrJn8+YvmTe1sOjb1JDhjqtbzZxjDDlm0jHGNKUCrRbMbOYLiy9j6543C+YVxVvIXe8v5Kax6wdmoObxTXQeOMFj2/bx2pG3CjP4XF65beGMQvK1OVMmxObdN4yhwgTfGNOkCbQKBHvD9h4e3bqX/rwnyo939LDcL2pSSpzdxdzVm3fTn9cB5pp50yYWauFu2N7Dd+5ZxOc+esnQfgCG4WCCb4xp0ubDcQX7YX8hN5c7O0DsPHSyMIMvNRsvCu7KCNfNuZALJoxjR8iObwu0xnBjgm+MaRbMbGblTXN5esdBlsybWlJgl81vYf32nqIB4oGnXuFrP+0C4GevvgF4GTGTrumaeIJZfUM2Q0a8IihZJx9+2API/PKNocIE3xjTdHT3FgR36543mTNlQskAqJU3nU2KBrDmZ11F7Z7ecTBR8OHsE8PqzbsLs/r+XB4vY6x65bIoXmNoyGb4yGUX8JNdR+nPpUvuZhjlYIJvjGnSJkuLW9xdvXk3+ZBf/ZJ5U1Nf3zXviAg537YfmIuAQv/6+vM88/LhwrFm9jGqjbllGmOaQHADd8g4G36c62Xz+Kaidr/1gWnMmTIhdY6cwLzzB//PHFYtnUdjg/+V84Oxgv5JxLHZrPnlG9XFZvjGmCaNhw3EL+72nuor+OhnBN49rqGifPr5vJdeOZdX/uSJHfzp0nkFD6HHtu2jP6fk8Tx7cGpVmF3fqAYm+MaYI6qiVSmRTFv1SiG1p03YTLR8fgu5/Nn9ubyy8okdPPqZa/jz376SZfNbeHDTLn6++w3y6u0PnjTKGWQMIw4TfGNMUW5FK5e4soTuQABehsw0ZQ/DZiIFGrPFeXfyqoVBww0Cc89fbtEWw4ijKoIvIjcCf4VX8Wqtqj4Q2v8HwD1AP3AU+D1V7a7GtQ3DZSjEMTwQpI2QDT8dLJ/fwvL5LXztJ6/x438+UqiW5Q4acU8ag62taxguFQu+iGSB1cANQA+wVUQ2qurLTrP/A7Sp6ikR+TfAXwC3VXptwwhTSeHxgGray5fPb0H938G5HvrdtsRrhAeYtOsQhlGKaszwFwK7VbULQEQeAZYCBcFX1c1O+3bgd6pwXcMYQLniGOSkD0QZGJDG2C1UHmUyAgZcL8p+H+5nOcJdbnvDiKIagj8d2Oe87wGuTmj/KeDpqB0isgJYAdDamhzYYhhxlBLHYHbdPL6J+3/QWShMHlSscv3iw1Wrwiaj9dt72LC9Z8AAsf/46aqaltw+uwOQuy94mrEnASOOYV20FZHfAdqAj0TtV9U1wBqAtra2cM1ow6gYd+adEaHfiao6k/OCoqICpQLRDpuMBCdw6szZAaIhm6Eh4x3vmpaSTDlx+9w+59Vz2RzXePbpohCpmxEQsShdI5ZqCP5+wK0A0eJvK0JEFgP3AR9R1XeqcF3DKBt3ho4qWV+UwfOgCRZWw3lwxA+UivLaCXLviAh59cod5nJ5bl/YyrRJ59A8von2rmPsPHSycL6wICd5FxX1GYoGIHDcRHNecudg/4btPTbbN4qohuBvBS4VkYvxhP524E63gYh8EPg74EZVPVKFaxpG2XR09/LLfccRvBDzpkbPBNN54MSAhdXg995jv2bNz7rozyv3b9xB54ETLJvfUpTWOJwoLZj9LwutCWQinhiiRD28L6i/G7h2ZiiOGnYzcyJCLpcnm81YrVxjABULvqr2i8jngR/iuWV+XVU7RWQVsE1VNwL/DXgP8JiXPIq9qnpzpdc2xg5DbYfu6O7ljjUvFHzgsxlh5U1zByRBC/dj7fOvny1GnlMe3rKX9X4ue3dwCF7PmTKhqO9u8jTwnihUdYAHUWAqcp8mAnYeOknOj9hqyAqfaJtRNDiFnzjau46x//hpHnlxr/nuG0VUxYavqk8BT4W2rXReL67GdYyxSVHGyArs0En28fauY77JwyOfV3pP9Q04/o41L3AmpzRmhVvbZpDX4qUkxbPVBwXPS7lUhm3+Ya8f97iVN81l5RM7yOWVVU92MmfKBABvm57t9/RJ5ww4NtyP9dt7aMhmyOXMd984i0XaGsNKlLdJkTkjZIdOOzMtFWG7aPbkoijXKBFcv72nsL8vpxw5+U5h1h2gerbg+dY9b5YckEq5ibqfx9M7DpJXHWCjzzkLyxmRRPEOD563L2xlmfM0YDl56hsTfGPYiPM2WXnT3Eg7dDkz0yQbeCBy9988j807j3DkV29z21WeKWf15t0F8QtnrLxwwriixGZnckqQzr4cU0l4Bu6K/KonOwu1dIN6t5lQZs9xjRn6zuTJZIRVS+clXs/9HHJ5ZZrzNPDwlr0D3ExN9OsLE3xj2IjzNuk91Rdphy5nFhoXYRsuMIIq/XnllUOdhdeB+C2b38JjHT1Fi64LZjbT3nWsUKNW/GpVgR2+eXxT0aBRirBbaLCIG3weGeDaS84vMheVE0iW9DmsfGJHwQ21z+z6dYkJvjEkRJkOAm8T10QSlPqLSidQDnGmk/DMHyhkvAwIZuqf++glrPt06YyZgR3eddtMO2MuHvS8Rdx83kuJnBHP4ya8NlBOlG3S51DKNGTmnrGPCb5RdeLs6QtmNnPLghbW+UXCBbi1bUZV7csHjp9m/fYewBO/ooLi2Qx5Vfqdxduwi2OajJnBftcDJ615x/XGyYhwz4cuZsI5jZERtGHSfj5R97Bo9uRE01AlWUaN2sEE36g6Sfb05fNbitILB77qlQpOR3cvdzzUXpQmYd2KawaI9d/95DV+5JQRvLJlIis/PrdsOzyUl6jNFWvXG+cbL+xJda/lfj5RNQGSTEOWgrk+MME3qk6SEKY1vZQrOO1dxwomG/C8fdw888G5zp8wrui4udMnDlrYAlfKp3ccZMm8qZHnCZKzuUFQy+a3FLxx+vrjXTzD9zfYwivuE1aQAC687lCNLKPG6McE36gqwcwyzt8cKp8tR7Fo9mQaHRfKxmy0++Ly+S087nvcZLNS6HOaGXZURszAhr91z5vMmTKhIKhRXjjgibVAwayTV3j+1dIunmk/n47uXh7ctKtwzSiPpbjBwI0YDlxCbZY/tjDBN6rGYM0ywSz4w5dewPkTxhVFkZZjt1736UWFVMfzpk2MFK0FM5tZt+Kawqz7kRf3siEUOZv2vuIKnwdtBc99MxB7gYIZKyhn+Pyrb6SKOShlkgn3MyoFAyQ/KQS/zZY/djHBN6pGWrNDOH2Ba3tv8hOYBe3KER/XZJF0nOtqmcZEsmF7T+SMOWrWHXY9DYhKiRBVzjCJUt467rUzUuzeGQyqR06+kxiBa7b8sY0JvlE10pgdAjF+54wXZPWx911YZHvvyykbtvfEzqDTLFSmyUVfjonksW37CrP0bLbYoyeuHKFrxhHgtrYZ/PlvXzmgr0mmr3IJ35Mr9u6g2pBhQARuuZ9LHObaOboxwTeqRhqzQ3vXMd4+4wlPf1559pXDZJwUxQCPbt3LsvktZXvBuAFWaXLRh4O9ogKogicB8IT7lgXx5qbg/d3XzOKFrmO8fPBX5P0+LHMqXoVz9qxbcU3s9cshaUHcHVRzeZg26ZzIa6b5G4Y/d/czNHPQ6MYE36gqrsi47+FsemKXvML7p5/LL3tOFLbl8hQCocLiEyey7qy+vz/PlS0TmTt9YsGEkrRYGbevo7uXA8dPFw0ey+e3RHrerLxpLvdv3FEooiJ4C8dRM+lwzp6v/eQ1fvbq0YqE0v1c3NTNEL2g3Ty+KVack0xH7nV2HjpZlKph+fwWMweNckzw65S0j97lpi1OEs/AlOOSzQi3XdXKK4c6S3rYhM+98qa5BS+ZYFbfn/OiVl/af4Kdh08W1gOSzEOlFl8bshluWzijqOata7LpO5Pnr5/dVRBx8BZqw7lsgnt40T9/wOtH36pIKMMJ025tm1E0yIQXtIMCL+VeM/wUlcvlC1k8g4XiuLQOZuYZHZjg1yGlFjXDboXhPDRJs9CwkKz3qy4d8Gfgrl0760d83nl1K3OmTCgSpKiZd3gG+fSOg2cThfkVpva+eYqf735jgHjvP366aLHSzYFTavE1l8sXUhIH0bXummweOPSr4iJugUeOO3AF9/N2aNA7791N7Dl2qlCB68Dx0zy8ZW+kbT9KPN2+psnZH+BGIO8/frqke2r4b+t+BhnxFtvnTZtYiEsA+PL3XuLxjp6S6a5tUBgeTPDrkFKZJeMqNAElXQiLUhlkpPBlb8hIQXCz2Qy3LGgZUGGqVPRneAa5ZN7UIi+XwE7ubnNNF0G64LnTJg7IgRO3+BqerYZTNVw+ZQIv7T9R5JUT5NJfHjLlBPfj0pgVfrHvOHlVJCPkgXV+4RLhbO3aJNNT8/gmMuIldNMUf6Pg83YzgaZxTw3fO36aiiBVA1D4XLe8/ib5fB73dt0BOGyWu+Oh9sJnve7TyS6yNjAMnqoIvojcCPwVXsWrtar6QGj/OOBbwALgGHCbqu6pxrXHKkP5jx2e0bqzXVdkAwEBT0QasoLmNXFG6C76HTh+uiBeZ3LKDVdcyL+YMSn1PYX76dabDc4RrjAFxdklo9IF957qGzDghU1IcYuXUYu9d61tLwxwYXNK3P0EbRUKlanE+byDz9wV7jjT06onO8nllYxAxk/GltbNsxz31AUziyOLgcLrO69uLcotFB7YgieeqLWDDdt7Cu37+vMFL60wlu+ncioWfBHJAquBG4AeYKuIbFTVl51mnwJ6VfUSEbkd+ApwW6XXHqukNbkMdjAIR1W6s103Nz1S7D3zsTkXcsGEcSVnhO5i6GMdPQUTyHM7j/CZj7w3dZ+TRDd8rajrB0TN1MMDXvjzhoHrFa6pq1Qf095PsAB8pt9LauZW5QqbhZJMT8FRt7bNYPqkcwY9qJbyhCrM4LuOFSqTBVHG7rnC/zvv93MWRQ1aoZCFAe8D0qw72BNAMtWY4S8EdqtqF4CIPAIsBVzBXwrc779+HPgbERFVjfvb1jVpTS6VzHICUQxnfHRz0/9i33GecRKNnT9hHNMmnVPWjNDNjpnLa+QjfZp+xlHqCx4nyHFPAWf8GeZ6f9bpDgDuYm2Qytj19klDuK3bP7cOrQAfurQ4L36S37/7FFTO/0M5A1bR5xRRmcz1qmoe38T9P+gsDGS3XdWaaCoL0l00OoF3YZKeTJNMXsZZqiH404F9zvse4Oq4Nn7R8xPAZOANt5GIrABWALS2tlKvJM264hZFBzujibqWO0P/yc4jA76I5QTmhLNjJrkDBpTjQZTmC17uU0CQ1CxsOnE9c9ImMUs78w/P9t3AqaR+lxLsNH2IG7A6untZv70HgQFxESKgeCa+4P8g6lrhGr1x/V234ppU/Yx7Mg22m1toMqNq0VZV1wBrANra2up29p/0JY5bFB3sjCbJzKDA/TfPG+AtUk5gTvj8pb6U5czSqvEFj4p4BYqEd9HsyTzTeajI1BDlhRM+b9x9xIlwObNtlyTBLnfG65qt7t+4o+Bq+lhHD+s+7S1ur9/eU/i/y2SElTfNLfjk9+eVBn8Rt/dU34AavUFfSw1kpe41qhZBpVHC9UA1BH8/MMN53+Jvi2rTIyINwES8xVsjhrgvQOyiaILguQIe9cjvXischt/kR4IGA0FccE8595L0pSxHxCv9gifVeA0vzK59/vWiYwNzC8RH6EbdR9I1q2F/ds9R7oAY9tDqd2zwrtmmvesY/bmgLrHSeeAEj27dV2jfn1dWPrGDVUvnJZpgKiHuyXQwA2Y9UQ3B3wpcKiIX4wn77cCdoTYbgU8CLwC3AD82+/3gcU0A60Mz0TBJhUGiCIfhB3nloTph86W+lOXkuAnPzMvpT6kar+4gtXrzbvLOv2tDRgpiH/eZRN1H3DWBolnzYAeCqMC0cgZEd4AAz+sn0PykxWOFos8H/727HjSYcpBJJC3om9DHU7Hg+zb5zwM/xHPL/LqqdorIKmCbqm4E/h74tojsBt7EGxSMCkkzo4kT8KRZc1Re+bSzxUpsxmnvKc5U8fCWvUVugkm0d5Wu8ep+JkH+ehHho++7sHCOpFTD4ftYvXn3gGsGaxrhnPlpBoKoe4pbgE8zIIb97POq5HNeMNj9TlWw8L1BcUZRARp8112Az330kkGVgyyFiXv5VMWGr6pPAU+Ftq10Xr8N3FqNaxnFlPqnjxPwpPOFw/CTvCtcomaYg5l9l7qnKKHdeegkX/7eSwD87FXPFyBJ9BfNTq7xGu6PG6T07CuH+dmrR7n7mllkRECj/d7D9xF1zSAmwI1Adhe3owaCckxc5YiiK+SBxxAA/mw96d6C406ePsMLXcfoPHCiyHXX7Oujg1G1aGtUnyQBTzomqs3y+S2J53CFuK8/H2urrpQo8Xhw066iNk/vOJgo+ElPElFPKQtmFgcp9fXnWfv86+RVCwuX5S5gB2a5cDDWMifXTXggKDVYV2rDds2F4YXrUsfBwDxDUS6bZl8fOUzw64BKH32jctoE290vcLHbXnFahg2++2jz+KaK879HCduSeVMLM3ugEAla6jzhPoQXLoNcP0Ds/QkDZ8Bpr5kk0lEDwWAH6yiSzG+DGTxKDVJmghl5TPCNksSF9EfZ0YPQ+7lTz+UbL+zxTEkiPLp1L7k8BRuvmyOmFHEzbvfYQJSDa/ee6ktVqzZ8nVU/6CwkN8ur520S+I+7IhgsQlbDRBHnpjgUM2K38tVPdh1NXBsoV6DDawDhfEnGyGOCb5QkKaQ/Kq9LX78Xbn/3NbNY+/zrRe59QCGwKbBHhwN8wuaVtN5BQdbNwaSlCIqSuCmOAfL54kVuVwTnTJlQ6PdQUO0ZcdhjK6DvTPUWUc1sM7oxwTdKEvdFThoE+vrz/GPnoQHuegHie8WEhTYI8AmuEZX2IGwDd98X9eFMngc37Soq9Rc3GLR3HSvKYwPek0hTY/LsPUj8FU5HPByUiq8IE/bYCshDUX6gSjCzzejGBN9IRVrbc+C+mFe8HO94eWeyGSHnFycBEH9aHBbasCdK2EwQrjIV9u0uuFCeyZMHfr77DbbueTMx9D6obJXNUEjnG5fi2MU93ztnPNGvhv087fHlxFfAQI+tAMEzhQWmK2PsYoJvDJq4QeDBTbt4/tU3vCRjwLWXeFGpG7b38LCfSE19U8mi2ZNpzEphhh+2h0e5CkYVQAl7gzy4adeAQijuYCC+D3y4itMN77uACyeMi10gdYV60ezJNGS8vivweEdPqpl21JMGpE8oF7QtJ74i+CwDj62CDb9/4MBooj92McE3qsqCmc18YfFlRUVIApPKzkMnyUbka7+1bQavHj7JO/35oqyK7jmDmXiQVTGbkQEFUFxvkHAfAiFdedPcgrvo/Rt3cPnUc4sqW31gxqTYtBFRpQSvm3Mhz7x82MsGmktnC48yU4Wzc5Y6R7nxFeHPMrifqIHRBH/sYoJvVJ04f/NVT3YW+a3DWU+fICXwzsOdyaYFEUBBhDlTJsQuEroeQ0vmTS3sCxJ6BeUAg+LpGSnt5160PpBTvrNlL9mMkM0WZ40sRdhMtWP/iaK1D3fdIekzLhVfkWQ2CvbFDZrG2MQE3xgSwuYeVywDv/Xi3C2lS/MFSbsU6M+dFcaoGblbrCMo0OHGCrjBQa7ZKUlkm8c3EXI4IpdXsgK3L2xN5ScffDZu5O4/9ZwouKvm1YsU3tJ1rKRNPmmBNCnqGahKRLRRe5jgG8NCXGi9u8ibIXmW7ea0ySs8/2q83TlugTaY+T+6dS8vH/xVwbxUSuwBOg+ciNyuCtP8IuflsPfNU5zx7f8Z4MJzxxWKofflNHYROM2Cb1LU87JQMfjeU31lZT81ahcTfGNYiPPqcQOZSs0yoxaF454I4rJVBrPqIG972pn5w1v2ns0tA2QzgHpiX8p1M4w7+w7Evqkxw/tbJnHIqTAW5d+fNi4hKepZKK+IjTF2MME3ho24iNJyZsZxC7IuQdbMu6+ZxYRzGovMGK4pJyhqHuXP7xKkNQ68RwW47arWAQXV0+LOvjNy1pwE8Nyuo4X7WuaU+gv6d+D46VRZJ5OigpfNbynk6zEzTn1hgm/UHEkRnQ9v2VuUNfO//PaVADy4aVdsnpegKEkur5EpH8KplLMZKSySDkYsw08frjlp3acH3lfYO6ghmyGXKz07D0cFR6WnMOoLE3yjJgjPwOPE9ukdB4vef/35Ll4/dqog2BnxcrUHeV6A4qIkZ87WCQ7MTM3jm1KnUo7qa9T2JO+ipPWIXF65beEMpk86p6zZuUXAGmCCb9QA5eTTCWfN7Hrj10WeNVdOn8hKp5hHuJqVCIWCI4Gr6LjG9J4scX2N2p5moTSIAm7IeHb4Rj9bqYm3MRgqEnwROQ94FJgF7AE+oaq9oTYfAP4WOBfIAX+uqo9Wcl2jviinNqubNfOcxiw/chZBAeZOn1h0rOv5kxHhY++7kE2vHB7gKprWkyWur0lpHaJMOOGygA3ZDLctLE71UI0auEZ9UekM/17gWVV9QETu9d9/KdTmFPC7qvqqiEwDOkTkh6p6vMJrG3VCudWS7ry6lTuvbqWju5fndh0tRKM2ZCiYcQLC6wEAP331aGpX0bR9jfMaikqxECwu+8W0ClG80x3Xz3KeegwjoFLBXwpc57/+JvAcIcFX1V3O6wMicgS4ADhe4bWNOiFpkbbUcWmqfYXt2+W4iqbta9T2qDqvQMGTKLA0BcnnDhw/XcjxH/azt5QIRhoqFfyLVDVYJTsEXJTUWEQWAk3AazH7VwArAFpbk4tQG/XFYBcdSx3nmk92HDhRyMlfbiBS1KJyqb7EPQ1kM1JUQ+DK6RN55dBJ1r24t5CG2Y36zWv10hsbY5uSgi8im4ApEbvuc9+oqopIdPJz7zxTgW8Dn1TVgUm5vXOsAdYAtLW1xZ7LMKIIi24pG7drFnEXdsM5+dNcdzDmlbingXs+dDFf+2lXod1F576Ll/afGPAk4GcVIgOpSywa9U1JwVfVxXH7ROSwiExV1YO+oB+JaXcu8L+A+1S1fdC9NYwYAtF954xXA/aeD11cKLEYJ8LhXD4BabJGuoNJOYvKYaKeBiac01gk5udPGBf5JDCu0aJljfKo1KSzEfgk8ID/+4lwAxFpAr4HfEtVH6/weoZRICy6ge27P6/8nT9DVjyb+IZQXpqO7l72Hz9NQzZTcMEMiBPQKO+ZIPlYNVMVLJo9uUjMl89viYzqtXKCRrlUKvgPAN8VkU8B3cAnAESkDfisqt7jb/swMFlE7vaPu1tVf1HhtY06JiobpGv7diftCjy2bV8hZ044cvWOha3MnTaRzTuPcORXb3PbVd760erNu4vMQ8ExGZFCiuXAZTNJfMt1nwxMPW693OC4wJyTFHxmGHFUJPiqegy4PmL7NuAe//U/AP9QyXUMI0zYjNJ7qo9VS+fxx99/aYCJBrwI1Sif+CCfzpwpEwqz9lcOdYJqoZRiuDwifk5/QYuKq5RaJyinqDoU18uNKudoYm+Ui0XaGjVJlIdLIIBBKuCGbAZUCxGqST7x4QEEivPzh48JR97GCXeSfT9uMHh4y17W/PS1gokqrpyjCb5RLib4Rk0S5+Fy59WtRYnCYGCt2HAmyeC3W4UqPFAkxQIkpVMI1gmikp2t395TJOrtXcfYeehkIfkbnE3yZpWpjGpggm/ULGn93ePaQHzlJ4geKKLOFTWL33no5NknjVDe/Y7uXtb7efkD61M249WkfXDTrqJzz5w8nq9+4gMsmNkcmfHSMMrBBN+oW6LWAYKAq45uLyXUzkMnS0bdhs09zeObijJw9ofy7ofz8gteIfcFM5sHJH9b8eH3lhxwDCMtJvjGmCZpUTQu0jUckBX4xAeZM90F0+D87tNBOH9+RoTm8U2s3ry7UMDEFftxjWeLnbjJ35bMm1p4bxjVwATfGLOU8pCJs8u7tnWg6Hc402Xc+d38+fd86OKzWS+dAiYiMHfaRG67qrWoX0HyN8OoNib4xpglTQRs2EzS0d3L4x09RX78btSr+ySQVCjdHUiiCpiAl3f/pf0n2Hm4kzlTJpi5xhhyTPCNMUuatMphk0971zH6c2dTPWUEbl/YyrxpEwfY8ONSHgfncxOwue2CqNkgwtfcLI3hwgTfGLOUSqscZZJZNHtyIZIWAIXpk86JNLFE5dKPMvHE9SNpsDDxN4YCE3xjTJPk2RJlkvncRy9h1dJ5nktlXmlqTF8o3M1vH85RH+UqmmawMIxqYoJv1C1xJp9w8FaaVMtA2Tnq4wYLM/EYQ4UJvlF3uOIdZ/JxxThNvvuO7l6e3nGw8L7cHPXllnE0jMFggm/UFVHiHVXdqpx89x3dvdzxUHuhdm4GBpiCSj0hDLaMo2GUgwm+UVekcdWMSr0cF6DV3nWMX+47XhB7gCtbJrLy43PLekIAi6Q1hh4TfKOuiEqD4Oa9h+iUC+HZtyviqsXXmDt9YuqMmYYxnJjgG3VFOFNmVI75pNTLQWK0p3ccLIrGBT+zZVZY7qdJCGge30RGBFTNPm+MKBUJvoicBzwKzAL2AJ9Q1d6YtucCLwPfV9XPV3Jdw6iEQLwf3LQrVaRsVNKzIPo2QIAPXXo+X1h82QD7/qonO8n7RVNW3jTXZvfGiFHpDP9e4FlVfUBE7vXffymm7Z8CP63weoZRMWHxzsjAOrZhe3pglnHz6gggvvI3NWYixd4dVAQty3PHMKpNpYK/FLjOf/1N4DkiBF9EFgAXAf8ItFV4TcOoCFe8M8C1lwycmYcJzDx9Z/Lk8QaJcA79qKIoSYOKYQw3lQr+RaoaOB8fwhP1IkQkA3wV+B1gcdLJRGQFsAKgtdWyBRpDQ9hGX0rsYaDtPy43fsBgBhXDGGpKCr6IbAKmROy6z32jqioiEeWj+X3gKVXtEZHEa6nqGmANQFtbW9S5DKNiBuvzXo7b5GAGFcMYakoKvqrGzspF5LCITFXVgyIyFTgS0ewa4DdE5PeB9wBNIvKWqt476F4bRoVEiXc1k5dZIJUxGqnUpLMR+CTwgP/7iXADVb0reC0idwNtJvbGaCNtcFQ5WCCVMdrIVHj8A8ANIvIqnn3+AQARaRORtZV2zjCGi6jgqICO7l5Wb95dqHNrGLVKRTN8VT0GXB+xfRtwT8T2bwDfqOSahjEUxAVHDcXM3zBGikpn+IZR8yQFRyXN/A2j1jDBN+oeV9RVi4OjAm+brPnRG2MAy6Vj1D2uC2U2Ixw4fpqO7t7ENAtJ3jfB/jT++oYxnIiGU/2NEtra2nTbtm0j3Q2jTujo7mX99h4e7+ihP5dc6CTJpu/u99IpwLhGs/0bw4eIdKhqZEYDM+kYBp4L5fRJ59CfS7bXl7Lpu/vBy7ljtn9jtGCCbxg+aez1SW06uns5cPw0DRkh4weVZzDbvzF6MJOOYTikibaNauOachqyGW5Z0MK8aRPNhm8MO0kmHVu0NQyHNNGxUW1cU04ul2f6pHO482pLAGiMLsykYxhVwNw3jVrAZviGkZIkc48lSzNqARN8w4jBFXigZIoFS5ZmjHZM8A0jgrC//bL5LZH1bw2jljAbvmFEEPa3FzAbvVHz2AzfMCIIV6xaNr+FZfNbWL+9h+S6bYYxejHBN4wI4nLobNjeQ19/nvXbeyxdglFzVCT4InIe8CgwC9gDfEJVB1SJEJFWYC0wAy/a/DdVdU8l1zaMoSa8CBuVVsEE36glKrXh3ws8q6qXAs/676P4FvDfVPVyYCHRtW8NY1RjvvZGrVOpSWcpcJ3/+pvAc8CX3AYicgXQoKrPAKjqWxVe0zBGBPO1N2qdSgX/IlU96L8+BFwU0eYy4LiIbAAuBjYB96pqrsJrG8aQEw62cithASb6Rk1RUvBFZBMwJWLXfe4bVVURicrE1gD8BvBBYC+ezf9u4O8jrrUCWAHQ2mp5SIyRJSr3PZQOwDKM0UpJwVfVxXH7ROSwiExV1YMiMpVo23wP8AtV7fKP+T6wiAjBV9U1wBrwsmWmugPDGCLict/bwq1Rq1S6aLsR+KT/+pPAExFttgKTROQC//3HgJcrvK5hDDlRi7TuNhH4UechHt6yd6S7ahipqCgfvohMBr4LtALdeG6Zb4pIG/BZVb3Hb3cD8FW8im8dwApV7Ys5LWD58I3RQVzu+6/95DWeeflwod0lF7yb3/vQbEuJbIw4Q5YPX1WPAddHbN8G3OO8fwZ4fyXXMoyRICoh2oKZzbx9ptjnYPfRX/Pl770EYKJvjFosl45hlEFHdy+rN+9m7tRzI/c/veNg5HbDGA1YagXDSEnYa+ezH57Npn8+wu4jZ0NLlsybOoI9NIxkTPANIyWu105ff57Og7/iK8vfz85DJ3l6x0GWzJtq5hxjVGOCbxgpCTx0AtF//tU32LrnTb5zzyITeqMmMBu+YaQkSK1w7SXnI3hZAF3/fMMY7ZjgG0YZLJjZzBcWX8a4RkuiZtQeZtIxjDKxJGpGrWKCbxhl4AZife6jl4x0dwyjLEzwDSMlrltmRoRVS+fZYq1RU5gN3zBS4rpl9ueVlU/soKN7QIE3wxi1mOAbRkoWzZ5MRs6WMM/n1Tx0jJrCBN8wUrJgZjOrls6jISNkgKZG89Axaguz4RtGGdx5dStzpkxg/fYepHRzwxhV2AzfMAbBhu09rHtxL3etbTc7vlEzmOAbRpnEVcIyjNGOCb5hlEFHdy/7j5+mIWuRtkbtUZENX0TOwytKPgvYg1fxasDzrYj8BfCv8QaYZ4B/r5WU2jKMEcD1w2/ICLcvbGXZ/BaLtDVqhkpn+PcCz6rqpcCz/vsiRORfAtfiVbyaB1wFfKTC6xrGsOOacnJ5Zdqkc0zsjZqiUsFfCnzTf/1N4Lci2ijwLqAJGAc0Aocj2hnGqCaqqLlh1BKVumVepKpBTbdDwEXhBqr6gohsBg7iFTH/G1V9JepkIrICWAHQ2moh68bowpKmGbVOScEXkU3AlIhd97lvVFVFZIBdXkQuAS4HWvxNz4jIb6jqz8JtVXUNsAagra3NbPzGqCOqqLlh1AolBV9VF8ftE5HDIjJVVQ+KyFTgSESz3wbaVfUt/5ingWuAAYJvGIZhDB2V2vA3Ap/0X38SeCKizV7gIyLSICKNeAu2kSYdwzAMY+ioVPAfAG4QkVeBxf57RKRNRNb6bR4HXgNeAn4J/FJVf1DhdQ3DMIwyqWjRVlWPAddHbN8G3OO/zgGfqeQ6hmEYRuVYpK1hGEadYIJvGIZRJ8hozXAgIkeB7pHuRxmcD7wx0p0YAey+6wu779HPTFW9IGrHqBX8WkNEtqlq20j3Y7ix+64v7L5rGzPpGIZh1Akm+IZhGHWCCX71WDPSHRgh7L7rC7vvGsZs+IZhGHWCzfANwzDqBBN8wzCMOsEEf5CIyK0i0ikieRGJddcSkRtFZKeI7BaRARXBag0ROU9EnhGRV/3fkbmCRSQnIr/wfzYOdz+rRam/n4iME5FH/f1bRGTWCHSz6qS477tF5KjzN75nJPpZTUTk6yJyRER2xOwXEflr/zP5JxGZP9x9rBQT/MGzA1gG/DSugYhkgdXAEuAK4A4RuWJ4ujdklCxr6XNaVT/g/9w8fN2rHin/fp8CelX1EuAvga8Mby+rTxn/t486f+O1EftrjW8ANybsXwJc6v+sAP52GPpUVUzwB4mqvqKqO0s0WwjsVtUuVe0DHsErC1nLpClrOVZI8/dzP4/HgetFRIaxj0PBWPy/LYmq/hR4M6HJUuBb6tEOTPLrgNQMJvhDy3Rgn/O+x99Wy5Qsa+nzLhHZJiLtIvJbw9O1qpPm71doo6r9wAmg1ovdpv2/Xe6bNh4XkRnD07URpea/z5XWtB3TJJV3VNWoYi9jgkrLWvrMVNX9IjIb+LGIvKSqr1W7r8aI8QNgnaq+IyKfwXvK+dgI98kogQl+AknlHVOyH3BnPi3+tlFNFcpaoqr7/d9dIvIc8EG8Qji1RJq/X9CmR0QagInAseHp3pBR8r79WhgBa4G/GIZ+jTQ1+X12MZPO0LIVuFRELhaRJuB2vLKQtUzJspYi0iwi4/zX5wPXAi8PWw+rR5q/n/t53AL8WGs/mrHkfYds1zdTH2VLNwK/63vrLAJOOObN2kBV7WcQP3jF2XuAd4DDwA/97dOAp5x2vwnswpvd3jfS/a7CfU/G8855FdgEnOdvbwPW+q//JWdLWr4EfGqk+13B/Q74+wGrgJv91+8CHgN2Ay8Cs0e6z8N03/8V6PT/xpuB9410n6twz+uAg8AZ/7v9KeCzwGf9/YLnvRSUbG0b6T6X+2OpFQzDMOoEM+kYhmHUCSb4hmEYdYIJvmEYRp1ggm8YhlEnmOAbhmHUCSb4hmEYdYIJvmEYRp3wfwEo/Cwav23OCwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "p = 0.5 # Which norm do we use\n",
    "M = 4000 # Number of sampling points\n",
    "b = []\n",
    "for i in range(M):\n",
    "    a = jax.random.normal(jax.random.PRNGKey(i), (1, 2))\n",
    "    if jnp.linalg.norm(a[i, :], p) <= 1:\n",
    "        b.append(a[i, :])\n",
    "b = jnp.array(b)\n",
    "plt.plot(b[:, 0], b[:, 1], '.')\n",
    "plt.axis('equal')\n",
    "plt.title('Unit disk in the p-th norm, $p={0:}$'.format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Why $L_1$-norm can be important?\n",
    "\n",
    "$L_1$ norm, as it was discovered quite recently, plays an important role in **compressed sensing**. \n",
    "\n",
    "The simplest formulation of the considered problem is as follows:\n",
    "\n",
    "- You have some observations $f$ \n",
    "- You have a linear model $Ax = f$, where $A$ is an $n \\times m$ matrix, $A$ is **known**\n",
    "- The number of equations, $n$, is less than the number of unknowns, $m$\n",
    "\n",
    "The question: can we find the solution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The solution is obviously non-unique, so a natural approach is to find the solution that is minimal in the certain sense:\n",
    "\n",
    "\\begin{align*}\n",
    "& \\Vert x \\Vert \\rightarrow \\min_x \\\\\n",
    "\\mbox{subject to } & Ax = f\n",
    "\\end{align*}\n",
    "\n",
    "- Typical choice of $\\Vert x \\Vert = \\Vert x \\Vert_2$ leads to the **linear least squares problem** (and has been used for ages).  \n",
    "\n",
    "- The choice $\\Vert x \\Vert = \\Vert x \\Vert_1$ leads to the [**compressed sensing**](https://en.wikipedia.org/wiki/Compressed_sensing)\n",
    "- It typically yields the **sparsest solution**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What is a stable algorithm?\n",
    "\n",
    "And we finalize the lecture by the concept of **stability**.\n",
    "\n",
    "- Let $x$ be an object (for example, a vector) \n",
    "- Let $f(x)$ be the function (functional) you want to evaluate \n",
    "\n",
    "You also have a **numerical algorithm** ``alg(x)`` that actually computes **approximation** to $f(x)$.  \n",
    "\n",
    "The algorithm is called **forward stable**, if $$\\Vert \\text{alg}(x) - f(x) \\Vert  \\leq \\varepsilon $$  \n",
    "\n",
    "The algorithm is called **backward stable**, if for any $x$ there is a close vector $x + \\delta x$ such that\n",
    "\n",
    "$$\\text{alg}(x) = f(x + \\delta x)$$\n",
    "\n",
    "and $\\Vert \\delta x \\Vert$ is small."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Classical example\n",
    "A classical example is the **solution of linear systems of equations** using Gaussian elimination which is similar to LU factorization (more details later)\n",
    "\n",
    "We consider the **Hilbert matrix** with the elements\n",
    "\n",
    "$$A = \\{a_{ij}\\}, \\quad a_{ij} = \\frac{1}{i + j + 1}, \\quad i,j = 0, \\ldots, n-1.$$\n",
    "\n",
    "And consider a linear system\n",
    "\n",
    "$$Ax = f.$$\n",
    "\n",
    "We will look into matrices in more details in the next lecture, and for linear systems in the upcoming weeks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.2877871109130684e-08\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2f927eeb0>]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0E0lEQVR4nO3deXyb1Z3v8c9PkuXdjuUtibPvdpIGQsgClDVpgbYEaIcC7UCnnUs35nbmdm4HpjNdaDvTdu5MmWm59DKUoTBd6UYotEBYSgIhkEBCVseOyUrs2JYd27Jj2dK5f0iPcRx50/NoifV7v15+WZYe6zlR5K+Pz3PO74gxBqWUUhOfK9UNUEoplRwa+EoplSE08JVSKkNo4CulVIbQwFdKqQyhga+UUhki7QNfRB4SkZMisnsMx84UkedE5C0ReVFEpiWjjUopdS5I+8AHHgauHuOx/wd4xBjzHuAe4J8T1SillDrXpH3gG2NeAvyD7xORuSLyRxHZLiKbRGRR9KEa4Pno7ReA9UlsqlJKpbW0D/xhPAD8lTHmAuBvgf8bvX8ncGP09g1AoYiUpqB9SimVdjypbsB4iUgBcBHwmIhYd2dHP/8t8AMR+QTwEnAcCCW7jUoplY7OucAn8ldJuzHmvKEPGGPeIdrDj/5i+LAxpj2prVNKqTR1zg3pGGM6gLdF5M8AJGJZ9HaZiFj/pruBh1LUTKWUSjtpH/gi8jNgC7BQRI6JyKeAjwGfEpGdwB7evTh7OVArIgeASuBbKWiyUkqlJdHyyEoplRnSvoevlFLKGWl70basrMzMmjUr1c1QSqlzyvbt21uMMeWxHkvbwJ81axbbtm1LdTOUUuqcIiKHh3tMh3SUUipDaOArpVSG0MBXSqkMoYGvlFIZQgNfKaUyhAa+UkplCA18pZTKEBr4Sqm0sv2wnz3vnEp1MyYkDXylVFr5h9/t4dt/2J/qZkxIabvSVimVmVq6elPdhAlLA18plTaMMbQFgrhk9GPV+OmQjlIqbXT09NMfNvgDQbR0u/M08JVSacPfHQSgL2To7O1PcWsmHg18pVTa8AfeHb/3dwVT2JKJSQNfKZU2WgeFfGtAA99pGvhKqbThHxTyfg18x2ngK6XShjWGD9Cmge84DXylVNrwdwXxROdk6pCO8zTwlVJpwx8IUlmUQ06W64wLuMoZuvBKKZU2WgNBfPnegdvKWRr4Sqm00db9buDrRVvnaeArpdJGa1eQeeUFgAZ+ImjgK6XShj8QpCTaw68/2ZXi1kw8GvhKqbTQEwzR0xfCl+9F0B5+ImjgK6XSgjUHvzTfiwj09IXoCYbI9bpT3LKJQ6dlKqXSglU7x5fvpXRgpo5OzXSSBr5SKi1Y4e7L9+LLzwZ0WMdpOqSjlEoLVrj78r2I6GrbRNDAV0qlBSvwS/OzcUUDX0skO8uRIR0ReUhETorI7mEev1xETonIjujHV5w4r1Jq4vAHgrhdQlGuB1+BLr5KBKd6+A8DPwAeGeGYTcaYDzp0PqXUBOMPBCnJiwznFGZ7yHKLDuk4zJEevjHmJcDvxHMppTJTayA4MDtHRPDle7WAmsOSOUtnjYjsFJE/iMjiJJ5XKXUOaBtUOA3Al5+tQzoOS1bgvwHMNMYsA74P/C7WQSJyh4hsE5Ftzc3NSWqaUiod+APBgbF7iCzA0iEdZyUl8I0xHcaYrujtp4AsESmLcdwDxpgVxpgV5eXlyWiaUipNtAaC+PIG9/C92sN3WFICX0QmS3RirYisjJ63NRnnVkqlv75QmFM9fUOGdLw6LdNhjszSEZGfAZcDZSJyDPgqkAVgjPkh8BHgsyLSD/QANxtjjBPnVkqd+9q7+wAoHTKk09nbT29/iGyP1tNxgiOBb4y5ZZTHf0Bk2qZSSp1l8CpbizWe3xboY3KxBr4TtJaOUirlBuro5J3Zwx/8mLJPA18plXIDPfyCM6dlDn5M2aeBr5RKubZYQzq6t63jNPCVUilnzbcviTWkozN1HKOBr5RKOX8gSFGOhyz3u5FUnJuF2yXaw3eQBr5SKuX8gSClBdln3OdyCSV5Wbra1kEa+EqplPMPqaNj0QJqztLAV0ql3MiBrz18p2jgK6VSbmgdHUtpfrYO6ThIA18plVLGmEhp5ALt4SeaBr5SKqU6TvfTHzYD0zAH8+V7ae/uoz8UTkHLJh4NfKVUSsWqo2Oxiqm1RYurKXs08JVSKWXNwikZpocfOUaHdZygga+USil/IFoaeYTA1wJqztDAV0qllNXDjzmkowXUHKWBr5RKKWvapRXug+mQjrM08JVSKeXvCpKT5SLXe/YmJyV5WYAWUHOKBr5SKqX83cGYvXsAj9vFpLws7eE7RANfKZVSw5VVsOjiK+do4CulUmq0wC/N9+osHYdo4CulUqq1S3v4yaKBr5RKqbbu0QI/WwPfIRr4SqmUOd0XojsYGnVIp627j3DYJLFlE5MGvlIqZd6dgz984JfkewmFDad6tJ6OXRr4SqmU8Ufn18eqo2MZ2Mxch3Vs08BXSqWMv3v0Hr6utnWOI4EvIg+JyEkR2T3M4yIi/yEi9SLylogsd+K8Sqlz20h1dCzvBr5OzbTLqR7+w8DVIzx+DTA/+nEHcL9D51VKncOskgnDrbSFd2vi65COfY4EvjHmJcA/wiHrgUdMxKvAJBGZ4sS5lVLnLn8giNslFOZ4hj1moIev9XRsS9YYfhVwdNDXx6L3nUFE7hCRbSKyrbm5OUlNU0qlSlt3kJI8Ly6XDHtMtsdNQbZHe/gOSKuLtsaYB4wxK4wxK8rLy1PdHKVUgrV2BUe8YGvR1bbOSFbgHwemD/p6WvQ+pVQGG62OjkUD3xnJCvwNwG3R2TqrgVPGmBNJOrdSKk2NNfBLNfAdMfyVknEQkZ8BlwNlInIM+CqQBWCM+SHwFHAtUA90A3/hxHmVUuc2/yh1dCy+fC973ulIQosmNkcC3xhzyyiPG+DzTpxLKTUx9IfCtHf3jS3wCyI9fGMMIsNf4FUjS6uLtkqpzNHWHamNY82zH0lpvpdgKExXb3+imzWhaeArpVLCGpMvyRvLkE72Gd+j4qOBr5RKCf8YKmVatICaMzTwlVIpYQW+bwxDOrra1hka+EqplBgonDamIR2tmOkEDXylVEpYwzMj1cK3aAE1Z2jgK6VSoi0QpCjHQ5Z79BjK83rIyXJpiWSbNPCVUinRGghSWjB8WeShSvOztYdvkwa+Uiol/IEgJXlZYz5e6+nYp4GvlEqJSB2dsffwNfDt08BXSqWEPzC20siW0nzvwA5ZKj4a+EqppDPG0NYdHNMcfIv28O3TwFdKJV3H6X76QmZMc/AtvgIvPX0heoKhBLZsYtPAV0ol3cAq23EO6QC06tTMuGngK6WSbjxlFSxaQM0+DXylVNKNp3CaxacF1GzTwFdKJZ21YnYspZEtpVpAzTYNfKVU0vkDY9/8xGIN/+iQTvw08JVSSecP9JKT5SLPO/ZdVguzPWS5RYd0bNDAV0olXWsgSOk4VtkCiEh0Lr7O0omXBr5SKun8gSAl+WOvo2Px5WfrkI4NGvhKOezTj27jwU0NqW5GWmsbZx0dS2m+V4d0bNDAV8pB3cF+ntnbxCNbDmOMSXVz0lbrOOvoWLS8gj0a+Eo5qK6pC2PgiL+bupNdqW5O2opUyowz8HVaZtw08JVyUG1j58DtZ/c2pbAl6et0X4juYCiuwC/N99LZ209vv9bTiYcjgS8iV4tIrYjUi8hdMR7/hIg0i8iO6MdfOnFepdLN/sZOcrJcLKkq0sAfRjx1dCzWXPy26Dx+NT62A19E3MB9wDVADXCLiNTEOPQXxpjzoh8P2j2vUumotqmDBZWFvL9mMjuOtnOy43Sqm5R27AS+FlCzx4ke/kqg3hjTYIwJAj8H1jvwvEqdc2obO1lYWci6xZUAbNx3MsUtSj+tcdTRsWgBNXucCPwq4Oigr49F7xvqwyLyloj8SkSmO3BepdJKS1cvLV1BFk4uZGFlIdN9uWzcp8M6Qw3U0Ynzom3kOTTw45Gsi7ZPALOMMe8BngV+HOsgEblDRLaJyLbm5uYkNU0pZ1gXbBdNLkJEWFtdyeb6FgK9/SluWXoZqKNjZ0hHZ+rExYnAPw4M7rFPi943wBjTaoyxBt0eBC6I9UTGmAeMMSuMMSvKy8sdaJpSybM/GvgLJxcCsK6mkmB/mE11LalsVtrxB3pxu4SinPGvtC3OzcIl2sOPlxOB/zowX0Rmi4gXuBnYMPgAEZky6MvrgH0OnFeptFLb2EFpvpfywsg484WzfBTleHS2zhD+QJCSPC8ul4z7e10uoSRPV9vGa+yl6oZhjOkXkTuBpwE38JAxZo+I3ANsM8ZsAP6niFwH9AN+4BN2z6tUuqlt7Bzo3QNkuV1cuaiC5/c30R8K43HrsheIDMf44qijY9ECavFz5B1ojHnKGLPAGDPXGPOt6H1fiYY9xpi7jTGLjTHLjDFXGGP2O3FepdJFOGw40NR1RuADrKuZTFt3H28caU9Nw9JQW3d8q2wtWl4hftrlUMoBR/zd9PSFWDQk8C9dUEaWW3h2b2OKWpZ+4imNPFhpgQ7pxEsDXykH7G/sAGDh5KIz7i/MyWLN3DKe3dukxdSi4q2jY9Eefvw08JVywP7GTkRgQWXBWY+tq6nkUGs3B5u1mFp/KEx7d19cc/Atvvxs2rv76A+FHWxZZtDAV8oBtY2dzPTlxdyyb211BQDP6Gwd2nvin4Nvsb63rVvr6YyXBr5SDhg6Q2ewKcW5LK0qZqMGvq06OhZdbRs/DXylbDrdF+JQa+Cs8fvB1lZX8ubRdpo7M3s6obVC1okevhZQGz8NfKVsqmvqImw4a4bOYOtqKjEGnt+f2b18q1duawy/QHv48dLAV8qmd2foDB/41VMKqZqUm/BVt8YYPv/TN7jvhfqEnide/m77PXwd0omfBr5SNtU2dpLtcTGrNH/YY0SEdTWVbKproTuYuGJqL9W18ORbJ/jhiwfpSsOibdb2hHZ6+CV5WkAtXhr4StlU29TJ/MoC3KPUhllXU0lvf5jNCSqmZozh3o0HKMz20Nnbz2/fPD76NyWZP9BLUY6HLBtlJrLcLopzs9K2h99xuo9TPek5g0gDXymb9jd2srBy+Au2lpWzfRQmsJjaproW3jzSzpeuWcSSqiIe3XIo7RZ7tdpcdGUpTdPFV8H+MDf+31f4wH9sovN0+oW+Br5SNvgDQZo7e0e8YGvJcru4YmEFz+8/SSjsbBBbvfupxTnctGIat62exYGmLra+7Xf0PHbZraNj8eV703KWzoObG6g/2cXx9h7ueWJvqptzFg18NaH0hcLseedU0jYdGcsF28HW1VTSGgjy5pE2R9uxqa6FN46089kr5pHtcfOhZVMpzs3i0S2HHT2PXZFKmfHX0bGkY3mF4+09fP+5et6/uJLPXT6Xx7Yf45k96VVDyXZ5ZKVSqS0Q5I0jbWw/3Ma2w228dayd031hbl8zk6+vX5Lw87+7y9XYAv+yheXRYmpNrJjlc6QNxhj+/bk6pkR79wC5Xjc3rZjGf718iKaO01QW5ThyLrv8gSDLpk2y/TylBV7ecPiXpl3feGIvBsNXPrSY8oJsXqxt5u7f7OL8GSUDeySkmvbw1TnDGEP9yS5+8foRvvSrnVz1ry9y/jee5VM/3sYDLzXQ2xfi1pUzWVJVxOb65OwyVdvYSUle1ph/oItyslg9p5RnHdzrdnN9C9sPt/G5aO/e8vHVMwkZw0+3HnHsXHYYY2jrDtqaoWPx5Xtp6+4j7PDQWLxerD3JH/c08ldXzqdqUi5ej4vvffQ8Onv7ufs3b6XNtRTt4au01RMMsfNYO9sPR3rwbxxpoz1aP2VSXhYXzCjhxuXTWDGzhPdMm0SuNxJ2P/zTQb79h/00d/YmvGe1L1pSQWTsuzetra7kqxv2cLC5i7nlZxdbG4/I2P2ZvXvLzNJ8LltQzs9eO8KdV86zNTPGCZ29/fSFjK05+BZffjahsOFUj71CbE7o7Q/xtQ17mFOWz1++d/bA/QsqC/nS+xfyzSf38cttR/nohTNS2MoIDXyVdh546SBPvnWCPe900B/twc0tz+d9NZWsmOlj+cwS5pbnDxuyq2ZHhkpee9vPB94zJeYxTgiHDXVNndy0YvroBw+ytiYS+Bv3NjH3MnuB/3J9K9sPt/GN9YvP6N1bblszk08+vI2n9zTywfdMtXUuu6w5+E7N0oHIrJ9UB/4Df2rgUGs3j35q5Vn/B5+8eDbP7TvJPU/sZc2cMmaU5qWolRE6pKPSyh93n+CfntoPItxx6Rx+dPsK3vzHdTz3xcv57keWcdOF05lXUTBij3pJVTF5Xjdb325NaFuPtnXTHTx705PRVE3KZfHUItvTM62ZOZOLcrjpwti/dC5bUMEMXx6PpMHFW2vTEqs0gh2+gYqZqb1we9TfzQ9eqOcDS6fw3vnlZz3ucgn/56ZluET44mM7HJ+dNV4a+CptnOrp4yuP76FmShG/+swavnT1Iq6qrhx3Dy7L7eKCmSVsbUjslMT90Qu2Y52hM9ja6kq2H2mjpSv+qYUv17ey7XAbn7tibszePYDbJXx89Qxee9s/MKMoVQYqZeY5F/ipXm379Sf24nYJ//DB6mGPqZqUy9euW8zrh9r4z00NSWzd2TTwVdr49h/209LVy3c+/B7b482rZvuobeqkLYFT96wZOgsqxx/47xZTOxnXuSMzc6K9+1GGlG5aMZ1sjyvlvfw2B0ojW0rToIDac/ua2LiviS9cNZ8pxbkjHnvj8iquWTKZf32mlr3vpO4Xrwa+SgtbG1r52WtH+NQls1k6rdj2862aUwrAa4cS18uvbexkhi+P/OzxXwpbPLWIqcU5cQ/rvHKwldcPRXr3OVmxe/eWSXlerls2ld+9eZyOFK7+tIZ0Sh0c0vGnaPHV6b4QX3tiD/MqCviLi2ePeryI8K0bllKc6+V//XIHvf2hJLTybBr4KuVO94W4+ze7mO7L5W/WLXDkOd8zrZhsjyuhwzr7GzviGs6BSACsralkU10zp/vG98Nvjd1XFmWP+YLxbWtm0R0M8evtx+JpriP8gV5yslwxdwUbr2yPm4JsT8o2M7//xYMc9fdwz/rFeD1ji1FfvpfvfmQp+xs7+bdnDiS4hbFp4KuzvNPeQ18S9wu974V6GloC/NMNSx0JA4gEwvIZJQm7cBvZ9KR73BdsB1tXU8npvvEXU9ti9e4vnzdq796ydFox502fxKOvHk7ZnPDWQNCR8XtLqlbbHm4NcP+fDnLdsqlcNLdsXN975aJKblk5gwc2NbC1IbGTCmLRwFcAhMKGp/c0cut/vspF336eTz+6PSmbRO9v7OD+Fw9y4/KqmLMc7Fg528feEx0JqVxYf7KLUNjE3cMHWDW7lMLs8RVTs+bdVxZl89FhZuYM57Y1M2loDvByffKDBiJj+E7M0LGkIvCNMXxtwx6yXMKXPzD8hdqR/MMHqpnhy+OLj+1MeoE1DfwM5w8Euf/Fg1z63Rf49KPbOdQS4Mbzq3h+/0m++eS+hJ47FDb83a93UZybxT9+oMbx5181x4cxsP2w88M64y2pEIvX4+KyheU8t79pzCtGtxxs5bVDfj572ehj90Ndu3QKvnwvj2w5FEdr7fMHnKmjYynN9yZ9ls6ze5t4obaZv1m3IO5yFfnZHv7tpvN4p72Hb/w+uQXWdOFVhtp9/BQ/fuUQj+98h2B/mIvmlvKPH6xhbXUFHrcLX76XBze/zczSvDFdlIrHI1sOsfNoO/9+83kJWTyzfEYJXndkHP/KRZWOPndtUyfeUTY9GYt1NZX8/q0TvHm0nQtmlox4rDGGe5+ro6Iwm5tXjn/VZk6Wm49eOJ3/96eDHG/voWrSyDNLnNYaCDLH5sriwXz5XvYkccZLTzDE15/Yy8LKQm6/aJat57pgZgmfvXwu971wkLXVlbxv8WRnGjkKR3r4InK1iNSKSL2I3BXj8WwR+UX08a0iMsuJ86rxCfaHeXzHcT58/yt88PubeXLXCW5aMY1n/uZSfvo/VnP1ksl4otMh7762mvfVVHLP7/eyMQH124+1dfMvT9dy+cJyrluWmBWgOVlulk0v5tUElAje39jJvPKCgdcrXpcvrMDjEjaOobbOloZWXnvbz+cuH3/v3vKxVZFfFD/dmvwpmm2B4MBuVU7wFUSGdJJ1TeK+F+ojZY/XL3akTMUXrlrA4qlF3P2bXbbWY4yH7VaLiBu4D7gGqAFuEZGhf59/CmgzxswDvgd8x+551dg1dZzme88e4OLvPM8Xfr4DfyDIVz5Yw5a7r+Kb1y+NOY/c7RLuvfk8llYV81c/e5Pdx0851h5jDP/wu90AfPP6JeOqQzNeq2aXsvv4Kce3+6tt7LA1nGMpzs1i1RzfmMbx790Yf+/eMq0kjysXVfLz144mdWrg6b4QgWDIkSmZltJ8L8FQOClbOTY0d/HASw3ccH7VwJRfuwYXWLvr17uS8ovLiSGdlUC9MaYBQER+DqwHBg9OrQe+Fr39K+AHIiImAf/CjtN9/O/Hdjr9tAOEs8MpVl6FjSEUjoRbyBhCYYMxkXHrkDGR+8OGkInUZAlHvw4bQ9hAYY6HsoLs6Id34HZp9HZ5QTZFuZ5hw9IYw/bDbTz8yiH+uLuRkDFcsbCC29bM5NL55bhG2Y4PIM/r4cHbV3DDfa/wyYdf53efv5ipDgwDbNj5Di/WNvOVD9YwrSSxtUVWzvbxgxfq2X64jcsWOHNRuL07SFNHr60LtoOtra7k60/s5e2WALPLYg8RbTkY6d1/9UM1cffuLbetmcnGfU08tesEN5w/bfRvcIDfwUVXFut6gD8QpDAny7HnHcoYw1c37CHb4+Luaxc5+tyDC6w9tu3YsCUynOJE4FcBRwd9fQxYNdwxxph+ETkFlAJnzEcTkTuAOwBmzIivFxMOGw63dsf1vaOJ9evJcPadxoBLBJdLcLuit0VwuwSXRL72uF1ke6LHSKRHLSK4RRCJ/OI60trNm0fa8AeCxLqml+UWSvOzKSv0Rj4XRG4XeD38YXcje090UJjj4RMXzeLjq2cya5gwGUlFYQ4PfeJCPnx/JPR/9dmLKIhjoZGlLRDknif2smz6JNvjoGNxwcwS3C7htbdbHQt8OyUVYllXEwn8jXub+B+Xzol5zL0bD1BRmM0tNnr3lkvmlTGnLJ9Hthw+pwN/cAG1mTavpYzkj7sb2VTXwlc/VENFofP7ClgF1r7+xB7WzC1lui9xnaC0umhrjHkAeABgxYoVcfX+J+V5+eNfX+pou1ItFI7UEW/p6qWlM0hroJfmzl5aA0FarM9dvdQ1ddLSFSQYCrOwspB/umEp158/1fbc9oWTC7nvY8v55MOvc+dP3+DB21bEPXb9zSf3caqnj598eOmom347IT/bw9KqYkcXYO0/EblQWD1l9H1sx2JaSR7VUyLF1GIF/paDrWx1qHcPkYJeH189k3t+v5fdx0+xpMr+yubRJKaHH11tm8CZOt3Bfu75/V6qpxTx56tnJuQcVoG1q7/3El/85U5+dsfqhP1sOBH4x4HBf4dMi94X65hjIuIBioHUTAY+B7ldMjCkwygX840xdAdD5Hndjo6NX7agnG+sX8Lf/3YXX39iL/esXzzu599U18yv3zjGnVfMY9FkZ8JyLFbN8fHQ5rfpCYYGaubbUdvUyaS8LCocrLW/rrqCH7xQH526eGYo/vtzByh3qHdv+fAF0/iXp2t5ZMshvvuRZY4973ASGvgJnIv/H8/Vc+LUab5/y/m2L9CPxCqw9sXHdvLgpgY+fdnchJzHiX/B68B8EZktIl7gZmDDkGM2ALdHb38EeD4R4/cqsmQ/P3v4sX07bl01g09fOodHXz3MQy8fGtf3dgf7+fvf7mJOWT53XjnP8baNZPXsUvpCxrF9ZPc3drKwcnybnoxmXc1kwjGKqW052MqrDfHNux9JcW4W159fxeM73qE9CSWGB+roODmkU/DukE4i1J/s5MFNDXzkgmmObUc5khuXV3H14sn86zMH2HciMdNNbQe+MaYfuBN4GtgH/NIYs0dE7hGR66KH/QgoFZF64H8BZ03dVOeGv7t6Edcsmcw3n9w7rg2a791Yx1F/D/9841JHg2ssLphVgktgqwPTM8Nhw4HGTkdm6Ay2pKqIyUU5PLv3zNfU6t3fusr53ZJuWzOT3v4wj21LfH0df6AXt0socvDiap7XQ06WKyEF1KwLtXleN3dd4+yF2uGICP9041KKcrP44i93JmT7Rkf+RjHGPGWMWWCMmWuM+Vb0vq8YYzZEb582xvyZMWaeMWalNaNHnXtcLuHfbjqP90ybxBd+voNdx0afrrnr2Cke3NTALStnODalbTyKcrKomVrkSF2d4+09BIIhFjo8JBUpplbBSwdaBoqpvdoQ6d1/xuHevaV6ShEXzirhv7ceTvjesP5AHyV5WWOaHTYepfnZCenh//6tE7xc38r/fv/CyFBqkvjyvdz70fP46odqHH+tQEsrqDjket38520X4Mv38skfv87x9p5hj+0Lhfm7X79FWUF20npKsayaXcqbR9ptzz13eobOYOtqJtPTF+KVg5HJa/++sY7ywuyBxVKJ8OdrZnG4tZs/1TUn7BwQ6eE7OX5vSUQ9na7efr755F6WVBVx66rEXKgdySXzyxLWMdLAV3GpKMzhv/7iQk4HQ3zq4deHLQL1o81vs/dEB/esX0JxbuLmSo9m1Wwfvf1hdh61t4CsNrprVCICf/UcHwXRYmpbG1rZ0tCasN695erFkykryObRBG+OEutitBMSEfjff76Opo5evrF+SVJmkiWTBr6K24LKQu7/+AXUn+zi8z9986zqmodaAnzv2QO8f3ElVy9JTq2Q4ayc7UME2yVp9zd2Mq0k19ZahOFke9xctqCcjftOcu/GOsoKEtu7h8hqz1tXTueF2pMcSdD6FYhcWC11sHCaxekCasYYfvPGca5ePJnzZ4xc2+hcpIGvbLlkfhnfvH4JLx1o5qsb9gwsDzfG8Pe/3YXX7eKe9UtS3MrI+oyFlYW2d8CqTcAF28HW1lTQ3Nkb7d3PScoF7ltXzcQlwn8nsL5OWyBISb7zf+E53cM/0NRFc2cvVy6qcOw504kGvrLt5pUz+Mxlc/nJ1iP8aPPbADy2/RivHGzl7mur4y4j67RVs31sP9wW9+Yuvf0hGloCCRnOsVyxsGJg3cXHkjR+PLk4h/fVVPLLbUfHvfvWWITChvaePkdLI1t8BV56+kL0BJ1p96botYyL549vY5NzhQa+csSX3r+Qa5dO5ltP7eMnWw/zrSf3sXKWj5sTXBtkPFbNKaU7GGJXnIXgDp4MRDc9SdyisUl5Xv7u6oX8841LHVkkNlZ/vmYm7d19bNj5juPP3dYdxBhn5+BbrB20Wh2amrm5voU5ZflJLx2dLBr4yhHWdM1l0ybx5d/upicY4p8/vDQhU8vitXJ2ZPFMvGUWapsiF2wTOaQDcMelc1lX42z9/tGsmVPK/IoCHt3i/BaIiVhla3FytW2wP8zWBj+XTNDePWjgKwflZLl58PYVLJs+iS9/oJq5Dm524YSygmzmVRTEPR9/f2MnWW4ZtqLluUxEuG3NTHYdP8WOo+2OPnciA9/J1bZvHGmjpy/EJfM08JUak7KCbB7//MVJqYQZj1WzfWw71EYojoVGtY2dzC0vcGTzi3R0w/JpFGR7HJ+imdgefrREsgMzdTbXteB2CavnJn9xYLJMzHeuUsNYOdtHV28/e+PYGq+2sdOxCpnpqCDbw43Lq/j9WydodXAHpkTU0bE4OaSzqb6FZdOKHS3/kG408FVGWR1dwTjeYZ1T3X2cOHU6oTN00sGfr55JMBTmF9uOjn7wGFm970TsW1yU4yHLLbaHdE5197HrWDuXzHdmz4R0pYGvMkplUQ6zSvN4dZwXbvcncIVtOplfWciFs0p4YucJx56zrTtIYY4nIUNhIkJJntd2AbUtDS2EDbx3Al+wBQ18lYFWzS7l9UP+cRUMq22K1NBJ9AyddLCuppJ9JzpGrJE0HpFVts737i1OLL7aVNdCvtfNedMnOdOoNKWBrzLOqjk+TvX0DYT4WOxv7KQox8PkNFlElkhrqyNTQp/bN/rG6mORqMJpltICr+0hnc31LayeUzphL8hbJva/TqkY3p2PP/Zx/EhJhaKEbCyTbuaUFzCnPJ9n9zoT+K1dwYSssrX48rNt9fCP+rs53No9oeffWzTwVcaZVpJH1aTcMW+IYkxk05OJPn4/2NrqSl5taB22Cup4tHUH8SWgjo6lNN9ra1rm5vpIOeqJPn4PGvgqQ62a4+O1t/1jWlV6vL2Hzt7+jAv8vpBhU12LrecxxkRLIyeyh++ls7c/7r0ONte1UFmUnXYLBRNBA19lpNWzS2kNBKk/2TXqsbWNmXPB1rJ8xiRK8rLYaHNYp7O3n76QSfhFW4C2wPj/GgmHDS8fbOGSeeUZMVynga8y0qo50XH8MQzrWLtcLcigwPe4XVyxqILna0+etc/BeFhDLQm9aJsffwG1Pe900N7dlxHDOaCBrzLUDF8elUXZYwr82sZOqiblTugVmLGsra6kvbuPN460x/0c/u7EB76d1bab6qPlkCdw/ZzBNPBVRhIRVs0uZWtD66jj+LUZdsHWcumCcrxuFxttTM9MSg+/IP7A31zXwqLJhZQXJm+j8lTSwFcZa9UcHyc7ezk0wtZ+wf4wB5u7MjLwC7I9rJ5bamscP5GF0yzWBeHxbnXYEwyx7VDbhK6OOZQGvspYq2ZH6uq8NkJdnYaWLvrDJqMu2A62trqChpYAB5tHv7gdy0DhtILEBf6k3CxcMv4e/uuH/ARD4YyYf2/RwFcZa255PmUF3hE3RHl3hs7ErZI5kqtsrrpt6w6S7XGRm8C9eV2uSD2d8a623VzfgtftGliIlwk08FXGEhFWzvaNeOF234nIpidzyifepidjUTUpl5opRWzcezKu72/titTRSfSUx0g9nfHN0tlU18LymZPI83oS1Kr0YyvwRcQnIs+KSF30c8kwx4VEZEf0Y4OdcyrlpFWzSzne3sNRf+xx/NrGjgm96clYrK2uYNthP21xXBT1B3rxJXA4xzLeAmotXb3sO9HBeyd4OeSh7L6L7wKeM8bMB56Lfh1LjzHmvOjHdTbPqZRjRpuPn6kzdAZbW1NJ2MALtePv5Sd6la1lvAXUXo6WU8ikC7ZgP/DXAz+O3v4xcL3N51MqqRZUFDIpLyvmhdtTPX28kwGbnoxmydRiKouy45qe6e8O4stL/PqF8fbwN9e1UJybxZKq4gS2Kv3YDfxKY4y1U0IjUDnMcTkisk1EXhWR622eUynHuFzCylmxx/EPZFAN/JG4XMKViyr5U23zuOvV+BNcKdPiy8+mvbtvTKuCjTFsrm/hormluF0Tv5zCYKMGvohsFJHdMT7WDz7ORFavDLeCZaYxZgVwK3CviMwd5lx3RH8xbGtubh7vv0WpuKyc7eNwazeNp06fcb9VUmFhhs7QGWxdTQWBYGjEGU1Dne4LEQiGEjol02KVV2jrHr2eTkNLgBOnTmfUdEzLqIFvjFlrjFkS4+NxoElEpgBEP8cc5DPGHI9+bgBeBM4f5rgHjDErjDEryssz62KKSp3h9rmtbeygMMfD1OKJv+nJaC6aW0ZulntcwzrJWHRlGU95hc3RCqDvnZd5GWN3SGcDcHv09u3A40MPEJESEcmO3i4DLgb22jyvUo6pnlJEYY7nrH1uaxs7WVhZmBFVFEeTk+XmkvllbNzbNKaS0vBu+JbkJa+HP5YCapvqWpjuy2VGaV6im5V27Ab+t4F1IlIHrI1+jYisEJEHo8dUA9tEZCfwAvBtY4wGvkobbpdw4SzfGRdujTHs1xk6Z1hXXck7p06z78TYtob0J2GVrcU3xno6faEwrza0ckkG9u4BbK04MMa0AlfFuH8b8JfR268AS+2cR6lEWzXbx/P7T9Lc2Ut5YTYnTp2m83R/xl+wHeyKRRWIwMZ9TdRMHf26RjoO6bx1rJ2u3v6MKYc8VOauJlFqEGt5/WvR2Tq1esH2LOWF2Zw3fdKYx/EH6ugkIfCtYaPRCqhtqmtBBC6aW5rwNqUjDXylgCVVxeR53QMXbgdm6FRqD3+wtdWVvHXsFE0dp0c9ti0QxO2SpOwjkOV2UZybNWoPf3NdC0uripmUhOsK6UgDXykigXHBzJKBaYe1jR1MLc6hOAmLhs4l62qsYmqjr7ptDQQpycvClaS57qX53oENV2LpPN3Hm0fbM2517WAa+EpFrZ5TSm1TJ22BoF6wHcb8igKm+3LHNKzjD/QmZfze4sv3Dmy4EsvWBj+hsMnI+fcWDXylolZFx/FfOdga3fREx++HEhHWVlfycn0L3cH+EY9tC/QlP/BHGNLZXN9CTlbkL7lMpYGvVNTSacVke1z8/PUj9IUyd9OT0ayrrqS3PzywgGk4rUnu4Y9WQG1TXTMrZ5eS7Ulcbf50p4GvVFS2x83yGSVsigaZDunEduFsH4U5nlGHdSKVMpPbw2/rDhIOn70w7MSpHg42B3hvBo/fgwa+UmewyiV7XMLc8oIUtyY9ZbldXL6wguf3n4wZrgChsKG9py8phdMsvvxsQmFDx+mz6+lYf41k8vg9aOArdQZrn9s55fl4PfrjMZy11RW0dAXZcaw95uPt3UGMSc4cfMu75RXOHtbZXN9CWYE344fp9B2t1CDnz5iE1+3SC7ajuHxBBR6XsHFv7GGdgTo6SR7SGXxuSzhseLm+hYvnlWV8XSQNfKUGycly8/1bz+d/Xjkv1U1Ja8V5WVw4yzfsOH4yV9larMAfutq2tqmTlq5gRs+/t2jgKzXE+xdPZr6usB3V2ppKDjR1caT17P2Ak1lHx1I6TAG1gXLIGbZ/bSwa+EqpuKytrgCI2cv3p7CH7x9SInlTfQvzKgqYrPsaaOArpeIzszSf+RUFIwZ+MmvWZHvcFGR7zrhoe7ovxGtvt+pwTpQGvlIqbmtrKnntbT+nes6cCukPBCnM8SR9plNJ/pkF1N440sbpvrAGfpQGvlIqbmurK+kPG/504Mw9qFsDwaQO51h8+dlnBP7muhY8LmF1hpZDHkoDXykVt/OmT6I033vW9My2JK+ytZTme8+YpbO5voXzZ0yiINvWXk8Thga+Uipubpdw5aIKXqw9SV8oPHB/a4oCf3ABtbZAkF3HT3GxDucM0MBXStmytqaSjtP9vH7o3U3gk10a2VIaDXxjDFsaWjGGjN3OMBYNfKWULe+dX4bX42Lj3simKMaYaOG05NXRsfjyvQRDYbp6+9lU10Jhtodl0yYlvR3pSgNfKWVLntfDxXNLeW5/E8YYunr76QuZFF20fXfx1eb6ZlbPLcXj1piz6CuhlLJtbU0lh1u7qT/ZlZI6OhZrte2Oo+0c9ffodMwhNPCVUrZdtSiy1+2z+5pSUkfHYg0jPb7jHUDLIQ+lga+Usm1ycQ5Lq4rZuLdpYF/ZVF20BXjpQDNTi3OYU5af9DakMw18pZQj1lZX8ubRdupOdgGpCXzrnP3RzcozvRzyUBr4SilHXFVdgTHwq+1HgdQEfp7XTXa0nIPOvz+brcAXkT8TkT0iEhaRFSMcd7WI1IpIvYjcZeecSqn0tHhqEVOKczjYHCDb4yLPm/zNwkVkYFhHA/9sdnv4u4EbgZeGO0BE3MB9wDVADXCLiNTYPK9SKs2ICGurIxdvS/O9KRtOKS/KoWZKEWUFyV8HkO5sFZgwxuwDRvuPXQnUG2Maosf+HFgP7LVzbqVU+llbU8mjrx5OyZRMy7euX4LHrWP3sSSjolAVcHTQ18eAVbEOFJE7gDsAZsyYkfiWKaUctXqOj3yvOyXj95YlVcUpO3e6GzXwRWQjMDnGQ182xjzuZGOMMQ8ADwCsWLHCOPncSqnEy/a4+fr6JQMLoFR6GTXwjTFrbZ7jODB90NfTovcppSagj1wwLdVNUMNIxrTM14H5IjJbRLzAzcCGJJxXKaXUIHanZd4gIseANcCTIvJ09P6pIvIUgDGmH7gTeBrYB/zSGLPHXrOVUkqNl91ZOr8Ffhvj/neAawd9/RTwlJ1zKaWUskdX2iqlVIbQwFdKqQyhga+UUhlCA18ppTKEBr5SSmUIMSY9F7SKSDNw2MZTlAEtDjUnEbR99mj77NH22ZPO7ZtpjCmP9UDaBr5dIrLNGDNsyeZU0/bZo+2zR9tnT7q3bzg6pKOUUhlCA18ppTLERA78B1LdgFFo++zR9tmj7bMn3dsX04Qdw1dKKXWmidzDV0opNYgGvlJKZYhzOvBF5GoRqRWRehG5K8bj2SLyi+jjW0VkVhLbNl1EXhCRvSKyR0S+EOOYy0XklIjsiH58JVntG9SGQyKyK3r+bTEeFxH5j+hr+JaILE9i2xYOem12iEiHiPz1kGOS+hqKyEMiclJEdg+6zyciz4pIXfRzyTDfe3v0mDoRuT2J7fsXEdkf/f/7rYhMGuZ7R3wvJLB9XxOR44P+D68d5ntH/HlPYPt+Mahth0RkxzDfm/DXzzZjzDn5AbiBg8AcwAvsBGqGHPM54IfR2zcDv0hi+6YAy6O3C4EDMdp3OfD7FL+Oh4CyER6/FvgDIMBqYGsK/78biSwqSdlrCFwKLAd2D7rvu8Bd0dt3Ad+J8X0+oCH6uSR6uyRJ7Xsf4Ine/k6s9o3lvZDA9n0N+Nsx/P+P+POeqPYNefxfga+k6vWz+3Eu9/BXAvXGmAZjTBD4ObB+yDHrgR9Hb/8KuEpEkrKdvTHmhDHmjejtTiKbv1Ql49wOWw88YiJeBSaJyJQUtOMq4KAxxs7qa9uMMS8B/iF3D36f/Ri4Psa3vh941hjjN8a0Ac8CVyejfcaYZ0xkIyKAV4lsM5oSw7x+YzGWn3fbRmpfNDtuAn7m9HmT5VwO/Crg6KCvj3F2oA4cE33DnwJKk9K6QaJDSecDW2M8vEZEdorIH0RkcXJbBoABnhGR7SJyR4zHx/I6J8PNDP+DlurXsNIYcyJ6uxGojHFMuryOnyTyF1sso70XEunO6JDTQ8MMiaXD6/deoMkYUzfM46l8/cbkXA78c4KIFAC/Bv7aGNMx5OE3iAxRLAO+D/wuyc0DuMQYsxy4Bvi8iFyagjaMSCJ7IV8HPBbj4XR4DQeYyN/2aTnXWUS+DPQDPxnmkFS9F+4H5gLnASeIDJuko1sYuXef9j9L53LgHwemD/p6WvS+mMeIiAcoBlqT0rrIObOIhP1PjDG/Gfq4MabDGNMVvf0UkCUiZclqX/S8x6OfTxLZrnLlkEPG8jon2jXAG8aYpqEPpMNrCDRZw1zRzydjHJPS11FEPgF8EPhY9JfSWcbwXkgIY0yTMSZkjAkD/znMeVP9+nmAG4FfDHdMql6/8TiXA/91YL6IzI72AG8GNgw5ZgNgzYb4CPD8cG92p0XH+34E7DPG/Nswx0y2rimIyEoi/x/J/IWULyKF1m0iF/d2DzlsA3BbdLbOauDUoOGLZBm2Z5Xq1zBq8PvsduDxGMc8DbxPREqiQxbvi96XcCJyNfAl4DpjTPcwx4zlvZCo9g2+JnTDMOcdy897Iq0F9htjjsV6MJWv37ik+qqxnQ8iM0gOELl6/+XoffcQeWMD5BAZBqgHXgPmJLFtlxD50/4tYEf041rgM8BnosfcCewhMuPgVeCiJL9+c6Ln3hlth/UaDm6jAPdFX+NdwIoktzGfSIAXD7ovZa8hkV88J4A+IuPInyJyXeg5oA7YCPiix64AHhz0vZ+Mvhfrgb9IYvvqiYx/W+9Da+baVOCpkd4LSWrfo9H31ltEQnzK0PZFvz7r5z0Z7Yve/7D1nht0bNJfP7sfWlpBKaUyxLk8pKOUUmocNPCVUipDaOArpVSG0MBXSqkMoYGvlFIZQgNfKaUyhAa+UkpliP8PR0pNnje5VWwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 20\n",
    "a = [[1.0/(i + j + 1) for i in range(n)] for j in range(n)] # Hilbert matrix\n",
    "A = jnp.array(a)\n",
    "#rhs =  jax.random.normal(jax.random.PRNGKey(0), (n,))\n",
    "rhs = jnp.ones(n)\n",
    "sol = jnp.linalg.solve(A, rhs)\n",
    "print(jnp.linalg.norm(A @ sol - rhs)/jnp.linalg.norm(rhs))\n",
    "plt.plot(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "rhs = jnp.ones(n)\n",
    "sol = jnp.linalg.solve(A, rhs)\n",
    "print(jnp.linalg.norm(A @ sol - rhs)/jnp.linalg.norm(rhs))\n",
    "#plt.plot(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More examples of instability\n",
    "\n",
    "How to compute the following functions in numerically stable manner?\n",
    "\n",
    "- $\\log(1 - \\tanh^2(x))$\n",
    "- $\\text{SoftMax}(x)_j = \\dfrac{e^{x_j}}{\\sum\\limits_{i=1}^n e^{x_i}}$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "u = 300\n",
    "eps = 1e-6\n",
    "print(\"Original function:\", jnp.log(1 - jnp.tanh(u)**2))\n",
    "eps_add = jnp.log(1 - jnp.tanh(u)**2 + eps)\n",
    "print(\"Attempt to improve stability by adding a small constant:\", eps_add)\n",
    "print(\"Use more numerically stable form:\", jnp.log(4) - 2 * jnp.log(jnp.exp(-u) + jnp.exp(u)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 5\n",
    "x = jax.random.normal(jax.random.PRNGKey(0), (n, ))\n",
    "x = jax.ops.index_update(x, [0], 1000)\n",
    "print(jnp.exp(x) / jnp.sum(jnp.exp(x)))\n",
    "print(jnp.exp(x - jnp.max(x)) / jnp.sum(jnp.exp(x - jnp.max(x))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take home message\n",
    "- Floating point  (double, single, number of bytes), rounding error\n",
    "- Norms are measures of smallness, used to compute the accuracy\n",
    "- $1$, $p$ and Euclidean norms \n",
    "- $L_1$ is used in compressed sensing as a surrogate for sparsity (later lectures) \n",
    "- Forward/backward error (and stability of algorithms)  (later lectures)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "\n",
    "- Matrix norms: what is the difference between matrix and vector norms\n",
    "- Unitary matrices, including elementary unitary matrices."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Слайд-шоу",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  },
  "nav_menu": {},
  "rise": {
   "scroll": true,
   "theme": "sky",
   "transition": "zoom"
  },
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}