Image by Editor | Midjourney

A stochastic matrix is a square matrix used in probability theory to represent state transitions in a Markov chain. Each element in the matrix represents a probability, and the rows (or columns) sum to 1.

Stochastic matrices are used in Markov chains, Google PageRank, and reinforcement learning. NumPy makes it easy to perform operations like multiplication, eigenvalue computation, and inverse.

This article will guide you through using NumPy to perform operations on stochastic matrices.

Creating a Stochastic Matrix

Right Stochastic Matrix

A right stochastic matrix is a square matrix in which each row sums to 1. Each row shows the probability of moving from one state to another.

To generate a right stochastic matrix, we can create a random matrix and normalize its rows so that each row sums to 1.

import numpy as np

def create_right_stochastic_matrix(size):
    matrix = np.random.rand(size, size)  # Generate a random matrix
    matrix /= matrix.sum(axis=1, keepdims=True)  # Normalize rows to sum to 1
    return matrix

size = 3  # Example size
stochastic_matrix = create_right_stochastic_matrix(size)
print("Right Stochastic Matrix:")
print(stochastic_matrix)

Output:

Right Stochastic Matrix:
[[0.37930486 0.12934182 0.49135332]
 [0.06001697 0.42928354 0.51069948]
 [0.60501106 0.37961044 0.0153785 ]]

Left Stochastic Matrix

A left stochastic matrix is a square matrix in which each column sums to 1. This means each column shows a probability distribution. This type of matrix is used in scenarios where transitions are described in terms of incoming probabilities rather than outgoing probabilities.

We can create a left stochastic matrix by normalizing the columns of a randomly generated matrix:

def create_left_stochastic_matrix(size):
    matrix = np.random.rand(size, size)
    matrix /= matrix.sum(axis=0, keepdims=True)  # Normalize columns to sum to 1
    return matrix

left_stochastic_matrix = create_left_stochastic_matrix(size)
print("Left Stochastic Matrix:")
print(left_stochastic_matrix)

Output:

Left Stochastic Matrix:
[[0.40377905 0.3474779  0.40288629]
 [0.4770592  0.33867685 0.44348314]
 [0.11916175 0.31384525 0.15363057]]

Doubly Stochastic Matrix

A doubly stochastic matrix is a matrix where both the rows and columns sum to 1.

One method to generate a doubly stochastic matrix is the Sinkhorn-Knopp algorithm. This iteratively normalizes the rows and columns until both sum to 1.

def create_doubly_stochastic_matrix(size, iterations=10):
    matrix = np.random.rand(size, size)
    for _ in range(iterations):
        matrix /= matrix.sum(axis=1, keepdims=True)  # Normalize rows
        matrix /= matrix.sum(axis=0, keepdims=True)  # Normalize columns
    return matrix

doubly_stochastic_matrix = create_doubly_stochastic_matrix(size)
print("Doubly Stochastic Matrix:")
print(doubly_stochastic_matrix)

Output:

Doubly Stochastic Matrix:
[[0.22041348 0.56070786 0.2188787 ]
 [0.50684479 0.0434798  0.44967536]
 [0.27274173 0.39581234 0.33144595]]

Verifying a Stochastic Matrix

Before performing operations, it is useful to verify whether a given matrix is stochastic. We check if all elements are non-negative and whether the rows (or columns) sum to 1.

def is_right_stochastic(matrix):
    return np.allclose(matrix.sum(axis=1), 1) and np.all(matrix >= 0)

def is_left_stochastic(matrix):
    return np.allclose(matrix.sum(axis=0), 1) and np.all(matrix >= 0)

def is_doubly_stochastic(matrix):
    return is_right_stochastic(matrix) and is_left_stochastic(matrix)

print("Is Right Stochastic:", is_right_stochastic(stochastic_matrix))
print("Is Left Stochastic:", is_left_stochastic(left_stochastic_matrix))
print("Is Doubly Stochastic:", is_doubly_stochastic(doubly_stochastic_matrix))

Output:

Is Right Stochastic: True
Is Left Stochastic: True
Is Doubly Stochastic: True

Operations on Stochastic Matrices

Matrix Multiplication

Multiplying a stochastic matrix with itself simulates multiple transitions in a Markov process. It helps analyze long-term behavior. This operation predicts how probabilities change over time.

The function np.dot() computes the matrix product.

matrix_product = np.dot(stochastic_matrix, stochastic_matrix)
print("Matrix Multiplication Result:")
print(matrix_product)

Output:

Matrix Multiplication Result:
[[0.44890907 0.29110715 0.25998378]
 [0.35750786 0.38591392 0.25657822]
 [0.26157087 0.24705159 0.49137755]]

Eigenvalues and Eigenvectors

The dominant eigenvector of a stochastic matrix gives the steady-state distribution. It helps find long-term probabilities. The largest eigenvalue is always 1. Its corresponding eigenvector represents the equilibrium state.

The function np.linalg.eig() computes eigenvalues and eigenvectors.

eigenvalues, eigenvectors = np.linalg.eig(stochastic_matrix.T)
steady_vector = eigenvectors[:, np.argmax(np.abs(eigenvalues))].real
steady_vector /= steady_vector.sum()
print("Steady-State Vector from Eigenvalues:")
print(steady_vector)

Output:

Steady-State Vector from Eigenvalues:
[0.35789897 0.30519978 0.33690125]

Inverse of a Stochastic Matrix

The inverse of a stochastic matrix gives insights into reversibility in Markov processes. It helps solve linear systems with probability distributions.

The function np.linalg.inv() calculates the inverse of the given stochastic matrix.

inverse_matrix = np.linalg.inv(stochastic_matrix)
print("Inverse of Stochastic Matrix:")
print(inverse_matrix)

Output:

Inverse of Stochastic Matrix:
[[ 1.26867847 -1.25017419  0.98149573]
 [-2.08700824  1.97444668  1.11256156]
 [ 1.60520243  0.44533862 -1.05054105]]

Computing the n-th Power of a Stochastic Matrix

Raising a stochastic matrix to a power simulates multiple steps in a Markov process. As n increases, the matrix may converge to a steady-state distribution. This depends on the matrix’s properties.

The function np.linalg.matrix_power() computes the power of the matrix.

def matrix_power(matrix, n):
    return np.linalg.matrix_power(matrix, n)

n = 5  # Example power
nth_power_matrix = matrix_power(stochastic_matrix, n)
print(f"{n}-th Power of Stochastic Matrix:")
print(nth_power_matrix)

Output:

5-th Power of Stochastic Matrix:
[[0.35323727 0.30137311 0.34538962]
 [0.35006877 0.30383606 0.34609518]
 [0.36994461 0.31050036 0.31955504]]

Conclusion

Stochastic matrices model probability transitions. Right stochastic matrices have row sums of 1, while left stochastic matrices have column sums of 1. Doubly stochastic matrices satisfy both conditions. These matrices are essential in Markov processes and steady-state analysis. Operations like multiplication, eigen decomposition, and exponentiation help study long-term behavior. Verifying stochastic properties ensures correctness. Stochastic matrices play an important role in probability, machine learning, and economics.