Skip to article frontmatterSkip to article content
Contents
and

38. The Perron-Frobenius Theorem

In addition to what’s in Anaconda, this lecture will need the following libraries:

!pip install quantecon
Output
Requirement already satisfied: quantecon in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (0.8.0)
Requirement already satisfied: numba>=0.49.0 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from quantecon) (0.61.0)
Requirement already satisfied: numpy>=1.17.0 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from quantecon) (2.1.3)
Requirement already satisfied: requests in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from quantecon) (2.32.3)

Requirement already satisfied: scipy>=1.5.0 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from quantecon) (1.15.2)
Requirement already satisfied: sympy in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from quantecon) (1.13.3)
Requirement already satisfied: llvmlite<0.45,>=0.44.0dev0 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from numba>=0.49.0->quantecon) (0.44.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from requests->quantecon) (3.4.1)
Requirement already satisfied: idna<4,>=2.5 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from requests->quantecon) (3.10)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from requests->quantecon) (2.3.0)
Requirement already satisfied: certifi>=2017.4.17 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from requests->quantecon) (2025.1.31)
Requirement already satisfied: mpmath<1.4,>=1.1.0 in /opt/hostedtoolcache/Python/3.12.9/x64/lib/python3.12/site-packages (from sympy->quantecon) (1.3.0)

In this lecture we will begin with the foundational concepts in spectral theory.

Then we will explore the Perron-Frobenius theorem and connect it to applications in Markov chains and networks.

We will use the following imports:

import numpy as np
from numpy.linalg import eig
import scipy as sp
import quantecon as qe

38.1Nonnegative matrices

Often, in economics, the matrix that we are dealing with is nonnegative.

Nonnegative matrices have several special and useful properties.

In this section we will discuss some of them --- in particular, the connection between nonnegativity and eigenvalues.

An n×mn \times m matrix AA is called nonnegative if every element of AA is nonnegative, i.e., aij0a_{ij} \geq 0 for every i,ji,j.

We denote this as A0A \geq 0.

38.1.1Irreducible matrices

We introduced irreducible matrices in the Markov chain lecture.

Here we generalize this concept:

Let aijka^{k}_{ij} be element (i,j)(i,j) of AkA^k.

An n×nn \times n nonnegative matrix AA is called irreducible if A+A2+A3+0A + A^2 + A^3 + \cdots \gg 0, where 0\gg 0 indicates that every element in AA is strictly positive.

In other words, for each i,ji,j with 1i,jn1 \leq i, j \leq n, there exists a k0k \geq 0 such that aijk>0a^{k}_{ij} > 0.

38.1.2Left eigenvectors

Recall that we previously discussed eigenvectors in Eigenvalues and Eigenvectors.

In particular, λ is an eigenvalue of AA and vv is an eigenvector of AA if vv is nonzero and satisfy

Av=λv.Av = \lambda v.

In this section we introduce left eigenvectors.

To avoid confusion, what we previously referred to as “eigenvectors” will be called “right eigenvectors”.

Left eigenvectors will play important roles in what follows, including that of stochastic steady states for dynamic models under a Markov assumption.

A vector ww is called a left eigenvector of AA if ww is a right eigenvector of AA^\top.

In other words, if ww is a left eigenvector of matrix AA, then Aw=λwA^\top w = \lambda w, where λ is the eigenvalue associated with the left eigenvector vv.

This hints at how to compute left eigenvectors

A = np.array([[3, 2],
              [1, 4]])

# Compute eigenvalues and right eigenvectors
λ, v = eig(A)

# Compute eigenvalues and left eigenvectors
λ, w = eig(A.T)

# Keep 5 decimals
np.set_printoptions(precision=5)

print(f"The eigenvalues of A are:\n {λ}\n")
print(f"The corresponding right eigenvectors are: \n {v[:,0]} and {-v[:,1]}\n")
print(f"The corresponding left eigenvectors are: \n {w[:,0]} and {-w[:,1]}\n")
The eigenvalues of A are:
 [2. 5.]

The corresponding right eigenvectors are: 
 [-0.89443  0.44721] and [0.70711 0.70711]

The corresponding left eigenvectors are: 
 [-0.70711  0.70711] and [0.44721 0.89443]

We can also use scipy.linalg.eig with argument left=True to find left eigenvectors directly

eigenvals, ε, e = sp.linalg.eig(A, left=True)

