Constant-time Fibonacci

This is the second part in a 2-part series on the “Fibonacci” interview problem. We are building off of a previous post, so take a look at Part I if you haven’t seen it.

Previously, we examined the problem and constructed a logarithmic-time solution based on computing the power of a matrix. Now we will derive a constant time solution using some more linear algebra. If you had trouble with the linear algebra in part I, it may help to read up on matrices, matrix multiplicaiton, and special matrix operations (specifically determinants and inverses) before moving on.

It’s time to pull out the eigenvalues.

Eigenvalues and Eigenvectors

Eigenvalues and eigenvectors are a bit of a magical concept in linear algebra. The definition of eigenvalues is: “A matrix times an eigenvector equals an eigenvalue times an eigenvector.” As an equation, the eigenvalues ($\lambda$) and eigenvectors ($\textbf{v}$) of a matrix ($\textbf{A}$) can be computed by solving:

$$ \textbf{A}\textbf{v} = \textbf{v}\lambda $$

In more intuitive terms, each matrix has a set of vectors (its eigenvectors) which are scaled by a constant factor (the eigenvalues) when multiplied by the matrix. If we think of matrix multiplication as a transformation of a shape in space, the eigenvectors are the directions along which the matrix stretches the shape, and the eigenvalues are the scaling factors.

A square $m \times m$ matrix can have up to $m$ eigenvalues. In the case of recurrence relations, the matrices we are manipulating should always have $m$ eigenvalues. As long as the matrix has enough eigenvalues, we can make a matrix out of all of the eigenvectors, which we will call $\textbf{Q}$, and we can make a matrix which has the eigenvalues along the diagonal, $\Lambda$. Now our equation looks like:

$$ \textbf{A}\textbf{Q} = \textbf{Q}\Lambda $$

We can multiply both sides by the inverse of the $\textbf{Q}$ matrix, and we will get:

$$ \textbf{A} = \textbf{Q}\Lambda\textbf{Q}^{-1} $$

Now we can do something interesting:

$$ \textbf{A}^2 = \textbf{A} \cdot \textbf{A} = \textbf{Q}\Lambda\textbf{Q}^{-1} \cdot \textbf{Q}\Lambda\textbf{Q}^{-1} = \textbf{Q}\Lambda^2\textbf{Q}^{-1} $$

The multiplication of $\textbf{Q}^{-1}$ and $\textbf{Q}$ in the middle allows us to cancel them out, so we can reduce the problem of squaring a matrix to the problem of squaring its eigenvalues. There is no need for us to stop at squaring. If we multiply by $\textbf{A}$ again, we get another pair of $\textbf{Q}^{-1}$ and $\textbf{Q}$ that cancel:

$$ \textbf{A}^3 = \textbf{A}^2 \cdot \textbf{A} = \textbf{Q}\Lambda^2\textbf{Q}^{-1} \cdot \textbf{Q}\Lambda\textbf{Q}^{-1} = \textbf{Q}\Lambda^3\textbf{Q}^{-1} $$

More generally, for every power of the matrix, we can compute it by taking a power of its eigenvalue matrix:

$$ \textbf{A}^n = \textbf{Q} \Lambda^n \textbf{Q}^{-1} $$

The “stretching” idea hints at this: stretching a shape multiple times using the same transformation should stretch it more along the same axes. However, formula won’t apply unless $\textbf{A}$ (which is an $m \times m$ matrix) has exactly $m$ eigenvalues.

At first glance, it’s not exactly clear why this helps: we have reduced computing the power of the $\textbf{A}$ matrix to the power of a different matrix. However, $\Lambda$ is a diagonal matrix, so computing its power is equivalent to computing the power of each of the elements along the diagonal:

$$ \Lambda^n = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}^n = \begin{bmatrix} \lambda_1^n & 0 \\ 0 & \lambda_2^n \end{bmatrix} $$

Through eigendecomposition, we have reduced the problem of an $m \times m$ matrix power to computing $m$ scalar powers. The scalar powers can be computed in constant time (using the venerable pow function), whereas the matrix power cannot otherwise be computed in constant time. Now we just have to apply this to the Fibonacci problem.

Computing our Eigenvalues and Eigenvectors

