numeric-linalg

 {
  "cells": [
   {
    "cell_type": "markdown",
    "id": "e1791697-e733-4721-9e10-fcdcfbda9064",
    "metadata": {},
    "source": [
     "# Solving Linear systems in SciPy\n",
     "\n",
     "A linear system is a system of equations of the form\n",
     "$$ \\left\\{ \\begin{aligned} a_{1 1} x_1 + a_{1 2} x_2 + \\cdots + a_{1 n} x_n &= b_1 \\\\ a_{2 1} x_1 + a_{2 2} x_2 + \\cdots + a_{2 n} x_n &= b_2 \\\\ & \\vdots \\\\ a_{n 1} x_1 + a_{n 2} x_2 + \\cdots + a_{n n} x_n &= b_n\\end{aligned} \\right.$$\n",
     "on the variables $x_1, \\ldots, x_n$.\n",
     "\n",
     "Solving this system is equivalent to solving the equation $A x = b$ on $x$ where $A = (a_{ij})_{ij}$ is a $n\\times n$ matrix and $x = (x_1, \\ldots, x_n) \\; \\& \\; b = (b_1, \\ldots, b_n)$ are vectors, which can always be done provided $A$ is invertible."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "26bc6a30-7853-4c78-adfb-320f0a65dd10",
    "metadata": {},
    "source": [
     "In SciPy, we can solve linear systems using the `la.solve` function."
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 1,
    "id": "b1ced47f-6783-48ed-a4e4-4f1f4e2b8835",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
        "array([[-0.12672176],\n",
        "       [ 0.1046832 ],\n",
        "       [ 1.19008264]])"
       ]
      },
      "execution_count": 1,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "import numpy as np\n",
     "import scipy.linalg as la\n",
     "\n",
     "A, b = np.array([[6,15,1],[8,7,12],[2,7,8]]), np.array([[2], [14], [10]])\n",
     "la.solve(A, b)"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "bfdacdbf-150f-4b4d-8072-ee18758d3b60",
    "metadata": {},
    "source": [
     "## But how does `la.solve` work???\n",
     "\n",
     "Internally, the `la.solve` function uses the the [LAPACK library](https://netlib.org/lapack), a Fortran package for numerical linear algebra. The LAPACK generic linear solver algorithm goes something like the following:\n",
     "\n",
     "1. Decompose $A$ as $A = PL U$, where $P$ is permutation matrix, $L$ is a lower triangular matrix with $1$ in the diagonal and $U$ is an upper triangular matrix.\n",
     "3. Solve $P b' = b$ for $b'$, i.e. compute $b' = P^{-1} b$. Since $P$ is a permutation matrix, this operation is $O(n)$.\n",
     "4. Solve $L b'' = b'$ for $b''$, i.e. compute $b'' = L^{-1} P^{-1} b$. Since $L$ is known to be lower triangular, this operation is $O(n^2)$.\n",
     "3. Solve $U x = b''$ for $x$, i.e. compute $x = U^{-1} L^{-1} P^{-1} b = A^{-1} b$. Since $U$ is known to be upper triangular, this operation is $O(n^2)$.\n",
     "\n",
     "This is implemented in the `GETRS` family of subroutines."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "5b7ad45f-ea81-4d46-9097-735cc159cf1a",
    "metadata": {},
    "source": [
     "As for the decomposition of $A$ in the first step, LAPACK uses a method called [_partial pivoting_](https://en.wikipedia.org/wiki/LU_decomposition#LU_factorization_with_partial_pivoting). A simple simple recurssive algorithm using such method might look something like the following:\n",
     "\n",
     "1. If $A = a_{11}$ is $1 \\times 1$ then take $P = L = 1$ and $U = a_{11}$.\n",
     "2. If $A$ is $n \\times n$ for $n > 1$, choose $i_0$ that maximizes $|a_{i_0, 1}|$ and consider the $n \\times n$ permutation matrix $S_{i_0}$ that swaps the first and $i_0$-th basis vectors. Searching for $i_0$ is an $O(n)$ operation.\n",
     "3. Write\n",
     "   $$S_{i_0} A = \\left( \\begin{array}{c|c} a_{i_0} & A_{12}' \\\\ \\hline A_{21}' & A_{22}' \\end{array} \\right), $$\n",
     "   where $A_{22}'$ is $(n - 1) \\times (n - 1)$ and $a_{i_0} \\ne 0$ — given $A$ is invertible. Since $S_{i_0}$ acts on $A$ by swaping the first and $i_0$-th rows, computing $S_{i_0} A$ is an $O(n)$ operation.\n",
     "4. We want to solve the equation\n",
     "   $$S_{i_0} A = \\left( \\begin{array}{c|c} 1 & 0 \\\\ \\hline 0 & P_{22} \\end{array} \\right) \\left( \\begin{array}{c|c} 1 & 0 \\\\ \\hline L_{21} & L_{22} \\end{array} \\right) \\cdot \\left( \\begin{array}{c|c} u_{11} & U_{12} \\\\ \\hline 0 & U_{22} \\end{array} \\right),$$\n",
     "   where $P_{22}$ is a permutation matrix, $L_{22}$ is lower triangular with $1$ in the diagonal entries and $U_{22}$ is upper triangular. In other words, we want to solve the equations\n",
     "   $$\n",
     "   \\begin{aligned}\n",
     "       a_{i_0} &= u_{11} & A_{12}' &= U_{12} \\\\\n",
     "       A_{21}' &= u_{11} P_{22} L_{21} & A_{22}' &= P_{22} L_{21} U_{12} + P_{22} L_{22} U_{22}.\n",
     "   \\end{aligned}\n",
     "   $$\n",
     "   We must take $u_{11} = a_{i_0}$, $U_{12} = A_{12}'$ and $L_{21} = a_{i_0}^{-1} P_{22}^{-1} A_{12}'$, so it remains to solve the bottom-right equation.\n",
     "5. Write $(A_{22}' - a_{i_0}^{-1} A_{21}' A_{12}') = P_{22} L_{22} U_{22}$, where $P_{22}$ is a permutation matrix, $L_{22}$ is lower triangular with $1$ in the diagonals and $U_{22}$ is upper triangular. Computing $A_{21}' A_{12}'$ (and thus $A_{22}' - a_{i_0}^{-1} A_{21}' A_{12}'$) is, of course, a $O(n^3)$ operation. Now since $P_{22}$ is a permutation matrix, computing $L_{21} = a_{i_0}^{-1} P_{22}^{-1} A_{12}$ is an $O(n^2)$ operation.\n",
     "7. Take\n",
     "   $$\n",
     "   \\begin{aligned}\n",
     "       L &= \\begin{pmatrix}      1 & 0      \\\\ L_{21} & L_{22} \\end{pmatrix} &\n",
     "       U &= \\begin{pmatrix} u_{11} & U_{12} \\\\      0 & U_{22} \\end{pmatrix}\n",
     "   \\end{aligned}\n",
     "   $$\n",
     "   for $L_{21}, L_{22}, u_{11}, U_{12}, U_{22}$ as above, so that\n",
     "   $$\n",
     "   S_{i_0} A = \\begin{pmatrix} 1 & 0 \\\\ 0 & P_{22} \\end{pmatrix} L U.\n",
     "   $$\n",
     "8. Hence by taking\n",
     "   $$\n",
     "   P = S_{i_0} \\begin{pmatrix} 1 & 0 \\\\ 0 & P_{22} \\end{pmatrix}\n",
     "   $$\n",
     "   we get $A = P L U$ as desired!\n",
     "\n",
     "In total, this algorithm takes $n$ recursive steps to solve $A = P L U$. Since each step involves $O(n^3)$ operations, the total complexity of our facotization algorithm is $O(n^4)$."
    ]
   },
   {
    "cell_type": "markdown",
    "id": "7c6678f2-9ae6-4daa-922c-828771c1a796",
    "metadata": {},
    "source": [
     "The LAPACK algorithm for factorizing a $m \\times n$ matrix $A$ improves on this simple concept by instead taking the decomposition\n",
     "$$\n",
     "A =\n",
     "\\left(\n",
     "\\begin{array}{c|c}\n",
     "    A_{11} & A_{12} \\\\ \\hline\n",
     "    A_{21} & A_{22}\n",
     "\\end{array}\n",
     "\\right),\n",
     "$$\n",
     "where $A_{11}$ is a $\\left\\lfloor \\frac{\\min \\{m, n\\}}{2} \\right\\rfloor \\times \\left\\lfloor \\frac{\\min \\{m, n\\}}{2} \\right\\rfloor$ matrix.\n",
     "This is implemented in the `GETRF` family of subroutines. A SciPy version of their precise algorithm might look something like the following:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": 3,
    "id": "f942e1ba-0449-4f35-b894-4dc5cad385a6",
    "metadata": {},
    "outputs": [],
    "source": [
     "import numpy as np\n",
     "import scipy.linalg as la\n",
     "\n",
     "def getrf(A: np.ndarray) -> (np.ndarray, np.ndarray, np.ndarray):\n",
     "    \"\"\"Returns the P, L, U\n",
     "\n",
     "    * A is m by n\n",
     "    * P is m by m\n",
     "    * L is m by n if m >= n and m by m if m <= n\n",
     "    * U is n by n if m >= n and m by n if m <= n\n",
     "    \"\"\"\n",
     "    m, n = A.shape\n",
     "    \n",
     "    # A is a row\n",
     "    if m == 1:\n",
     "        return np.eye(1), np.eye(1), A\n",
     "\n",
     "    # A is a column\n",
     "    elif n == 1:\n",
     "        i0 = 0\n",
     "\n",
     "        for i in range(m):\n",
     "            if abs(A[i, 0]) > abs(A[i0, 0]): i0 = i\n",
     "\n",
     "        P = np.eye(m)\n",
     "        P[0,0],  P[i0,i0] = 0, 0\n",
     "        P[i0,0], P[0,i0]  = 1, 1\n",
     "\n",
     "        if A[i0, 0] != 0:\n",
     "            L = P@A / A[i0, 0]\n",
     "            U = A[i0, 0] * np.eye(1)\n",
     "        else:\n",
     "            L = A\n",
     "            U = np.zeros((1, 1))\n",
     "\n",
     "        return P, L, U\n",
     "    else:\n",
     "        n1 = min(m, n)//2\n",
     "        n2 = n - n1\n",
     "\n",
     "        # Write\n",
     "        #\n",
     "        #   A = [[A11, A12],\n",
     "        #        [A21, A22]],\n",
     "        #\n",
     "        #   A1 = [[A11, \n",
     "        #          A21]],\n",
     "        #\n",
     "        #   A2 = [[A12, \n",
     "        #          A22]]\n",
     "        #\n",
     "        # where A11 is n1 by n1 and A22 is n2 by n2\n",
     "        A11, A12 = A[:n1,:n1], A[:n1,n1:]\n",
     "        A21, A22 = A[n1:,:n1], A[n1:,n1:]\n",
     "        A1, A2   = A[:,:n1],   A[:,n1:]\n",
     "\n",
     "        # Solve the A1 block\n",
     "        P1, L1, U11 = getrf(A1)\n",
     "\n",
     "        # Apply pivots\n",
     "        A2 = la.inv(P1) @ A2\n",
     "        A12, A22 = A2[:n1,:], A2[n1:,:]\n",
     "        \n",
     "        # Solve A12 \n",
     "        L11, L21 = L1[:n1,:], L1[n1:,:]\n",
     "        A12 = la.inv(L11) @ A12\n",
     "\n",
     "        # Update A22\n",
     "        A22 = -L21@A12 + A22\n",
     "        \n",
     "        # Solve the A22 block\n",
     "        P2, L22, U22 = getrf(A22)\n",
     "\n",
     "        # Take P = P1 @ P2_ for\n",
     "        #\n",
     "        # P2_ = [[1, 0,\n",
     "        #         0, P2]]\n",
     "        P2_ = np.eye(m)\n",
     "        P2_[n1:,n1:] = P2\n",
     "        P = P1 @ P2_\n",
     "\n",
     "        # Apply interchanges to L21\n",
     "        L21 = la.inv(P2) @ L21\n",
     "\n",
     "        # Take\n",
     "        # \n",
     "        # L = [[L11, 0],\n",
     "        #      [L21, L22]],\n",
     "        #\n",
     "        # U = [[U11, A12],\n",
     "        #      [  0, U22]]\n",
     "        if m >= n:\n",
     "            L, U = np.zeros((m, n)), np.zeros((n, n))\n",
     "        else:\n",
     "            L, U = np.zeros((m, m)), np.zeros((m, n))\n",
     "        L[:n1,:n1] = L11\n",
     "        L[n1:,:n1] = L21\n",
     "        L[n1:,n1:] = L22\n",
     "        U[:n1,:n1] = U11\n",
     "        U[:n1,n1:] = A12\n",
     "        U[n1:,n1:] = U22\n",
     "\n",
     "        return P, L, U"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "229ed6e3-493e-49c3-835a-0b8582aa6586",
    "metadata": {},
    "outputs": [],
    "source": []
   }
  ],
  "metadata": {
   "kernelspec": {
    "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
     "version": 3
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.12.6"
   }
  },
  "nbformat": 4,
  "nbformat_minor": 5
 }