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",
    "- Как измерять точность: векторные нормы\n",
    "- Матричные нормы\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"
   ]
  },
  {
   "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": [
    "## Альтернатива стандарту 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": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error in np sum: 1.9e-04\n",
      "Error in Kahan sum: -1.3e-07\n",
      "Error in dumb sum: -1.0e-02\n",
      "Error in math fsum: 1.3e-10\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from numba import jit\n",
    "\n",
    "n = 10 ** 8\n",
    "sm = 1e-10\n",
    "x = np.ones(n, dtype=np.float32) * sm\n",
    "x[0] = 1.0\n",
    "true_sum = 1.0 + (n - 1)*sm\n",
    "approx_sum = np.sum(x)\n",
    "math_fsum = math.fsum(x)\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def dumb_sum(x):\n",
    "    s = np.float32(0.0)\n",
    "    for i in range(len(x)):\n",
    "        s = s + x[i]\n",
    "    return s\n",
    "\n",
    "@jit(nopython=True)\n",
    "def kahan_sum(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",
    "k_sum = kahan_sum(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: {0:3.1e}'.format(k_sum - 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: 1%\">\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": [
    "## Матрицы и матричные нормы"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- Векторные нормы позволяли измерить размер вектора и расстояние между векторами\n",
    "\n",
    "- Как обобщить эту концепцию на матрицы?\n",
    "\n",
    "- Если считать матрицу вектором, элементы которого переставлены в другом порядке, то мы по аналогии с евклидовой нормой вектора получим фробениусову норму матрицы:\n",
    "\n",
    "$$\n",
    "  \\Vert A \\Vert_F \\stackrel{\\mathrm{def}}{=} \\Big(\\sum_{i=1}^n \\sum_{j=1}^m |a_{ij}|^2\\Big)^{1/2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матричные нормы\n",
    "\n",
    "$\\Vert \\cdot \\Vert$ называется матричной нормой, если это векторная норма на векторном пространстве матриц размера $n \\times m$:\n",
    "1. $\\|A\\| \\geq 0$ и если $\\|A\\| = 0$ то $A = O$\n",
    "3. $\\|\\alpha A\\| = |\\alpha| \\|A\\|$\n",
    "4. $\\|A+B\\| \\leq \\|A\\| + \\|B\\|$ (неравенство треугольника)\n",
    "\n",
    "Дополнительно некоторые нормы обладают свойством субмультипликативности\n",
    "\n",
    "* <font color='red'> $\\Vert A B \\Vert \\leq \\Vert A \\Vert \\Vert B \\Vert$ </font>\n",
    "\n",
    "Эти нормы также называются **субмультипликативными нормами**.\n",
    "\n",
    "Это свойство необходимо во многих случаях, например для оценки ошибки решения линейной системы (позднее обсудим эту тему). \n",
    "\n",
    "Примером не субмультипликативной нормы является чебышёвская норма\n",
    "\n",
    "$$\n",
    "\\|A\\|_C = \\displaystyle{\\max_{i,j}}\\, |a_{ij}|\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Операторные нормы\n",
    "\n",
    "Наиболее важным классом матричных норм являются **операторные нормы**, которые определяются следующим образом\n",
    "\n",
    "$$\n",
    "    \\Vert A \\Vert_{*,**} = \\sup_{x \\ne 0} \\frac{\\Vert A x \\Vert_*}{\\Vert x \\Vert_{**}},\n",
    "$$\n",
    "\n",
    "где $\\Vert \\cdot \\Vert_*$ и $\\| \\cdot \\|_{**}$ **векторные нормы**.\n",
    "\n",
    "Легко проверить, что операторные нормы являются субмультипликативными."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матричные $p$-нормы\n",
    "\n",
    "Важный класс операторных норм – матричные $p$-нормы, которые определены для $\\|\\cdot\\|_* = \\|\\cdot\\|_{**} = \\|\\cdot\\|_p$. <br>\n",
    "\n",
    "Среди всех $p$-норм три нормы наиболее часто используются:  \n",
    "\n",
    "- $p = 1, \\quad \\Vert A \\Vert_{1} = \\displaystyle{\\max_j \\sum_{i=1}^n} |a_{ij}|$.\n",
    "\n",
    "- $p = 2, \\quad$ спектральная норма, обозначается как $\\Vert A \\Vert_2$.\n",
    "\n",
    "- $p = \\infty, \\quad \\Vert A \\Vert_{\\infty} = \\displaystyle{\\max_i \\sum_{j=1}^m} |a_{ij}|$.\n",
    "\n",
    "Покажем справедливость выражения для $p=\\infty$ на доске."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Спектральная норма\n",
    "\n",
    "- Спектральная норма, $\\Vert A \\Vert_2$, одна из наиболее часто используемых, наряду с нормой Фробениуса\n",
    "- Она не может быть вычислена напрямую через элементы матрицы, но для её вычисления существуют эффективные алгоритмы\n",
    "- Она напрямую связана с сингулярным разложением матрицы (SVD) поскольку выполнено следующее равенство\n",
    "\n",
    "$$\n",
    "   \\Vert A \\Vert_2 = \\sigma_1(A) = \\sqrt{\\lambda_\\max(A^*A)}\n",
    "$$\n",
    "\n",
    "где $\\sigma_1(A)$ – максимальное сингулярное число матрицы $A$. \n",
    "\n",
    "- Подробнее про SVD будет рассказано далее, а сейчас давайте вычислим все вышеупомянутые нормы"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Spectral: 55.06036768932167 \n",
      "Frobenius: 446.9316116203868 \n",
      "1-norm: 102.20049588176245 \n",
      "infinity: 1649.7258667392314\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "n = 100\n",
    "m = 2000\n",
    "a = np.random.randn(n, m) # Random n x n matrix\n",
    "s1 = np.linalg.norm(a, 2) # Spectral\n",
    "s2 = np.linalg.norm(a, 'fro') # Frobenius\n",
    "s3 = np.linalg.norm(a, 1) # 1-norm\n",
    "s4 = np.linalg.norm(a, np.inf) \n",
    "print('Spectral: {0:} \\nFrobenius: {1:} \\n1-norm: {2:} \\ninfinity: {3:}'.format(s1, s2, s3, s4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Примеры\n",
    "\n",
    "Несколько примеров задач оптимизации, в которых необходимы матричные нормы:\n",
    "* <font color='red'> $\\displaystyle{\\min_{\\mathrm{rank}(A_r)=r}}\\| A - A_r\\|$ </font> –– задача нахождения наилучшей малоранговой аппроксимации матрицы $A$. SVD помогает решить эту задачу для спектральной нормы $\\|\\cdot\\|_2$ и фробениусовой нормы $\\|\\cdot\\|_F$.\n",
    "\n",
    "\n",
    "* $\\displaystyle{\\min_B}\\| P_\\Omega \\odot(A - B)\\| + \\lambda \\cdot \\mathrm{rank}(B)$ –– задача заполнения матрицы (matrix completion problem). \n",
    "$$\n",
    "(P_\\Omega)_{ij} = \\begin{cases} 1 & i,j\\in\\Omega \\\\ 0 & \\text{otherwise} \\end{cases}\n",
    "$$\n",
    "$\\odot$ обозначает произведение Адамара (поэлементное)\n",
    "\n",
    "\n",
    "* $\\displaystyle{\\min_{B,C\\geq 0}} \\|A - BC\\|$ –– неотрицательное матричное разложение (nonnegative matrix factorization). Запись $B\\geq0$ означает, что все элементы $B$ неотрицательны."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Скалярное произведение\n",
    "\n",
    "Если норма помогает измерять расстояние, то **скалярное произведение** позволяет учесть угол между объектами.  \n",
    "\n",
    "Скалярное произведение определено\n",
    "\n",
    "* **для векторов:**\n",
    "$$\n",
    "   (x, y) =  x^{\\top} y = \\sum_{i=1}^n x_i y_i.\n",
    "$$\n",
    "Тогда евклидова норма записывается как\n",
    "$$\n",
    "   \\Vert x \\Vert_2 = \\sqrt{(x, x)};\n",
    "$$  \n",
    "\n",
    "\n",
    "* **для матриц**:\n",
    "$$\n",
    "    (A, B)_F = \\displaystyle{\\sum_{i=1}^{n}\\sum_{j=1}^{m}} a_{ij} b_{ij} \\equiv \\mathrm{trace}(A^{\\top} B),\n",
    "$$\n",
    "где $\\mathrm{trace}(A)$ обозначает след матрицы, то есть сумму диагональных элементов. \n",
    "\n",
    "**Упражнение:** покажите, что $\\|A\\|_F = \\sqrt{(A, A)_F}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Замечание**. Угол между векторами или матрицами определяется как\n",
    "\n",
    "$$\n",
    "   \\cos \\phi = \\frac{(x, y)}{\\Vert x \\Vert_2 \\Vert y \\Vert_2}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матрицы, сохраняющие норму\n",
    "\n",
    "- Для устойчивости вычислений необходимо, чтобы ошибка не возрослала после применения некоторого преобразования. Найдём эти преобразования. \n",
    "\n",
    "- Пусть дан вектор $\\widehat{x}$ – аппроксимация вектора $x$ такая что,  \n",
    "\n",
    "$$\n",
    "  \\frac{\\Vert x - \\widehat{x} \\Vert}{\\Vert x \\Vert} \\leq \\varepsilon.\n",
    "$$\n",
    "\n",
    "- Вычислим образ векторов $x$ и $\\widehat{x}$ после применения линейного преобразования $U$:  \n",
    "\n",
    "$$\n",
    "   y = U x, \\quad \\widehat{y} = U \\widehat{x}.\n",
    "$$\n",
    "\n",
    "- Для построения алгоритмов необходимо использовать преобразования, которые не увеличивают (или даже сохраняют) ошибку:\n",
    "\n",
    "$$\n",
    "   \\frac{\\Vert y - \\widehat{y} \\Vert}{\\Vert y \\Vert } = \\frac{\\Vert U ( x - \\widehat{x}) \\Vert}{\\Vert U  x\\Vert}  \\leq \\varepsilon.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Вопрос состоит в том, какой класс матриц не меняет норму вектора после умножения матрицы из этого класса на вектор?\n",
    "\n",
    "$$\n",
    "\\frac{\\Vert U ( x - \\widehat{x}) \\Vert}{\\Vert U  x\\Vert} = \\frac{ \\|x - \\widehat{x}\\|}{\\|x\\|}.\n",
    "$$\n",
    "\n",
    "Для евклидовой нормы $\\|\\cdot\\|_2$ таким свойством обладают **унитарные** (или ортогональные) матрицы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Унитарные (ортогональные) матрицы\n",
    "\n",
    "Пусть $U$ комплексная $n \\times n$ матрица, и $\\Vert U z \\Vert_2 = \\Vert z \\Vert_2$ для вех $z$. \n",
    "\n",
    "Это выполнено тогда и только тогда, когда\n",
    "\n",
    "$$\n",
    "   U^* U = I_n,\n",
    "$$\n",
    "\n",
    "где $I_n$ единичная матрица $n\\times n$.\n",
    "\n",
    "Комплексная $n\\times n$ квадратная матрица называется **унитарной** если\n",
    "\n",
    "$$\n",
    "    U^*U = UU^* = I_n,\n",
    "$$\n",
    "\n",
    "что означает, что столбцы и строки матрицы образуют базис в $\\mathbb{C}^{n}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Свойство ортогональных матриц\n",
    "\n",
    "Произведение двух унитарных матриц – унитарная матрица:  \n",
    "\n",
    "$$(UV)^* UV = V^* (U^* U) V = V^* V = I,$$\n",
    "\n",
    "- Позже мы покажем, что существуют классы унитарных матриц, произведение которых может дать произвольную унитарную матрицу \n",
    "- Эта идея лежит в основе некоторых алгоритмов, например вычисления QR разложения"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Примеры  унитарных матриц\n",
    "\n",
    "Два важных класса унитарных матриц, произведение которых может дать любую унитарную матрицу:\n",
    "1. Матрицы Хаусхолдера\n",
    "2. Матрицы Гивенса\n",
    "\n",
    "Этот факт станет очевидным после того, как мы рассмотрим QR разложение и способы его вычисления.\n",
    "\n",
    "Другие важные примеры\n",
    "* **матрица перестановки** $P$ строки (столбцы) которой получены перестановкой строк (столбцов) единичной матрицы.\n",
    "* **матрица Фурье** $F_n = \\frac{1}{\\sqrt{n}} \\left\\{ e^{-i\\frac{2\\pi kl}{n}}\\right\\}_{k,l=0}^{n-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матрица Хаусхолдера\n",
    "\n",
    "Матрица Хаусхолдера – это матрица вида \n",
    "\n",
    "$$H \\equiv H(v) = I - 2 vv^*,$$\n",
    "\n",
    "где $v$ – $n \\times 1$ вектор и $v^* v = 1$. Покажите, что $H$ унитарна и Эрмитова. Это также матрица отражения относительно плоскости с нормалью $v$:\n",
    "\n",
    "$$ Hx = x - 2(v^* x) v$$\n",
    "\n",
    "<img src=\"./householder.jpeg\">  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Важное свойство матрицы Хаусхолдера\n",
    "\n",
    "Преобразование Хаусхолдера может занулить все элементы в столбце матрицы, кроме первого:\n",
    "\n",
    "$$\n",
    "    H \\begin{bmatrix} \\times \\\\ \\times \\\\ \\times \\\\ \\times  \\end{bmatrix} = \n",
    "      \\begin{bmatrix} \\times \\\\ 0 \\\\ 0 \\\\ 0  \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Доказательство на семинаре."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Матрица Гивенса\n",
    "\n",
    "Матрица Гивенса имеет вид\n",
    "\n",
    "$$\n",
    "    G = \\begin{bmatrix}\n",
    "          \\cos \\alpha & -\\sin \\alpha \\\\\n",
    "          \\sin \\alpha & \\cos \\alpha\n",
    "        \\end{bmatrix},\n",
    "$$\n",
    "\n",
    "что является поворотом на угол $\\alpha$ на плоскости. В общем случае, выберем две плоскости $(i, j)$ и повернём вектор $x$  \n",
    "\n",
    "$$\n",
    "x' = G x,\n",
    "$$\n",
    "\n",
    "только в $i$-ой и $j$-ой координатах:\n",
    "\n",
    "$$\n",
    "    x'_i =  x_i\\cos \\alpha - x_j\\sin \\alpha , \\quad x'_j = x_i \\sin \\alpha  +  x_j\\cos\\alpha,\n",
    "$$\n",
    "\n",
    "и оставим остальные координаты $x_i$ неизменными.\n",
    "\n",
    "Поэтому мы можем занулить элемент в позиции $j$, если выберем $\\alpha$ так что\n",
    "$$\n",
    "     \\cos \\alpha = \\frac{x_i}{\\sqrt{x_i^2 + x_j^2}}, \\quad \\sin \\alpha = -\\frac{x_j}{\\sqrt{x_i^2 + x_j^2}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0, 1.0)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAD8CAYAAACLrvgBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl4VPXZ//H3za6oNcgqIkJBJewygqCoRUAEEawUQ11QUUrripUKD48bbpRHBbeqUXChFhd6qah4IasrKFEJhFAkYltTqKAoLiDr/ftjDvzGMElOmJlMEj6v65przvl+v+ecO4chn5wzZ86YuyMiIhJGtXQXICIilYdCQ0REQlNoiIhIaAoNEREJTaEhIiKhKTRERCS0pISGmU0zsw1mlldMv5nZA2ZWYGbLzeyEmL7hZrYmeAxPRj0iIpIayTrSeAroV0L/WUDr4DESeATAzOoBtwDdgK7ALWaWkaSaREQkyZISGu7+NrCphCGDgGc8aglwuJk1Ac4E5rr7Jnf/BphLyeEjIiJpVKOcttMU+CJmvjBoK659H2Y2kuhRCnXr1u1y/PHHp6ZSEZEq6qOPPvrK3Rskso7yCg2L0+YltO/b6J4NZANEIhHPyclJXnUiIgcAM/tXousor6unCoFmMfNHAetKaBcRkQqovEJjFnBxcBXVScBmd18PzAH6mllG8AZ436BNREQqoKScnjKzGcDpQH0zKyR6RVRNAHd/FJgN9AcKgC3ApUHfJjO7HVgarGqCu5f0hrqIiKRRUkLD3YeV0u/AlcX0TQOmJaMOERFJLX0iXEREQlNoiIhIaAoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEppCQ0REQlNoiIhIaAoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEppCQ0REQlNoiIhIaEkJDTPrZ2arzazAzMbG6Z9sZsuCx6dm9m1M366YvlnJqEdERFIj4e8IN7PqwMNAH6AQWGpms9w9f88Ydx8dM/5qoHPMKra6e6dE6xARkdRLxpFGV6DA3de6+3bgOWBQCeOHATOSsF0RESlnyQiNpsAXMfOFQds+zKw50AJYENNcx8xyzGyJmQ1OQj0iIpIiCZ+eAixOmxczNguY6e67YtqOdvd1ZtYSWGBmK9z9s302YjYSGAlw9NFHJ1qziIjsh2QcaRQCzWLmjwLWFTM2iyKnptx9XfC8FljEz9/viB2X7e4Rd480aNAg0ZpFRGQ/JCM0lgKtzayFmdUiGgz7XAVlZscBGcDimLYMM6sdTNcHTgbyiy4rIiIVQ8Knp9x9p5ldBcwBqgPT3H2lmU0Actx9T4AMA55z99hTV22Ax8xsN9EAmxh71ZWIiFQs9vPf4ZVDJBLxnJycdJex33bt2sWUKVNo0aIFffr04dBDD013SSJyADCzj9w9ksg6kvFGuJRR9erVOffcc+nWrRubN2/m1FNPZcCAAQwYMIBjjz023eWJiBRLRxpp9N5779GrVy+2b9++t61Vq1YMGDCA/v37c9ppp1G7du00VigiVUkyjjQUGmn27LPPcuGFF8btq1u3Lr1792bAgAGcd9551KtXr5yrE5GqJBmhoRsWptkFF1zAzTffHLfvxx9/ZNGiRWzdulXve4hIhaDQqABuvfVWsrKy4vaNHTuWq6++mpo1a5ZzVSIi+1JoVABmxpNPPslJJ520T9+4cePo27cvn3/+eRoqExH5OYVGBVGnTh1efvllmjdvDkDHjh1p3LgxAPPmzaNdu3bcf//97Nq1q6TViIiklEKjAmnUqBGvvfYahx56KP379yc/P5/LLrsMgC1btnDdddfRs2dP8vP1+UcRSQ+FRgXTrl07XnjhBVq0aEFGRgZTp05l7ty5HHPMMQAsXryYzp07c8cdd7Bjx470FisiBxxdcltB7dix42dvfv/www/cdNNN3H///ez5N+vQoQPTpk2jS5cu6SpTRCoRXXJbhRW9WuqQQw5h8uTJvPfee7Rp0waA5cuX07VrV2688Ua2bt2ajjJF5ACj0KhkunfvzieffMJNN91EjRo12L17N5MmTaJDhw689dZb6S5PRKo4hUYlVLt2bSZMmEBOTs7eU1MFBQWcfvrp/P73v+e7775Lc4UiUlUpNCqxjh07smTJEiZNmkSdOnUAePTRR2nbti2zZ89Oc3UiUhUpNCq5GjVqMGbMGHJzczn11FMBKCwsZMCAAVx44YV89dVXaa5QRKoShUYVceyxx7Jw4UIeeeSRvfepevbZZ8nMzOT555+nMl4lJyIVj0KjCqlWrRqjRo1i5cqV9O/fH4CNGzeSlZXF4MGD+c9//pPmCkWkslNoVEHNmjXjtdde469//StHHHEEALNmzSIzM5PHH39cRx0ist8UGlWUmXHBBReQn5+/9w663333HSNHjuSMM87gs88+S3OFIlIZJSU0zKyfma02swIzGxun/xIz22hmy4LH5TF9w81sTfAYnox65P9r2LAhM2bM4JVXXuHII48EYOHChbRv35777rtPN0AUkTJJODTMrDrwMHAWkAkMM7PMOEOfd/dOweOJYNl6wC1AN6ArcIuZZSRak+zrnHPOYeXKlVxxxRUAbN26lT/+8Y/06NGDvLy8NFcnIpVFMo40ugIF7r7W3bcDzwGDQi57JjDX3Te5+zfAXKBfEmqSOA4//HCys7OZP38+LVu2BODDDz/khBNO4LbbbvvZd5WLiMSTjNBoCnwRM18YtBV1npktN7OZZtasjMtiZiPNLMfMcjZu3JiEsg9cvXr1YsWKFVx//fVUq1aNHTt2cOutt9KlSxc+/PDDdJcnIhVYMkLD4rQVvTznVeAYd+8AzAOeLsOy0Ub3bHePuHukQYMG+12sRB188MHce++9vP/++7Rt2xaAvLw8unfvzg033MCWLVvSXKGIVETJCI1CoFnM/FHAutgB7v61u28LZh8HuoRdVlKrW7dufPzxx9xyyy3UrFmT3bt3c++999K+fXsWLlyY7vJEpIJJRmgsBVqbWQszqwVkAbNiB5hZk5jZc4BVwfQcoK+ZZQRvgPcN2qQc1apVi1tvvZWPPvqIE088EYC1a9fSq1cvfve737F58+Y0VygiFUXCoeHuO4GriP6yXwW84O4rzWyCmZ0TDLvGzFaaWS5wDXBJsOwm4HaiwbMUmBC0SRq0b9+exYsXc88993DQQQcBkJ2dTWZmJq+++mqaqxORikDf3CdxFRQUcMUVV7Bo0aK9bcOGDeP+++9H7ymJVE765j5JmVatWrFgwQKys7M57LDDAJgxYwZt2rThb3/7m25FInKAUmhIscyMK664gvz8fAYOHAjA119/zQUXXMDAgQP54osvSlmDiFQ1Cg0pVdOmTXnllVeYMWMG9evXB+D111+nbdu2PPbYY+zevTvNFYpIeVFoSChmRlZWFqtWreK3v/0tAN9//z2jRo2iV69erFmzJs0Vikh5UGhImdSvX59nn32WV199laZNox/ef+utt+jQoQP33HMPO3fuTHOFIpJKCg3ZL2effTb5+fmMGjUKgJ9++okxY8bQvXt3li9fnubqRCRVFBqy3w477DAeeeQRFi1aRKtWrQDIycmhS5cu3HzzzWzbtq2UNYhIZaPQkISddtpp5ObmMmbMGKpVq8bOnTu5/fbbOeGEE1iyZEm6yxORJFJoSFIcfPDBTJo0iQ8++ID27dsDkJ+fT48ePRg9ejQ//vhjmisUkWRQaEhSRSIRcnJymDBhAjVr1sTdmTJlCu3bt2f+/PnpLk9EEqTQkKSrVasWN910E5988gknnXQSAJ9//jm9e/fm8ssv59tvv01zhSKyvxQakjJt27bl3XffZcqUKRx88MEATJ06lczMTF5++eU0Vyci+0OhISlVvXp1rr32WlasWMEZZ5wBwPr16zn33HMZOnQoX375ZZorFJGyUGhIuWjZsiVz585l6tSp/OIXvwDgxRdfJDMzk+nTp+sGiCKVhEJDyo2Zcdlll5Gfn8+gQYMA2LRpExdffDEDBgzg3//+d5orFJHSKDSk3B155JG89NJLvPDCCzRs2BCAN954g7Zt2/KXv/xFN0AUqcAUGpIWZsZvfvMb8vPzueiiiwD44YcfuPLKKzn99NNZvXp1misUkXgUGpJWRxxxBM888wyzZ8+mWbNmALzzzjt07NiRiRMn6gaIIhVMUkLDzPqZ2WozKzCzsXH6rzezfDNbbmbzzax5TN8uM1sWPGYlox6pfM466yxWrlzJlVdeCcC2bdsYN24c3bp1Y9myZWmuTkT2SDg0zKw68DBwFpAJDDOzzCLDPgEi7t4BmAlMiunb6u6dgsc5idYjldehhx7KQw89xNtvv03r1q0B+Pjjj4lEIowfP56ffvopzRWKSDKONLoCBe6+1t23A88Bg2IHuPtCd98SzC4BjkrCdqWK6tmzJ7m5uYwdO5bq1auza9cu7rrrLjp37sz777+f7vJEDmjJCI2mQOyXRRcGbcUZAbwRM1/HzHLMbImZDS5uITMbGYzL2bhxY2IVS4V30EEHcffdd/PBBx/QsWNHAP7xj39wyimncM011/DDDz+kuUKRA1MyQsPitMX9pJaZXQhEgP+LaT7a3SPAb4EpZvbLeMu6e7a7R9w90qBBg0RrlkqiS5cuLF26lDvvvJNatWrh7jz44IO0a9eON998M93liRxwkhEahUCzmPmjgHVFB5lZb2A8cI677/12HndfFzyvBRYBnZNQk1QhNWvW5H/+53/Izc2lR48eAPzrX//izDPP5NJLL2XTpk1prlDkwJGM0FgKtDazFmZWC8gCfnYVlJl1Bh4jGhgbYtozzKx2MF0fOBnIT0JNUgUdf/zxvPPOOzzwwAPUrVsXgKeeeorMzEz+/ve/p7k6kQNDwqHh7juBq4A5wCrgBXdfaWYTzGzP1VD/BxwCvFjk0to2QI6Z5QILgYnurtCQYlWrVo2rr76avLw8+vTpA8CXX37JkCFDGDJkCP/973/TXKFI1WaV8UZxkUjEc3Jy0l2GpJm78/TTTzN69Oi939GRkZHBfffdx/DhwzGL93abyIHLzD4K3kPeb/pEuFRaZsYll1zCqlWrOO+88wD45ptvuPTSS+nXrx///Oc/01ugSBWk0JBKr3HjxsycOZOZM2fSqFEjAN58803atWvHgw8+qBsgiiSRQkOqjPPOO4/8/HwuueQSAH788UeuueYaevbsyapVq9JbnEgVodCQKqVevXo8+eSTzJkzh+bNo7c4e//99+nUqRN33XUXO3bsSHOFIpWbQkOqpL59+5KXl8fVV1+NmbF9+3bGjx9P165d+fjjj9NdnkilpdCQKuuQQw7hgQce4J133uG4444DYNmyZXTt2pVx48axdevWNFcoUvkoNKTKO/nkk1m2bBnjx4/fewPEiRMn0qlTJ9599910lydSqSg05IBQp04d7rjjDnJycujcOXqnmk8//ZSePXty1VVX8f3336e5QpHKQaEhB5ROnTrx4YcfMnHiRGrXrg3Aww8/TNu2bXnjjTdKWVpEFBpywKlRowY33ngjubm5nHLKKQB88cUX9O/fn4svvpivv/46zRWKVFwKDTlgHXfccbz11ls8/PDDHHLIIQBMnz6dzMxMXnzxRSrjLXZEUk2hIQe0atWq8Yc//IG8vDz69esHwIYNGxg6dCi//vWvWb9+fZorFKlYFBoiQPPmzZk9ezbPPPMM9erVA+Dll1+mTZs2TJs2TUcdIgGFhkjAzLjooovIz89n6NChAGzevJkRI0bQp08f1q5dm+YKRdJPoSFSRKNGjXj++ed56aWXaNy4MQDz58+nffv2TJkyhV27dqW5QpH0UWiIFGPw4MHk5+czYsQIALZs2cLo0aM55ZRTyM/Xd4XJgUmhIVKCjIwMnnjiCebOncsxxxwDwJIlS+jcuTO3334727dvT2+BIuVMoSESQu/evcnLy+O6667bewPEm2++mRNPPBF9i6Skw+OPP878+fPL/Q+XpISGmfUzs9VmVmBmY+P01zaz54P+D8zsmJi+cUH7ajM7Mxn1iKRC3bp1mTx5Mu+99x6ZmZkALF++nG7duvGnP/2JLVu2pLlCOZAcdthh9O7dmwYNGnD++eczffr0cvlgasLfEW5m1YFPgT5AIbAUGObu+TFj/gB0cPdRZpYFnOvu55tZJjAD6AocCcwDjnX3Et9p1HeES7pt27aNO++8k7vvvpudO3cC0KpVK2677Tbat2+f5urkQLBr1y769+//s88SVatWjR49ejBw4EAGDhzI8ccfj5nt7U/Gd4Tj7gk9gO7AnJj5ccC4ImPmAN2D6RrAV4AVHRs7rqRHly5dXKQiyM3N9S5dujjggDds2HDvtB56pPvRsmVLv/baa33evHm+bds2B3I8wd/5yTg91RT4Ima+MGiLO8bddwKbgSNCLguAmY00sxwzy9m4cWMSyhZJXIcOHViyZAmTJk3il7/85d5vCxSpCAoLC1m9ejWrV69O2qmrGklYh8Vp85BjwiwbbXTPBrIhenqqLAWKpFKNGjUYM2YM1113HW+//bZueCjlwt0ZPXr0Pre6adiwIQMGDGDgwIH06dNn733VkiUZoVEINIuZPwpYV8yYQjOrAfwC2BRyWZFKoWbNmpxxxhnpLkMOEAsWLNgbGB06dNj7PsaJJ55ItWqpuzA2GaGxFGhtZi2A/wBZwG+LjJkFDAcWA0OABe7uZjYL+JuZ3Uf0jfDWwIdJqElEpEpbvHgxDz30EGeffXa5nhZNODTcfaeZXUX0TezqwDR3X2lmE4i+6TILmApMN7MCokcYWcGyK83sBSAf2Alc6aVcOSUiIjB+/Pi0bDfhS27TQZfcioiUXTIuudUnwkVEJDSFhoiIhKbQEBGR0BQaIiISmkJDRERCU2iIiEhoCg0REQlNoSEiIqEpNEREJDSFhoiIhKbQEBGR0BQaIiISmkJDRERCU2iIiEhoCg0REQlNoSEiIqEpNEREJDSFhoiIhJZQaJhZPTOba2ZrgueMOGM6mdliM1tpZsvN7PyYvqfM7HMzWxY8OiVSj4iIpFaiRxpjgfnu3hqYH8wXtQW42N3bAv2AKWZ2eEz/GHfvFDyWJViPiIikUKKhMQh4Oph+GhhcdIC7f+rua4LpdcAGoEGC2xURkTRINDQauft6gOC5YUmDzawrUAv4LKb5zuC01WQzq13CsiPNLMfMcjZu3Jhg2SIisj9KDQ0zm2dmeXEeg8qyITNrAkwHLnX33UHzOOB44ESgHnBjccu7e7a7R9w90qCBDlRERNKhRmkD3L13cX1m9qWZNXH39UEobChm3GHA68D/uvuSmHWvDya3mdmTwA1lql5ERMpVoqenZgHDg+nhwCtFB5hZLeAl4Bl3f7FIX5Pg2Yi+H5KXYD0iIpJCiYbGRKCPma0B+gTzmFnEzJ4IxgwFTgUuiXNp7bNmtgJYAdQH7kiwHhERSSFz93TXUGaRSMRzcnLSXYaISKViZh+5eySRdegT4SIiEppCQ0REQlNoiIhIaAoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEppCQ0REQlNoiIhIaAoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEppCQ0REQksoNMysnpnNNbM1wXNGMeN2xXw/+KyY9hZm9kGw/PNmViuRekREJLUSPdIYC8x399bA/GA+nq3u3il4nBPT/mdgcrD8N8CIBOsREZEUSjQ0BgFPB9NPA4PDLmhmBvQCZu7P8iIiUv4SDY1G7r4eIHhuWMy4OmaWY2ZLzGxPMBwBfOvuO4P5QqBpcRsys5HBOnI2btyYYNkiIrI/apQ2wMzmAY3jdI0vw3aOdvd1ZtYSWGBmK4Dv4ozz4lbg7tlANkAkEil2nIiIpE6poeHuvYvrM7MvzayJu683sybAhmLWsS54Xmtmi4DOwN+Bw82sRnC0cRSwbj9+BhERKSeJnp6aBQwPpocDrxQdYGYZZlY7mK4PnAzku7sDC4EhJS0vIiIVR6KhMRHoY2ZrgD7BPGYWMbMngjFtgBwzyyUaEhPdPT/ouxG43swKiL7HMTXBekREJIUs+gd/5RKJRDwnJyfdZYiIVCpm9pG7RxJZhz4RLiIioSk0REQkNIWGiIiEptAQEZHQFBoiIhKaQkNEREJTaIiISGgKDRERCU2hISIioSk0REQkNIWGiIiEptAQEZHQFBoiIhKaQkNEREJTaIiISGgKDRERCU2hISIioSk0REQktIRCw8zqmdlcM1sTPGfEGfMrM1sW8/jJzAYHfU+Z2ecxfZ0SqUdERFIr0SONscB8d28NzA/mf8bdF7p7J3fvBPQCtgBvxgwZs6ff3ZclWI+IiKRQoqExCHg6mH4aGFzK+CHAG+6+JcHtiohIGiQaGo3cfT1A8NywlPFZwIwibXea2XIzm2xmtROsR0REUqhGaQPMbB7QOE7X+LJsyMyaAO2BOTHN44D/ArWAbOBGYEIxy48ERgIcffTRZdm0iIgkSamh4e69i+szsy/NrIm7rw9CYUMJqxoKvOTuO2LWvT6Y3GZmTwI3lFBHNtFgIRKJeGl1i4hI8iV6emoWMDyYHg68UsLYYRQ5NRUEDWZmRN8PyUuwHhERSaFEQ2Mi0MfM1gB9gnnMLGJmT+wZZGbHAM2At4os/6yZrQBWAPWBOxKsR0REUqjU01MlcfevgTPitOcAl8fM/xNoGmdcr0S2LyIi5UufCBcRkdAUGiIiEppCQ0REQlNoiIhIaAoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEppCQ0REQlNoiIhIaAoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEppCQ0REQksoNMzsN2a20sx2m1mkhHH9zGy1mRWY2diY9hZm9oGZrTGz582sViL1iIhIaiV6pJEH/Bp4u7gBZlYdeBg4C8gEhplZZtD9Z2Cyu7cGvgFGJFiPiIikUEKh4e6r3H11KcO6AgXuvtbdtwPPAYPMzIBewMxg3NPA4ETqERGR1KpRDttoCnwRM18IdAOOAL51950x7U2LW4mZjQRGBrPbzCwvBbUmW33gq3QXEUJlqLMy1AiqM9lUZ3Idl+gKSg0NM5sHNI7TNd7dXwmxDYvT5iW0x+Xu2UB2UFOOuxf7HkpFoTqTpzLUCKoz2VRncplZTqLrKDU03L13gtsoBJrFzB8FrCOayoebWY3gaGNPu4iIVFDlccntUqB1cKVULSALmOXuDiwEhgTjhgNhjlxERCRNEr3k9lwzKwS6A6+b2Zyg/Ugzmw0QHEVcBcwBVgEvuPvKYBU3AtebWQHR9zimhtx0diJ1lyPVmTyVoUZQncmmOpMr4Tot+ge/iIhI6fSJcBERCU2hISIioVXY0KgMtygxs3pmNjfYxlwzy4gz5ldmtizm8ZOZDQ76njKzz2P6OiW7xrB1BuN2xdQyK6a9XG73EnJ/djKzxcFrY7mZnR/Tl9L9WdxrLaa/drB/CoL9dUxM37igfbWZnZnMuvajzuvNLD/Yf/PNrHlMX9zXQBpqvMTMNsbUcnlM3/DgNbLGzIanqsaQdU6OqfFTM/s2pq9c9mWwrWlmtsGK+fyaRT0Q/BzLzeyEmL6y7U93r5APoA3RD6IsAiLFjKkOfAa0BGoBuUBm0PcCkBVMPwr8PgU1TgLGBtNjgT+XMr4esAk4OJh/ChhSDvsyVJ3AD8W0p3xfhq0TOBZoHUwfCawHDk/1/izptRYz5g/Ao8F0FvB8MJ0ZjK8NtAjWUz2Ndf4q5jX4+z11lvQaSEONlwAPxVm2HrA2eM4IpjPSVWeR8VcD08pzX8Zs61TgBCCvmP7+wBtEPx93EvDB/u7PCnuk4ZXjFiWDgnWH3cYQ4A1335KCWkpS1jr3Ksd9CSHqdPdP3X1NML0O2AA0SFE9seK+1oqMia1/JnBGsP8GAc+5+zZ3/xwoCNaXljrdfWHMa3AJ0c9Ilacw+7I4ZwJz3X2Tu38DzAX6VZA6hwEzUlRLidz9baJ/kBZnEPCMRy0h+hm5JuzH/qywoRFSvFuUNKWMtyhJQCN3Xw8QPDcsZXwW+76o7gwOFyebWe0U1Ajh66xjZjlmtmTPKTTKb1+WpU4AzKwr0b8AP4tpTtX+LO61FndMsL82E91/YZYtzzpjjSD6F+ge8V4DyRa2xvOCf8uZZrbnA8IVcl8Gp/haAAtimstjX4ZV3M9S5v1ZHveeKpZVkFuUlLiBEmos43qaAO2Jfl5lj3HAf4n+4ssm+rmVCWms82h3X2dmLYEFZrYC+C7OuP2+TjvJ+3M6MNzddwfNSduf8TYZp63ofkj56zGE0NsyswuBCHBaTPM+rwF3/yze8imu8VVghrtvM7NRRI/geoVcNlnKsq0sYKa774ppK499GVbSXptpDQ2vBLcoKalGM/vSzJq4+/rgl9iGElY1FHjJ3XfErHt9MLnNzJ4EbtifGpNVZ3C6B3dfa2aLgM7A30ni7V6SUaeZHQa8DvxvcKi9Z91J259xFPdaizem0MxqAL8gesogzLLlWSdm1ptoUJ/m7tv2tBfzGkj2L7pSa3T3r2NmHyf6NQp7lj29yLKLklzfHmX5d8sCroxtKKd9GVZxP0uZ92dlPz2V7luUzArWHWYb+5zvDH4x7nnfYDDR7ydJhVLrNLOMPadzzKw+cDKQX477MmydtYCXiJ6ffbFIXyr3Z9zXWgn1DwEWBPtvFpBl0aurWgCtgQ+TWFuZ6jSzzsBjwDnuviGmPe5rIE01NomZPYfo3SQgeqTeN6g1A+jLz4/ey7XOoNbjiL6JvDimrbz2ZVizgIuDq6hOAjYHf2SVfX+W17v7ZX0A5xJNwW3Al8CcoP1IYHbMuP7Ap0QTfHxMe0ui/zELgBeB2imo8QhgPrAmeK4XtEeAJ2LGHQP8B6hWZPkFwAqiv9z+ChySon1Zap1Aj6CW3OB5RHnuyzLUeSGwA1gW8+hUHvsz3muN6Omvc4LpOsH+KQj2V8uYZccHy60Gzkrx/53S6pwX/J/as/9mlfYaSEONdwMrg1oWAsfHLHtZsI8LgEvTuS+D+VuBiUWWK7d9GWxvBtErCXe4OIZfAAAATklEQVQQ/b05AhgFjAr6jeiX4X0W1BOJWbZM+1O3ERERkdAq++kpEREpRwoNEREJTaEhIiKhKTRERCQ0hYaIiISm0BARkdAUGiIiEtr/A62sBYUb1b13AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "alpha = -3*np.pi / 4\n",
    "G = np.array([\n",
    "    [np.cos(alpha), -np.sin(alpha)],\n",
    "    [np.sin(alpha), np.cos(alpha)]\n",
    "])\n",
    "x = np.array([-1./np.sqrt(2), 1./np.sqrt(2)])\n",
    "y = G.dot(x)\n",
    "\n",
    "plt.quiver([0, 0], [0, 0], [x[0], y[0]], [x[1], y[1]], angles='xy', scale_units='xy', scale=1)\n",
    "plt.xlim(-1., 1.)\n",
    "plt.ylim(-1., 1.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Сингулярное разложение (SVD)\n",
    "\n",
    "SVD расмотрим более подробно далее в курсе. Пока же приведём теорему о его существовании.\n",
    "\n",
    "**Теорема.** Любая матрица $A\\in \\mathbb{C}^{n\\times m}$ может быть записана как произведение трёх матриц:  \n",
    "\n",
    "$$\n",
    "   A = U \\Sigma V^*,\n",
    "$$\n",
    "\n",
    "где \n",
    "- $U$ унитарная матрица $n \\times n$\n",
    "- $V$ унитарная матрица $m \\times m$\n",
    "- $\\Sigma$ диагональная матрица с неотрицательными диагональными элементами $\\sigma_1 \\geq  \\ldots, \\geq \\sigma_{\\min (m,n)}$.\n",
    "\n",
    "Более того, если $\\text{rank}(A) = r$, тогда $\\sigma_{r+1} = \\dots = \\sigma_{\\min (m,n)} = 0$.\n",
    "\n",
    "Если заменить все сингулярные числа кроме первых $r$ нулями, то результат будет наилучшим приближением ранга $r$ исходной матрицы в нормах $\\|\\cdot\\|_2$ и $\\|\\cdot\\|_F$. Это утверждение называется теоремой Eckart'a-Young'a и будет доказано позже в нашем курсе."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Краткое резюме\n",
    "\n",
    "- Форматы представления действительных чисел и ошибка округления\n",
    "- Векторные нормы: $1$, $p$ и евклидова норма \n",
    "- $L_1$ норма как спопособ получить разреженное решение, задача compressed sensing\n",
    "- Прямая (forward) и обратная (backward) устойчивость алгоритмов\n",
    "- Матричные нормы\n",
    "- Унитарные матрицы: примеры и свойства"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Анонс на следующее занятие (уже в понедельник!)\n",
    "\n",
    "- Две лекции\n",
    "- Задача о малоранговом приближении: SVD и приложения\n",
    "- Линейные системы и как их решать: знакомство с базовыми методами"
   ]
  },
  {
   "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 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  },
  "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
}