A nice interviewer might let you use Wolfram Alpha to do an eigendecomposition, but just in case, let’s do it by hand. We determine the eigenvalues using the characteristic polynomial of the matrix. The polynomial is:

$$ p(\lambda) = \det (\textbf{A} - \lambda \textbf{I}) $$

Then, for each eigenvalue ($\lambda_i$), we solve the eigenvalue equation to get the corresponding eigenvector ($\textbf{v}_i$):

$$ (\textbf{A} - \lambda_i \textbf{I}) \textbf{v}_i = 0 $$

Let’s do it for the Fibonacci recurrence relation. We will start with the eigenvalues:

$$ 0 = \begin{vmatrix} 1 - \lambda & 1 \\ 1 & -\lambda \end{vmatrix} = -\lambda (1 - \lambda) - 1 = \lambda^2 - \lambda - 1 $$

Solving this with the quadratic formula gives:

$$ \lambda = \frac{1 \pm \sqrt{5}}{2} $$

That’s cool. That’s the golden ratio (and its conjugate)! Our eigenvalues are the golden ratio and its conjugate. We will use $\phi = \frac{1 + \sqrt{5}}{2}$ for the golden ratio from now on, and $\psi = \frac{1 - \sqrt{5}}{2}$ for its conjugate. Also, there is an interesting property of the golden ratio that helps us simplify the algebra:

$$ -\frac{1}{\phi} = 1 - \phi = \psi $$

Let’s solve for the eigenvectors, starting with the eigenvector for $\phi$:

$$ \begin{bmatrix} 1 - \phi & 1 \\ 1 & -\phi \end{bmatrix}\begin{bmatrix} x \\ y \end{bmatrix} = 0 $$

This simplifies, using some properties of the golden ratio, to the equations:

$$ 0 = \psi x + y $$ $$ 0 = x + \phi y $$

So we can use $x = \phi$ and $y = 1$ for this eigenvector. Note that any constant factor of these numbers also solves this set of equations, so $x = 2\phi$ and $y = 2$ also works. Conventionally, eigenvectors are usually scaled so that their magnitude is 1, but we can use any scaling factor that makes the math easy. With that one out of the way, let’s do the other eigenvector (for the eigenvalue $\psi$):

$$ \begin{bmatrix} 1 - \psi & 1 \\ 1 & -\psi \end{bmatrix}\begin{bmatrix} x \\ y \end{bmatrix} = 0 $$

Which gives $x = \psi$ and $y = 1$ as its solution. The resulting $\textbf{Q}$ matrix, which has the eigenvectors as its columns, is:

$$ \textbf{Q} = \begin{bmatrix} \phi & \psi \\ 1 & 1 \end{bmatrix} $$

Finishing Constant-time Fibonacci

So we have $\textbf{Q}$ and $\Lambda$, and now we just need the inverse matrix of $\textbf{Q}$. There is a pretty straightforward formula for the inverse of a $2 \times 2$ matrix, and we will apply it to get the inverse matrix: negate the off-diagonal elements (bottom left and top right), swap the diagonal elements, and divide by the determinant. For a $3 \times 3$ matrix, you are better off using a computer than trying to memorize a formula. Fortunately, linear algebra packages (including LAPACK and NumPy) can do this. The inverse of the eigenvector matrix is:

$$ \textbf{Q}^{-1} = \frac{1}{\phi - \psi}\begin{bmatrix} 1 & -\psi \\ -1 & \phi \end{bmatrix} = \frac{1}{\sqrt{5}}\begin{bmatrix} 1 & -\psi \\ -1 & \phi \end{bmatrix} $$

Putting it all together, we have our formula for $\textbf{A}^n$, based on constant-time calculations:

$$ \textbf{A}^n = \frac{1}{\sqrt{5}} \begin{bmatrix} \phi & \psi \\ 1 & 1 \end{bmatrix} \begin{bmatrix} \phi^n & 0 \\ 0 & \psi^n \end{bmatrix} \begin{bmatrix} 1 & -\psi \\ -1 & \phi \end{bmatrix} $$

