Phase 19 - Lesson 19.13

Matrix Exponentiation and Linear Recurrence Acceleration

Any linear recurrence is a matrix product; square-and-multiply then evaluates term \(10^{18}\) in \(O(k^3\log n)\) operations.

⏱ 60 min● Advanced🔗 Prereqs: 19.3 fast exponentiation; 19.6 recurrences; Phase 4 linear algebra
↖ Phase 19 hub
Builds on: 19.3 exponentiated numbers; 19.6 wrote the recurrences. This lesson multiplies them together.
Leads to: 19.14 uses this as the canonical example of ‘change the exponent, not the constant’. The Euler Lab’s Number Theory and Dynamic Programming paths both have huge-index recurrence problems.

Learning Objectives

Click a status chip to cycle: Not started → In progress → Studied → Practiced → Needs review → Mastered.

Key Vocabulary

Companion matrix
The \(k\times k\) matrix that shifts a state vector of the last \(k\) terms and appends the new one; its characteristic polynomial IS the recurrence’s.
State vector
\(v_n=(a_n,a_{n-1},\dots,a_{n-k+1})^\top\); the recurrence is exactly \(v_n=Av_{n-1}\), hence \(v_n=A^{n-1}v_1\).
Matrix power by squaring
Same square-and-multiply as 19.3 with matrix multiplication as the operation: \(O(\log n)\) multiplications, each \(O(k^3)\).
Adjacency matrix
\(A_{ij}=1\) iff there is an edge \(i\to j\); then \((A^n)_{ij}\) counts walks of length exactly \(n\) from \(i\) to \(j\).
Eigenvalue growth
If \(A\) is diagonalisable with eigenvalues \(\lambda_i\), then \(a_n=\sum_i c_i\lambda_i^n\) - the characteristic-root formula of 19.6, in matrix clothing.
Associativity
The reason this works at all: \(A^{2m}=(A^m)^2\) requires only that matrix multiplication be associative, which it is.

Intuition & Motivation

Intuition
A linear recurrence is a machine that turns the last \(k\) values into the next one - a linear map. Applying it \(n\) times is applying a single matrix \(n\) times, and a matrix raised to the \(n\)-th power can be computed by repeated squaring, exactly as with numbers. The cost collapses from \(n\) steps to \(\log_2 n\) matrix multiplications. This is the cleanest example in the phase of the deepest habit: change the exponent, not the constant.

The companion matrix

Take Fibonacci. Stack the last two values in a vector and write the recurrence as a matrix acting on it:

\[\begin{pmatrix}F_{n+1}\\ F_n\end{pmatrix}=\begin{pmatrix}1&1\\ 1&0\end{pmatrix}\begin{pmatrix}F_{n}\\ F_{n-1}\end{pmatrix}\quad\Longrightarrow\quad \begin{pmatrix}F_{n+1}&F_n\\ F_n&F_{n-1}\end{pmatrix}=\begin{pmatrix}1&1\\ 1&0\end{pmatrix}^{\!n}.\] (19.29)

The right-hand identity follows by induction, and taking determinants of both sides gives Cassini’s identity for free: \(F_{n+1}F_{n-1}-F_n^2=(-1)^n\) (because \(\det\begin{pmatrix}1&1\\1&0\end{pmatrix}=-1\)). One matrix identity, two theorems.

Definition - General companion form
For \(a_n=c_1a_{n-1}+c_2a_{n-2}+\cdots+c_ka_{n-k}\), let \(v_n=(a_n,a_{n-1},\dots,a_{n-k+1})^\top\) and
\[A=\begin{pmatrix} c_1&c_2&\cdots&c_{k-1}&c_k\\ 1&0&\cdots&0&0\\ 0&1&\cdots&0&0\\ \vdots& &\ddots& &\vdots\\ 0&0&\cdots&1&0\end{pmatrix},\qquad v_n=Av_{n-1}=A^{n-1}v_1 .\] (19.30)

The top row computes the new term; the subdiagonal shifts the window. The characteristic polynomial of \(A\) is exactly \(x^k-c_1x^{k-1}-\cdots-c_k\) - the same polynomial as 19.6. So the eigenvalues of the companion matrix ARE the characteristic roots, and \(A^n\) growing like \(\lambda_{\max}^n\) is the matrix version of \(a_n\sim\alpha\varphi^n\). Two lessons, one theorem.

Cost accounting

\[\text{iterate: } O(n)\ \text{steps}\qquad\text{vs}\qquad \text{matrix power: } O(k^3\log n)\ \text{scalar multiplications (mod }p).\] (19.31)
TaskIterationMatrix powerWinner
\(F_{10^6}\) exact\(10^6\) big-int adds\(8\cdot20\) big-int muls on huge numbersiteration (the numbers themselves are the cost)
\(F_{10^{18}}\bmod p\)\(10^{18}\) steps - impossible\(8\cdot60=480\) modmulsmatrix power, decisively
Walks of length \(10^9\) in a 50-node graphimpossible\(50^3\cdot30\approx4\times10^6\)matrix power
\(F_{100})\) exact100 adds20 matrix mulsiteration (simpler, and no faster to be clever)

