part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: Power Method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Implement the power method with normalization.\n",
    "Here we've implemented the following algorithm:\n",
    "- Take a random vector $x^{(0)}$.\n",
    "- For $k = 1, 2, ..., i$ \n",
    "  - Max:\n",
    "    - $z^{(k+1)} = Ax^{(k)}$\n",
    "    - $x^{(k+1)} = \\frac{z^{(k+1)}}{||z^{(k+1)}||}$\n",
    "    - Look for $i$ such that $z_i^{(k+1)}$ is the largest element in the vector $z^{(k+1)}$.\n",
    "    - $\\lambda^{(k+1)} = \\frac{z_i^{(k+1)}}{x_i^{(k)}}$\n",
    "\n",
    "This algorithm will be applied on the following matrix, after which the calculated eigenvalues with their corresponding eigenvectors will be printed.\n",
    "\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    "5 & -1 & -2 & 0\\\\\n",
    "-1 & 3 & -2 & -1 \\\\\n",
    "-2 & -2 & 5 & 0 \\\\\n",
    "0 & -1 & 0 & 4\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Largest (absolute) eigenvalue: 7.179335040371525, associated eigenvector: [-0.59084255 -0.24322828  0.76543608  0.07650288]\n",
      "Ax = [-4.24185663 -1.74621733  5.49532206  0.5492398 ] = lambda x = [-4.24185663 -1.74621733  5.49532206  0.5492398 ]\n"
     ]
    }
   ],
   "source": [
    "def power_method(matrix: np.ndarray, maximum_iterations: int):\n",
    "    x_k = np.random.rand(matrix.shape[1])\n",
    "    lmbda = 0\n",
    "    for _ in range(maximum_iterations):\n",
    "        z_k1 = np.dot(matrix, x_k)\n",
    "        x_k = z_k1 / np.linalg.norm(z_k1)\n",
    "        i = np.argmax(z_k1)\n",
    "        lmbda = z_k1[i]/x_k[i]\n",
    "\n",
    "    return lmbda, x_k\n",
    "A = np.array([[5, -1, -2, 0], [-1, 3, -2, -1], [-2, -2, 5, 0], [0, -1, 0, 4]])\n",
    "eigenvalue, eigenvector = power_method(A, 100)\n",
    "print(f\"Largest (absolute) eigenvalue: {eigenvalue}, associated eigenvector: {eigenvector}\")\n",
    "print(f\"Ax = {A @ eigenvector} = lambda x = {eigenvalue * eigenvector}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Extend your function, to a new one, such that it takes an extra argument $\\epsilon = 10^{−5}$. Keep track on the number of iterations that you need to reach the precision $|λ(k+1) − λ(k)| ≤ ε$.\n",
    "Here we've used the same function as the previous excercise, with the added ability to specify an argument $\\epsilon$ to dictate the wanted precision. The for loop which contains our algorithm breaks (stops) when the difference between eigenvalue $k-1$ and eigenvalue $k$ is smaller than $\\epsilon$. We also added a new variable ```itr``` to keep track of the iterations.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Largest (absolute) eigenvalue: 7.179327213465006, associated eigenvector: [ 0.59181702  0.24219264 -0.76508741 -0.07574036], iterations: 18\n",
      "Ax = [ 4.24706727  1.74067607 -5.49345634 -0.54515409] = lambda x = [ 4.24884803  1.73878018 -5.49281284 -0.54376486]\n"
     ]
    }
   ],
   "source": [
    "def power_method_limit(matrix: np.ndarray, maximum_iterations: int, epsln: float):\n",
    "    x_k = np.random.rand(matrix.shape[1])\n",
    "    oldlmbda = 0\n",
    "    lmbda = 0\n",
    "    itr = 0\n",
    "    \n",
    "    for _ in range(maximum_iterations):\n",
    "        z_k1 = np.dot(matrix, x_k)\n",
    "        x_k = z_k1 / np.linalg.norm(z_k1)\n",
    "        i = np.argmax(z_k1)\n",
    "        oldlmbda = lmbda\n",
    "        lmbda = z_k1[i]/x_k[i]\n",
    "        itr += 1\n",
    "        if (np.abs(lmbda - oldlmbda) <= epsln):\n",
    "            break\n",
    "\n",
    "    return lmbda, x_k, itr\n",
    "A = np.array([[5, -1, -2, 0], [-1, 3, -2, -1], [-2, -2, 5, 0], [0, -1, 0, 4]])\n",
    "eigenvalue_lim, eigenvector_lim, iterations_lim = power_method_limit(A, 100, 1e-5)\n",
    "print(f\"Largest (absolute) eigenvalue: {eigenvalue_lim}, associated eigenvector: {eigenvector_lim}, iterations: {iterations_lim}\")\n",
    "print(f\"Ax = {A @ eigenvector_lim} = lambda x = {eigenvalue_lim * eigenvector_lim}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Apply your code to find the eigenvalue of $A$ with the smallest absolute value.\n",
    "We've used the inverse power method to find the smallest eigenvalue. For this we first found the inverse of $A$: $A^{-1}$. We then applied the standard power method on $A^{-1}$ to find the smallest eigenvalue and it's corresponding eigenvector.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Smalles (absolute) eigenvalue: 0.6800134632155518, associated eigenvector: [0.40746162 0.71808884 0.52108295 0.21632384], iterations: 6\n",
      "Ax = [0.27705338 0.48831518 0.35431379 0.14720651] = lambda x = [0.27707939 0.48831008 0.35434342 0.14710312]\n"
     ]
    }
   ],
   "source": [
    "def power_method_inverse(matrix: np.ndarray, maximum_iterations: int, epsln: float):\n",
    "    inverse_matrix = np.linalg.inv(matrix)\n",
    "    x_k = np.random.rand(inverse_matrix.shape[1])\n",
    "    oldlmbda = 0\n",
    "    lmbda = 0\n",
    "    itr = 0\n",
    "    \n",
    "    for _ in range(maximum_iterations):\n",
    "        z_k1 = np.dot(inverse_matrix, x_k)\n",
    "        x_k = z_k1 / np.linalg.norm(z_k1)\n",
    "        i = np.argmax(z_k1)\n",
    "        oldlmbda = lmbda\n",
    "        lmbda = z_k1[i]/x_k[i]\n",
    "        itr += 1\n",
    "        if (np.abs(lmbda - oldlmbda) <= epsln):\n",
    "            break\n",
    "    smallest_lbmda = 1/lmbda\n",
    "    return smallest_lbmda, x_k, itr\n",
    "A = np.array([[5, -1, -2, 0], [-1, 3, -2, -1], [-2, -2, 5, 0], [0, -1, 0, 4]])\n",
    "eigenvalue_inv, eigenvector_inv, iterations_inv = power_method_inverse(A, 100, 1e-5)\n",
    "print(f\"Smalles (absolute) eigenvalue: {eigenvalue_inv}, associated eigenvector: {eigenvector_inv}, iterations: {iterations_inv}\")\n",
    "print(f\"Ax = {A @ eigenvector_inv} = lambda x = {eigenvalue_inv * eigenvector_inv}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4. Apply your code to find the eigenvalue of A that is closest to 3.\n",
    "\n",
    "We've used the shifted inverse power method to find the eigenvalue closest to 3. We first shifted the matrix by $s$, which is a new variable that can be specified. We then applied the inverse power method to find the eigenvalue closest to 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Closest eigenvalue to 3: 3.690000971986961\n"
     ]
    }
   ],
   "source": [
    "def power_method_shifted(matrix: np.ndarray, maximum_iterations: int, epsln: float, shift: float):\n",
    "    matrix = matrix - shift\n",
    "    inverse_matrix = np.linalg.inv(matrix)\n",
    "    x_k = np.random.rand(inverse_matrix.shape[1])\n",
    "    oldlmbda = 0\n",
    "    lmbda = 0\n",
    "    itr = 0\n",
    "    \n",
    "    for _ in range(maximum_iterations):\n",
    "        z_k1 = np.dot(inverse_matrix, x_k)\n",
    "        x_k = z_k1 / np.linalg.norm(z_k1)\n",
    "        i = np.argmax(z_k1)\n",
    "        oldlmbda = lmbda\n",
    "        lmbda = z_k1[i]/x_k[i]\n",
    "        itr += 1\n",
    "        if (np.abs(lmbda - oldlmbda) <= epsln):\n",
    "            break\n",
    "    smallest_lbmda = 1/lmbda\n",
    "    return smallest_lbmda, x_k, itr\n",
    "\n",
    "A = np.array([[5, -1, -2, 0], [-1, 3, -2, -1], [-2, -2, 5, 0], [0, -1, 0, 4]])\n",
    "eigenvalue_shift, eigenvector_shift, iterations_shift = power_method_shifted(A, 100, 1e-10, 3)\n",
    "print(f\"Closest eigenvalue to 3: {eigenvalue_shift}\")"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "1c7434731747d88815ea98a4470429bf1ebbf9287f8433fda1e0ecccc300fe0c"
  },
  "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}