print(f"The eigenvalues of A are:\n {eigenvals.real}\n")
print(f"The corresponding right eigenvectors are: \n {e[:,0]} and {-e[:,1]}\n")
print(f"The corresponding left eigenvectors are: \n {ε[:,0]} and {-ε[:,1]}\n")
The eigenvalues of A are:
 [2. 5.]

The corresponding right eigenvectors are: 
 [-0.89443  0.44721] and [0.70711 0.70711]

The corresponding left eigenvectors are: 
 [-0.70711  0.70711] and [0.44721 0.89443]

The eigenvalues are the same while the eigenvectors themselves are different.

(Also note that we are taking the nonnegative value of the eigenvector of dominant eigenvalue, this is because eig automatically normalizes the eigenvectors.)

We can then take transpose to obtain Aw=λwA^\top w = \lambda w and obtain wA=λww^\top A= \lambda w^\top.

This is a more common expression and where the name left eigenvectors originates.

38.1.3The Perron-Frobenius theorem

For a square nonnegative matrix AA, the behavior of AkA^k as kk \to \infty is controlled by the eigenvalue with the largest absolute value, often called the dominant eigenvalue.

For any such matrix AA, the Perron-Frobenius theorem characterizes certain properties of the dominant eigenvalue and its corresponding eigenvector.

(This is a relatively simple version of the theorem --- for more details see here).

We will see applications of the theorem below.

Let’s build our intuition for the theorem using a simple example we have seen before.

Now let’s consider examples for each case.

38.1.3.1Example: irreducible matrix

Consider the following irreducible matrix AA:

A = np.array([[0, 1, 0],
              [.5, 0, .5],
              [0, 1, 0]])

We can compute the dominant eigenvalue and the corresponding eigenvector

eig(A)
EigResult(eigenvalues=array([-1.00000e+00, 2.90566e-17, 1.00000e+00]), eigenvectors=array([[ 5.77350e-01, 7.07107e-01, 5.77350e-01], [-5.77350e-01, 1.35087e-16, 5.77350e-01], [ 5.77350e-01, -7.07107e-01, 5.77350e-01]]))

Now we can see the claims of the Perron-Frobenius theorem holds for the irreducible matrix AA:

  1. The dominant eigenvalue is real-valued and non-negative.
  2. All other eigenvalues have absolute values less than or equal to the dominant eigenvalue.
  3. A non-negative and nonzero eigenvector is associated with the dominant eigenvalue.
  4. As the matrix is irreducible, the eigenvector associated with the dominant eigenvalue is strictly positive.
  5. There exists no other positive eigenvector associated with the dominant eigenvalue.

38.1.4Primitive matrices

We know that in real world situations it’s hard for a matrix to be everywhere positive (although they have nice properties).

The primitive matrices, however, can still give us helpful properties with looser definitions.

Let AA be a square nonnegative matrix and let AkA^k be the kthk^{th} power of AA.

A matrix is called primitive if there exists a kNk \in \mathbb{N} such that AkA^k is everywhere positive.

We can see that if a matrix is primitive, then it implies the matrix is irreducible but not vice versa.

Now let’s step back to the primitive matrices part of the Perron-Frobenius theorem

38.1.4.1Example 1: primitive matrix

Consider the following primitive matrix BB:

B = np.array([[0, 1, 1],
              [1, 0, 1],
              [1, 1, 0]])

np.linalg.matrix_power(B, 2)
array([[2, 1, 1], [1, 2, 1], [1, 1, 2]])

We compute the dominant eigenvalue and the corresponding eigenvector

eig(B)
EigResult(eigenvalues=array([-1., 2., -1.]), eigenvectors=array([[-0.8165 , 0.57735, 0.22646], [ 0.40825, 0.57735, -0.79259], [ 0.40825, 0.57735, 0.56614]]))

Now let’s give some examples to see if the claims of the Perron-Frobenius theorem hold for the primitive matrix BB:

  1. The dominant eigenvalue is real-valued and non-negative.
  2. All other eigenvalues have absolute values strictly less than the dominant eigenvalue.
  3. A non-negative and nonzero eigenvector is associated with the dominant eigenvalue.
  4. The eigenvector associated with the dominant eigenvalue is strictly positive.
  5. There exists no other positive eigenvector associated with the dominant eigenvalue.
  6. The inequality λ<r(B)|\lambda| < r(B) holds for all eigenvalues λ of BB distinct from the dominant eigenvalue.

