{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arnoldi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defined as Algorithm 2 in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Arnoldi(A,v0,m):\n",
    "    v = v0\n",
    "    beta = sp.linalg.norm(v)\n",
    "    v = v/beta\n",
    "    H = sp.sparse.lil_matrix((m+1,m)) \n",
    "    V = sp.sparse.lil_matrix((A.shape[0],m+1))\n",
    "    V[:,0] = v # each column of V is a vector v\n",
    "\n",
    "    for j in range(m):\n",
    "        w = A @ v  \n",
    "        for i in range(j):\n",
    "            H[i,j] = v.T @ w # tmp is a 1x1 matrix, so it's O(1) in memory\n",
    "            w = w - H[i,j]*v \n",
    "            \n",
    "        H[j+1,j] = np.linalg.norm(w)\n",
    "\n",
    "        if H[j+1,j] == 0:\n",
    "            print(\"Arnoldi breakdown\")\n",
    "            m = j\n",
    "            v = 0\n",
    "            break\n",
    "        else:\n",
    "            if j < m-1:\n",
    "                v = w/H[j+1,j]\n",
    "                V[:,j+1] = v\n",
    "\n",
    "    return V, H, beta, j  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Small test case\n",
    "\n",
    "The final implementation will be using all sparse arrays and matrices, no numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 110\n",
    "\n",
    "# start with a random sparse matrix\n",
    "A = sp.sparse.rand(n,n, density=0.1, format='lil')\n",
    "\n",
    "# Starting vector\n",
    "v = np.repeat(1/n,n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the Arnoldi Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of iterations is:  109\n",
      "The matrix H is:  (111, 110)\n",
      "The matrix V is:  (110, 111)\n"
     ]
    }
   ],
   "source": [
    "V, H, beta, j = Arnoldi(A,v,110)\n",
    "print(\"The number of iterations is: \", j)\n",
    "print(\"The matrix H is: \", H.shape)\n",
    "print(\"The matrix V is: \", V.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, it returns $H_{m+1}$ that is an upper-hassemberg matrix, with one extra row."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.8 64-bit",
   "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.10.9"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}