Note the honest caveat in row 1: when you need the exact \(F_{10^6}\) (a 208 988-digit number), the bottleneck is big-integer arithmetic, not the number of steps, and the matrix method’s \(O(\log n)\) advantage is eaten by multiplying enormous integers. Matrix exponentiation is a modular weapon: it shines when every entry stays bounded, which is exactly the ‘give the answer mod \(10^9+7\)’ setting.

Worked Example - Compute \(F_{10}\) by squaring, and then \(F_{10^{18}}\bmod (10^9+7)\)
1
Let \(M=\begin{pmatrix}1&1\\1&0\end{pmatrix}\). Then \(M^2=\begin{pmatrix}2&1\\1&1\end{pmatrix}\) and \(M^4=(M^2)^2=\begin{pmatrix}5&3\\3&2\end{pmatrix}\).
2
\(M^8=(M^4)^2=\begin{pmatrix}34&21\\21&13\end{pmatrix}\). Now \(10=8+2\), so \(M^{10}=M^8M^2=\begin{pmatrix}34&21\\21&13\end{pmatrix}\begin{pmatrix}2&1\\1&1\end{pmatrix}=\begin{pmatrix}89&55\\55&34\end{pmatrix}\).
3
By (19.29), \(F_{10}=55\) (the top-right entry) and \(F_{11}=89\). ✓ Four matrix multiplications instead of ten additions - not a win at \(n=10\), but the ALGORITHM is now independent of \(n\) except through \(\log n\).
4
For \(n=10^{18}\): the exponent has 60 bits, so at most 60 squarings and 60 multiplies of \(2\times2\) matrices, each 8 modular multiplications: fewer than 1000 modmuls total. Runtime: microseconds.
5
Sanity check the whole pipeline on a value you can verify independently - e.g. confirm the routine returns \(F_{100}=354{,}224{,}848{,}179{,}261{,}915{,}075\) exactly before you trust it at \(10^{18}\).

Beyond recurrences: counting walks

If \(A\) is the adjacency matrix of a digraph, then \((A^n)_{ij}\) is the number of walks of length exactly \(n\) from \(i\) to \(j\) - because matrix multiplication sums over intermediate vertices, which is precisely the ‘first step, then a walk of length \(n-1\)’ decomposition. So ‘how many length-\(10^9\) paths?’ is a matrix power, not a DP. Any DP whose transition is linear and time-invariant can be accelerated this way - that is the general statement, and recurrences are just the \(k\)-state special case.

Coding workbench

Whenever an Euler Lab problem gives you a recurrence and then asks for term \(10^{15}\) (mod something), the intended solution is this lesson. Recognising a linear, time-invariant transition is the whole skill; the code is 20 lines you will reuse forever.

Common Mistakes to Avoid
  • Using floating-point matrices (or NumPy’s default float64) for exact counting. Entries overflow 53 bits almost immediately - stay in Python ints mod \(p\) (19.4).
  • Forgetting to reduce mod \(p\) inside the multiplication loop; entries then grow to thousands of digits and the \(O(k^3)\) cost becomes a lie.
  • Getting the index off by one: \((M^n)_{01}=F_n\), not \(F_{n+1}\). Always verify against a known small value before scaling up.
  • Applying matrix exponentiation to a NON-linear or time-varying recurrence (e.g. \(a_n=a_{n-1}^2\) or coefficients that depend on \(n\)). The method requires a fixed linear map.
  • Using it when \(n\) is small. For \(n\le10^7\), plain iteration is simpler, has a better constant, and cannot be off by one.
Quant Practitioner Tips
  • Any DP with a fixed linear transition can be accelerated: write the transition as a matrix and exponentiate. Path counting, restricted-string counting, and Markov-chain steps are all this.
  • The companion matrix’s characteristic polynomial is the recurrence’s characteristic polynomial - so eigenvalues give you the growth rate (19.6) and the matrix gives you the exact values. Use both.
  • Determinant tricks are free: \(\det(M^n)=(\det M)^n\) gave Cassini’s identity in one line.
  • For Pell (19.12), the map \((x,y)\mapsto(x_1x+Dy_1y,\ x_1y+y_1x)\) is linear, so the \(n\)-th solution is a matrix power - \(O(\log n)\) instead of \(n\) multiplications.

Practice this in the Euler Lab

Computational problems that exercise exactly this technique. Each opens in the Euler Lab with a Python workbench, a progressive hint ladder, and answer checking. Tier A/B run at full scale in the browser.

Warm-up:
#2 Even Fibonacci Numbers (1%, tier A) #25 $1000$-digit Fibonacci Number (4%, tier A) #81 Path Sum: Two Ways (7%, tier A) #104 Pandigital Fibonacci Ends (8%, tier A)

