Matrix Exponentiation: Fast Linear Recurrence

Key Insights

Matrix exponentiation transforms O(n) linear recurrence computations into O(log n) by representing the recurrence as matrix multiplication and applying binary exponentiation.
Any k-th order linear recurrence can be solved using a k×k companion matrix, making this technique universally applicable to sequences like Fibonacci, Tribonacci, and custom recurrences.
The technique extends beyond sequences to graph problems, dynamic programming optimization, and any scenario involving repeated linear transformations.

The Fibonacci Problem at Scale

Computing the nth Fibonacci number seems trivial. Loop n times, track two variables, done. But what happens when n equals 10^18?

At one billion iterations per second, that naive O(n) loop would take over 31 years. This isn’t a theoretical concern—competitive programming problems regularly demand computing sequence values at astronomical indices, and cryptographic applications need fast modular exponentiation of recurrence relations.

Matrix exponentiation solves this by reducing the problem to O(log n) matrix multiplications. For n = 10^18, that’s roughly 60 operations instead of a quintillion. The technique leverages two key insights: linear recurrences can be expressed as matrix multiplication, and matrices can be exponentiated using the same binary exponentiation trick that works for integers.

From Recurrence to Matrix Form

The Fibonacci recurrence F(n) = F(n-1) + F(n-2) involves two previous terms. We can express the transition from one state to the next as a matrix equation:

[F(n)  ]   [1 1] [F(n-1)]
[F(n-1)] = [1 0] [F(n-2)]

Verify this mentally: the top row computes 1×F(n-1) + 1×F(n-2) = F(n), and the bottom row computes 1×F(n-1) + 0×F(n-2) = F(n-1). The matrix preserves F(n-1) while computing the next term.

The critical observation: applying this transformation n-1 times starting from the base case gives us F(n):

[F(n)  ]   [1 1]^(n-1) [F(1)]   [1 1]^(n-1) [1]
[F(n-1)] = [1 0]       [F(0)] = [1 0]       [0]

Here’s the basic implementation:

def matrix_mult(A, B):
    """Multiply two 2x2 matrices."""
    return [
        [A[0][0]*B[0][0] + A[0][1]*B[1][0], A[0][0]*B[0][1] + A[0][1]*B[1][1]],
        [A[1][0]*B[0][0] + A[1][1]*B[1][0], A[1][0]*B[0][1] + A[1][1]*B[1][1]]
    ]

def fibonacci_naive_matrix(n):
    """Compute F(n) using O(n) matrix multiplications."""
    if n <= 1:
        return n
    
    result = [[1, 0], [0, 1]]  # Identity matrix
    transform = [[1, 1], [1, 0]]
    
    for _ in range(n - 1):
        result = matrix_mult(result, transform)
    
    return result[0][0]

This is still O(n), but we’ve set up the framework for the real optimization.

Fast Exponentiation: The Core Algorithm

Binary exponentiation computes x^n in O(log n) multiplications by exploiting the binary representation of n. The same principle applies to matrices.

Consider computing M^13. Since 13 = 1101 in binary (8 + 4 + 1), we have M^13 = M^8 × M^4 × M^1. We compute M^1, M^2, M^4, M^8 by repeated squaring, then multiply together only the powers corresponding to set bits.

def matrix_mult(A, B, size=2):
    """Multiply two square matrices of given size."""
    result = [[0] * size for _ in range(size)]
    for i in range(size):
        for j in range(size):
            for k in range(size):
                result[i][j] += A[i][k] * B[k][j]
    return result

def matrix_pow(M, n, size=2):
    """Compute M^n using binary exponentiation."""
    # Start with identity matrix
    result = [[1 if i == j else 0 for j in range(size)] for i in range(size)]
    
    base = [row[:] for row in M]  # Copy M
    
    while n > 0:
        if n & 1:  # If current bit is set
            result = matrix_mult(result, base, size)
        base = matrix_mult(base, base, size)  # Square the base
        n >>= 1
    
    return result

def fibonacci_fast(n):
    """Compute F(n) in O(log n) time."""
    if n <= 1:
        return n
    
    transform = [[1, 1], [1, 0]]
    result = matrix_pow(transform, n - 1)
    return result[0][0]

The matrix_pow function performs O(log n) matrix multiplications. Each multiplication of k×k matrices takes O(k³) time, giving overall complexity O(k³ log n). For Fibonacci with k=2, this is effectively O(log n).

Generalizing to k-th Order Recurrences

The Tribonacci sequence follows T(n) = T(n-1) + T(n-2) + T(n-3). The companion matrix becomes 3×3:

[T(n)  ]   [1 1 1] [T(n-1)]
[T(n-1)] = [1 0 0] [T(n-2)]
[T(n-2)]   [0 1 0] [T(n-3)]

For a general k-th order recurrence f(n) = c₁f(n-1) + c₂f(n-2) + … + cₖf(n-k), the companion matrix has coefficients in the first row and a shifted identity below:

