A matrix based way to calculate recurrence relations.

matrix_fibonacci.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "wireless-slope",
   "metadata": {},
   "source": [
    "# Fibonacci\n",
    "\n",
    "I came across a new way to calculate fibonacci numbers quickly with matrices. This method has been known for quite some time, but I thought it was an interesting technique to calculate recurrence relations.\n",
    "\n",
    "The fibonacci numbers can be defined with a recurrence relation as follows:\n",
    "\n",
    "$$\n",
    "F_n = \\begin{cases}\n",
    "    F_{n-1} + F_{n-2} &\\text{If } n > 1 \\\\\n",
    "    1 &\\text{If } n=0 \\text{ or } n=1\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "We can represent this relationship in matrix form as well:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{bmatrix}\n",
    "F_n \\\\\n",
    "F_{n-1}\n",
    "\\end{bmatrix}\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "F_{n-1} \\\\\n",
    "F_{n-2}\n",
    "\\end{bmatrix}\\\\\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}^{n-1}\n",
    "\\begin{bmatrix}\n",
    "F_1 \\\\\n",
    "F_0\n",
    "\\end{bmatrix}\\\\\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}^n\n",
    "\\begin{bmatrix}\n",
    "1 \\\\\n",
    "0\n",
    "\\end{bmatrix}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "To change the power from $n-1$ to $n$, we used \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{bmatrix}\n",
    "F_1 \\\\\n",
    "F_0\n",
    "\\end{bmatrix}\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}^{-1}\n",
    "\\begin{bmatrix}\n",
    "F_1 \\\\\n",
    "F_0\n",
    "\\end{bmatrix}\\\\\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "0 &1 \\\\\n",
    "1 &-1\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "F_1 \\\\\n",
    "F_0\n",
    "\\end{bmatrix}\\\\\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "0 &1 \\\\\n",
    "1 &-1\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "1 \\\\\n",
    "1\n",
    "\\end{bmatrix}\\\\\n",
    "&=\n",
    "\\begin{bmatrix}\n",
    "1 &1 \\\\\n",
    "1 &0\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "1 \\\\\n",
    "0\n",
    "\\end{bmatrix}\\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "*The logic we applied to change the power to $n$ could be used to calculate negative terms of the fibonacci sequence too!*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "instructional-helmet",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "fib_matrix = np.array([[1, 1],\n",
    "                       [1, 0]])\n",
    "\n",
    "fib_inv = np.array([[0, 1],\n",
    "                    [1, -1]])\n",
    "\n",
    "def fibonacci(N):\n",
    "    \"\"\"\n",
    "    Calculates the Nth fibonacci number.\n",
    "    \"\"\"\n",
    "    # np.linalg.matrix_power uses the square and multiply method to\n",
    "    # calculate this in O(log N) time and O(1) space.\n",
    "    if N > 0:\n",
    "        matrix = np.linalg.matrix_power(fib_matrix, N)\n",
    "    elif N < 0:\n",
    "        matrix = np.linalg.matrix_power(fib_inv, -N)\n",
    "    else:\n",
    "        matrix = np.identity(2, dtype=int)\n",
    "    \n",
    "    return (matrix @ [1, 0])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "rocky-yemen",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5 -3\n",
      "-4 2\n",
      "-3 -1\n",
      "-2 1\n",
      "-1 0\n",
      "0 1\n",
      "1 1\n",
      "2 2\n",
      "3 3\n",
      "4 5\n",
      "5 8\n"
     ]
    }
   ],
   "source": [
    "for i in range(-5, 6):\n",
    "    print(i, fibonacci(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "south-utilization",
   "metadata": {},
   "source": [
    "## Other Recurrence Relations\n",
    "\n",
    "Something great about this approach is that it translates well across a lot of similar problems. If we can describe the recurrence relation in the form of a matrix, then we can write a simple and fast implementation!\n",
    "\n",
    "### Lucas Numbers\n",
    "\n",
    "The [Lucas Number Sequence](https://en.wikipedia.org/wiki/Lucas_number) is a number sequence very similar to the fibonacci sequence. It has the same recurrence relation, but different initial conditions.\n",
    "\n",
    "$$\n",
    "L_n = \\begin{cases}\n",
    "    L_{n-1} + L_{n-2} &\\text{If } n > 1 \\\\\n",
    "    2 &\\text{If } n=0 \\\\\n",
    "    1 &\\text{If } n=1\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "We can use the same technique for any such sequence, given the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "annual-oriental",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fibonacci_like(N, f_0, f_1):\n",
    "    \"\"\"\n",
    "    Calculates the Nth term of the fibonacci-like sequence:\n",
    "    \n",
    "    f(N) = f(N - 1) + f(N - 2)\n",
    "    f(0) = f_0\n",
    "    f(1) = f_1\n",
    "    \"\"\"\n",
    "    \n",
    "    # This accounts for the offset from before\n",
    "    N -= 1\n",
    "    \n",
    "    if N > 0:\n",
    "        matrix = np.linalg.matrix_power(fib_matrix, N)\n",
    "    elif N < 0:\n",
    "        matrix = np.linalg.matrix_power(fib_inv, -N)\n",
    "    else:\n",
    "        matrix = np.identity(2, dtype=int)\n",
    "    \n",
    "    return (matrix @ [f_1, f_0])[0]\n",
    "\n",
    "def lucas(N):\n",
    "    return fibonacci_like(N, 2, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "chemical-vienna",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5 -11\n",
      "-4 7\n",
      "-3 -4\n",
      "-2 3\n",
      "-1 -1\n",
      "0 2\n",
      "1 1\n",
      "2 3\n",
      "3 4\n",
      "4 7\n",
      "5 11\n"
     ]
    }
   ],
   "source": [
    "for i in range(-5, 6):\n",
    "    print(i, lucas(i))"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}