Matrix Exponentiation: Fast Computation of Linear Recurrences

Have you ever needed to compute the 10^18th Fibonacci number? Doing it iteratively would take billions of years. But with matrix exponentiation, you can do it in milliseconds. This technique is a game-changer for problems involving linear recurrences, such as Fibonacci, Tribonacci, or any sequence defined by a linear combination of previous terms.

Matrix exponentiation works by representing the recurrence as a matrix multiplication. For Fibonacci, the recurrence F(n) = F(n-1) + F(n-2) can be expressed as: [F(n), F(n-1)]^T = [[1,1],[1,0]] * [F(n-1), F(n-2)]^T.

By raising the matrix to the power n-1, you can directly compute F(n) in O(log n) time using exponentiation by squaring. This is not just a theoretical exercise; it's used in real-world applications like cryptography, random number generation, and dynamic programming optimization.

In this tutorial, you'll learn the mathematics behind matrix exponentiation, how to implement it in Python, and how to debug it in production. We'll also explore a real incident where a naive implementation caused a system outage. By the end, you'll be able to apply this technique to any linear recurrence with confidence.

Understanding Linear Recurrences and Their Matrix Form

A linear recurrence is a sequence where each term is a linear combination of previous terms. For example, Fibonacci: F(n) = F(n-1) + F(n-2). More generally, a k-th order linear recurrence is: a_n = c1a_{n-1} + c2a_{n-2} + ... + ck*a_{n-k}.

We can represent this using matrices. Define a state vector S(n) = [a_n, a_{n-1}, ..., a_{n-k+1}]^T. Then the recurrence can be written as S(n) = M * S(n-1), where M is a k x k matrix called the transformation matrix. For Fibonacci (k=2): M = [[1, 1], [1, 0]] S(n) = [F(n), F(n-1)]^T S(n-1) = [F(n-1), F(n-2)]^T

Then S(n) = M S(n-1). By induction, S(n) = M^(n-1) S(1). So to compute a_n for large n, we just need to compute M^(n-1) efficiently.

import numpy as np

# Define transformation matrix for Fibonacci
M = np.array([[1, 1], [1, 0]])

# Initial state: F(1)=1, F(0)=0
S1 = np.array([1, 0])

# Compute S(10) = M^9 * S1
n = 10
M_pow = np.linalg.matrix_power(M, n-1)
S_n = M_pow @ S1
print(f"F({n}) = {S_n[0]}")

Matrix Multiplication and Exponentiation by Squaring

Matrix multiplication is the core operation. For two k x k matrices A and B, the product C = A B is defined as C[i][j] = sum_{r=0}^{k-1} A[i][r] B[r][j]. This is O(k^3) for naive multiplication.

Exponentiation by squaring (binary exponentiation) computes M^n in O(log n) multiplications. The idea: if n is even, M^n = (M^(n/2))^2; if n is odd, M^n = M * M^(n-1). This is similar to fast exponentiation for numbers.

Here's an iterative implementation:

def matrix_mult(A, B, mod=None):
    n = len(A)
    C = [[0]*n for _ in range(n)]
    for i in range(n):
        for k in range(n):
            aik = A[i][k]
            if aik:
                for j in range(n):
                    C[i][j] = (C[i][j] + aik * B[k][j]) % mod if mod else C[i][j] + aik * B[k][j]
    return C

def matrix_pow(M, power, mod=None):
    n = len(M)
    # Initialize result as identity matrix
    result = [[1 if i==j else 0 for j in range(n)] for i in range(n)]
    base = M
    while power:
        if power & 1:
            result = matrix_mult(result, base, mod)
        base = matrix_mult(base, base, mod)
        power >>= 1
    return result

Implementing Fibonacci with Matrix Exponentiation

Let's implement a complete function to compute the n-th Fibonacci number using matrix exponentiation. We'll use the transformation matrix M = [[1,1],[1,0]]. The initial state vector for F(0)=0, F(1)=1 is [1,0]^T. Then F(n) = (M^n * [1,0])[1] (the second element). Alternatively, we can compute M^(n-1) and multiply by [1,0] to get [F(n), F(n-1)].

Here's a clean implementation:

def fib_matrix(n, mod=None):
    if n == 0:
        return 0
    M = [[1, 1], [1, 0]]
    # Compute M^(n-1)
    result = matrix_pow(M, n-1, mod)
    # Multiply by initial vector [1,0]
    # F(n) = result[0][0]*1 + result[0][1]*0 = result[0][0]
    return result[0][0] % mod if mod else result[0][0]

# Example
print(fib_matrix(10))  # 55
print(fib_matrix(100, mod=10**9+7))  # 354224848179261915075 % mod

Generalizing to Arbitrary Linear Recurrences

Any k-th order linear recurrence a_n = c1a_{n-1} + c2a_{n-2} + ... + ck*a_{n-k} can be represented with a k x k companion matrix: M = [[c1, c2, ..., ck-1, ck], [1, 0, ..., 0, 0], [0, 1, ..., 0, 0], ... [0, 0, ..., 1, 0]]

The state vector S(n) = [a_n, a_{n-1}, ..., a_{n-k+1}]^T. Then S(n) = M * S(n-1).

For example, Tribonacci: a_n = a_{n-1} + a_{n-2} + a_{n-3} with initial a_0=0, a_1=0, a_2=1. Then M = [[1,1,1],[1,0,0],[0,1,0]]. S(2) = [1,0,0]^T. To compute a_n, compute M^(n-2) * S(2) and take the first element.

def tribonacci(n, mod=None):
    if n == 0:
        return 0
    if n == 1 or n == 2:
        return 0 if n==0 else 1 if n==2 else 0
    M = [[1, 1, 1],
         [1, 0, 0],
         [0, 1, 0]]
    # Compute M^(n-2)
    result = matrix_pow(M, n-2, mod)
    # Initial vector S(2) = [1,0,0]
    # a_n = result[0][0]*1 + result[0][1]*0 + result[0][2]*0 = result[0][0]
    return result[0][0] % mod if mod else result[0][0]

print(tribonacci(10))  # 81

Performance Optimization and Modulo Arithmetic

Matrix multiplication is O(k^3) per multiplication, and we do O(log n) multiplications, so total O(k^3 log n). For k up to 100, this is feasible. But for larger k, consider using Strassen's algorithm (O(k^2.81)) or other optimizations.

Modulo arithmetic is crucial to keep numbers small. Use modulo after each addition and multiplication to avoid overflow. In Python, integers are big but modulo prevents them from growing unboundedly.

Another optimization: use list comprehensions and local variable bindings to speed up loops. Also, consider using NumPy for large matrices, but be careful with integer overflow in NumPy (use dtype=object for arbitrary precision).

def matrix_mult_opt(A, B, mod):
    n = len(A)
    C = [[0]*n for _ in range(n)]
    for i in range(n):
        Ai = A[i]
        Ci = C[i]
        for k in range(n):
            aik = Ai[k]
            if aik:
                Bk = B[k]
                for j in range(n):
                    Ci[j] = (Ci[j] + aik * Bk[j]) % mod
    return C

Common Pitfalls and How to Avoid Them

1. Off-by-one errors: Ensure the exponent is correct. For Fibonacci, F(0)=0, F(1)=1, so to compute F(n), we need M^(n-1). For n=0, handle separately.
2. Modulo mistakes: Apply modulo consistently. If mod is None, skip modulo but be aware of huge numbers.
3. Matrix dimensions: Ensure the transformation matrix is k x k and the state vector is k x 1.
4. Recursive exponentiation: Use iterative to avoid stack overflow.
5. Integer overflow in other languages: In Python, not an issue, but in C++/Java, use long long and modulo.
6. Incorrect initial state: Double-check initial values and their order in the state vector.

# Mistake: Using wrong exponent for n=0
def fib_wrong(n):
    if n == 0:
        return 0
    M = [[1,1],[1,0]]
    result = matrix_pow(M, n)  # Should be n-1
    return result[0][0]

print(fib_wrong(1))  # Should be 1, but gives 2?

Real-World Applications Beyond Fibonacci

For example, counting the number of ways to tile a 2xn board with dominoes follows a Fibonacci-like recurrence. Matrix exponentiation can compute the number for large n quickly.

# Number of ways to tile 2xn board with 2x1 dominoes: F(n) = F(n-1) + F(n-2)
# So same as Fibonacci.
def tiling_ways(n, mod):
    return fib_matrix(n+1, mod)  # Because F(1)=1, F(2)=2, etc.

print(tiling_ways(100, 10**9+7))