$$ = \frac{1}{\sqrt{5}}\begin{bmatrix} \phi^{n+1}-\psi^{n+1} & -\phi^{n + 1}\psi + \psi^{n+1}\phi \\ \phi^n-\psi^n & -\phi^n\psi + \psi^n\phi \end{bmatrix} $$

Simplifying based on our earlier identity that $\psi = -\frac{1}{\phi}$:

$$ A^n = \frac{1}{\sqrt{5}}\begin{bmatrix} \phi^{n+1} - \psi^{n+1} & \phi^n - \psi^{n} \\ \phi^n - \psi^{n} & \phi^{n-1} - \psi^{n-1} \end{bmatrix} $$

The last step is to multiply by the initial conditions:

$$ \begin{bmatrix} F_{n+1} \\ F_{n} \end{bmatrix} = \frac{1}{\sqrt{5}}\begin{bmatrix} \phi^{n+1} - \psi^{n+1} & \phi^n - \psi^{n} \\ \phi^n - \psi^{n} & \phi^{n-1} - \psi^{n-1} \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} $$

Simplifying, and only taking the bottom row since that is the formula for $F_n$, we get:

$$ F_n = \frac{\phi^n - \psi^n}{\sqrt{5}} $$

For the Fibonacci recurrence relation, his formula is fairly well known under the name “Binet’s formula.”

Writing the Code

In theory, if this is an interview, we are going to be expected to produce some code. We are now going to pretend that we don’t know the formula we just derived in the previous section, and we are going to write a more generic solution rather than attempting to re-derive the formula during a high-pressure situation. I don’t know about you, but I’m pretty bad at doing algebra by hand, so it would be better to have the comptuer do everything. With that in mind, let’s do this in Python this time–Python is generally not my choice of language, but NumPy is very user-friendly and I would recommend it for linear algebra like this. Honestly, if I had to use C or C++ (without LAPACK), I would probably do the logarithmic-time method from part I rather than risking the algebra to derive a closed-form $O(1)$ solution.

import numpy as np
from numpy import linalg as LA

def constant_fibonacci(n):
    # Set up the A matrix and the initial conditions
    a = np.array([[1, 1],
                  [1, 0]])
    initial_conditions = np.array([1, 0])

    # Do eigendecomposition and invert the eigenvector matrix
    eigenvalues, eigenvectors = LA.eig(a)
    vector_inverse = np.inv(eigenvectors)

    # Compute the power of the eigenvalues, done elementwise
    lambda_pow = np.diag(eigenvalues ** (n - 1))

    # Put it all together - compute A^n and apply initial conditions
    a_pow = eigenvectors @ lambda_pow @ vector_inverse
    return round((a_pow @ initial_conditions)[0])

Note the rounding here at the end: we are working with transcendental numbers as intermediate values, so we cannot avoid predominantly using floating-point math for this form of solution. Additionally, eigenvalues and eigenvectors may end up being computed using numerical methods, and they may not be exact. Even with rounding, the use of floats will cause divergence from the acutal Fibonacci numbers when $Fib(n)$ is big: Double-precision floating point numbers have 53 bits of precision, and once you are past $2^{53}$, they cannot represent integers exactly. Unsigned ints have the same problem at $2^64$. This happens pretty fast: the Fibonacci numbers grow exponentially, so the numpy code above (which uses doubles as its underlying numeric type) has its first error at n = 71. Using a bignum or arbitrary precision numeric library will bring this solution back to at least $O(log n)$, if not $O(n)$, depending on the algorithms used for exponentiation and dvision.

In response to my last post on reddit, a number of commenters suggested that you just use NumPy’s matrix power function. This is how that matrix power calculation works.

Conclusions

Over the last two posts, we have discussed how to do better than O(n) on Fibonacci-style problems. We have an easier logarithmic solution and a tricky constant-time solution straight out of a differential equations course.

The Fibonacci problem is a nice example of the improvements you can get from recursion with memoization and one-dimensional dynamic programming, but can’t we find another nice example where there isn’t a better solution than recursion? Teaching students that they should be using recursion to find Fibonacci numbers teaches them to develop a blind spot - instead, they should be learning to examine the problem carefully for a better solution rather than reaching into a bag of algorithmic tricks.

Subscribe for Email Notifications