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