Turn the Fibonacci rule into one matrix step on a pair of numbers, diagonalise the Fibonacci matrix, derive Binet’s formula from the eigen-expansion, and find the exact place where the formula fails in floating point.
M03.2.3L2 Intermediate25 minIn review (what this means)
The maths and code on this page are checked by automated tests on every change, but an independent reader has not signed it off yet. Spot a mistake or a confusing step? Tell us.
Rewrite a second-order linear recurrence as a first-order matrix system u_(k+1) = F u_k.
Find the eigenvalues phi and -1/phi and the eigenvectors of the Fibonacci matrix, and derive Binet’s formula from the eigen-expansion of the start.
Explain why F_(k+1)/F_k tends to phi and why its error shrinks by a factor 1/phi^2 each step.
Compute Fibonacci numbers exactly with integer matrix powers and locate where the floating-point formula first fails.
The question
The Fibonacci numbers \(0, 1, 1, 2, 3, 5, 8, 13, \dots\) obey one rule: each is the sum of the two before it. Is there a formula for the hundredth one that does not walk through the first ninety-nine? The answer comes from a 2 by 2 matrix and its eigenvalues, and the same trick analyses every update rule that remembers two steps, including gradient descent with momentum.
From the last lesson
1. If \(\vect{x} = c_1\vect{v}_1 + c_2\vect{v}_2\) in eigenvectors, what is \(\mat{A}^k\vect{x}\)?
Show answer
\(c_1\lambda_1^k\vect{v}_1 + c_2\lambda_2^k\vect{v}_2\): each eigen-coordinate is multiplied by its own eigenvalue \(k\) times.
2. What do the eigenvalues say about the long run of \(\mat{A}^k\vect{x}\)?
Show answer
All of size below 1: decay to zero. One above 1 (with a non-zero coefficient): growth along its eigenvector. Largest exactly 1 and the others smaller: a steady state along the eigenvector of 1.
Intuition
The rule \(F_{k+2} = {F_{k+1} + F_k}\) needs two numbers to take a step. So keep two numbers: the state\(\vect{u}_k = {(F_{k+1}, F_k)}\), the newest Fibonacci number on top and the one before it underneath. One step of the rule turns \(\vect{u}_k\) into \(\vect{u}_{k+1}\), and that step is a linear map of the plane: a matrix.
Plot the states as points. They start at \(\vect{u}_0 = (1, 0)\) and run out through \((1, 1), (2, 1), (3, 2), (5, 3), \dots\) They do not wander. They hug one line, ever more tightly, and along that line the ratio of the two coordinates is a fixed number, about 1.618.
M03.2.3-F1
The Fibonacci states line up along an eigen-line of the Fibonacci matrix, and the ratio of neighbours zigzags onto the golden ratio \(\varphi\).
That line is an eigen-line of the step matrix, and 1.618 is the golden ratio \(\varphi\). Lesson M03.2.2 explains why: after many steps only the term with the largest eigenvalue matters.
Build it up
Step 1. Write the step for both entries of the state. The new top entry follows the rule, and the new bottom entry is the old top entry, copied down:
Step 2. Read off the matrix. The coefficients, row by row, are the Fibonacci matrix:
M03.2.3-F2
One step of the recurrence as a matrix. The first row adds the two entries; the second row only copies, so that the next step has the two numbers it needs.
Step 3. Many steps are a power. From \(\vect{u}_0 = (F_1, F_0) = (1, 0)\), \(\vect{u}_k = \mat{F}^k\vect{u}_0\). A question about the hundredth Fibonacci number is now a question about \(\mat{F}^{100}\), and lesson M03.2.2 answers questions about powers through eigenvalues.
Step 4. The eigenvalues. The trace of \(\mat{F}\) is 1 and its determinant \({0 - 1} = -1\), so its characteristic polynomial is \(\lambda^2 - \lambda - 1\) (lesson M03.1.2). The quadratic formula gives
The middle step uses \(\lambda^2 = \lambda + 1\), which is the characteristic equation. So \(\vect{v}_1 = (\varphi, 1)\) and \(\vect{v}_2 = (\psi, 1)\), with slopes \(1/\varphi\) and \(1/\psi = -\varphi\): the two lines of the anchor figure.
The start in eigen-coordinates. We need \(c_1, c_2\) with \({c_1(\varphi, 1) + c_2(\psi, 1)} = (1, 0)\). The second entries give \(c_1 + c_2 = 0\), so \(c_2 = -c_1\). The first entries then give \({c_1(\varphi - \psi)} = 1\), and \(\varphi - \psi = \sqrt5\), so
Left: \(\mat{F}\) stretches \(\vect{v}_1\) by \(\varphi\) and flips and shrinks \(\vect{v}_2\) by \(\psi\). Right: the start \(\vect{u}_0 = (1, 0)\) is the difference of the two eigenvectors, divided by \(\sqrt5\).
NoteResult: Binet’s formula
For every \(k \ge 0\),
\[
F_k = \frac{\varphi^k - \psi^k}{\sqrt5},
\]
with \(\varphi = {(1 + \sqrt5)/2}\) and \(\psi = {(1 - \sqrt5)/2}\).
Proof. By the vector form of lesson M03.2.2, each eigen-coordinate is multiplied by its eigenvalue at every step:
The second entry of \(\vect{u}_k\) is \(F_k\), and the second entry of each eigenvector is 1. So \(F_k = {(\varphi^k - \psi^k)/\sqrt5}\). \(\blacksquare\)
Rounding is enough.\(\abs{\psi} < 1\), so \(\abs{\psi^k/\sqrt5} \le 1/\sqrt5 < \tfrac12\) for every \(k \ge 0\). The whole number \(F_k\) is therefore within \(\tfrac12\) of \(\varphi^k/\sqrt5\): \(F_k\) is \(\varphi^k/\sqrt5\) rounded to the nearest whole number.
The ratio of neighbours. Divide \(F_{k+1}\) by \(F_k\) and subtract \(\varphi\). Over the common denominator \({\varphi^k - \psi^k}\), the numerator is \({\varphi^{k+1} - \psi^{k+1} - \varphi^{k+1} + \varphi\psi^k} = {\psi^k(\varphi - \psi)}\), so
As \(k\) grows the denominator behaves like \(\varphi^k\), so the error behaves like \({\sqrt5\,(\psi/\varphi)^k}\). It changes sign at every step, because \(\psi < 0\), and shrinks by the factor \(\abs{\psi/\varphi} = 1/\varphi^2 \approx 0.382\), the zigzag of the anchor figure.
Any linear recurrence of order two. The same steps work for \(a_{k+2} = {p\,a_{k+1} + q\,a_k}\). The state \((a_{k+1}, a_k)\) moves by the companion matrix\(\begin{bmatrix} p & q \\ 1 & 0 \end{bmatrix}\), whose characteristic polynomial is \({\lambda^2 - p\lambda - q}\): exactly the equation you get by guessing \(a_k = \lambda^k\) in the recurrence. With two distinct roots, \(a_k = {c_1\lambda_1^k + c_2\lambda_2^k}\), and the two starting values fix \(c_1\) and \(c_2\) (exercise 2).
Frontier note: momentum is a 2 by 2 recurrence too
Gradient descent with momentum (the heavy-ball method, module M08.4) updates a weight using the current step and the previous one, \(w_{k+1} = {w_k - \eta\,\nabla f(w_k) + \beta\,(w_k - w_{k-1})}\). On a quadratic \(f(w) = \tfrac12\lambda w^2\) the gradient is \(\lambda w\), so
a companion matrix exactly like \(\mat{F}\). The error shrinks at the rate of its largest eigenvalue in size. Its determinant is \(\beta\), so when the two eigenvalues are complex, they are conjugates with product \(\beta\) and both have size \(\sqrt\beta\), whatever \(\eta\) and \(\lambda\) are. Gabriel Goh’s article analyses momentum in exactly this way, coordinate by coordinate in the eigenbasis of the quadratic. With curvatures from 0.1 to 10 and his tuned step and momentum, every direction shrinks by \(9/11 \approx 0.818\) per step, against about \(0.980\) for plain gradient descent at its best step: about 69 steps instead of 691 to cut the error a millionfold.
Worked example
Compute \(\mat{F}^{10}\), and \(F_{10}\) from Binet’s formula.
The matrix. Column 1 of \(\mat{F}^k\) is \(\mat{F}^k(1, 0) = \mat{F}^k\vect{u}_0 = \vect{u}_k\), which is \({(F_{k+1}, F_k)}\). Column 2 is \(\mat{F}^k(0, 1)\). One step sends \((0, 1)\) to \({(0 + 1,\; 0)} = (1, 0) = \vect{u}_0\), so \((0, 1)\) is one step behind \(\vect{u}_0\), and \(\mat{F}^k(0, 1) = \mat{F}^{k-1}\vect{u}_0 = \vect{u}_{k-1} = {(F_k, F_{k-1})}\). Together
reading \(F_9 = 34\), \(F_{10} = 55\) and \(F_{11} = 89\) off the sequence.
Binet. Now the formula, with numbers.
M03.2.3-F3
Binet’s formula for \(F_{10}\), from the start to the number. The \(\psi\) term is tiny, which is why rounding \(\varphi^{10}/\sqrt5\) already gives 55.
\(\varphi^{10} \approx 122.99187\) and \(\psi^{10} \approx 0.00813\), so
Dropping \(\psi^{10}\) gives \(122.99187 / 2.23607 \approx 55.004\), which rounds to 55, as the rounding rule promised.
The ratio.\(F_{11}/F_{10} = 89/55 \approx 1.61818\), against \(\varphi \approx 1.61803\). The error is about \({1.5 \times 10^{-4}}\), and the previous one, \(55/34 - \varphi \approx {-3.9 \times 10^{-4}}\), was about \(1/0.382\) times as large and of the other sign.
Code
From scratch: exact integer powers of \(\mat{F}\) by repeated squaring (lesson M03.2.2), in Python’s unlimited integers.
def matmul(X, Y):return [[sum(X[i][t] * Y[t][j] for t inrange(2)) for j inrange(2)] for i inrange(2)]def mat_power(M, k):"""M^k by repeated squaring, exact for integer matrices.""" result, square = [[1, 0], [0, 1]], Mwhile k:if k &1: result = matmul(result, square) square = matmul(square, square) k >>=1return resultF = [[1, 1], [1, 0]]print("F^10 =", mat_power(F, 10))print("F_100 =", mat_power(F, 100)[0][1])
Now Binet’s formula in floating point, compared with the exact numbers until the first mistake. To make the result the same on every machine, \(\varphi^n\) is built by repeated multiplication: IEEE arithmetic rounds each product and math.sqrt the same way everywhere.
import mathphi = (1+ math.sqrt(5)) /2psi = (1- math.sqrt(5)) /2print("phi, psi:", phi, psi, " phi * psi:", phi * psi)def binet(n): power =1.0for _ inrange(n): power *= phireturnround(power / math.sqrt(5))first_wrong =next(n for n inrange(200) if binet(n) != fib[n])print("first wrong index:", first_wrong)print("exact:", fib[first_wrong], " formula:", binet(first_wrong))
The formula is right up to \(F_{70}\) and wrong from \(F_{71} = 308061521170129\) on, a 15-digit number. A double holds about 16 significant digits, and the 71 roundings in \(\varphi^{71}\) eat into the last of them, so the error passes \(\tfrac12\) already here. Beyond about \({9 \times 10^{15}}\) (that is, \(2^{53}\)) a double cannot even hold every whole number, so no amount of care in the formula could save the last digit. NumPy’s int64 has a limit too: np.linalg.matrix_power on an integer \(\mat{F}\) is exact up to \(F_{92}\) and wraps round silently at \(F_{93}\), which exceeds \(2^{63} - 1\).
import numpy as npF_np = np.array([[1, 1], [1, 0]])print("F_92 via int64:", np.linalg.matrix_power(F_np, 92)[0, 1], " exact:", fib[92])print("F_93 via int64:", np.linalg.matrix_power(F_np, 93)[0, 1], " exact:", fib[93])w, V = np.linalg.eig(F_np.astype(float))print("eigenvalues:", np.sort(w))
F_92 via int64: 7540113804746346429 exact: 7540113804746346429
F_93 via int64: -6246583658587674878 exact: 12200160415121876738
eigenvalues: [-0.61803399+0.j 1.61803399+0.j]
Explore
This module’s interactive, the power stepper, lives in lesson M03.2.2. Choose the Fibonacci matrix and put the start at \((1, 0)\). Each press of step produces the next state \((F_{k+1}, F_k)\). Watch the second bar, \(c_2\psi^k\): it changes sign at every step and shrinks, which is the zigzag of the ratio.
Pitfalls
M03.2.3-F4
Wrong beside right. Rounded Binet in floating point keeps up with the true digits until about 15 of them are needed, then fails for good from \(n = 71\). Integer matrix powers stay exact.
“Binet’s formula is exact, so its computer version is exact too.”
The formula is exact in real numbers. A floating-point number carries about 16 significant digits, so \(\varphi^n/\sqrt5\) is computed with a tiny relative error that becomes an absolute error above \(\tfrac12\) once the numbers are about \(10^{15}\). From then on rounding picks the wrong whole number. Use the formula for insight and for the growth rate; use integer powers for exact values.
“The state is just \(F_k\), so the matrix should be 1 by 1.”
The rule looks back two steps, so one number cannot carry enough information to take a step. The state must hold everything the rule needs, here \(F_{k+1}\) and \(F_k\). A recurrence that looks back \(m\) steps needs a state of \(m\) numbers and an \(m\) by \(m\) matrix.
“Since \(\abs{\psi} < 1\), the \(\psi^k\) term can be dropped: \(F_k = \varphi^k/\sqrt5\).”
Dropping it gives a number that is close to \(F_k\) but never equal to it, because \(\varphi^k/\sqrt5\) is irrational for \(k \ge 1\). For \(k = 1\) it gives \(0.7236\), not 1. Drop it only together with rounding to the nearest whole number, which is exactly what the bound \(\abs{\psi^k/\sqrt5} < \tfrac12\) allows.
Check yourself
Recall
1. What is the state vector of the Fibonacci recurrence, and which matrix moves it one step?
Show answer
\(\vect{u}_k = (F_{k+1}, F_k)\), and \(\vect{u}_{k+1} = \mat{F}\vect{u}_k\) with \(\mat{F} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\).
2. What are the eigenvalues and eigenvectors of \(\mat{F}\)?
Show answer
\(\varphi = {(1 + \sqrt5)/2}\) and \(\psi = {(1 - \sqrt5)/2} = -1/\varphi\), with eigenvectors \((\varphi, 1)\) and \((\psi, 1)\).
3. State Binet’s formula.
Show answer
\(F_k = {(\varphi^k - \psi^k)/\sqrt5}\).
4. At what rate does \(F_{k+1}/F_k\) approach \(\varphi\)?
Show answer
The error shrinks by the factor \(\abs{\psi/\varphi} = 1/\varphi^2 \approx 0.382\) at each step, changing sign each time.
5. What is \(\mat{F}^{10}\)?
Show answer
\(\begin{bmatrix} 89 & 55 \\ 55 & 34 \end{bmatrix}\), that is, \(\begin{bmatrix} F_{11} & F_{10} \\ F_{10} & F_9 \end{bmatrix}\).
Exercises
1.★ Use \(\mat{F}^k = \begin{bmatrix} F_{k+1} & F_k \\ F_k & F_{k-1} \end{bmatrix}\) to write down \(\mat{F}^5\), then check it by squaring \(\mat{F}^2\) and multiplying by \(\mat{F}\).
by the Fibonacci rule, which is the claim for \(k + 1\). For Cassini, take determinants. The determinant of a product is the product of the determinants (M02.8), so \(\det\mat{F}^k = (\det\mat{F})^k = (-1)^k\). The determinant of the right side is \({F_{k+1}F_{k-1} - F_k^2}\). For \(k = 10\): \({89 \cdot 34} - 55^2 = 3026 - 3025 = 1\).
Summary card
A recurrence that looks back two steps is a 2 by 2 matrix power.
Binet.\(F_k = {(\varphi^k - \psi^k)/\sqrt5}\), and \(F_k\) is \(\varphi^k/\sqrt5\) rounded.
Ratio.\(F_{k+1}/F_k \to \varphi\), error times about \(-0.382\) per step.
In floats. Rounded Binet first fails at \(F_{71}\); integer powers stay exact.
Sources and next steps
Sources. Each was opened and checked against this lesson on 15 November 2026.
Gilbert Strang, Diagonalization and powers of A, MIT OpenCourseWare 18.06SC, lecture 22 and its lecture summary, section Fibonacci sequence. It converts \(F_{k+2} = F_{k+1} + F_k\) into \(u_{k+1} = Au_k\) with \(u_k = (F_{k+1}, F_k)\), finds the eigenvalues \({\tfrac12(1 \pm \sqrt5)}\) and eigenvectors \((\lambda, 1)\), solves \(u_0 = (1, 0)\) for \(c_1 = -c_2 = 1/\sqrt5\) and arrives at the closed form for \(F_k\).
Gabriel Goh, Why Momentum Really Works, Distill, 2017. On a quadratic, each eigen-coordinate of momentum follows a 2 by 2 linear recurrence, the rate of convergence is the larger size of that matrix’s two eigenvalues, and the step size and momentum are tuned from the smallest and largest curvature. The \(\sqrt\beta\) argument through the determinant is this lesson’s own.
NumPy reference, numpy.linalg.matrix_power. Integer powers by repeated squaring, with the element type of the input kept for non-negative exponents (so an int64 input stays int64, and can overflow).
Next.M03.2.4 The power method runs the same repeated multiplication on purpose, to find the dominant eigenvector of a matrix whose eigenvectors we do not yet know.