Furthermore, we can verify the convergence property (7) of the theorem on the following examples:

def compute_perron_projection(M):

    eigval, v = eig(M)
    eigval, w = eig(M.T)

    r = np.max(eigval)

    # Find the index of the dominant (Perron) eigenvalue
    i = np.argmax(eigval)

    # Get the Perron eigenvectors
    v_P = v[:, i].reshape(-1, 1)
    w_P = w[:, i].reshape(-1, 1)

    # Normalize the left and right eigenvectors
    norm_factor = w_P.T @ v_P
    v_norm = v_P / norm_factor

    # Compute the Perron projection matrix
    P = v_norm @ w_P.T
    return P, r

def check_convergence(M):
    P, r = compute_perron_projection(M)
    print("Perron projection:")
    print(P)

    # Define a list of values for n
    n_list = [1, 10, 100, 1000, 10000]

    for n in n_list:

        # Compute (A/r)^n
        M_n = np.linalg.matrix_power(M/r, n)

        # Compute the difference between A^n / r^n and the Perron projection
        diff = np.abs(M_n - P)

        # Calculate the norm of the difference matrix
        diff_norm = np.linalg.norm(diff, 'fro')
        print(f"n = {n}, error = {diff_norm:.10f}")


A1 = np.array([[1, 2],
               [1, 4]])

A2 = np.array([[0, 1, 1],
               [1, 0, 1],
               [1, 1, 0]])

A3 = np.array([[0.971, 0.029, 0.1, 1],
               [0.145, 0.778, 0.077, 0.59],
               [0.1, 0.508, 0.492, 1.12],
               [0.2, 0.8, 0.71, 0.95]])

for M in A1, A2, A3:
    print("Matrix:")
    print(M)
    check_convergence(M)
    print()
    print("-"*36)
    print()
Matrix:
[[1 2]
 [1 4]]
Perron projection:
[[0.1362  0.48507]
 [0.24254 0.8638 ]]
n = 1, error = 0.0989045731
n = 10, error = 0.0000000001
n = 100, error = 0.0000000000
n = 1000, error = 0.0000000000
n = 10000, error = 0.0000000000

------------------------------------

Matrix:
[[0 1 1]
 [1 0 1]
 [1 1 0]]
Perron projection:
[[0.33333 0.33333 0.33333]
 [0.33333 0.33333 0.33333]
 [0.33333 0.33333 0.33333]]
n = 1, error = 0.7071067812
n = 10, error = 0.0013810679
n = 100, error = 0.0000000000
n = 1000, error = 0.0000000000
n = 10000, error = 0.0000000000

------------------------------------

Matrix:
[[0.971 0.029 0.1   1.   ]
 [0.145 0.778 0.077 0.59 ]
 [0.1   0.508 0.492 1.12 ]
 [0.2   0.8   0.71  0.95 ]]
Perron projection:
[[0.12506 0.31949 0.20233 0.43341]
 [0.07714 0.19707 0.1248  0.26735]
 [0.12158 0.31058 0.19669 0.42133]
 [0.13885 0.3547  0.22463 0.48118]]
n = 1, error = 0.5361031549
n = 10, error = 0.0000434043
n = 100, error = 0.0000000000
n = 1000, error = 0.0000000000
n = 10000, error = 0.0000000000

------------------------------------

The convergence is not observed in cases of non-primitive matrices.

Let’s go through an example

B = np.array([[0, 1, 1],
              [1, 0, 0],
              [1, 0, 0]])

# This shows that the matrix is not primitive
print("Matrix:")
print(B)
print("100th power of matrix B:")
print(np.linalg.matrix_power(B, 100))

check_convergence(B)
Matrix:
[[0 1 1]
 [1 0 0]
 [1 0 0]]
100th power of matrix B:
[[1125899906842624                0                0]
 [               0  562949953421312  562949953421312]
 [               0  562949953421312  562949953421312]]
Perron projection:
[[0.5     0.35355 0.35355]
 [0.35355 0.25    0.25   ]
 [0.35355 0.25    0.25   ]]
n = 1, error = 1.0000000000
n = 10, error = 1.0000000000
n = 100, error = 1.0000000000
n = 1000, error = 1.0000000000
n = 10000, error = 1.0000000000

The result shows that the matrix is not primitive as it is not everywhere positive.

