Home Fibonacci with linear algebra
Post
Cancel

Fibonacci with linear algebra

Fibonacci sequence

The Fibonacci sequence is one of the few sequences in mathematics that are widely known among people outside of this field. Let’s remind how it works: You start with 0 and 1, then add the two previous terms to generate the next one.

$$0 \hspace{1em} 1 \hspace{1em} 1 \hspace{1em} 2 \hspace{1em} 3 \hspace{1em} 5 \hspace{1em} 8 \hspace{1em} 13 \hspace{1em} 21 \hspace{1em} 34\dots$$

If we define \(f_n\) as the \(n\)-th element in this sequence, the following holds:

$$ \begin{cases} \begin{align*} f_0 &= 0 \\ f_1 &= 1 \\ f_{n} &= f_{n-1} + f_{n-2}, \hspace{1em} \text{for all } n \ge 2 \end{align*} \end{cases} $$

We can use this definition to compute any Fibonacci number, but it would take \(\mathcal{O}(n)\) additions. What if we want the \(10^{16}\)-th Fibonacci number? (of course, it’s too large, but maybe we want to compute some of the last digits, or use some extended but finite precision float numbers for the computation). I don’t think we want to wait some months until the computation is done. Let’s try to find a different approach.

Linear recurrences are related to linear algebra, as usually happens with linear things. Let’s use a linear algebra approach.

Matrix form

Since every new term requires the last two terms, let’s define the \(n\)-th Fibonacci vector as:

\[F_n = \begin{pmatrix} f_{n} \\\ f_{n-1} \end{pmatrix}\]

In the regular sequence form, we were using \((f_{n-1}, f_{n-2})\) to get \(f_{n}.\) It’s reasonable to think that now we can get \(F_n\) from \(F_{n-1} \sim (f_{n-1}, f_{n-2}).\) Let’s try to find a linear application \(M\) such that \(F_n = M F_{n-1}.\)

\[F_n = M F_{n-1} \iff \begin{pmatrix}f_{n} \\\ f_{n-1}\end{pmatrix} = \begin{pmatrix}a & b \\\ c & d\end{pmatrix} \begin{pmatrix}f_{n-1} \\\ f_{n-2}\end{pmatrix} = \begin{pmatrix}a f_{n-1} + b f_{n-2} \\\ c f_{n-1} + d f_{n-2}\end{pmatrix}\]

Let’s find the constants \(a,b,c,d\) that satisfy the equation above. Starting with the bottom row, it’s pretty clear that \(d = 0\) and \(c = 1.\) Then looking at the top row, we use the definition of \(f_n\) and get \(a = b = 1.\) There are no other possible constants that satisfy that equation, so we found the unique linear application that gives us the next Fibonacci vector.

\[\begin{pmatrix}f_{n} \\\ f_{n-1}\end{pmatrix} = \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \begin{pmatrix}f_{n-1} \\\ f_{n-2}\end{pmatrix}\]

Moreover, we can repeat substitution with this equation to get an interesting result:

\[\begin{align*} \begin{pmatrix}f_{n} \\\ f_{n-1}\end{pmatrix} &= \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \begin{pmatrix}f_{n-1} \\\ f_{n-2}\end{pmatrix} = \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \left[ \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \begin{pmatrix}f_{n-2} \\\ f_{n-3}\end{pmatrix} \right] \\ &= \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \left[ \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \left\{ \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix} \begin{pmatrix}f_{n-3} \\\ f_{n-4}\end{pmatrix} \right\} \right] = \\ &= \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix}^3\begin{pmatrix}f_{n-3} \\\ f_{n-4}\end{pmatrix} = \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix}^4\begin{pmatrix}f_{n-4} \\\ f_{n-5}\end{pmatrix} = \dots\\ \dots &= \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix}^{n-1}\begin{pmatrix}f_{1} \\\ f_{0}\end{pmatrix} = \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix}^{n-1}\begin{pmatrix}1 \\\ 0\end{pmatrix} \end{align*}\]

This is a matrix closed form for the \(n\)-th Fibonacci number, which may be written in the following nicer form:

The matricial expression

\[\begin{pmatrix}f_{n+1} \\\ f_{n}\end{pmatrix} = \begin{pmatrix}1 & 1 \\\ 1 & 0\end{pmatrix}^{n}\begin{pmatrix}1 \\\ 0\end{pmatrix}\]