def solve_linear_recurrence(coefficients, initial_values, n):
    """
    Solve a linear recurrence relation.
    
    coefficients: [c1, c2, ..., ck] where f(n) = c1*f(n-1) + c2*f(n-2) + ... + ck*f(n-k)
    initial_values: [f(0), f(1), ..., f(k-1)]
    n: index to compute
    """
    k = len(coefficients)
    
    if n < k:
        return initial_values[n]
    
    # Build companion matrix
    matrix = [[0] * k for _ in range(k)]
    
    # First row contains coefficients
    for j in range(k):
        matrix[0][j] = coefficients[j]
    
    # Subdiagonal is 1s (shifted identity)
    for i in range(1, k):
        matrix[i][i-1] = 1
    
    # Compute matrix^(n-k+1)
    result = matrix_pow(matrix, n - k + 1, k)
    
    # Multiply by initial state vector (reversed order)
    answer = 0
    for j in range(k):
        answer += result[0][j] * initial_values[k - 1 - j]
    
    return answer

# Tribonacci: T(n) = T(n-1) + T(n-2) + T(n-3), with T(0)=0, T(1)=1, T(2)=1
def tribonacci(n):
    return solve_linear_recurrence([1, 1, 1], [0, 1, 1], n)

Handling Modular Arithmetic

Large sequence values overflow quickly. Competitive programming problems typically require answers modulo 10^9 + 7. Apply the modulus after every arithmetic operation:

MOD = 10**9 + 7

def matrix_mult_mod(A, B, size, mod):
    """Multiply matrices with modular arithmetic."""
    result = [[0] * size for _ in range(size)]
    for i in range(size):
        for j in range(size):
            for k in range(size):
                result[i][j] = (result[i][j] + A[i][k] * B[k][j]) % mod
    return result

def matrix_pow_mod(M, n, size, mod):
    """Matrix exponentiation with modulus."""
    result = [[1 if i == j else 0 for j in range(size)] for i in range(size)]
    base = [row[:] for row in M]
    
    while n > 0:
        if n & 1:
            result = matrix_mult_mod(result, base, size, mod)
        base = matrix_mult_mod(base, base, size, mod)
        n >>= 1
    
    return result

def fibonacci_mod(n, mod=MOD):
    """Compute F(n) mod m in O(log n) time."""
    if n <= 1:
        return n
    
    transform = [[1, 1], [1, 0]]
    result = matrix_pow_mod(transform, n - 1, 2, mod)
    return result[0][0]

# F(10^18) mod (10^9 + 7) computes in milliseconds
print(fibonacci_mod(10**18))

For C++ implementations where performance matters, use fixed-size arrays and avoid dynamic allocation:

#include <vector>
using namespace std;

const long long MOD = 1e9 + 7;

vector<vector<long long>> matmul(const vector<vector<long long>>& A,
                                  const vector<vector<long long>>& B) {
    int n = A.size();
    vector<vector<long long>> C(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++)
            for (int j = 0; j < n; j++)
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
    return C;
}

vector<vector<long long>> matpow(vector<vector<long long>> M, long long p) {
    int n = M.size();
    vector<vector<long long>> result(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++) result[i][i] = 1;
    
    while (p > 0) {
        if (p & 1) result = matmul(result, M);
        M = matmul(M, M);
        p >>= 1;
    }
    return result;
}

Practical Applications & Variations

Beyond sequences, matrix exponentiation shines in graph problems. The number of paths of exactly length n between vertices i and j in a directed graph equals (A^n)[i][j], where A is the adjacency matrix.

def count_paths(adj_matrix, length, start, end, mod=MOD):
    """
    Count paths of exact length from start to end vertex.
    adj_matrix[i][j] = 1 if edge exists from i to j, else 0.
    """
    n = len(adj_matrix)
    result = matrix_pow_mod(adj_matrix, length, n, mod)
    return result[start][end]

# Example: 4-vertex graph
# 0 -> 1, 0 -> 2, 1 -> 2, 1 -> 3, 2 -> 3
adj = [
    [0, 1, 1, 0],
    [0, 0, 1, 1],
    [0, 0, 0, 1],
    [0, 0, 0, 0]
]

# Paths of length 2 from vertex 0 to vertex 3
print(count_paths(adj, 2, 0, 3))  # Output: 2 (0->1->3 and 0->2->3)

This technique also optimizes certain dynamic programming problems where transitions form a fixed linear pattern, transforming O(n×k) DP into O(k³ log n).

Complexity Analysis & When to Use

Time complexity: O(k³ log n) where k is the matrix dimension (recurrence order).

Space complexity: O(k²) for storing matrices.

Use matrix exponentiation when:

Computing the nth term of a linear recurrence for very large n (typically n > 10^6)
The recurrence order k is small (k ≤ 100 is practical; k³ grows fast)
You need exact values or modular results, not approximations

Avoid it when:

n is small enough that O(n) iteration is acceptable
The recurrence is non-linear (matrix exponentiation won’t help)
You only need approximate values (closed-form solutions like Binet’s formula work)
The recurrence order k is large relative to n

For Fibonacci specifically, the crossover point where matrix exponentiation beats iteration is around n = 100-1000 depending on implementation. But when n reaches millions or billions, there’s no alternative.

Matrix exponentiation is a fundamental technique that transforms intractable computations into millisecond operations. Master the pattern of converting recurrences to matrices, and you’ll recognize opportunities to apply it across domains from competitive programming to scientific computing.