Applied:
#345 Matrix Sum (17%, tier B) #186 Connectedness of a Network (19%, tier B) #437 Fibonacci Primitive Roots (23%, tier B)

Challenge:
#157 Base-10 Diophantine Reciprocal (42%, tier C) #304 Primonacci (42%, tier C)

57 Project Euler problems in total are mapped to this lesson. Open the Euler Lab to filter them all.

Knowledge Check

Q1 Medium
Why does \(\begin{pmatrix}1&1\\1&0\end{pmatrix}^n=\begin{pmatrix}F_{n+1}&F_n\\F_n&F_{n-1}\end{pmatrix}\) let you compute \(F_n\) fast?
Because matrices are faster than integers
Because a matrix power can be computed by repeated squaring in \(O(\log n)\) multiplications
Because the matrix is symmetric
Because the determinant is \(-1\)
Q2 Medium
Matrix exponentiation is the wrong tool when:
The recurrence is linear with constant coefficients
The recurrence is non-linear (e.g. \(a_n=a_{n-1}^2+1\)) or its coefficients depend on \(n\)
You need the answer mod a prime
n is large
Q3 Medium
The cost of computing term \(n\) of a \(k\)-term linear recurrence mod \(p\) by matrix power is:
\(O(n)\)
\(O(k^3\log n)\)
\(O(k\log n)\)
\(O(n\log k)\)

Practical Exercise

(a) Write the companion matrix for the tribonacci recurrence \(T_n=T_{n-1}+T_{n-2}+T_{n-3}\) and state the cost of computing \(T_{10^{18}}\bmod p\). (b) Prove Cassini’s identity \(F_{n+1}F_{n-1}-F_n^2=(-1)^n\) from (19.29) in one line. (c) A frog hops on a 5-vertex cycle. In how many ways can it return to its start after exactly \(10^9\) hops? Give the method and the complexity, not the number.

▶ Show full solution

(a) \(A=\begin{pmatrix}1&1&1\\1&0&0\\0&1&0\end{pmatrix}\) acting on \(v_n=(T_n,T_{n-1},T_{n-2})^\top\). Cost: \(\lceil\log_2 10^{18}\rceil\approx60\) squarings plus up to 60 multiplies, each a \(3\times3\) matmul at \(27\) modular multiplications - roughly \(120\times27\approx3{,}200\) modmuls. Microseconds.

(b) Take determinants of both sides of (19.29). The left is \(\det(M^n)=(\det M)^n=(-1)^n\); the right is \(F_{n+1}F_{n-1}-F_n^2\). Equate. □ (This is the sort of identity that takes half a page by induction and one line by structure - which is the argument for finding the right representation.)

(c) The frog’s moves are a linear, time-invariant transition on 5 states, so let \(A\) be the adjacency matrix of the 5-cycle (\(A_{ij}=1\) iff \(j=i\pm1\bmod 5\)). The number of closed walks of length \(10^9\) from vertex 0 is \((A^{10^9})_{00}\). Compute it by square-and-multiply: \(\lceil\log_2 10^9\rceil=30\) squarings, each a \(5\times5\) matmul at \(125\) multiplications - about \(7{,}500\) operations, versus \(10^9\) for a step-by-step DP. Complexity \(O(k^3\log n)=O(125\cdot30)\). (The counts are astronomically large, so work mod a prime unless you want a 300-million-digit integer.)

After the reveal, answer for yourself: In (b) the identity fell out of the determinant. What does that suggest about looking for structure BEFORE reaching for induction?

Lesson Summary

Any linear, time-invariant recurrence is a single matrix applied repeatedly: \(v_n=Av_{n-1}=A^{n-1}v_1\) with \(A\) the companion matrix, whose characteristic polynomial is exactly the recurrence’s (so its eigenvalues are 19.6’s characteristic roots). Since matrix multiplication is associative, square-and-multiply computes \(A^n\) in \(O(k^3\log n)\) scalar operations - turning \(F_{10^{18}}\bmod p\) from impossible into microseconds. The same trick counts length-\(n\) walks in a graph and generates Pell solutions.

Retrieval Practice

Close the lesson and answer from memory before checking. This is deliberate, effortful recall - the single highest-yield study action.

▶ Show retrieval prompts & answers
Q: Why can a linear recurrence be evaluated in \(O(\log n)\) matrix multiplications?
A: Because the recurrence is a fixed linear map \(A\) on the state vector, so \(v_n=A^{n-1}v_1\), and matrix multiplication is associative - so \(A^n\) yields to square-and-multiply. Cost \(O(k^3\log n)\).
Q: When does matrix exponentiation NOT apply, and when is it not worth it?
A: It does not apply to non-linear or time-varying recurrences (there is no single \(A\)). It is not worth it for small \(n\) (iteration is simpler and has a better constant), nor for exact huge values, where big-integer cost dominates.

Completion Checklist

Confidence / mastery rating
Personal notes

Source References

This lesson synthesizes and paraphrases concepts from the sources below. No copyrighted text, problem sets, or solutions are reproduced. Return to the originals for full depth.