Is this a real improvement somehow? If the exponent is computed in a naive way, it will cost \(\mathcal{O}(n)\) additions and multiplications, but we can use fast exponentiation, which is simple and computes the result in \(\mathcal{O}(\log n).\) What an improvement! Thank the linear algebra.

But linear algebra can do more for us! Now we can also find a closed-form expression for \(f_n.\) We will use the power of the Jordan form! The idea is that we decompose \(M\) in a base that makes \(M\) as diagonal as possible. Since \(M\) is a real symmetric matrix, the Spectral Theorem ensures we can diagonalize it! And the Jordan form will do the job.

We’ll find matrices \(J, P\) such that \(J\) is diagonal and \(M = PJP^{-1}.\) Note that now exponentiation is pretty easy since \(M^n = P J^n P^{-1},\) and \(J^n\) can be computed by simply exponentiating the diagonal elements to the \(n\)-th power.

Jordan form using SymPy

Since I am lazy, I will use SymPy, a Python symbolic algebra library. I run it on a Jupyter notebook:

1
2
3
4
5
6
from sympy import *
M = Matrix([
  [1,1],
  [1,0]
])
M
\[\left[\begin{matrix}1 & 1\\1 & 0\end{matrix}\right]\]
1
2
3
P, J = M.jordan_form()  # SymPy does all the Jordan work!
P_ = simplify(P.inv())  # P_ = P^-1
simplify(J)
\[\left[\begin{matrix}\frac{1}{2} - \frac{\sqrt{5}}{2} & 0\\0 & \frac{1}{2} + \frac{\sqrt{5}}{2}\end{matrix}\right]\]

\(J\) it’s the diagonal matrix of the eigenvalues of \(M,\) which are the roots of the characteristic polynomial of \(M\):

\[0 = \det(M-\lambda I) = \begin{vmatrix}1-\lambda & 1 \\\ 1 & -\lambda\end{vmatrix} = (1-\lambda)(-\lambda) - 1 = \lambda^2 - \lambda - 1\]

The solutions are:

\[\lambda_{-} = \frac{1-\sqrt{5}}{2}, \hspace{1em} \lambda_{+} = \frac{1+\sqrt{5}}{2}\]
1
simplify(P)
\[\left[\begin{matrix}\frac{1}{2} - \frac{\sqrt{5}}{2} & \frac{1}{2} + \frac{\sqrt{5}}{2}\\1 & 1\end{matrix}\right]\]
1
simplify(P_)
\[\left[\begin{matrix}- \frac{\sqrt{5}}{5} & \frac{\sqrt{5}}{10} + \frac{1}{2}\\\frac{\sqrt{5}}{5} & \frac{1}{2} - \frac{\sqrt{5}}{10}\end{matrix}\right]\]
1
2
3
4
5
6
7
8
F_1 = Matrix([
  [1],
  [0]
])
n = Symbol('n')
M_power = P @ J**(n-1) @ P_  # @ is the matrix product operator, M_power = M^(n-1)
F_n = M_power @ F_1
simplify(F_n)
\[\left[\begin{matrix}\frac{2^{- n} \sqrt{5} \left(- \left(1 - \sqrt{5}\right)^{n} + \left(1 + \sqrt{5}\right)^{n}\right)}{5}\\\frac{2^{1 - n} \sqrt{5} \left(- \left(1 - \sqrt{5}\right)^{n - 1} + \left(1 + \sqrt{5}\right)^{n - 1}\right)}{5}\end{matrix}\right]\]

And that’s it! All computations are doable by hand, but using SymPy we get in a few lines of code our desired result:

\[F_n = \begin{pmatrix}f_{n} \\\ f_{n-1}\end{pmatrix} = \left(\begin{matrix}\frac{2^{- n} \sqrt{5} \left(- \left(1 - \sqrt{5}\right)^{n} + \left(1 + \sqrt{5}\right)^{n}\right)}{5}\\\frac{2^{1 - n} \sqrt{5} \left(- \left(1 - \sqrt{5}\right)^{n - 1} + \left(1 + \sqrt{5}\right)^{n - 1}\right)}{5}\end{matrix}\right)\]

The closed-form expression

Comparing the first entries of the previous vectors we have:

