text: In this post we solve the Fibonacci sequence using linear algebra. The Fibonacci equation is a second-order difference equation which is a particular type of sequence.
  
  ## Sequences
  
  A sequence is a (possibly infinite) set of numbers with a defined order.
  
  $a_n = \frac{1}{n}, \textit{ for } n = 0, 1, 2, ...,$
  
  ### Difference Equations
  
  A difference equation is a type of sequence. A difference equation expresses the value of the sequence relative to the value of another point in the sequence.
  
  $x_{n+1} = 1.5x_n, \textit{ for } n = 0, 1, 2, ...,$
  
  To evaluate this equation at $n=3$ we need to know the initial value $x_0$, and evaluate $n=1, 2, 3$ in succession. For example, if $x_0=100$:
  
  \begin{align} x_0 & && = 100 \\ x_1 & = 1.5 \times 100 && = 150\\ x_2 & = 1.5 \times 150 && = 225\\ x_3 & = 1.5 \times 225 && = 337.5 \end{align}
  
  In the example above, we evaluate the equation recursively. It would be better to be able to evaluate the equation at $n = 100$ without performing $100$ calculations. Note that we can also write the equation above in this form:
  
  $x_3 = 100 \times 1.5 \times 1.5 \times 1.5 = 337.5$
  
  Solving a difference equation means writing the same equation such that it can be evaluated at a given $n$ without first evaluating $n=1,2,etc$. The solution to the equation above is:
  
  $x_{n} = 1.5^{n} \times x_0, \textit{ for } n = 0, 1, 2, ...,$
  
  ### Second Order Difference Equations
  
  The previous example is a first order difference equation because it expresses the value of $x_n$ as a function of one other sample, $x_{n-1}$. A second order difference equation is a function of two other samples, for example
  
  $x_{n+2} = 1.5x_n + 1.4x_{n+1}, \textit{ for } n = 0, 1, 2, ...,$
  
  ## Fibonacci Numbers
  
  The easiest way to describe the Fibonacci sequence is with a second order difference equation. Assume $F_0 = 0$ and $F_1 = 1$
  
  $F_{k+2} = F_{k+1} + F_{k}$
  
  The first eight Fibonacci numbers are:
  
  $F_k=0,1,1,2,3,5,8,13,...$
  
  We can solve the sequence with linear algebra. Let us start our sequence with a vector, $u_0=\begin{bmatrix}1\\0\end{bmatrix}$. Find the next vector in our sequence, $u_1$ like this:
  
  $u_{1} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$
  
  More generally:
  
  $u_{k+1} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} u_k$
  
  Note that $u_0$ contains $\begin{bmatrix}F_1\\F_0\end{bmatrix}$.
  
  The key to solving the sequence is in the eigenvalues and eigenvectors of $\begin{bmatrix}1&1\\1&0\end{bmatrix}$. We can write $u_3$ as
  
  $u_3 = \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix} \begin{bmatrix}1 \\ 0\end{bmatrix} = \begin{bmatrix}3 \\ 2\end{bmatrix} = \begin{bmatrix}F_4 \\ F_3\end{bmatrix}$
  
  import numpy as np
  fib = np.matrix([[1,1],[1,0]])
  u_0 = np.matrix([[1],[0]])
  
  fib * fib * fib * u_0 # The * operator is dot product for np.matrix (but not np.array)
  matrix([[3],
  [2]])
  
  ### Eigenvector Solution
  
  The general case of the previous example is:
  
  \begin{align} u_k &= \begin{bmatrix}1 & 1 \\ 1 & 0\end{bmatrix}^k\begin{bmatrix}1 \\ 0\end{bmatrix}\\ u_k &= A^ku_0 \end{align}
  
  We can use eigenvectors and eigenvalues to solve the matrix power (see Eigenvalues and Eigenvectors). First, find the eigenvalues of $A$:
  
  \begin{align} det\begin{vmatrix} 1 - \lambda & 1\\ 1 & - \lambda \end{vmatrix} & = 0\\ -\lambda(1-\lambda) - 1 & = 0\\ \lambda{}^2 - \lambda - 1 & = 0 \end{align}
  
  Solve using the quadratic formula $\frac{-b \pm \sqrt{b^2-4ac}}{2a} \implies \frac{1 \pm \sqrt{5}}{2} \approx \frac{1 \pm 2.236}{2}\\$, which gives the eigenvalues for our matrix: \begin{align} \lambda_1 & \approx &1.618 \\ \lambda_2 & \approx &-0.618 \end{align}
  
  Now we plug our constant eigenvalues ($\lambda_1$ and $\lambda_2$) into $(A-\lambda I)\vec{x}=0$ to find the eigen vectors.
  
  $\begin{bmatrix}1-\lambda_1 & 1 \\ 1 & -\lambda_1 \end{bmatrix}\vec{x}_1 = 0 \implies \vec{x}_1 = \begin{bmatrix}\lambda_1 \\ 1\end{bmatrix}$
  
  Note that this solution is not necessarily intuitive. This is not a valid solution for any arbitrary value of $\lambda$: our values, $\lambda_1$ and $\lambda_2$ are valid solutions. We can use the same process to solve for $\lambda_2$, revealing eigenvectors $\begin{bmatrix}\lambda_1 \\1\end{bmatrix}$ and $\begin{bmatrix}\lambda_2 \\1\end{bmatrix}$. Now create the eigenvector matrix, $S$ and eigenvalue matrix, $\Lambda$.
  
  \begin{align} S &= \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix}\\ \Lambda &= \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{bmatrix} \end{align}
  
  Remeber that $A^k = S^{-1}\Lambda^kS$. To proceed, invert $S$. Remember this shortcut for inverting a $2 \times 2$ matrix.
  
  \begin{align} \begin{bmatrix} a & b \\ c & d \end{bmatrix}^{-1} =&& \frac{1}{ad-cb} &\begin{bmatrix} d & -b \\ -c & a \end{bmatrix} \\ S^{-1} =&& &\begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix}^{-1} \\ S^{-1} =&& \frac{1}{\lambda_1-\lambda_2} &\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix}\\ S^{-1} =&& \tfrac{1}{\sqrt{5}} &\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix} \end{align}
  
  Using the components $(A, \Lambda, S, S^{-1},u_0)$:
  
  \begin{align} u_k &= A^ku_0 \\ u_k &= S\Lambda^kS^{-1}u_o \\ u_k &= S\Lambda^k\frac{1}{\sqrt{5}}\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix}u_0 \\ u_k &= \frac{1}{\sqrt{5}}S\Lambda^k\begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix}u_0 \\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1 & 0 \\ 0 & \lambda_2\end{bmatrix}^k \begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix} u_0 \\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1^k & 0 \\ 0 & \lambda_2^k\end{bmatrix} \begin{bmatrix}1 & -\lambda_2\\ -1 & \lambda_1 \end{bmatrix} \begin{bmatrix}1\\0\end{bmatrix} \\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1^k & 0 \\ 0 & \lambda_2^k\end{bmatrix} \begin{bmatrix}1 \\ -1 \end{bmatrix}\\ u_k &= \tfrac{1}{\sqrt{5}} \begin{bmatrix}\lambda_1 & \lambda_2 \\ 1 & 1\end{bmatrix} \begin{bmatrix}\lambda_1^k \\ -\lambda_2^k\end{bmatrix} \end{align}
  
  ## Verify Solution in Python
  
  import numpy as np
  
  lmbd_1 = (1+np.sqrt(5))/2
  lmbd_2 = (1-np.sqrt(5))/2
  
  a       = np.matrix([[1,      1],      [1,  0]])
  s       = np.matrix([[lmbd_1, lmbd_2], [1,  1]])
  s_inv   = np.matrix([[1,     -lmbd_2], [-1, lmbd_1]]) * (1/np.sqrt(5))
  eig_val = np.matrix([[lmbd_1, 0],      [0,  lmbd_2]])
  u_0     = np.matrix([[1],              [0]])
  
  def eig_k(k):
  return eig_val ** k
  
  def u_k(k):
  return s * eig_k(k) * s_inv * u_0
  
  def fib(n):
  f_n = u_k(n)[1,0]
  return int(round(f_n))
  
  [fib(n) for n in range(15)]
  [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377]