codes/14_linear-algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "📝 **Author:** Amirhossein Heydari - 📧 **Email:** AmirhosseinHeydari78@gmail.com - 📍 **Linktree:** [linktr.ee/mr_pylin](https://linktr.ee/mr_pylin)\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NumPy - Linear Algebra\n",
    "📝 Doc:\n",
    "   - Linear algebra (`numpy.linalg`): [numpy.org/doc/stable/reference/routines.linalg.html](https://numpy.org/doc/stable/reference/routines.linalg.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix and vector products\n",
    "\n",
    "<table style=\"margin: 0 auto;\">\n",
    "   <thead>\n",
    "      <tr>\n",
    "         <th style=\"text-align: center;\">Function</th>\n",
    "         <th style=\"text-align: center;\">Operator</th>\n",
    "         <th style=\"text-align: center;\">Description</th>\n",
    "         <th style=\"text-align: center;\">Details</th>\n",
    "      </tr>\n",
    "   </thead>\n",
    "   <tbody>\n",
    "      <tr>\n",
    "         <td><code>np.dot</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Dot product of two arrays</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.dot.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.vdot</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Return the dot product of two vectors</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.vdot.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.vecdot</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Vector dot product of two arrays</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.vecdot.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.inner</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Inner product of two arrays</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.inner.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.outer</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Compute the outer product of two vectors</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.outer.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.matmul</code></td>\n",
    "         <td style=\"text-align: center;\"><code>@</code></td>\n",
    "         <td>Matrix product of two arrays</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.matmul.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.tensordot</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Compute tensor dot product along specified axes</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.tensordot.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.einsum</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Evaluates the Einstein summation convention on the operands</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.einsum.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.kron</code></td>\n",
    "         <td style=\"text-align: center;\">-</td>\n",
    "         <td>Kronecker product of two arrays</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.kron.html\">link</a></td>\n",
    "      </tr>\n",
    "   </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dot/Scalar Product\n",
    "Let $\\mathbf{u}$ and $\\mathbf{v}$ be two vectors:\n",
    "$$\\mathbf{u} = \\begin{bmatrix} u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_n \\end{bmatrix}, \\quad \\mathbf{v} = \\begin{bmatrix} v_1 \\\\ v_2 \\\\ \\vdots \\\\ v_n \\end{bmatrix}$$\n",
    "\n",
    "The dot product of $\\mathbf{u}$ and $\\mathbf{v}$ is:\n",
    "$$\\mathbf{u} \\cdot \\mathbf{v} = u_1 v_1 + u_2 v_2 + \\cdots + u_n v_n$$\n",
    "\n",
    "This can be written as $\\mathbf{u}^T \\mathbf{v}$ (the multiplication of two matrices):\n",
    "$$\\mathbf{u}^T \\mathbf{v} = \\begin{bmatrix} u_1 & u_2 & \\cdots & u_n \\end{bmatrix} \\begin{bmatrix} v_1 \\\\ v_2 \\\\ \\vdots \\\\ v_n \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dot_1: 32\n",
      "dot_2: 32\n"
     ]
    }
   ],
   "source": [
    "arr_1d_1 = np.array([1, 2, 3])\n",
    "arr_1d_2 = np.array([4, 5, 6])\n",
    "\n",
    "# np.dot on 1D arrays performs the dot product\n",
    "dot_1 = np.dot(arr_1d_1, arr_1d_2)\n",
    "dot_2 = np.dot(arr_1d_2, arr_1d_1)\n",
    "\n",
    "# log\n",
    "print(f\"dot_1: {dot_1}\")\n",
    "print(f\"dot_2: {dot_2}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dot Product for 2D Arrays\n",
    "For matrices, the dot product is usually referred to as matrix multiplication."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dot_3:\n",
      "[[ 7  9]\n",
      " [17 19]]\n",
      "--------------------------------------------------\n",
      "dot_4:\n",
      "[[ 6 10]\n",
      " [14 20]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_1 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_2 = np.array([[3, 1], [2, 4]])\n",
    "\n",
    "# np.dot on 2D arrays performs the matrix multiplication [using <np.matmul(a, b)> or <a @ b> is preferred]\n",
    "dot_3 = np.dot(arr_2d_1, arr_2d_2)\n",
    "dot_4 = np.dot(arr_2d_2, arr_2d_1)\n",
    "\n",
    "# log\n",
    "print(f\"dot_3:\\n{dot_3}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"dot_4:\\n{dot_4}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### np.vdot vs np.dot\n",
    "   - `np.vdot` conjugates the first argument if it is complex. Unlike `np.dot`, it always returns a scalar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dot_5: [2 4 6]\n",
      "dot_6: [2 4 6]\n"
     ]
    }
   ],
   "source": [
    "arr_0d_1 = np.array(2)\n",
    "arr_1d_3 = np.array([1, 2, 3])\n",
    "\n",
    "# np.vdot on arrays which one of them is scalar, performs the multiply [using <np.multiply(a, b)> or <a * b> is preferred]\n",
    "dot_5 = np.dot(arr_0d_1, arr_1d_3)\n",
    "dot_6 = np.dot(arr_1d_3, arr_0d_1)\n",
    "\n",
    "# log\n",
    "print(f\"dot_5: {dot_5}\")\n",
    "print(f\"dot_6: {dot_6}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vdot_1: (9+4j)\n",
      "dot_7 : (11+0j)\n",
      "--------------------------------------------------\n",
      "vdot_2: 27\n",
      "--------------------------------------------------\n",
      "dot_8 :\n",
      "[[ 7  9]\n",
      " [17 19]]\n"
     ]
    }
   ],
   "source": [
    "arr_1d_4 = np.array([1, 2-1j, 3])\n",
    "arr_1d_5 = np.array([3, 2+1j, 1])\n",
    "arr_2d_3 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_4 = np.array([[3, 1], [2, 4]])\n",
    "\n",
    "# np.vdot vs np.dot\n",
    "vdot_1 = np.vdot(arr_1d_4, arr_1d_5)  # complex conjugate of the first argument is used\n",
    "dot_7 = np.dot(arr_1d_4, arr_1d_5)\n",
    "\n",
    "vdot_2 = np.vdot(arr_2d_3, arr_2d_4)\n",
    "dot_8 = np.dot(arr_2d_3, arr_2d_4)\n",
    "\n",
    "# log\n",
    "print(f\"vdot_1: {vdot_1}\")\n",
    "print(f\"dot_7 : {dot_7}\")\n",
    "print('-' * 50)\n",
    "print(f\"vdot_2: {vdot_2}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"dot_8 :\\n{dot_8}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inner Product\n",
    "Inner Product of 1D arrays is the same concept as dot/scalar product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner_1: 32\n",
      "dot_9  : 32\n"
     ]
    }
   ],
   "source": [
    "arr_1d_6 = np.array([1, 2, 3])\n",
    "arr_1d_7 = np.array([4, 5, 6])\n",
    "\n",
    "# np.inner on 1D arrays does the same calculation as np.dot\n",
    "inner_1 = np.inner(arr_1d_6, arr_1d_7)\n",
    "dot_9 = np.dot(arr_1d_6, arr_1d_7)\n",
    "\n",
    "# log\n",
    "print(f\"inner_1: {inner_1}\")\n",
    "print(f\"dot_9  : {dot_9}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inner Product for 2d Arrays\n",
    "Computes the inner product of each pair of rows from the first matrix with each row from the second matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner_2:\n",
      "[[ 5 10]\n",
      " [13 22]]\n",
      "--------------------------------------------------\n",
      "dot_10 :\n",
      "[[ 7  9]\n",
      " [17 19]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_5 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_6 = np.array([[3, 1], [2, 4]])\n",
    "\n",
    "# np.inner performs different operations on 2D arrays than np.dot\n",
    "inner_2 = np.inner(arr_2d_5, arr_2d_6)\n",
    "dot_10 = np.dot(arr_2d_5, arr_2d_6)\n",
    "\n",
    "# log\n",
    "print(f\"inner_2:\\n{inner_2}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"dot_10 :\\n{dot_10}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outer Product\n",
    "Let $\\mathbf{u}$ and $\\mathbf{v}$ be two vectors:\n",
    "$$\\mathbf{u} = \\begin{bmatrix} u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_n \\end{bmatrix}, \\quad \\mathbf{v} = \\begin{bmatrix} v_1 \\\\ v_2 \\\\ \\vdots \\\\ v_n \\end{bmatrix}$$\n",
    "The outer product of $\\mathbf{u}$ and $\\mathbf{v}$ is:\n",
    "$$\\mathbf{u} \\otimes \\mathbf{v} = \\begin{bmatrix} u_1 v_1 & u_1 v_2 & \\cdots & u_1 v_n \\\\ u_2 v_1 & u_2 v_2 & \\cdots & u_2 v_n \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ u_n v_1 & u_n v_2 & \\cdots & u_n v_n \\end{bmatrix}$$\n",
    "This can be written as $\\mathbf{u} \\mathbf{v}^T$ (the multiplication of two matrices):\n",
    "$$\\mathbf{u} \\mathbf{v}^T = \\begin{bmatrix} u_1 \\\\ u_2 \\\\ \\vdots \\\\ u_n \\end{bmatrix} \\begin{bmatrix} v_1 & v_2 & \\cdots & v_n \\end{bmatrix}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "outer_1:\n",
      "[[ 4  5  6]\n",
      " [ 8 10 12]\n",
      " [12 15 18]]\n"
     ]
    }
   ],
   "source": [
    "arr_1d_8 = np.array([1, 2, 3])\n",
    "arr_1d_9 = np.array([4, 5, 6])\n",
    "\n",
    "# outer\n",
    "outer_1 = np.outer(arr_1d_8, arr_1d_9)\n",
    "\n",
    "# log\n",
    "print(f\"outer_1:\\n{outer_1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outer Product for 2D Arrays\n",
    "Flattens the matrices into 1D arrays and then computes the outer product of these flattened arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "outer_2:\n",
      "[[ 3  1  2  4]\n",
      " [ 6  2  4  8]\n",
      " [ 9  3  6 12]\n",
      " [12  4  8 16]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_7 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_8 = np.array([[3, 1], [2, 4]])\n",
    "\n",
    "# outer\n",
    "outer_2 = np.outer(arr_2d_7, arr_2d_8)\n",
    "\n",
    "# log\n",
    "print(f\"outer_2:\\n{outer_2}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Multiplication\n",
    "Let $\\mathbf{A}$ and $\\mathbf{B}$ be two matrices:\n",
    "$$\\mathbf{A} = \\begin{bmatrix} a_{11} & a_{12} & \\cdots & a_{1n} \\\\ a_{21} & a_{22} & \\cdots & a_{2n} \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ a_{m1} & a_{m2} & \\cdots & a_{mn} \\end{bmatrix}, \\quad \\mathbf{B} = \\begin{bmatrix} b_{11} & b_{12} & \\cdots & b_{1p} \\\\ b_{21} & b_{22} & \\cdots & b_{2p} \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ b_{n1} & b_{n2} & \\cdots & b_{np} \\end{bmatrix}$$\n",
    "\n",
    "The matrix multiplication of $\\mathbf{A}$ and $\\mathbf{B}$ is:\n",
    "$$\\mathbf{C} = \\mathbf{A} \\cdot \\mathbf{B} = \\begin{bmatrix} c_{11} & c_{12} & \\cdots & c_{1p} \\\\ c_{21} & c_{22} & \\cdots & c_{2p} \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ c_{m1} & c_{m2} & \\cdots & c_{mp} \\end{bmatrix}$$\n",
    "\n",
    "where each element $c_{ij}$ of the resulting matrix $\\mathbf{C}$ is computed as:\n",
    "$$c_{ij} = \\sum_{k=1}^{n} a_{ik} b_{kj}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matmul_1:\n",
      "[[19 22]\n",
      " [43 50]]\n",
      "--------------------------------------------------\n",
      "matmul_2:\n",
      "[[19 22]\n",
      " [43 50]]\n",
      "--------------------------------------------------\n",
      "matmul_3:\n",
      "[[23 34]\n",
      " [31 46]]\n",
      "--------------------------------------------------\n",
      "matmul_4:\n",
      "[[23 34]\n",
      " [31 46]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_9 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_10 = np.array([[5, 6], [7, 8]])\n",
    "\n",
    "# matmul\n",
    "matmul_1 = np.matmul(arr_2d_9, arr_2d_10)\n",
    "matmul_2 = arr_2d_9 @ arr_2d_10\n",
    "\n",
    "matmul_3 = np.matmul(arr_2d_10, arr_2d_9)\n",
    "matmul_4 = arr_2d_10 @ arr_2d_9\n",
    "\n",
    "# log\n",
    "print(f\"matmul_1:\\n{matmul_1}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"matmul_2:\\n{matmul_2}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"matmul_3:\\n{matmul_3}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"matmul_4:\\n{matmul_4}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fourier_coefficients_1:\n",
      "[[10 -2]\n",
      " [-4  0]]\n",
      "--------------------------------------------------\n",
      "fourier_coefficients_2:\n",
      "[[10.+0.j -2.+0.j]\n",
      " [-4.+0.j  0.+0.j]]\n",
      "--------------------------------------------------\n",
      "reconstructed_image_1:\n",
      "[[1. 2.]\n",
      " [3. 4.]]\n",
      "--------------------------------------------------\n",
      "reconstructed_image_2:\n",
      "[[1. 2.]\n",
      " [3. 4.]]\n"
     ]
    }
   ],
   "source": [
    "image_1 = np.array([[1, 2], [3, 4]])\n",
    "fourier_2x2_forward_basis_images_1 = np.array([[[[1, 1], [1, 1]], [[1, -1], [1, -1]]], [[[1, 1], [-1, -1]], [[1, -1], [-1, 1]]]])\n",
    "fourier_2x2_inverse_basis_images_1 = np.conj(fourier_2x2_forward_basis_images_1).transpose(0, 1, 3, 2)\n",
    "\n",
    "# tensordot\n",
    "fourier_coefficients_1 = np.tensordot(fourier_2x2_forward_basis_images_1, image_1, axes=([2, 3], [0, 1]))\n",
    "reconstructed_image_1 = np.tensordot(fourier_2x2_forward_basis_images_1, fourier_coefficients_1, axes=([2, 3], [0, 1])) / np.prod(image_1.shape)\n",
    "\n",
    "fourier_coefficients_2 = np.fft.fft2(image_1)\n",
    "reconstructed_image_2 = np.fft.ifft2(fourier_coefficients_2).real.clip(0, 4)\n",
    "\n",
    "# log\n",
    "print(f\"fourier_coefficients_1:\\n{fourier_coefficients_1}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"fourier_coefficients_2:\\n{fourier_coefficients_2}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"reconstructed_image_1:\\n{reconstructed_image_1}\", end=f\"\\n{'-' * 50}\\n\")\n",
    "print(f\"reconstructed_image_2:\\n{reconstructed_image_2}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Einstein Summation\n",
    "Provides a way to perform tensor operations efficiently using Einstein summation notation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "einsum_1:\n",
      "[[19 22]\n",
      " [43 50]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_11 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_12 = np.array([[5, 6], [7, 8]])\n",
    "\n",
    "# einsum\n",
    "einsum_1 = np.einsum('ij,jk->ik', arr_2d_11, arr_2d_12)\n",
    "\n",
    "# log\n",
    "print(f\"einsum_1:\\n{einsum_1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kronecker Product\n",
    "the Kronecker product, sometimes denoted by ⊗, is an operation on two matrices that produces a block matrix.\n",
    "\n",
    "Let $\\mathbf{A}$ and $\\mathbf{B}$ be two matrices:\n",
    "$$\\mathbf{A} = \\begin{bmatrix} a_{11} & a_{12} \\\\ a_{21} & a_{22} \\end{bmatrix}, \\quad \\mathbf{B} = \\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22} \\end{bmatrix}$$\n",
    "\n",
    "The Kronecker product of $\\mathbf{A}$ and $\\mathbf{B}$ is:\n",
    "$$\\mathbf{A} \\otimes \\mathbf{B} = \\begin{bmatrix}\n",
    "a_{11} \\mathbf{B} & a_{12} \\mathbf{B} \\\\\n",
    "a_{21} \\mathbf{B} & a_{22} \\mathbf{B}\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "Expanding this, we get:\n",
    "$$\\mathbf{A} \\otimes \\mathbf{B} = \\begin{bmatrix}\n",
    "a_{11} \\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22} \\end{bmatrix} &\n",
    "a_{12} \\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22} \\end{bmatrix} \\\\\n",
    "a_{21} \\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22} \\end{bmatrix} &\n",
    "a_{22} \\begin{bmatrix} b_{11} & b_{12} \\\\ b_{21} & b_{22} \\end{bmatrix}\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "a_{11} b_{11} & a_{11} b_{12} & a_{12} b_{11} & a_{12} b_{12} \\\\\n",
    "a_{11} b_{21} & a_{11} b_{22} & a_{12} b_{21} & a_{12} b_{22} \\\\\n",
    "a_{21} b_{11} & a_{21} b_{12} & a_{22} b_{11} & a_{22} b_{12} \\\\\n",
    "a_{21} b_{21} & a_{21} b_{22} & a_{22} b_{21} & a_{22} b_{22}\n",
    "\\end{bmatrix}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kron_1: [1 2 2 4]\n"
     ]
    }
   ],
   "source": [
    "arr_1d_10 = np.array([1, 2])\n",
    "arr_1d_11 = np.array([1, 2])\n",
    "\n",
    "# kron\n",
    "kron_1 = np.kron(arr_1d_10, arr_1d_11)\n",
    "\n",
    "# log\n",
    "print(f\"kron_1: {kron_1}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kron_2:\n",
      "[[1 1 2 2]\n",
      " [1 1 2 2]\n",
      " [3 3 4 4]\n",
      " [3 3 4 4]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_13 = np.array([[1, 2], [3, 4]])\n",
    "arr_2d_14 = np.array([[1, 1], [1, 1]])\n",
    "\n",
    "# kron\n",
    "kron_2 = np.kron(arr_2d_13, arr_2d_14)\n",
    "\n",
    "# log\n",
    "print(f\"kron_2:\\n{kron_2}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decompositions\n",
    "\n",
    "<table style=\"margin: 0 auto;\">\n",
    "   <thead>\n",
    "      <tr>\n",
    "         <th style=\"text-align: center;\">Function</th>\n",
    "         <th style=\"text-align: center;\">Description</th>\n",
    "         <th style=\"text-align: center;\">Details</th>\n",
    "      </tr>\n",
    "   </thead>\n",
    "   <tbody>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.cholesky</code></td>\n",
    "         <td>Cholesky decomposition</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.cholesky.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.qr</code></td>\n",
    "         <td>Compute the qr factorization of a matrix</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.qr.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.svd</code></td>\n",
    "         <td>Returns the singular values of a matrix (or a stack of matrices) <code>x</code></td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html\">link</a></td>\n",
    "      </tr>\n",
    "   </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cholesky_1:\n",
      "[[2.         0.        ]\n",
      " [1.         1.41421356]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_15 = np.array([[4, 2], [2, 3]])\n",
    "\n",
    "# Compute Cholesky decomposition of a positive-definite matrix\n",
    "cholesky_1 = np.linalg.cholesky(arr_2d_15)\n",
    "\n",
    "# log\n",
    "print(f\"cholesky_1:\\n{cholesky_1}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Q:\n",
      "[[-0.31622777 -0.9486833 ]\n",
      " [-0.9486833   0.31622777]]\n",
      "\n",
      "R:\n",
      "[[-3.16227766 -4.42718872]\n",
      " [ 0.         -0.63245553]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_16 = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "# Compute QR decomposition\n",
    "Q, R = np.linalg.qr(arr_2d_16)\n",
    "\n",
    "# log\n",
    "print(f\"Q:\\n{Q}\\n\")\n",
    "print(f\"R:\\n{R}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U:\n",
      "[[-0.40455358 -0.9145143 ]\n",
      " [-0.9145143   0.40455358]]\n",
      "\n",
      "S:\n",
      "[5.4649857  0.36596619]\n",
      "\n",
      "Vt:\n",
      "[[-0.57604844 -0.81741556]\n",
      " [ 0.81741556 -0.57604844]]\n"
     ]
    }
   ],
   "source": [
    "arr_2d_17 = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "# Compute SVD\n",
    "U, S, Vt = np.linalg.svd(arr_2d_17)\n",
    "\n",
    "# Log\n",
    "print(f\"U:\\n{U}\\n\")\n",
    "print(f\"S:\\n{S}\\n\")\n",
    "print(f\"Vt:\\n{Vt}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix eigenvalues\n",
    "\n",
    "<table style=\"margin: 0 auto;\">\n",
    "   <thead>\n",
    "      <tr>\n",
    "         <th style=\"text-align: center;\">Function</th>\n",
    "         <th style=\"text-align: center;\">Description</th>\n",
    "         <th style=\"text-align: center;\">Details</th>\n",
    "      </tr>\n",
    "   </thead>\n",
    "   <tbody>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.eig</code></td>\n",
    "         <td>Compute the eigenvalues and right eigenvectors of a square array</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.eigh</code></td>\n",
    "         <td>Return the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigh.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.eigvals</code></td>\n",
    "         <td>Compute the eigenvalues of a general matrix</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigvals.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.eigvalsh</code></td>\n",
    "         <td>Compute the eigenvalues of a complex Hermitian or real symmetric matrix</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigvalsh.html\">link</a></td>\n",
    "      </tr>\n",
    "   </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Norms and other numbers\n",
    "\n",
    "<table style=\"margin: 0 auto;\">\n",
    "   <thead>\n",
    "      <tr>\n",
    "         <th style=\"text-align: center;\">Function</th>\n",
    "         <th style=\"text-align: center;\">Description</th>\n",
    "         <th style=\"text-align: center;\">Details</th>\n",
    "      </tr>\n",
    "   </thead>\n",
    "   <tbody>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.norm</code></td>\n",
    "         <td>Matrix or vector norm</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.norm.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.det</code></td>\n",
    "         <td>Compute the determinant of an array</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.det.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.matrix_rank</code></td>\n",
    "         <td>Return matrix rank of array using SVD method</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.matrix_rank.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.trace</code></td>\n",
    "         <td>Return the sum along diagonals of the array</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.trace.html\">link</a></td>\n",
    "      </tr>\n",
    "   </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm_1 (L_1       norm): 7.0\n",
      "norm_2 (L_2       norm): 5.0\n",
      "norm_3 (L_inf     norm): 4.0\n",
      "norm_4 (Frobenius norm): 5.477225575051661\n"
     ]
    }
   ],
   "source": [
    "arr_1d_12 = np.array([3, 4])\n",
    "arr_2d_18 = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "# norm\n",
    "norm_1 = np.linalg.norm(arr_1d_12, ord=1)\n",
    "norm_2 = np.linalg.norm(arr_1d_12, ord=2)\n",
    "norm_3 = np.linalg.norm(arr_1d_12, ord=np.inf)\n",
    "norm_4 = np.linalg.norm(arr_2d_18, ord='fro')\n",
    "\n",
    "# log\n",
    "print(f\"norm_1 (L_1       norm): {norm_1}\")\n",
    "print(f\"norm_2 (L_2       norm): {norm_2}\")\n",
    "print(f\"norm_3 (L_inf     norm): {norm_3}\")\n",
    "print(f\"norm_4 (Frobenius norm): {norm_4}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "determinant_1: -2.0000000000000004\n"
     ]
    }
   ],
   "source": [
    "arr_2d_19 = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "# determinants\n",
    "determinant_1 = np.linalg.det(arr_2d_19)\n",
    "\n",
    "# log\n",
    "print(f\"determinant_1: {determinant_1}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "trace_1: 15\n"
     ]
    }
   ],
   "source": [
    "arr_2d_20 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
    "\n",
    "# trace\n",
    "trace_1 = np.linalg.trace(arr_2d_20)\n",
    "\n",
    "# log\n",
    "print(f\"trace_1: {trace_1}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rank_1: 2\n"
     ]
    }
   ],
   "source": [
    "arr_2d_21 = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])\n",
    "\n",
    "# rank\n",
    "rank_1 = np.linalg.matrix_rank(arr_2d_21)\n",
    "\n",
    "# log\n",
    "print(f\"rank_1: {rank_1}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving equations and inverting matrices\n",
    "\n",
    "<table style=\"margin: 0 auto;\">\n",
    "   <thead>\n",
    "      <tr>\n",
    "         <th style=\"text-align: center;\">Function</th>\n",
    "         <th style=\"text-align: center;\">Description</th>\n",
    "         <th style=\"text-align: center;\">Details</th>\n",
    "      </tr>\n",
    "   </thead>\n",
    "   <tbody>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.solve</code></td>\n",
    "         <td>Solve a linear matrix equation, or system of linear scalar equations</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.lstsq</code></td>\n",
    "         <td>Return the least-squares solution to a linear matrix equation</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.inv</code></td>\n",
    "         <td>Compute the inverse of a matrix</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html\">link</a></td>\n",
    "      </tr>\n",
    "      <tr>\n",
    "         <td><code>np.linalg.pinv</code></td>\n",
    "         <td>Compute the (Moore-Penrose) pseudo-inverse of a matrix</td>\n",
    "         <td style=\"text-align: center;\"><a href=\"https://numpy.org/doc/stable/reference/generated/numpy.linalg.pinv.html\">link</a></td>\n",
    "      </tr>\n",
    "   </tbody>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\mathbf{A}$ be an $m \\times n$ matrix, $\\mathbf{x}$ be an $n \\times 1$ column vector, and $\\mathbf{b}$ be an $m \\times 1$ column vector. The system of linear equations can be written in matrix form as:\n",
    "\n",
    "$$\\mathbf{A} \\mathbf{x} = \\mathbf{b}$$\n",
    "\n",
    "where:\n",
    "\n",
    "$$\\mathbf{A} = \\begin{bmatrix}\n",
    "a_{11} & a_{12} & \\cdots & a_{1n} \\\\\n",
    "a_{21} & a_{22} & \\cdots & a_{2n} \\\\\n",
    "\\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "a_{m1} & a_{m2} & \\cdots & a_{mn}\n",
    "\\end{bmatrix}, \\quad\n",
    "\\mathbf{x} = \\begin{bmatrix}\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "\\vdots \\\\\n",
    "x_n\n",
    "\\end{bmatrix}, \\quad\n",
    "\\mathbf{b} = \\begin{bmatrix}\n",
    "b_1 \\\\\n",
    "b_2 \\\\\n",
    "\\vdots \\\\\n",
    "b_m\n",
    "\\end{bmatrix}$$\n",
    "The $i$-th equation in the system is:\n",
    "$$\\sum_{j=1}^{n} a_{ij} x_j = b_i$$\n",
    "\n",
    "Let $\\mathbf{A}$ be an $m \\times n$ matrix, $\\mathbf{x}$ be an $n \\times 1$ column vector, and $\\mathbf{b}$ be an $m \\times 1$ column vector. The system of linear equations can be written in matrix form as:\n",
    "\n",
    "$$\\mathbf{A} \\mathbf{x} = \\mathbf{b}$$\n",
    "\n",
    "### Case 1: Square Matrix and Invertibility\n",
    "\n",
    "If $\\mathbf{A}$ is a square matrix (i.e., $m = n$) and invertible, the solution to the system can be obtained by multiplying both sides of the equation by the inverse of $\\mathbf{A}$:\n",
    "\n",
    "$$\\mathbf{x} = \\mathbf{A}^{-1} \\mathbf{b}$$\n",
    "\n",
    "where $\\mathbf{A}^{-1}$ denotes the inverse of $\\mathbf{A}$. The matrix $\\mathbf{A}$ is invertible if and only if its determinant is non-zero:\n",
    "\n",
    "$$\\det(\\mathbf{A}) \\neq 0$$\n",
    "\n",
    "Thus, the solution vector $\\mathbf{x}$ is:\n",
    "\n",
    "$$\\mathbf{x} = \\mathbf{A}^{-1} \\mathbf{b}$$\n",
    "\n",
    "### Case 2: Non-Square Matrix or Singular Matrix\n",
    "\n",
    "If $\\mathbf{A}$ is not a square matrix (i.e., $m \\neq n$) or is a square but singular matrix (i.e., $\\det(\\mathbf{A}) = 0$), the inverse does not exist. In such cases, we use the Moore-Penrose pseudo-inverse to find a least-squares solution.\n",
    "\n",
    "The Moore-Penrose pseudo-inverse of $\\mathbf{A}$, denoted $\\mathbf{A}^+$, is defined as:\n",
    "\n",
    "$$\\mathbf{A}^+ = (\\mathbf{A}^\\top \\mathbf{A})^{-1} \\mathbf{A}^\\top$$\n",
    "\n",
    "when $\\mathbf{A}^\\top \\mathbf{A}$ is invertible. The solution vector $\\mathbf{x}$ can then be approximated by:\n",
    "\n",
    "$$\\mathbf{x} = \\mathbf{A}^+ \\mathbf{b}$$\n",
    "\n",
    "For general cases, where $\\mathbf{A}^\\top \\mathbf{A}$ is not invertible, the pseudo-inverse can be computed using Singular Value Decomposition (SVD):\n",
    "\n",
    "$$\\mathbf{A} = \\mathbf{U} \\mathbf{\\Sigma} \\mathbf{V}^\\top$$\n",
    "\n",
    "where $\\mathbf{\\Sigma}$ is a diagonal matrix with singular values. The pseudo-inverse is given by:\n",
    "\n",
    "$$\\mathbf{A}^+ = \\mathbf{V} \\mathbf{\\Sigma}^+ \\mathbf{U}^\\top$$\n",
    "\n",
    "where $\\mathbf{\\Sigma}^+$ is the pseudo-inverse of $\\mathbf{\\Sigma}$.\n",
    "\n",
    "Thus, the least-squares solution vector $\\mathbf{x}$ is:\n",
    "\n",
    "$$\\mathbf{x} = \\mathbf{A}^+ \\mathbf{b}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x: [0.5  1.75]\n"
     ]
    }
   ],
   "source": [
    "A_1 = np.array([[3, 2], [1, 2]])\n",
    "b_1 = np.array([5, 4])\n",
    "\n",
    "# solve\n",
    "x = np.linalg.solve(A_1, b_1)\n",
    "\n",
    "# log\n",
    "print(f\"x: {x}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x: [1.80069003e-16 1.00000000e+00]\n"
     ]
    }
   ],
   "source": [
    "A_2 = np.array([[1, 1], [1, 2], [1, 3]])\n",
    "b_2 = np.array([1, 2, 3])\n",
    "\n",
    "# solve using least square error\n",
    "x, residuals, rank, s = np.linalg.lstsq(A_2, b_2, rcond=None)\n",
    "\n",
    "# log\n",
    "print(f\"x: {x}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_inv:\n",
      "[[-2.   1. ]\n",
      " [ 1.5 -0.5]]\n",
      "\n",
      "eye_1:\n",
      "[[1.0000000e+00 0.0000000e+00]\n",
      " [8.8817842e-16 1.0000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "A_3 = np.array([[1, 2], [3, 4]])\n",
    "\n",
    "# inverse\n",
    "A_inv = np.linalg.inv(A_3)\n",
    "eye_1 = A_3 @ A_inv\n",
    "\n",
    "# log\n",
    "print(f\"A_inv:\\n{A_inv}\\n\")\n",
    "print(f\"eye_1:\\n{eye_1}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A_pinv:\n",
      "[[-1.33333333 -0.33333333  0.66666667]\n",
      " [ 1.08333333  0.33333333 -0.41666667]]\n",
      "\n",
      "eye_2:\n",
      "[[ 1.0000000e+00  0.0000000e+00]\n",
      " [-4.4408921e-16  1.0000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "A_4 = np.array([[1, 2], [3, 4], [5, 6]])\n",
    "\n",
    "# pseudo-inverse\n",
    "A_pinv = np.linalg.pinv(A_4)\n",
    "eye_2 = A_pinv @ A_4\n",
    "\n",
    "# log\n",
    "print(f\"A_pinv:\\n{A_pinv}\\n\")\n",
    "print(f\"eye_2:\\n{eye_2}\")"
   ]
  }
 ],
 "metadata": {
  "author_email": "AmirhosseinHeydari78@gmail.com",
  "author_github": "https://github.com/mr-pylin",
  "author_name": "Amirhossein Heydari",
  "kernelspec": {
   "display_name": "venv",
   "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.3"
  },
  "origin_repo": "https://github.com/mr-pylin/numpy-workshop"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}