These examples show how the Perron-Frobenius theorem relates to the eigenvalues and eigenvectors of positive matrices and the convergence of the power of matrices.

In fact we have already seen the theorem in action before in the Markov chain lecture.

38.1.4.2Example 2: connection to Markov chains

We are now prepared to bridge the languages spoken in the two lectures.

A primitive matrix is both irreducible and aperiodic.

So Perron-Frobenius theorem explains why both Imam and Temple matrix and Hamilton matrix converge to a stationary distribution, which is the Perron projection of the two matrices

P = np.array([[0.68, 0.12, 0.20],
              [0.50, 0.24, 0.26],
              [0.36, 0.18, 0.46]])

print(compute_perron_projection(P)[0])
[[0.56146 0.15565 0.28289]
 [0.56146 0.15565 0.28289]
 [0.56146 0.15565 0.28289]]
mc = qe.MarkovChain(P)
ψ_star = mc.stationary_distributions[0]
ψ_star
array([0.56146, 0.15565, 0.28289])
P_hamilton = np.array([[0.971, 0.029, 0.000],
                       [0.145, 0.778, 0.077],
                       [0.000, 0.508, 0.492]])

print(compute_perron_projection(P_hamilton)[0])
[[0.8128  0.16256 0.02464]
 [0.8128  0.16256 0.02464]
 [0.8128  0.16256 0.02464]]
mc = qe.MarkovChain(P_hamilton)
ψ_star = mc.stationary_distributions[0]
ψ_star
array([0.8128 , 0.16256, 0.02464])

We can also verify other properties hinted by Perron-Frobenius in these stochastic matrices.

Another example is the relationship between convergence gap and convergence rate.

In the exercise, we stated that the convergence rate is determined by the spectral gap, the difference between the largest and the second largest eigenvalue.

This can be proven using what we have learned here.

Please note that we use 1\mathbb{1} for a vector of ones in this lecture.

With Markov model MM with state space SS and transition matrix PP, we can write PtP^t as

Pt=i=1n1λitviwi+1ψ,P^t=\sum_{i=1}^{n-1} \lambda_i^t v_i w_i^{\top}+\mathbb{1} \psi^*,

This is proven in Sargent & Stachurski (2023) and a nice discussion can be found here.

In this formula λi\lambda_i is an eigenvalue of PP with corresponding right and left eigenvectors viv_i and wiw_i .

Premultiplying PtP^t by arbitrary ψD(S)\psi \in \mathscr{D}(S) and rearranging now gives

ψPtψ=i=1n1λitψviwi\psi P^t-\psi^*=\sum_{i=1}^{n-1} \lambda_i^t \psi v_i w_i^{\top}

Recall that eigenvalues are ordered from smallest to largest from i=1...ni = 1 ... n.

As we have seen, the largest eigenvalue for a primitive stochastic matrix is one.

This can be proven using Gershgorin Circle Theorem, but it is out of the scope of this lecture.

So by the statement (6) of Perron-Frobenius theorem, λi<1\lambda_i<1 for all i<ni<n, and λn=1\lambda_n=1 when PP is primitive.

Hence, after taking the Euclidean norm deviation, we obtain

ψPtψ=O(ηt) where η:=λn1<1\left\|\psi P^t-\psi^*\right\|=O\left(\eta^t\right) \quad \text { where } \quad \eta:=\left|\lambda_{n-1}\right|<1

Thus, the rate of convergence is governed by the modulus of the second largest eigenvalue.

38.2Exercises

Solution to Exercise 1
A = np.array([[0.3, 0.2, 0.3],
              [0.2, 0.4, 0.3],
              [0.2, 0.5, 0.1]])

evals, evecs = eig(A)

r = max(abs(λ) for λ in evals)   #dominant eigenvalue/spectral radius
print(r)
0.8444086477164554

Since we have r(A)<1r(A) < 1 we can thus find the solution using the Neumann Series Lemma.

I = np.identity(3)
B = I - A

d = np.array([4, 5, 12])
d.shape = (3,1)

B_inv = np.linalg.inv(B)
x_star = B_inv @ d
print(x_star)
[[38.30189]
 [44.33962]
 [46.47799]]
References
  1. Sargent, T. J., & Stachurski, J. (2023). Economic Networks: Theory and Computation. arXiv Preprint arXiv:2203.11972.
CC-BY-SA-4.0

Creative Commons License – This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International.