\[\begin{align*} f_{n} &= \frac{2^{- n} \sqrt{5} \left[- \left(1 - \sqrt{5}\right)^{n} + \left(1 + \sqrt{5}\right)^{n}\right]}{5} \\ &= \frac{2^{-n}}{\sqrt{5}} \left[ \left(1 + \sqrt{5}\right)^{n} - \left(1 - \sqrt{5}\right)^{n} \right] \\ &= \frac{1}{\sqrt{5}} \left[ {\left( \frac{1+\sqrt{5}}{2} \right)}^n - {\left( \frac{1-\sqrt{5}}{2} \right)}^n \right] \end{align*}\]

Linear algebra did an amazing job here, allowing us to find both an efficient computation method and a closed-form formula by (almost) the same price :)

The Golden ratio

The Golden ratio is the proportion \(\varphi\) of any two lengths \(a, b \in \mathbb{R}^{+}\) satisfying:

\[\varphi := \frac{a}{b} = \frac{a+b}{a}\]

Notice that \(\varphi\) is uniquely defined since \(\varphi \in \mathbb{R}^{+}\):

\[\frac{a}{b} = \frac{a+b}{a} \iff \frac{a}{b} = 1 + \frac{b}{a} \iff \varphi = 1 + \frac{1}{\varphi} \iff \varphi^2 = \varphi + 1\]

This quadratic equation has only one positive solution:

\[\varphi = \frac{1+\sqrt{5}}{2}\]

Notice that the Golden ratio is defined as the positive solution to \(\varphi^2 = \varphi + 1 \iff \varphi^2 - \varphi - 1 = 0.\)
Sounds familiar? The characteristic polynomial of \(M\) is \(\lambda^2 - \lambda - 1\), which is the same expression as the Golden ratio. This is why the Golden ratio appears in the closed-form expression if we substitute \(\varphi\):

\[f_n = \frac{\varphi^n - (1-\varphi)^n}{\sqrt{5}}\]

Now that we have shown the relation between Fibonacci and the Golden ratio, we can make an interesting observation:

Notice that \(\varphi \approx 1.618,\) and \((1-\varphi) \approx -0.618.\) In the previous expression, the dominating term is \(\varphi^n,\) and \((1-\varphi)^n\) will become closer to \(0\) as \(n\) increases. This means that for big enough \(n\), \(\ f_n \approx \varphi^n / \sqrt{5}\).

The ratio of consecutive Fibonacci numbers is the Golden ratio!

We have seen that \(\ f_n \approx \varphi^n / \sqrt{5}.\) Intuitively, this implies that the ratio \(f_{n+1}/f_n \to \varphi\) as \(n\) increases. Let’s prove it by computing \(\lim_{n \to \infty} f_{n+1}/f_n\)

\[\begin{align*} \lim_{n \to \infty} \frac{f_{n+1}}{f_n} &= \lim_{n \to \infty} \frac{\varphi^{n+1} - (1-\varphi)^{n+1}}{\sqrt{5}} \bigg/ \frac{\varphi^n - (1-\varphi)^n}{\sqrt{5}} = \lim_{n \to \infty} \frac{\varphi^{n+1} - (1-\varphi)^{n+1}}{\varphi^{n} - (1-\varphi)^{n}} = \\ &= \lim_{n \to \infty} \frac{\varphi^n \left[\varphi - \frac{(1-\varphi)^{n+1}}{\varphi^n}\right]}{\varphi^n \left[ 1 - \frac{(1-\varphi)^{n}}{\varphi^n}\right] } = \lim_{n \to \infty} \frac{\varphi - \frac{(1-\varphi)^{n+1}}{\varphi^n}}{1 - \frac{(1-\varphi)^{n}}{\varphi^n}} = \frac{\lim_{n \to \infty} \varphi - \frac{(1-\varphi)^{n+1}}{\varphi^n}}{\lim_{n \to \infty} 1 - \frac{(1-\varphi)^{n}}{\varphi^n}} \end{align*}\]

But \(\|1-\varphi\| < 1 \implies (1-\varphi)^n \to 0,\) and \(\varphi^n \to \infty,\) so we have \(\lim_{n \to \infty} (1-\varphi)^n/\varphi^n = 0\). Finally, our desired result:

\[\lim_{n \to \infty} \frac{f_{n+1}}{f_n} = \varphi\]
This post is licensed under CC BY 4.0 by the author.