{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Линейная алгебра\n",
    "\n",
    "Подмодуль [scipy.linalg](https://docs.scipy.org/doc/scipy/reference/linalg.html#module-scipy.linalg) в целом является надмножеством над [numpy.linalg](https://numpy.org/doc/stable/reference/routines.linalg.html), т.е. кроме того, что он содержит все методы из `numpy.linalg`, в него включены еще ряд дополнительных методов.\n",
    "\n",
    "Ряд примеров приведен в [документации](https://docs.scipy.org/doc/scipy/reference/linalg.html#module-scipy.linalg). \n",
    "- [Обращение матриц](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#finding-the-inverse) (функция [linalg.inv](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.linalg.inv.html#scipy.linalg.inv)) и [решение СЛАУ](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#solving-a-linear-system) (функция [linalg.solve](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.linalg.solve.html));\n",
    "- [Вычисление норм векторов и матриц](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#computing-norms);\n",
    "- [Разложение матрицы](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#decompositions): собственные значения и вектора, сингулярное разложение, LU и QR разложения и другое;\n",
    "- [Псевдообращение матриц](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#generalized-inverse);\n",
    "- [Методы создания специальных матриц](https://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#special-matrices), [блочно-диагональных](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.linalg.block_diag.html#scipy.linalg.block_diag) и многих других;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Обращение матриц\n",
    "\n",
    "Функция [linalg.inv](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.linalg.inv.html#scipy.linalg.inv) позволяет обратить матрицу. \n",
    "\n",
    "Чтобы убедиться, что она возвращает обратную матрицу, проверим выполнение равенства $A A^{-1} = A^{-1} A = I$, где $I$ --- единичная матрица."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 2.]\n",
      " [3. 4.]]\n",
      "[[-2.   1. ]\n",
      " [ 1.5 -0.5]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg\n",
    "\n",
    "A = np.array([[1., 2.], [3., 4.]])\n",
    "A_inv = linalg.inv(A)\n",
    "\n",
    "print(A, A_inv, sep=\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.0000000e+00 0.0000000e+00]\n",
      " [8.8817842e-16 1.0000000e+00]]\n",
      "[[1.00000000e+00 0.00000000e+00]\n",
      " [1.11022302e-16 1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "print(A @ A_inv)\n",
    "print(A_inv @ A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Как видим, получилась не совсем единичная матрица, но т.к. все вычисления с числами с плавающей точкой производятся в компьютере не точно, большего ожидать нельзя. \n",
    "\n",
    "Сверять целый массив глазами не очень удобно. Давайте упростим себе задачу и посчитаем норму невязки:\n",
    "\n",
    "$$\n",
    "R_1 = || A A^{-1} - I ||\n",
    "$$\n",
    "\n",
    "и\n",
    "\n",
    "$$\n",
    "R_2 = || A^{-1} A - I ||.\n",
    "$$\n",
    "\n",
    "\n",
    "Для этого воспользуемся функцией [scipy.linalg.norm](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.linalg.norm.html#scipy.linalg.norm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_1 = 9.930136612989092e-16, R_2 = 2.482534153247273e-16\n"
     ]
    }
   ],
   "source": [
    "R_1 = linalg.norm(A @ A_inv - np.eye(2))\n",
    "R_2 = linalg.norm(A_inv @ A - np.eye(2))\n",
    "print(f\"R_1 = {R_1}, R_2 = {R_2}\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "Функция `linalg.norm` позволяет считать разные нормы, в том числе и [матричные](https://ru.wikipedia.org/wiki/%D0%9D%D0%BE%D1%80%D0%BC%D0%B0_%D0%BC%D0%B0%D1%82%D1%80%D0%B8%D1%86%D1%8B). \n",
    "\n",
    "Если эта функция встречает вектор, то по умолчанию вычисляется евклидова норма, если матрицу, то по умолчанию вычисляется норма Фробениуса.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Решение СЛАУ\n",
    "\n",
    "Функция [linalg.solve](https://docs.scipy.org/doc/scipy/reference/reference/generated/scipy.linalg.solve.html) используется для решения СЛАУ вида $Ax = b$, где $A$ --- квадратная матрица размера $N \\times N$. Правая часть уравнения $b$ может быть формы `(N, NRHS)`, т.е. если $b$ --- одномерный массив, то решается обычное СЛАУ, а если $b$ --- матрица из $N$ строк и `NRHS` столбцов, то решается матричное уравнение (или сразу решается `NRHS` СЛАУ с общей матрицей $A$, а в качестве векторов правых частей выступают столбцы матрицы $b$). \n",
    " \n",
    "Для иллюстрации создадим матрицу `A` и вектор `x`. Сгенерируем вектор правой части `b = A @ x` и попробуем восстановить значение `x`, решив СЛАУ `A x = b`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Погрешность: 4.002966042486721e-16\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import linalg\n",
    "\n",
    "A = np.array([[1., 2.], [3., 4.]])\n",
    "x = np.array([1., 1.])\n",
    "b = A @ x\n",
    "\n",
    "x_approx = linalg.solve(A, b)\n",
    "error = linalg.norm(x_approx - x)\n",
    "print(f\"Погрешность: {error}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Проиллюстрируем одновременное решение нескольких СЛАУ. Одновременно решим уравнение $A x = b$ и $A y = c$, где $A$ --- общая матрица, а правые части $b$ и $c$ --- разные. Для этого составим матрицу $B$, которая состоит из столбцов $b$ и $c$. Решение $X$ матричного уравнения $A X = B$ --- матрица, столбцы которой --- решения исходных СЛАУ $x$ и $y$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Точные решения: [1. 1.] и [42. 42.]\n",
      "Численные: [1. 1.] и [42. 42.]\n"
     ]
    }
   ],
   "source": [
    "y = np.array([42., 42.])\n",
    "c = A @ y\n",
    "print(f\"Точные решения: {x} и {y}\")\n",
    "\n",
    "B = np.vstack([b, c]).T\n",
    "X = linalg.solve(A, B)\n",
    "\n",
    "x, y = X.T\n",
    "print(f\"Численные: {x} и {y}\")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "8617202e12f254480e1fae3258716b685f1a56bcbf234a446366b4fd3345ed22"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}