{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 04 - Linear algebra (cont'd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterative and gradient methods for solving linear algebraic equation systems, Jacobi method, Gauss-Seidel method, successive over-relaxation method, conjugate gradient method, power iteration and eigensystems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterative methods for linear algebraic equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Jacobi method\n",
    "\n",
    "A system of $ n $ linear algebraic equations $ \\mathbb{A} x = b $ with [diagonally dominant](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix) coefficient matrix $ \\mathbb{A} $ can be solved iteratively using [Jacobi method](https://en.wikipedia.org/wiki/Jacobi_method):\n",
    "\n",
    "$$\n",
    "x_i^{(k+1)} = \\frac{1}{a_{ii}} \\left( b_i - \\sum_{j = 1}^{i-1} a_{ij} x_j^{(k)} - \\sum_{j = i + 1}^{n} a_{ij} x_j^{(k)} \\right), \\qquad i = 0, 1, \\dots, n-1,\n",
    "$$\n",
    "\n",
    "where $ x^{(k)} $ and $ x^{(k+1)} $ are, respectively, the k-th and k+1-th iterations of $ x $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "**04.1 Example:** Implement the Jacobi method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jacobi_method(A, b, error_tolerance):\n",
    "    \"\"\"\n",
    "    Solves system of linear equations iteratively using Jacobi's algorithm.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n diagonally dominant matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "        error_tolerance (float): Error tolerance\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3) + 10. * np.eye(3) # create diagonally dominant matrix to ensure convergence\n",
    "b = np.random.rand(3)\n",
    "np.testing.assert_almost_equal(jacobi_method(A, b, 1.0e-15), la.solve(A, b), decimal=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauss-Seidel method\n",
    "\n",
    "A system of $ n $ linear algebraic equations $ \\mathbb{A} x = b $ with [diagonally dominant](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix) coefficient matrix $ \\mathbb{A} $ can be solved iteratively using [Gauss-Seidel method](https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method):\n",
    "\n",
    "$$\n",
    "x_i^{(k+1)} = \\frac{1}{a_{ii}} \\left( b_i - \\sum_{j = 1}^{i-1} a_{ij} x_j^{(k+1)} - \\sum_{j = i + 1}^{n} a_{ij} x_j^{(k)} \\right), \\qquad i = 0, 1, \\dots, n-1,\n",
    "$$\n",
    "\n",
    "where $ x^{(k)} $ and $ x^{(k+1)} $ are, respectively, the k-th and k+1-th iterations of $ x $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "**Exercise 04.2:** Implement the Gauss-Seidel method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_seidel_method(A, b, error_tolerance):\n",
    "    \"\"\"\n",
    "    Solves system of linear equations iteratively using Gauss-Seidel's algorithm.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n diagonally dominant matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "        error_tolerance (float): Error tolerance\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3) + 10. * np.eye(3) # create diagonally dominant matrix to ensure convergence\n",
    "b = np.random.rand(3)\n",
    "np.testing.assert_almost_equal(gauss_seidel_method(A, b, 1.0e-15), la.solve(A, b), decimal=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Successive over-relaxation method\n",
    "\n",
    "A system of $ n $ linear algebraic equations $ \\mathbb{A} x = b $ with [diagonally dominant](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix) coefficient matrix $ \\mathbb{A} $ can be solved iteratively using the method of [successive over-relaxation](https://en.wikipedia.org/wiki/Successive_over-relaxation):\n",
    "\n",
    "$$\n",
    "x_i^{(k+1)} = (1 - \\omega) x_i^{(k)}  + \\frac{\\omega}{a_{ii}} \\left( b_i - \\sum_{j = 1}^{i-1} a_{ij} x_j^{(k+1)} - \\sum_{j = i + 1}^{n} a_{ij} x_j^{(k)} \\right), \\qquad i = 0, 1, \\dots, n-1,\n",
    "$$\n",
    "\n",
    "where $ x^{(k)} $ and $ x^{(k+1)} $ are, respectively, the k-th and k+1-th iterations of $ x $ and $ \\omega \\in (0, 2) $ is the relaxation factor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "    \n",
    "**Exercise 04.3:** Implement the successive over-relaxation (SOR) method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def successive_overrelaxation_method(A, b, omega, error_tolerance):\n",
    "    \"\"\"\n",
    "    Solves system of linear equations iteratively using successive over-relaxation (SOR) method.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "        omega (float): Relaxation factor\n",
    "        error_tolerance (float): Error tolerance\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3) + 10. * np.eye(3) # create diagonally dominant matrix to ensure convergence\n",
    "b = np.random.rand(3)\n",
    "omega = 2.0 * np.random.rand()\n",
    "np.testing.assert_almost_equal(successive_overrelaxation_method(A, b, 0.9, 1.0e-15), la.solve(A, b), decimal=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 04.4:** Compare the number of iterations required to solve the system of 100 linear algebraic equations using Jacobi, Gauss-Seidel, and SOR methods. Try to find an optimal value of relaxation factor $ \\omega $ in the case of SOR method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient methods for linear algebraic equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conjugate gradient method\n",
    "\n",
    "A system of $ n $ linear algebraic equations $ \\mathbb{A} x = b $ with real, [symmetric](https://en.wikipedia.org/wiki/Symmetric_matrix), and [positive-definite](https://en.wikipedia.org/wiki/Definite_symmetric_matrix) coefficient matrix $ \\mathbb{A} $ can be solved using the [conjugate gradient](https://en.wikipedia.org/wiki/Conjugate_gradient_method) method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 04.5:** Implement the conjugate gradient method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conjugate_gradient_method(A, b, error_tolerance):\n",
    "    \"\"\"\n",
    "    Solves system of linear equations using conjugate gradient method.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n real, symmetric, and positive-definite matrix\n",
    "        b (array_like): RHS vector of size n\n",
    "        error_tolerance (float): Error tolerance\n",
    "    Returns:\n",
    "        numpy.ndarray: Vector of solution\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.solve`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3)\n",
    "A = np.tril(A) + np.tril(A, -1).T + 10 * np.eye(3) # create symmetric and diagonally dominant matrix to ensure convergence\n",
    "b = np.random.rand(3)\n",
    "np.testing.assert_almost_equal(conjugate_gradient_method(A, b, 1.0e-15), la.solve(A, b), decimal=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenvalues and eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Power iteration method\n",
    "\n",
    "The [power iteration](https://en.wikipedia.org/wiki/Power_iteration) method finds the greatest (in absolute value) eigen value and its corresponding eigenvector of a [diagonalizable](https://en.wikipedia.org/wiki/Diagonalizable_matrix) matrix $ \\mathbb{A} $. Starting with the vector $ v_{0} $, the algorithm is described by the recurrence\n",
    "\n",
    "$$\n",
    "v_{k+1} = \\frac{\\mathbb{A} v_k}{\\| \\mathbb{A} v_k \\|}.\n",
    "$$\n",
    "\n",
    "The sequence above converges to an eigenvector associated with the dominant eigenvalue if the following two conditions are satisfied:\n",
    "1. The matrix $ \\mathbb{A} $ must have an eigenvalue that is strictly greater in magnitude than its other eigenvalues.\n",
    "2. The starting vector $ v_{0} $ has a nonzero component in the direction of an eigenvector associated with the dominant eigenvalue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 04.6:** Implement the power iteration method.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def power_iteration(A, max_it):\n",
    "    \"\"\"\n",
    "    Finds the greatest (in absolute value) eigen value of given matrix and its corresponding eigenvector.\n",
    "    Args:\n",
    "        A (array_like): A n-by-n diagonalizable matrix\n",
    "        max_it (int): Maximum number of iterations\n",
    "    Returns:\n",
    "        numpy.ndarray: Eigenvector corresponding to a greatest eigenvalue (in absolute value)\n",
    "        float: Greatest eigenvalue (in absolute value)\n",
    "    \"\"\"\n",
    "    \n",
    "    # add your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may evaluate the correctness of your implementation using the [`scipy.linalg.eigvals`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eigvals.html) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.random.rand(3, 3)\n",
    "e_vec, e_val = power_iteration(A, 50)\n",
    "np.testing.assert_almost_equal(e_val, np.max(np.abs(la.eigvals(A))), decimal=15)"
   ]
  }
 ],
 "metadata": {
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}