lectures/lecture1/lecture-1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'theme': 'sky',\n",
       " 'transition': 'zoom',\n",
       " 'start_slideshow_at': 'selected',\n",
       " 'scroll': True}"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from traitlets.config.manager import BaseJSONConfigManager\n",
    "cm = BaseJSONConfigManager()\n",
    "cm.update('livereveal', {\n",
    "              'theme': 'sky',\n",
    "              'transition': 'zoom',\n",
    "              'start_slideshow_at': 'selected',\n",
    "            'scroll': True\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Лекция 1. Арифметика чисел с плавающей точкой. Векторные нормы. Матричные нормы. Унитарные матрицы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## План на сегодня\n",
    "\n",
    "- Арифметика чисел с фиксированной и плавающей точкой\n",
    "- Концепция **обратной (backward)** и **прямой (forward)** устойчивости алгоритмов\n",
    "- Как измерять точность: векторные нормы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Представление чисел \n",
    "\n",
    "- Действительные числа обозначают различные величины: вероятности, массы, скорости, длины...\n",
    "\n",
    "- Важно знать, как они представляются в компьютере, который оперирует только битами."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Представление с фиксированной точкой\n",
    "\n",
    "- Наиболее очевидный формат представления чисел – это формат чисел с **фиксированной точкой**, также известный как **Qm.n** формат\n",
    "\n",
    "- Число **Qm.n** лежит в отрезке $[-(2^m), 2^m - 2^{-n}]$ и его точность $2^{-n}$.\n",
    "\n",
    "- Общий размер памяти для хранения такого числа $m + n + 1$ бит.\n",
    "\n",
    "- Диапазон чисел, представимых таким образом, ограничен."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Числа с плавающей точкой (floating point numbers)\n",
    "\n",
    "Числа в памяти компьютера обычно представляются в виде **чисел с плавающей точкой.**\n",
    "\n",
    "Число с плавающей точкой представляется в виде\n",
    "\n",
    "$$ \\textrm{number} = \\textrm{significand} \\times \\textrm{base}^{\\textrm{exponent}},$$\n",
    "\n",
    "где *significand* – целое число (aka мантисса), *base* – натуральное число (основание) и *exponent* – целое число (может быть отрицательным), например\n",
    "\n",
    "$$ 1.2 = 12 \\cdot 10^{-1}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Формат фиксированной точки vs формат плавающей точки\n",
    "\n",
    "**Q**: какие достоинства и недостатки у рассмотренных форматов представления действительных чисел?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**A**: они подходят для большинства случаев.\n",
    "\n",
    "- Однако, числа с фиксированной точкой представляют числа из фиксированного интервала и ограничивают **абсолютную** точность.\n",
    "\n",
    "- Числа с плавающей точкой представлют числа с **относительной** точностью и удобны для случаев, когда числа, которые участвуют в вычислениях, имеют различный порядок (например $10^{-1}$ и $10^{5}$).\n",
    "\n",
    "- На практике, если скорость не критически важна, стоит использовать float32 или float64."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## IEEE 754\n",
    "В современных компьютерах представление чисел в виде чисел с плавающей точкой регулируется стандартом [IEEE 754](https://en.wikipedia.org/wiki/IEEE_floating_point), который был опубликован в **1985 г.** До этого компьютеры обрабатывали числа с плавающей точкой по-разному!\n",
    "\n",
    "IEEE 754 содержит следующие элементы:\n",
    "- Представление чисел с плавающей точкой (как было описано выше), $(-1)^s \\times c \\times b^q$.\n",
    "- Две бесконечности $+\\infty$ and $-\\infty$\n",
    "- Два типа **NaN**: \"тихий\" NaN (**qNaN**) and сигнализирующий NaN (**sNaN**) \n",
    "    - qNaN не бросает исключение на уровне блока, производящего операции с плавающей точкой, (floating point unit – FPU), до того как вы проверите результат вычислений\n",
    "    - значение sNaN бросает исключение из FPU, если вы используете это значение в вычислениях. Этот тип NaN может быть полезен для инициализиции\n",
    "    - C++11 имеет [стандартный интерфейс](https://en.cppreference.com/w/cpp/numeric/math/nan) для создания разных типов NaN\n",
    "- Правила **округления**\n",
    "- Правила для операций типа $\\frac{0}{0}, \\frac{1}{-0}, \\ldots$\n",
    "\n",
    "Возможные значения определяются через\n",
    "- основание $b$\n",
    "- точность $p$ - число цифр в записи\n",
    "- максимально возможное значение $e_{\\max}$\n",
    "\n",
    "и имеет следующие ограничения\n",
    "\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": [
    "## Два наиболее используемых формата: single & double\n",
    "\n",
    "Наиболее часто используются следующие форматы: **binary32** и **binary64** (также известные как **single** и **double**).\n",
    "В последнее время популярность набирают вычисления с половинной точностью **binary16**.\n",
    "\n",
    "| Формальное название | Другое название | Основание | Число цифр в записи | 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",
    "\n",
    "<img src=\"./double64.png\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Примеры\n",
    "\n",
    "- Для числа +0\n",
    "    - *sign* - 0\n",
    "    - *exponent* - 00000000000\n",
    "    - *fraction* - все нули\n",
    "- Для числа -0\n",
    "    - *sign* - 1\n",
    "    - *exponent* - 00000000000\n",
    "    - *fraction* - все нули\n",
    "- Для +infinity\n",
    "    - *sign* - 0\n",
    "    - *exponent* - 11111111111\n",
    "    - *fraction* - все нули\n",
    "\n",
    "**Q**: Как будет выглядеть -infinity и NaN ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Точность и размер памяти \n",
    "\n",
    "**Относительная точность** для различных форматов\n",
    "\n",
    "- половинная точность или float16: $10^{-3} - 10^{-4}$,\n",
    "- одинарная точность или float32: $10^{-7}-10^{-8}$,\n",
    "- двойная точность или float64: $10^{-14}-10^{-16}$.\n",
    "\n",
    "<font color='red'> Crucial note 1: </font> **float16** занимает **2 байта**, **float32** занимает **4 байта**, **float64** занимает **8 байт**\n",
    "\n",
    "<font color='red'> Crucial note 2: </font> Обычно в \"железе\" поддерживается одинарная и двойная точность. Для обучения нейросетей становится популярным использовать половинную точность или даже меньше... Подробности [тут](https://arxiv.org/pdf/1905.12334.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Влияние формата представления чисел на обучение нейросетей\n",
    "\n",
    "- Веса в слоях (полносвязном, свёрточном, действие функций активации) могут храниться с различной точностью\n",
    "- Это важно для повышении энергоэффективности устройств, на которых запускается обученная нейросеть\n",
    "- Проект [DeepFloat](https://github.com/facebookresearch/deepfloat) от Facebook показывает, как переизобрести операции с плавающей точкой, чтобы они были эффективны для обучения нейросетей, подробности см. в [статье](https://arxiv.org/pdf/1811.01721.pdf)\n",
    "- Влияние формата представления чисел на значения градиентов активаций\n",
    "\n",
    "<img width=500, src=\"./grad_norm_fp16.png\">\n",
    "\n",
    "- И на кривые обучения\n",
    "\n",
    "<img width=500, src=\"./train_val_curves.png\">\n",
    "\n",
    "Графики взяты из [этой статьи](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",
    "- Этот формат требует 16 битов\n",
    "    - 1 бит для знака\n",
    "    - 8 битов для экспоненты\n",
    "    - 7 bits для мантиссы\n",
    "    <img src=\"./bfloat16.png\">\n",
    "- Усечённая одинарная точность из стандарта IEEE\n",
    "- Какая разница между float32 и float16?\n",
    "- Этот формат используется в Intel FPGA, Google TPU, Xeon CPUs и других платформах"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor Float от Nvidia ([блог об этом формате](https://blogs.nvidia.com/blog/2020/05/14/tensorfloat-32-precision-format/))\n",
    "\n",
    "- Сравнение с другими форматами\n",
    "\n",
    "<img src=\"./tensor_float_cf.png\">\n",
    "\n",
    "- Результаты\n",
    "\n",
    "<img src=\"./TF32-BERT.png\">\n",
    "\n",
    "- PyTorch и Tensorflow, поддерживающие этот формат, доступны в [Nvidia NGC](https://ngc.nvidia.com/catalog/all)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Смешанная точность ([документация от Nvidia](https://docs.nvidia.com/deeplearning/performance/mixed-precision-training/index.html))\n",
    "\n",
    "- Основная идея:\n",
    "    - Поддерживать копии весов в одинарной точности\n",
    "    - Тогда на каждой итерации\n",
    "        - Сделать копию весов в половинной точности\n",
    "        - Сделать проход вперёд с весами в половинной точности\n",
    "        - Умножить значение функции ошибки на множитель $S$\n",
    "        - Вычислить градиент в половинной точности\n",
    "        - Умножить градиенты весов на $1/S$\n",
    "        - Обновить веса \n",
    "    - Множитель $S$ это гиперпараметр\n",
    "    - Постоянный: такое значение, что его произведение на максимальное абсолютное значение градиента меньше чем 65504 (максимальное число, представимое в половинной точности).\n",
    "    - Динамическое обновление, основанное на статистиках текущих градиентов\n",
    "- Сравнение производительности\n",
    "<img src=\"./mixed_precision_res.png\" width=500>\n",
    "\n",
    "- Существуют расширения для автоматического включения этой опции, больше подробностей [тут](https://developer.nvidia.com/automatic-mixed-precision)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Альтернатива стандарту IEEE 754\n",
    "\n",
    "Недостатки IEEE 754:\n",
    "- переполнения до бесконечности или нуля\n",
    "- много разных NaN\n",
    "- невидимые ошибки округления (об этом ниже)\n",
    "- точность либо очень хорошая, либо очень плохая\n",
    "- субнормальные числа – числа между 0 и минимально возможным представимым числом, то есть мантисса субнормального числа начинается с ведущего 0\n",
    "\n",
    "Концепция posits может заменить числа с плавающей точкой, см [статью](http://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf)\n",
    "\n",
    "<img width=600 src=\"./posit.png\">\n",
    "\n",
    "- выражают числа с некоторой точностью, но указывают предел изменения\n",
    "- отсутствуют переполнения до бесконечности или нуля\n",
    "- пример хранения числа\n",
    "\n",
    "<img width=600 src=\"./posit_example.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Пример 1: потеря точности при делении"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9259246\n",
      "0.1040364727377892\n",
      "-5.9604645e-08\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "#c = random.random()\n",
    "#print(c)\n",
    "c = np.float32(0.925924589693)\n",
    "print(c)\n",
    "a = np.float32(8.9)\n",
    "b = np.float32(c / a)\n",
    "print('{0:10.16f}'.format(b))\n",
    "print(a * b - c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Пример 2: потеря точности при извлечении корня"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0000001468220603\n"
     ]
    }
   ],
   "source": [
    "#a = np.array(1.585858585887575775757575e-5, dtype=np.float)\n",
    "a = np.float32(5.0)\n",
    "b = np.sqrt(a)\n",
    "print('{0:10.16f}'.format(b ** 2 - a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Пример 3: потеря точности при вычислении экспоненты"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "a = np.array(2.28827272710, dtype=np.float32)\n",
    "b = np.exp(a)\n",
    "print(np.log(b) - a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Выводы\n",
    "\n",
    "- Для некоторых чисел обратные функции дают неточный ответ\n",
    "- Относительная точность должна сохраняться в соответствии со стандартом IEEE\n",
    "- Это требование не выполняется на многих современных GPU\n",
    "- Подробный анализ выполнимости стандарта IEEE 754 на GPU от NVIDIA можно найти [здесь](https://docs.nvidia.com/cuda/floating-point/index.html#considerations-for-heterogeneous-world) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Потеря значимых цифр\n",
    "\n",
    "- Многие операции приводят к потере значимых цифр в результате, этот эффект называется [loss of significance](https://en.wikipedia.org/wiki/Loss_of_significance)\n",
    "- Например, при вычитании двух близких больших чисел результат будет иметь меньше правильных цифр, чем исходные цифры\n",
    "- Это связано с алгоритмами и их свойствами (прямой/обратной устойчивости), которые мы обсудим далее"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "123456789999999.11\n",
      "123456789999998.1\n",
      "1.01562500000000000\n"
     ]
    }
   ],
   "source": [
    "a = np.float64(123456789999999.11)\n",
    "print(a)\n",
    "b = np.float64(123456789999998.1)\n",
    "print(b)\n",
    "c = a - b\n",
    "print('{0:10.17f}'.format(c))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Алгоритм суммирования\n",
    "\n",
    "Однако ошибка округления зависит от алгоритма вычисления.\n",
    "\n",
    "- Рассмотрим простейшую задачу: дано $n$ чисел с плавающей точкой $x_1, \\ldots, x_n$  \n",
    "\n",
    "- Необходимо вычислить их сумму\n",
    "\n",
    "$$ S = \\sum_{i=1}^n x_i = x_1 + \\ldots + x_n.$$\n",
    "\n",
    "- Простейший алгоритм: складывать числа $x_1, \\ldots, x_n$ одно за другим\n",
    "\n",
    "- Какова ошибка такого алгоритма при работе в неточной арифметике?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Наивный алгоритм\n",
    "\n",
    "Сложим числа одно за другим: \n",
    "\n",
    "$$y_1 = x_1, \\quad y_2 = y_1 + x_2, \\quad y_3 = y_2 + x_3, \\ldots.$$\n",
    "\n",
    "- В **худшем случае** ошибка пропорциональна $\\mathcal{O}(n)$, в то время как **средне-квадратичная** ошибка $\\mathcal{O}(\\sqrt{n})$.\n",
    "\n",
    "- **Алгоритм Кахана** (Kahan algorithm) даёт ошибку в худшем случае $\\mathcal{O}(1)$ (то есть она не зависит от $n$).  \n",
    "\n",
    "- <font color='red'> Можете ли вы найти алгоритм с ошибкой $\\mathcal{O}(\\log n)$? </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Алгоритм Кахана (Kahan summation)\n",
    "\n",
    "Следующий алгоритм даёт ошибку $2 \\varepsilon + \\mathcal{O}(n \\varepsilon^2)$, где $\\varepsilon$ – машинная точность.\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",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/envs/pytorch/lib/python3.6/site-packages/jax/lib/xla_bridge.py:122: UserWarning: No GPU/TPU found, falling back to CPU.\n",
      "  warnings.warn('No GPU/TPU found, falling back to CPU.')\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error in np sum: 6.0e-07\n",
      "Error in Kahan sum Numba: -1.3e-07\n",
      "Error in Kahan sum JAX: -1.0e-02\n",
      "Error in dumb sum: -1.0e-02\n",
      "Error in math fsum: 1.3e-10\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 ** 8\n",
    "sm = 1e-10\n",
    "x = jnp.ones(n, dtype=jnp.float32) * sm\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": [
    "## Ещё один пример"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0\n",
      "0.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "test_list = [1, 1e20, 1, -1e20]\n",
    "print(math.fsum(test_list))\n",
    "print(np.sum(test_list))\n",
    "print(1 + 1e20 + 1 - 1e20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Выводы по операциям с числами с плавающей точкой\n",
    "\n",
    "Необходимо быть особенно осторожными при работе с числами с плавающей точкой, поскольку операции с ними могут дать неверный ответ из-за ошибок округления.\n",
    "\n",
    "Для многих стандартных алгоритмов вопросы устойчивости хорошо изучены и проблемы, связанные с этим, легко обнаруживаются."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Векторы\n",
    "\n",
    "- В рамках курса, мы будем работать не с числами, а с векторами\n",
    "- Вектор в фиксированном базисе размера $n$ может быть представлен как одномерный массив из $n$ чисел.\n",
    "- Обычно вектор рассматривается, как матрица $n \\times 1$, то есть вектор-столбец\n",
    "\n",
    "**Примеры:** \n",
    "- Многочлены степени $\\leq n$ образуют линейное пространство.\n",
    "- Многочлен $ x^3 - 2x^2 + 1$ может быть представлен в виду вектора $\\begin{bmatrix}1 \\\\ -2 \\\\ 0 \\\\ 1\\end{bmatrix}$ в базисе $\\{x^3, x^2, x, 1\\}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Векторные нормы\n",
    "\n",
    "- Векторы обычно содержат приближённое описание физических или каких-нибудь других объектов\n",
    "\n",
    "- Один из главных вопросов – это насколько приближение точное (1%, 10%)\n",
    "\n",
    "- Что является достаточно точным представлением объектов, конечно, зависит от конкретного приложения. Например:\n",
    "    - При решении уравнений в частных производных чаще всего встречается точность порядка $10^{-5} - 10^{-10}$\n",
    "    - В приложениях связанных с обработкой данных иногда точность порядка $80\\%$ приемлема при серьёзном зашумлении исходных данных."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Расстояния и нормы\n",
    "\n",
    "- Норма – это количественная мера малости вектора, обычно обознается как $\\Vert x \\Vert$.\n",
    "\n",
    "Норма должна удовлетворять следующим свойствам:\n",
    "\n",
    "- $\\Vert \\alpha x \\Vert = |\\alpha| \\Vert x \\Vert$\n",
    "- $\\Vert x + y \\Vert \\leq \\Vert x \\Vert + \\Vert y \\Vert$ (неравенство треугольника)\n",
    "- Если $\\Vert x \\Vert = 0$, то $x = 0$\n",
    "\n",
    "Расстояние между двумя векторами можно определить как норму разности между ними\n",
    "\n",
    "$$\n",
    "   d(x, y) = \\Vert x - y \\Vert.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Стандартные нормы\n",
    "\n",
    "Наиболее известная и чаще всего используемая норма – это **евклидова норма**:\n",
    "\n",
    "$$\\Vert x \\Vert_2 = \\sqrt{\\sum_{i=1}^n |x_i|^2},$$\n",
    "\n",
    "расстояние по которой соответствует расстоянию в реальном мире. Если векторы состоят из комплексных чисел, используются их модули."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $p$-норма\n",
    "\n",
    "Евклидова норма или $2$-норма – это частный случай важного класса $p$-норм:\n",
    "\n",
    "$$\n",
    " \\Vert x \\Vert_p = \\Big(\\sum_{i=1}^n |x_i|^p\\Big)^{1/p}.\n",
    "$$\n",
    "\n",
    "Также важны следующие два частных случая :\n",
    "- Бесконечная норма или Чебышёвская норма определяется как максимальный модуль элемента вектора $x$: \n",
    "\n",
    "$$\n",
    "\\Vert x \\Vert_{\\infty} = \\max_i | x_i|\n",
    "$$\n",
    "\n",
    "<img src=\"chebyshev.jpeg\" style=\"height\">\n",
    "\n",
    "- $L_1$ норма определяется как сумма модулей элементов вектора $x$: \n",
    "\n",
    "$$\n",
    "\\Vert x \\Vert_1 = \\sum_i |x_i|\n",
    "$$\n",
    "\n",
    "Расстояние, определённое как $L_1$ норма разности между векторами, называется **Манхэттенским расстоянием**\n",
    "\n",
    "<img src=\"manhattan.jpeg\" style=\"height\">  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Эквивалентность норм\n",
    "\n",
    "Все нормы эквивалентны в следующем смысле\n",
    "\n",
    "$$\n",
    "   C_1 \\Vert x \\Vert_* \\leq  \\Vert x \\Vert_{**} \\leq C_2 \\Vert x \\Vert_*\n",
    "$$  \n",
    "\n",
    "для некоторых положительных констант $C_1(n), C_2(n)$, $x \\in \\mathbb{R}^n$ для любых пар норм $\\Vert \\cdot \\Vert_*$ и $\\Vert \\cdot \\Vert_{**}$. Эквивалентность норм означает, что если некоторый вектор мал в одной норме, то он мал и в другой норме. Однако константы могут быть большими."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Вычисление норм в Python\n",
    "\n",
    "Пакет NumPy содержит всё необходимое для вычисления норм: функция ```np.linalg.norm```.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative error in L1 norm: 0.0008498416786543509\n",
      "Relative error in L2 norm: 0.001056723991554265\n",
      "Relative error in Chebyshev norm: 0.0030252783065499726\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 100\n",
    "a = np.ones(n)\n",
    "b = a + 1e-3 * np.random.randn(n)\n",
    "print('Relative error in L1 norm:', np.linalg.norm(a - b, 1) / np.linalg.norm(b, 1))\n",
    "print('Relative error in L2 norm:', np.linalg.norm(a - b) / np.linalg.norm(b))\n",
    "print('Relative error in Chebyshev norm:', np.linalg.norm(a - b, np.inf) / np.linalg.norm(b, np.inf))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Единичные диски в различных нормах\n",
    "\n",
    "- Единичный диск – это множество точек такое что $\\Vert x \\Vert \\leq 1$\n",
    "- Для евклидовой нормы единичный диск сопадает с обычным диском\n",
    "- Для других норм единичный диск сильно отличается от привычного нам диска"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Unit disk in the p-th norm, $p=1$')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt8lPWd6PHPdyYJF0EI4SIQEkSoBaGiREG71daqVU8r3hW6q91Ti+zq6atnT3el2ro9aLvYs/uqp2fdRbT2spWL9YLU1ZdXFHc1SEJRblJiJCFyi2G4Q27zPX/MM+wkzEwmmWfmeWbm+369Bmae6/dJJvOd53cVVcUYY4yJCngdgDHGGH+xxGCMMaYLSwzGGGO6sMRgjDGmC0sMxhhjurDEYIwxpgtLDMYYY7qwxGCMMaYLSwwmLhHZLCJfTmP/X4vIQ705lojsEJHLsxFfX8+bDX6LxxQeSwx5SkRURCZ2W/ZjEfldKvur6jmq+pazX1ofVLHHcktfj+m3D12/xZMrROQeEakRkVYR+bXX8eSbIq8DMMZ4T0SKVLXD6zh6YRfwEPA1YIDHseQdu2MoUM431e+LyIciclBEVohI/27rLxeRfwMqgD+IyBER+bsExztPRNaLyGERWQGccqyY1/eKyKfOtttE5Ktxjvd5EflERG5LEv/lqVxLzD7JrmV6kp/FGBF5VkSanZi+28PP9QciskVEQiLyq3ixpBNPnPMl+z1OFpG3ROSAU/x2bbd97xWRD4GjIlLkLPtb53hHReSXIjJKRF52fl+vi0hpouvvFts3ReRdJ6Y9IrJTRK5OZd+eqOpzqroSaHHjeKYbVbVHHj4ABSZ2W/Zj4HfO8x3A+8AYYBiwFZgfs+0O4PLuzxOcqwRoAP4nUAzcBLQDD8U51tnATmCM83o8cFbsdsD5QCPw9STn7B5fwmtJtF+3ZXH3J/LlqRZ4wLnOCUA98LUkx98EjHOO9Z/Rn4Mb8STYP1HsxUAdcJ8T+2XAYeDsmH03OLEOiFlWDYwCxgL7gPXAeUA/4E3g71N8Dy4CjgO3OLF8H2iIs92LwIEEjxd7OMdDwK+9/nvLt4fdMRS2X6jqLlXdD/wBmN7H48wi8of/iKq2q+ozwLoE23YS+YCZIiLFqrpDVT+OWf8lYBVwh6q+2IsY0r2WRPtfAIxQ1YWq2qaq9cDjQNw7Gcc/q+pO51g/Aeb0MpZk8fRm21nAIGCRE/ubRD6E53Tbd6eqHo9Z9v9Uda+qfgq8A6xV1T+qaivwPJEkkYppwM9V9WlVbQd+C1R0v/tR1a+r6tAEj6+neC7jIksM+auTyId1rGIi3+Sj9sQ8P0bkQ6QvxgCfqvMVztEQb0NVrQO+R+TuZZ+ILBeRMTGbzAfeVdXVvYwh3WtJtH8lMMYpijkgIgeIfAMfleRYO2OeNxD5+USLVo44j5f7GE9vth0D7FTVcLd4xiaINWpvzPPjcV6n+rOdBjwT83okcERVT6S4v/GIJYb81UikmCbWmST4wO5BT5N27AbGiojELKtIeDDVpar6Z0Q+dBV4OGb1fCLfKn/ehzhT0dsJSHYCn3T7FjtYVa9Jss+4mOcVRCpKUdWnVHWQ84iWtWdyQpRdwDgRif07rwA+jXmdkfOLyFAiP4fmmMU3AackRKf+4kiCR08J1GSAJYb8tQL4oYiUi0jAqaj9Bl2/waVqL5Gy9UTeAzqA7zoVmDcAF8bbUETOFpHLRKQfcILIN9DOmE0OA1cBl4jIoj7E2pOerqW794FDTiXtABEJishUEbkgyT53Oz/3YUTuLla4GE9vrAWOAn8nIsUS6ffxDWC5WyeQSH+VX8dZNY3I73Wu8574b8BfE7lT7EJVr45JmN0fcSurnWP2B4JAUET6i4i1snSJJYb8tRB4F/gPIAT8DPimqm7qw7H+gUiSOSAi3+++UlXbgBuAbznnuhV4LsGx+hGplPyMSBHISCIfnrHHOwBcAVwtIg/2Id5kkl5Ld6raSeTDdDrwiRP3E8CQJLstBV4lUkldT6SC1JV4esP5vVwLXE0k7n8BblfVj1w8zTgiFezdTQOeAi4i8p7438B1qrrFpfP+kMiXigXAnzvPf+jSsQuedC0WNsakQ0R2AHeq6utex5JpIlICfAB8walcjl33r8CfVDVTRYImg+yOwRjTJ05Lp8ndk4JjGpGmsyYHWWIwxmTCVMDNIiuTRVaUZIwxpgu7YzDGGNNFTjbvGj58uI4fP97rMIwxJqfU1tZ+pqojetouJxPD+PHjqamp8ToMY4zJKSKSUgdXK0oyxhjThSUGY4wxXVhiMMYY04UlBmOMMV1YYjDGGNOFK4lBRJ4UkX0iEneANon4hYjUOVMGnh+z7g4R2e487nAjHmOMMX3n1h3Dr4kMlZzI1cAk5zEP+FcAZ1jivwdmEhmm+e9TnU/WGD+pbQjx6Oo6ahtCXodiTNpc6cegqmtEZHySTWYDv3Vm+KoWkaEiMhr4MvCaMyUhIvIakQSzzI24jMmG2oYQc5a8R3unUhwUls27iBmV9v3G5K5s1TGMpesUgk3OskTLTyEi80SkRkRqmpub421ijCeeXd9EW6eiQFun8uz6Jq9DMiYt2UoMEmeZJll+6kLVJapapapVI0b02KPbmKzp/iZ+/5P9VqRkclq2EkMTXefBLScyH22i5cbkhNqGEPsOt3ZZVrfvCLcuec+Sg8lZ2RoraRVwj4gsJ1LRfFBVd4vIK8BPYyqcrwR+kKWYjElLbUOIOY9X09YRPmVdh1OktG3PYV7etJurp45m7swKD6I0pvdcSQwisoxIRfJwEWki0tKoGEBVFwMvAdcAdcAx4C+ddfudOX3XOYdaGK2INsbvnlvfFDcpRG1oDLF0bSMA72z/DMCSg8kJbrVKmtPDegXuTrDuSeBJN+IwJhtqG0I8u76JNdv2Jd1u/9G2Lq+f/M9PLDGYnJCTw24b44XahhDPrW9ixbpGktwoABAMQMWwgew59F/1D3X7jjDvtzXcdelZ1pzV+JolBmNSUNsQ4ptPVNPaHo7fbC7GsIHFHGnrpCZO5fOrW/ayZnszT905y5KD8S0bK8mYFFTXt9DW0XNSAJg4chAdnWHCCTZuaw9TXd/ianzGuMkSgzE9qG0IsevA8QQ9bE61adehhEkBIAwcPt7uSmzGZIIVJRmTRLQIqa0jTA/VCicda+vscZvH1tRT/9lRq28wvmR3DMYkES1CSnYH0BdKpL5hzuPV1hHO+I4lBmOSKB1YQkAk7tgtbmjrCHPHk2tZ9NLWDJ3BmN6zoiRjuqltCFFd30LpwBIWvriZDud2YeiAYg6kWTdQOWwgDfuPdVl2pLWTxWvq2XPoBI/cdl5axzfGDZYYjHFEO649U9tER2eYgMjJpACknRQAOsOJaypWbtjFX1w03uocjOcsMRhD/HGPwupyxQLQdOBE0vXV9S2WGIznLDEYAzz88tak4x5ly9r6Fl7dvIeRp/dnvrVYMh6xxGAK3qKXtvL+Dn+0DFrjDLYHB1n90V5W3HWxJQeTddYqyRS02oYQi9fUex1GXB1hbDY44wlLDKagPefzD95l7zdaU1aTdVaUZApSdKTUFz/094SBqpy8o1lwzWSPozGFwhKDKSjRJqkr1u2k0+3uzBm0eE09FWWn2XwOJivcmsHtKuD/AkHgCVVd1G39z4GvOC8HAiNVdaizrhPY6KxrVNVr3YjJmO56M3S2H933/EYaW47anYPJuLQTg4gEgUeBK4AmYJ2IrFLVLdFtVPV/xmz/P4DY7p3HVXV6unEYk8zStY0sWfMxJ9q9b5KaDrtzMNngRuXzhUCdqtarahuwHJidZPs5wDIXzmtMSpaubeS+5zeyo+VY3PX9itL/Mxg2sDjtY6RqxbrGrJ3LFCY3EsNYYGfM6yZn2SlEpBI4E3gzZnF/EakRkWoRuS7RSURknrNdTXNzswthm0JQ2xDiH1/dlnSbVhc6tp0+IHuJocSFRGZMMm68w+INPJmoCPc24BlVjR2wvkJVq4C5wCMicla8HVV1iapWqWrViBEj0ovY5L3ahhD3Pb+RW5e8x/6jbRk/X6K7kUyoaQgx77c1Nly3yRg3Kp+bgHExr8uBRG0AbwPujl2gqruc/+tF5C0i9Q8fuxCXKVC5XsncE9XIXA5vfLSXqWOGcOsFFVbnYFzlxh3DOmCSiJwpIiVEPvxXdd9IRM4GSoH3YpaVikg/5/lw4IvAlu77GtMbvZmfOZd1huGDpoPc9/xG6wRnXJV2YlDVDuAe4BVgK/C0qm4WkYUiEtv0dA6wXLXLkJWTgRoR+QBYDSyKbc1kTF8cPt5OBgZG9bUl79Rb0ZJxjSv9GFT1JeClbsse6Pb6x3H2exeY5kYMxtQ2hFj89se8tmWv16FkXVhtyG7jHuv5bPLCope2suSdetfnZs4lh12YSMgYsMRg8sDStY2+HSE1mxavqefFD3cxecwQm8vBpMUaRJuc9/Km3acsGzogf77zFAfjtQiPr+nACV7bspdbH3vX6hxMn+XPX48pOEvXNvLypt1xO9IcON6R9Xgypb2z9+VjHWGrczB9Z4nB5KToMBe5qHLYQAaWBNm653BGz/P7mp18euA4N55fbgnC9IolBpOT4hUfCYm73PtJw/7s9JLe0XKMHS2NLH+/kcsnj+Iuq3cwKbI6BpNTahtCPLq6jrLTSk5ZlwtJwQthp6f0nCXvWb2DSYndMZicUdsQYs7j1bR3hJHU62Pz1hmD+9GpyoHj7SnVQ7R1qtU7mJRYYjA547n1TbQ5I6EWWs/mePYcbgUgGIi0XOroVIqLAhQJHEsw70TpwFPvtIzpzhKD8b3ahhDV9S2sL4BikL7Uk3SGQUVRoKMjTLJubm9t20foWBuzJpTZnYNJyBKD8bV8HCk1IMTtoZ1oeSqi+/U0s8QbH+3j9a17KSkK8NSdsyw5mLgsMRjfqm0I8cjrf8r56Ti783LYjk7n5Cfawzz29secO26o3T2YU1hiML60dG0jP3ph08kPskKQ7Ut9dcteu3swcVlzVeM7tQ0hfrRyY0ElBa+EFdo7wlTXt3gdivERu2MwvlHbEOKxtz/mw6YD9GEUCNNHnQpvb9tnRUrmJFfuGETkKhHZJiJ1IrIgzvpviUiziGxwHnfGrLtDRLY7jzvciMfkntqGELc89h6vbtnLnkOtXodTcN7fEeI26wBnHGnfMYhIEHgUuILI/M/rRGRVnJnYVqjqPd32HQb8PVBFpJVerbOvvTsLzGNvf2xFR2kY3C/IVyePYtOuQ9TtO9KnY7R3Knc/VcsXyofa8BkFzo07hguBOlWtV9U2YDkwO8V9vwa8pqr7nWTwGnCVCzGZHLJ0bSOvbS28WdfcdLi1k5UbdjF1zOkE0ugVvudQK69u2ctNi23Y7kLmRmIYC+yMed3kLOvuRhH5UESeEZFxvdzX5Kmlaxv54cqNBduT2e2RPVZu2OXKcCGq8KOVuTl6rUmfG4kh3tuw+5/5H4DxqvoF4HXgN73YN7KhyDwRqRGRmubm5j4Ha/yjtiHEAy9sKujpODNx6Z0udfuo23eER1fX2Z1DAXIjMTQB42JelwO7YjdQ1RZVjdYoPg7MSHXfmGMsUdUqVa0aMWKEC2Ebr1XXt1i9go+1dSr/9Oo2vvlEtSWHAuNGYlgHTBKRM0WkBLgNWBW7gYiMjnl5LbDVef4KcKWIlIpIKXCls8zksdqGEPc/v5G3tu3Lm2EuAIYPLmFQSdDrMFxl/RwKU9qtklS1Q0TuIfKBHgSeVNXNIrIQqFHVVcB3ReRaoAPYD3zL2Xe/iDxIJLkALFTV/enGZPwrn3s0n2jt5EhbZ6/28fvkQgGB4qIApQNLeHR1nfV1KBCiOVjrV1VVpTU1NV6HYXqptiHELYvftc5rGVIUEDpcSrjBgPCdPzuTw60d/GnvYWobQoQVSooCLPuODZ+Rq0SkVlWretrOhsQwWfPwy1stKWSQW0kBIBxWDrd2sGJdI+t2hE42EGjrCPPc+ibXzmP8yRKDyYrvLf8j7+/IzQrMYCH+lQjsO9xKR5wWTps+PWiV0XmuEN/yJssWvbSVlRviNjbLCW41/8wlqon7WHzQdNDmj85zNoieyajahhCPvVPvdRimD0LH2igOStz5pNs6le8uW8/wQf249YIK5s6s8CBCkymWGIzrolNxzppQxnPrmwq2V3Ouq2kIcdeXJvBefQsfNB08Zf2nB07w6YETfNAU6SFtySF/WGIwroqdilOAQf3tLZZtJUGhzYVaflV4bE09d10ygU27DiVtYvyzVz4CLDnkC6tjMK56dn0TJ5z5mcPAoRMdXodUcNxIClEKLF5Tz+dGDkq63YFj7dz3/EaWrm107dzGO5YYjGtqG0KsWGcfDH4RdHGEvo/2HE5pO/v95wdLDMYVtQ0h7n3mg5xuwZPOcNVumV4+xLURV93sM5LqoXaGjtldQx6wAmCTttqGEHOWvOdqEYYX/DBKx4Y4lbzZFk1Mfflx7D8aKVICq2/IZXbHYNJWXd/i66QwoNje5r1RWTaQ2dPHpHWMJ//zE5eiMV6wOwbTZ9Fmqdv3plb+7JXj7TlcvuWBhpZjNOw/ltYx6vcdobYhZGMq5ShLDKZPos1S2zrCvh4d1PSenvyn78LAwj9sZtTp/Rk+uB83nl9uSSKHWGIwfVJd30JbR9gX5fLGnyKd4iJ1Js/U7GTZvIssOeQIK3w1vRKdZOftbfu8DsXkkLZOZeEfNtv4SjnC7hhMymobQsx5PFJ8ZExUQFJr0RUdfM/uHPzPlTsGEblKRLaJSJ2ILIiz/m9EZIuIfCgib4hIZcy6ThHZ4DxWdd/X+ENtQ4hHXv+TJYU8JfSuQ1x004D0rv9HW6fyrM3n4HtpJwYRCQKPAlcDU4A5IjKl22Z/BKpU9QvAM8DPYtYdV9XpzuPadOMx7oveKbyz/TOvQzEZovSuQ1x007ASd86GZJ6pbbIiJZ9z447hQqBOVetVtQ1YDsyO3UBVV6tqtP1bNVDuwnlNhtU2hHh0dR2Pvf1xTtwp+KDj8ikG9wt6HYLvtHeEeeT1P1ly8DE36hjGAjtjXjcBM5Ns/23g5ZjX/UWkBugAFqnqShdiMmmKrU/w4wduPBeML/XdLHGHWzu9DsF3FPiP7Z+xbsd+nrrT5o/2IzcSQ7zPjbg3pSLy50AVcGnM4gpV3SUiE4A3RWSjqn4cZ995wDyAigrrap9pz61vOnmXEO+XKQmWe2l9o7+SgklMidw5VNe3WGLwITeKkpqAcTGvy4FT5nEUkcuB+4FrVbU1ulxVdzn/1wNvAefFO4mqLlHVKlWtGjFihAthm2R6+tAXH95G5PIAfoXq1c17bNA9H3IjMawDJonImSJSAtwGdGldJCLnAY8RSQr7YpaXikg/5/lw4IvAFhdiMmmaOmZI0jdHPnRs82Fu86URg0oyctxOjTRhtXkc/CftoiRV7RCRe4BXgCDwpKpuFpGFQI2qrgL+DzAI+L1Evmo2Oi2QJgOPiUiYSJJapKqWGDxS2xDi2fVN1O09TG1DyHdFRT3pbby5dn29Mb18CAeOt7OjJb0xjwCaj7S5EFFyS9Z8zNlnDLZiJZ8QzcEJeauqqrSmpsbrMPLK0rWN/GjlRlfH8Dfe8WMdUE9KgsLNVeO4wcZVyhgRqVXVqp62syExDLUNIR54YZMlhTySi7/Ktk7lqbWN3LrkPWvK6jFLDIbq+hbCOXjnmIsCQNlpxV6H4Wsdncpjb5/SMNFkkSUGw6wJZRT5YV7LAhAGWo62ex2G773/yX6vQyhoNohegYpOsjNrQhkAI07vz6eh4x5HZUzEgePtXPqz1Rw+0c7EkYO49+rJVu+QRZYYCtDStY3c//xGlPTm9zUmk6KzyL2/I8TNi9/l9/MvtuSQJVaUVGCWrm3kPicpQCQhWFIwfhdWrN4hiywxFJDahhA/fH6j12EY0ydvfLTPWitliSWGArL47Y+xUSO8YVX76esMq43KmiWWGArIJ58dPWVZ5bCBjB3a34NoCosV17njne2f8c0nqi05ZJhVPue5aOuj0oEl7Gg5NTHsDKU/ZIIx2WSjsmaeJYY8tuilrSx5p56wQlFA4nZiy4fB8AqBAGOG9mf3wRMF/zsLBoTSgSU8urqOWRPKLEFkgCWGPLV0bSOL19SffN0RVoIBAY96OOfi2D1+osCugye8+vVl1djSAUn71EwfN5SFL26mrSNMSVHAJvvJAKtjyFP/snr7KcsyNXxyKgrg8yzjspEURgwqwetO8Lt66GhZ2xCitT1MWOFEe5h7n/3Q6hxcZokhz9Q2hJj32xqaDpw4Zd2eQ61x9jB+NyiL80Y3H2nzvKiqp9OHtes2dfuOcMtjNvCemywx5JHahhDffKKaV7fs7fW+wQBMPmNwBqIy6TqSR/NGu3E3Eu8QnWEbeM9NlhjySHV9Cyfa+9ZToTMM25uPuByR6Umh9W8oO62EsaUDGDKw79Wbid7hr23Za3cNLnElMYjIVSKyTUTqRGRBnPX9RGSFs36tiIyPWfcDZ/k2EfmaG/EUotqGELsOpDcIXodNyJB1XpfnZ1vzkTY+DR3n4LEO14+tRDpxmvSl3SpJRILAo8AVQBOwTkRWdZui89tASFUnishtwMPArSIyhcgc0ecAY4DXReRzqpo/985ZkGz2tfKh/akaP4yVG3ZlPzDTI8vF7tq662CXkYOttVLfuNFc9UKgTlXrAURkOTAbiE0Ms4EfO8+fAf5ZIpM/zwaWq2or8ImI1DnHe8+FuArC0rWN3L9yY8IWK00HTrDrA0sKpjA0HTjBzYvfPdl3Z+HsqcydWeF1WDnHjaKkscDOmNdNzrK426hqB3AQKEtxX5NAdErOnpoxet3KxJhsir7fO8LKAy9ssnqHPnAjMcQrJe3+UZRom1T2jRxAZJ6I1IhITXNzcy9DzE/V9S10+vxTXwqsDN34S1iV6voWr8PIOW4khiZgXMzrcqB72cXJbUSkCBgC7E9xXwBUdYmqVqlq1YgRI1wIO/fNmlBGv2L/NiwLSG41extUEiy4VkL5TICSosDJWQpN6tz4u10HTBKRM0WkhEhl8qpu26wC7nCe3wS8qarqLL/NabV0JjAJeN+FmArCjMpSnrpzFl+aNNzrUOIKa25Vrh5p68xID+0BRe6nx1xKuF64YHwp3//a2TZcRh+lXfmsqh0icg/wChAEnlTVzSKyEKhR1VXAL4F/cyqX9xNJHjjbPU2koroDuNtaJPXOjMpSvnf551i3Yz+t7eGsDj1h4x+l5niH+7Ng5Nq8Gtl8r5QEhQU2R3RaRHNwVK6qqiqtqanxOgxfqW0I8ez6Jlas2+n7egdjMuVLk4bzvcs/Z0khARGpVdWqnraz0VXzxIzKUmZUlrL504N80HTQ63CMybqiAJYUXGJFlXnm1guszbYpTLdeUGFJwSWWGPLM3JkVXDd9jNdhmAKVzZFgowToXxzghvPLs37ufGWJIQ89ctt5zL9kgi+aXgaAS7LQamr4oBKGDrCSUa8NH9Qvq+cT4M8mDbfWRy6zxJCnFlwzmWf+6mLOLR+S0vaZaFIJkdYza7Z/lpFjx9p/tI0xQwdk/DwmuYaW7M0hLkC/4oDVK2SAfcXKYzMqS3ngG+fwzSeqaWsPJ23iKDk+zGdYYcvuw16HUfCy1R7uuuljmDRqsA2UlyGWGPJctBNcdX0La+tbEn57P9aWeveRoQOLOXCs3a0QjQuied3rlsoByU4Mk0YN5u6vTMz8iQqUJYYCEG3KCvDO9s/S/lbn16RQFIAM9CVzVaY6eqX6YdyvKEBrBn5I0fvNbHSLKgqIDXORYVbHUECiYysFJfLHlUhAcnNmsUkj/T81qdddDzORFCByXdFHJgUFFs6easVHGWZ3DAUktlipdGAJP3phU9xe0l4XR/TVR3usjiGfnVs+hAe+cY4lhSywxFBgYouVnv9jE+t2pD5W/WklQY72oi4i23I0n5kUFAXFkkIWWVFSAfvcqN4VvQzuZ98jTHYNG1jMlVNGsWLeRZYUssj+0gvYDeeX83TNTtpTHBt7z+HWDEeUmI3kWniKAvD4HRdYQvCA3TEUsBmVpSyfdxFXThnldSg9Sicp5HgXjYIk2NhHXrLEUOBmVJay5PaqjCQHvwxREbTMkFMCEunRbGMfeccSgwHgrkvPcvWbdUDg4IkO9w4YhwCVwwb2uF2qRWXGewJ8caKNfeS1tBKDiAwTkddEZLvz/ym/SRGZLiLvichmEflQRG6NWfdrEflERDY4j+npxGP6bkZlKQ9dN821440eOiDjnZ1mTx/DLReMy8k+Fya+YEBs7CMfSPeOYQHwhqpOAt5wXnd3DLhdVc8BrgIeEZGhMev/VlWnO48NacZj0jB3ZgU/vX5alzuH8qH9+3SsT0PHXYoqsRc27GLp2garlM4TAeu85hvpFgLPBr7sPP8N8BZwb+wGqvqnmOe7RGQfMAI4kOa5TQbMnVnB2WcMprq+hcPH21m8pt7rkBJS4NMDJ7wOI+8M7hfkcGt2+6sEBR68bhpzZ9pEU36Q7h3DKFXdDeD8PzLZxiJyIVACfByz+CdOEdPPRSS7g7mbuGZUlnL3Vyayefchr0NJixUx9U02k4IQ6dH89PyLLSn4SI+JQUReF5FNcR6ze3MiERkN/Bvwl6oaHbDlB8DngQuAYXS72+i2/zwRqRGRmubm5t6c2vTR1VNHex1CWqyIyf+KrUezL/VYlKSqlydaJyJ7RWS0qu52Pvj3JdjudODfgR+qanXMsXc7T1tF5FfA95PEsQRYAlBVVWV/81kwd2YFjS1HWfJOfc6On2T8IV4HRQFurhpnScGH0i1KWgXc4Ty/A3ih+wYiUgI8D/xWVX/fbd1o538BrgM2pRmPcdmCaybze+c2vyRDs7wZf5AMlr3FSwrWV8G/RNNoUygiZcDTQAXQCNysqvtFpAqYr6p3isifA78CNsfs+i1V3SAibxKpiBZgg7PPkZ7OW1VVpTU1NX2O2/RNbUOIe5/9kLp9XX9FA0sCHGvz+UQIHhNggP2cACgJCjdXjeOG88vtbiHLRKRWVat62i6tVkmq2gJ8Nc7yGuBO5/nvgN8l2P+ydM4HvvrmAAAQ0klEQVRvsmtGZSkXnjnslMRQVTksK/M6Q+QWNxc/WhUsKQBXTBnF/EvPsoTgc1Y2YHrlxvPLCXYrcshWUoDcTAomYuiAYh6/vcqSQg6wxGB6ZUZlKQ9eN80GpnNZIYzn9HdXfd7rEEyKLDGYXps7s4KHrpuWdHpQk7qAwOdGDvI6jIwZXzaQn15vnddyiT+GvzQ5J9pD+rG3P+bVLXu9DiehXJjHIaz5Oy1pUQD+6ZbpVnyUY+yOwfRZdMjuiT7+tuv3pBCVK3H2xtih/Vlx18WWFHKQJQaTtv/+xTO9DsH4TDAg/GLO+ZYUcpQlBpO2uTMrmH/JBKuQNkAkKTxoo6TmNKtjMK5YcM1krjjnDKrrWygdWMKmXQep23uYmh0ha2JaQM4tH2JjH+UBSwzGNTMqS09+INQ2hJjzeLUlhQISFCwp5AkrSjIZUV3fQnuHpYVC8uB10ywp5AlLDCYjZk0oozgHBt0bPqjE6xAy6vT+RZxxemanOQkK1k8hz1hRksmIGZWlLPvOLJ5b38S+w628tW0f7Z3+a5TZcqQtJ/o69NWhEx0cOtHh+nEnnzGYv7hoPKFjbcyaUGZ3CnnGEoPJmGidw6Or63hjqz87wSlwWkmQo22ZmbVMBNIYwNi3vn7uGLtDyGP+v9c3OW/WhLKkczmMHdqfsoHFWYyoq0wlhclnDOYn101j4ojT8mqa0f7FAWZNKPM6DJNBac3H4BWbjyH31DaE+NHKjWzZ7e+hH07vX+Ra0YsQGQfJhyVofTZx5CAevvELVnSUo1Kdj8HuGExWzKgspWxQZitB0xUMwOfPGOza8ZT8SgoBsKRQINJKDCIyTEReE5Htzv9x3zEi0ikiG5zHqpjlZ4rIWmf/Fc40oCZPXT11dK/3yeR0k93HeOoMw/s7Qpk7YY6bd8kESwoFIt07hgXAG6o6CXjDeR3PcVWd7jyujVn+MPBzZ/8Q8O004zE+NndmBT+9fhrnlg+hKEBK5e6ZLOmsb+5xFlkTY/AA7+qBTHal2yppNvBl5/lvgLeAe1PZUUQEuAyYG7P/j4F/TTMm42NzZ1Ywd2YFtQ0hqutbOHy8nde37mXPwRMcyVAlcCLhPCrm6at+RQFaU+iIWBQQq3AuIOkmhlGquhtAVXeLyMgE2/UXkRqgA1ikqiuBMuCAqkZr+pqAsWnGY3JE7PAZC66ZzKOr6/jHV7blbX8Cv2rtCCdtUitEBsVbaIPiFZQeE4OIvA6cEWfV/b04T4Wq7hKRCcCbIrIROBRnu4SfCyIyD5gHUFFh7afzzawJZfQrDtDeEUZE6LCv81mTKClcOL6US88eaR3YClCPiUFVL0+0TkT2isho525hNLAvwTF2Of/Xi8hbwHnAs8BQESly7hrKgV1J4lgCLIFIc9We4ja5ZUZlKU/dOevk6KwPvLARG2rJOwGBe6+ebAmhQKVb+bwKuMN5fgfwQvcNRKRURPo5z4cDXwS2aKQDxWrgpmT7m8Ixo7KUu78yEYi0EDLeCAaEh2xAvIKWbh3DIuBpEfk20AjcDCAiVcB8Vb0TmAw8JiJhIolokapucfa/F1guIg8BfwR+mWY8JsfVNoT40QubktY15MLYRrkQY3djh/Zn7sxKKzoy6SUGVW0BvhpneQ1wp/P8XWBagv3rgQvTicHkl+r6Fjp7qF/IhQ/cXIixu7u/MsnGPzKADaJnfCY6rlKbVTBkxZTRgykb1I+rp462pGBOssRgfCV2uO43tu5lz6FWr0PKa9MrSvnp9XFv6E0Bs8RgfCfax+GcMUO47/mNXofTRS7WHSRSEhRuPL/c6zCMD1liML4VLdpYsa6RD5oOehxNRC4mhXjJzEZJNcnY6KrG1+bOrODKc+L1rzSpUmDiiNO6LLvwzGGWFExClhiM782aUEYwkE9T3WTf1LFDKAkKghUhmZ5ZYjC+N6OylAdnTyUYiHywFQeFC8fbt93eePHD3fz42ql8/2tns2zeRXa3YJKyOgaTE+bOrODsMwZTXd/CrAllVNe3FMTcCaf3L2JgcZA9h9NrndUZVkLH2k72LDcmGUsMJmfEjsgKEMyzaTPjGTW4H9ubj6Z9nKANm216wYqSTE6aUVnKVyeP8jqMjHMjKQTAhs02vWKJweSsuy49i5IiewsnM3HEafz+ry62Xs2mV6woyeSsaC/p6Exwi9fUex1SRpUP7c+nB06k3JciGBAevulcu1MwvWZft0xOiw7VPXhAcUpzSOeywf2LU04KIvCgFR+ZPrLEYPJCdAa4oEQmmfHSgOLM/Flt3XM4pe0CAj+5bpoVH5k+s8Rg8kJ0Bri/ufJsLve4UnrU6f09O3dQ4CFLCiZNVsdg8ka0OWttQ4jXtu5NOJdxpjWFjnly3iumjGL+pWdZ8ZFJW1p3DCIyTEReE5Htzv+nvCNF5CsisiHmcUJErnPW/VpEPolZNz2deIyBSIK4wsO7Bi+mkrhyyigev73KkoJxRbpFSQuAN1R1EvCG87oLVV2tqtNVdTpwGXAMeDVmk7+NrlfVDWnGYwzgNGUN5nt1dEQwINx16Vleh2HySLqJYTbwG+f5b4Dretj+JuBlVfXmXtsUjBmVpSybdxFXTvFfJ7hB/YKuHStgrY9MBqSbGEap6m4A5/+RPWx/G7Cs27KfiMiHIvJzEemXaEcRmSciNSJS09zcnF7UpiDMqCxlye1VvksOR1o7XTlOMCBW0WwyosfKZxF5HYg3IP79vTmRiIwGpgGvxCz+AbAHKAGWAPcCC+Ptr6pLnG2oqqrK8xFyjJvuuvQs3vxoHx3h/HjblA/tz6Vnj+SG88vtTsFkRI+JQVUvT7RORPaKyGhV3e188O9LcqhbgOdVtT3m2Ludp60i8ivg+ynGbUzKZlSWsnD2VO5fudGzlkpuuuTskfzE5mk2GZRuUdIq4A7n+R3AC0m2nUO3YiQnmSAiQqR+YlOa8RgT19yZFTwz/2IuHF/KYBfL+LPNJtkx2ZBuYlgEXCEi24ErnNeISJWIPBHdSETGA+OAt7vt/5SIbAQ2AsOBh9KMx5iEZlSW8vT8i5lekZvFL+eWD7FJdkxWpNXBTVVbgK/GWV4D3BnzegcwNs52l6VzfmP64uqpo3ln+2cZPYdAyuMapaI4KDzwjXMsKZissCExTMGZO7OC+ZdMyOige24mhYkjB7Hc7hRMFlliMAVpwTWT+cn10yjK0oh7fT1NMAAP3/gFSwomq2ysJFOwovNI/6+nN7CjJbN9LsOaevGSAGeNHMSE4adxl419ZDxgicEUtBmVpcy75Czue35jxs+VavHSFVNGseT2qozGYkwylhhMwYv2HH55027KTith5YZdnsbz5bN7GkDAmMyyxGAMkeQQTRBrtjez/2h7D3tkTuhYm2fnNgas8tmYU9wyY5xn5y4JCrMmlHl2fmPA7hiMOcWCayaz59CJrBYpTRw5iAvPHMaNNv6R8QFLDMbE8cht53HhmWVZGV+ppChgTVKNr1hiMCaB0LE2d3uqxXHllFHWJNX4jiUGYxKYNaGM4qDQ1pmZ7GDNUo1fWWIwJoHoLHDPrm/is8OtfNh0gD2HWlPeP1mHtpKiAPNtOk7jU5YYjEliRmXpyWKeRS9tZfGa+pT3LT2thP1HT216em75EBsQz/iaNVc1JkWDBxT3avt4ScFGSTW5wO4YjEnRrAllFAWkT1OECnDB+FLuvXqyJQXje5YYjElRdIrQB17YRGdYe2ywJEBxUYCbZpRb/wSTU9IqShKRm0Vks4iERSRh8woRuUpEtolInYgsiFl+poisFZHtIrJCRErSiceYTJs7s4IVd13EnJkVlBQFks7pcMWUUSz7zix+ev00Swomp6Rbx7AJuAFYk2gDEQkCjwJXA1OAOSIyxVn9MPBzVZ0EhIBvpxmPMRk3o7KUn14/jWXfmcWcmRUUB+Onh+GD+1lCMDkprcSgqltVdVsPm10I1Klqvaq2AcuB2SIiwGXAM852vwGuSyceY7IpmiCWz7uIL00afsr67EwBZIz7stEqaSywM+Z1k7OsDDigqh3dlsclIvNEpEZEapqbmzMWrDG9NaOylO9d/jlKiv7rz6koKNxwfrmHURnTdz1WPovI68AZcVbdr6ovpHCOeF+cNMnyuFR1CbAEoKqqKsMDFRjTOzMqS1n2nVk8t74JBatsNjmtx8SgqpeneY4mIHYc43JgF/AZMFREipy7huhyY3JSbGc4Y3JZNoqS1gGTnBZIJcBtwCpVVWA1cJOz3R1AKncgxhhjMijd5qrXi0gTcBHw7yLyirN8jIi8BODcDdwDvAJsBZ5W1c3OIe4F/kZE6ojUOfwynXiMMcakTzTTg81nQFVVldbU1HgdhjHG5BQRqVXVHof0tbGSjDHGdGGJwRhjTBeWGIwxxnSRk3UMItIMNKR5mOFEmszmOrsOf7Hr8Be7jq4qVXVETxvlZGJwg4jUpFIJ43d2Hf5i1+Evdh19Y0VJxhhjurDEYIwxpotCTgxLvA7AJXYd/mLX4S92HX1QsHUMxhhj4ivkOwZjjDFxWGIwxhjTRcEkhl7MT71DRDaKyAYR8d2ATOnOs+0XIjJMRF5z5vt+TUTijlctIp3O72KDiKzKdpyJ9PTzFZF+zjzmdc685uOzH2VyKVzDt0SkOebnf6cXcfZERJ4UkX0isinBehGRXzjX+aGInJ/tGFORwnV8WUQOxvw+HshYMKpaEA9gMnA28BZQlWS7HcBwr+NN5zqAIPAxMAEoAT4Apngde7cYfwYscJ4vAB5OsN0Rr2Pty88X+GtgsfP8NmCF13H34Rq+Bfyz17GmcC2XAOcDmxKsvwZ4mcjkYLOAtV7H3Mfr+DLwYjZiKZg7Bk1tfmrfS/E64s6znfnoemU2kXm+Iffm+07l5xt7fc8AX3XmOfeLXHiPpERV1wD7k2wyG/itRlQTmSBsdHaiS10K15E1BZMYekGBV0WkVkTmeR1MHyWaZ9tPRqnqbgDn/5EJtuvvzPVdLSJ+SR6p/HxPbqOROUkOEplzxC9SfY/c6BS/PCMi4+KszwW58PeQqotE5AMReVlEzsnUSXqc2jOXuDA/NcAXVXWXiIwEXhORj5xMnjUZnGc7q5JdRy8OU+H8PiYAb4rIRlX92J0I+yyVn68vfgdJpBLfH4BlqtoqIvOJ3AFdlvHI3Of330Wq1hMZ6+iIiFwDrAQmZeJEeZUYNP35qVHVXc7/+0TkeSK33FlNDC5cR6J5trMq2XWIyF4RGa2qu53b+n0JjhH9fdSLyFvAeUTKxr2Uys83uk2TiBQBQ/BJMYGjx2tQ1ZaYl48DD2chrkzwxd9DulT1UMzzl0TkX0RkuKq6PkigFSXFEJHTRGRw9DlwJRC3hYDPxZ1n2+OYultFZJ5vSDDft4iUikg/5/lw4IvAlqxFmFgqP9/Y67sJeFOdGkSf6PEaupXDX0tkat5ctAq43WmdNAs4GC3GzCUicka0nkpELiTy+d2SfK8+8romPlsP4Hoi3xxagb3AK87yMcBLzvMJRFpnfABsJlJ043nsvb0O5/U1wJ+IfLv243WUAW8A253/hznLq4AnnOcXAxud38dG4Ntex53s5wssBK51nvcHfg/UAe8DE7yOuQ/X8A/O38EHwGrg817HnOA6lgG7gXbnb+PbwHxgvrNegEed69xIklaJPr+Oe2J+H9XAxZmKxYbEMMYY04UVJRljjOnCEoMxxpguLDEYY4zpwhKDMcaYLiwxGGOM6cISgzHGmC4sMRhjjOni/wPOuZTnvJ4jxwAAAABJRU5ErkJggg==\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",
    "import numpy as np\n",
    "p = 1 # Which norm do we use\n",
    "M = 40000 # Number of sampling points\n",
    "a = np.random.randn(M, 2)\n",
    "b = []\n",
    "for i in range(M):\n",
    "    if np.linalg.norm(a[i, :], p) <= 1:\n",
    "        b.append(a[i, :])\n",
    "b = np.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": [
    "## Почему $L_1$ норма может быть важна?\n",
    "\n",
    "$L_1$ норма играет важную роль в задаче **compressed sensing**. \n",
    "\n",
    "Простейшая формулировка этой задачи следующая:\n",
    "\n",
    "- Даны некоторые наблюдения $f$ \n",
    "- Известно, что модель получения наблюдений линейная $Ax = f$, где $A$ – это $n \\times m$ известная матрица\n",
    "- Число уравнений $n$ меньше, чем число неизвестных $m$\n",
    "\n",
    "**Q**: можем ли мы найти решение такой системы?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Решение, очевидно, не единственно, поэтому естественный подход – это искать решение минимальное в некотором смысле:\n",
    "\n",
    "\\begin{align*}\n",
    "& \\Vert x \\Vert \\rightarrow \\min_x \\\\\n",
    "\\mbox{subject to } & Ax = f\n",
    "\\end{align*}\n",
    "\n",
    "- Выбор стандартной евклидовой нормы $\\Vert x \\Vert = \\Vert x \\Vert_2$ приводит к **линейной задаче наименьших квадратов**  \n",
    "\n",
    "- Выбор первой нормы $\\Vert x \\Vert = \\Vert x \\Vert_1$ приводит к задаче [**compressed sensing**](https://en.wikipedia.org/wiki/Compressed_sensing)\n",
    "- Обычно решение этой задачи является наиболее **разреженным** (sparse) решением данной системы "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Что такое устойчивый алгоритм?\n",
    "\n",
    "Введём понятия **устойчивости** алгоритма.\n",
    "\n",
    "- Пусть $x$ некоторый объект, с которым производят операции (например, вектор)\n",
    "- Пусть $f(x)$ функция от этого объекта, которую необходимо вычислить \n",
    "\n",
    "Также имеется **вычислительный алгоритм** ``alg(x)``, который вычисляет некоторую аппроксимацию функции $f(x)$.  \n",
    "\n",
    "Будем называть алгоритм **forward stable**, если $$\\Vert alg(x) - f(x) \\Vert  \\leq \\varepsilon $$  \n",
    "\n",
    "Будем называть алгоритм **backward stable**, если для любого $x$ найдётся близкий вектор $x + \\delta x$ такой что\n",
    "\n",
    "$$alg(x) = f(x + \\delta x)$$\n",
    "\n",
    "и $\\Vert \\delta x \\Vert$ мало."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Классический пример\n",
    "\n",
    "Классическим примером для демонстрации (не)устойчивости алгоритма является задача решения системы линейных уравнений с помощью метода Гаусса, который связан с вычислением LU разложения матрицы (подробности будут далее в курсе).\n",
    "\n",
    "Рассмотрим **матрицу Гильберта** с элементами\n",
    "\n",
    "$$A = \\{a_{ij}\\}, \\quad a_{ij} = \\frac{1}{i + j + 1}, \\quad i,j = 0, \\ldots, n-1.$$\n",
    "\n",
    "и соответствующую линейную систему с этой матрицей\n",
    "\n",
    "$$Ax = f.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.860287202800179\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 500\n",
    "a = [[1.0/(i + j + 1) for i in range(n)] for j in range(n)] # Hilbert matrix\n",
    "A = np.array(a)\n",
    "rhs =  np.random.random(n)\n",
    "sol = np.linalg.solve(A, rhs)\n",
    "print(np.linalg.norm(A.dot(sol) - rhs)/np.linalg.norm(rhs))\n",
    "#plt.plot(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.287160152440011e-08\n"
     ]
    }
   ],
   "source": [
    "rhs =  np.ones(n)\n",
    "sol = np.linalg.solve(A, rhs)\n",
    "print(np.linalg.norm(A.dot(sol) - rhs)/np.linalg.norm(rhs))\n",
    "#plt.plot(sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Другие примеры неустойчивости\n",
    "\n",
    "Какой вычислительно устойчивый способ вычисления следующих функций?\n",
    "\n",
    "- $\\log(1 - \\tanh^2(x))$\n",
    "- $SoftMax(x)_j = \\dfrac{e^{x_j}}{\\sum\\limits_{i=1}^n e^{x_i}}$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original function: -inf\n",
      "Attempt imporove stability with add small constant: -13.815510557964274\n",
      "Use more numerically stable form: -598.6137056388801\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    }
   ],
   "source": [
    "u = 300\n",
    "eps = 1e-6\n",
    "print(\"Original function:\", np.log(1 - np.tanh(u)**2))\n",
    "eps_add = np.log(1 - np.tanh(u)**2 + eps)\n",
    "print(\"Attempt imporove stability with add small constant:\", eps_add)\n",
    "print(\"Use more numerically stable form:\", np.log(4) - 2 * np.log(np.exp(-u) + np.exp(u)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[nan  0.  0.  0.  0.]\n",
      "[1. 0. 0. 0. 0.]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/alex/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in true_divide\n",
      "  after removing the cwd from sys.path.\n"
     ]
    }
   ],
   "source": [
    "n = 5\n",
    "x = np.random.randn(n)\n",
    "x[0] = 1000\n",
    "print(np.exp(x) / np.sum(np.exp(x)))\n",
    "print(np.exp(x - np.max(x)) / np.sum(np.exp(x - np.max(x))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Краткое резюме\n",
    "\n",
    "- Форматы представления действительных чисел и ошибка округления\n",
    "- Векторные нормы: $1$, $p$ и евклидова норма \n",
    "- $L_1$ норма как спопособ получить разреженное решение, задача compressed sensing\n",
    "- Прямая (forward) и обратная (backward) устойчивость алгоритмов"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"../styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python (pytorch)",
   "language": "python",
   "name": "pytorch"
  },
  "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.6.10"
  },
  "nav_menu": {},
  "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
}