Program works correcly only with certain initial conditions

user_2297237 0
Viewed 22

A square matrix is called stochastic if its entries are non-negative and the sum of each row is 1. For example, the following matrix is stochastic:

[0.3, 0.7]

[0.9, 0.1]

It is well known that any power of a stochastic matrix is still stochastic. I wrote a program to verify this fact, but I've ran into an issue: the program works for some stochastic matrices, but not others. Here is my code:

import numpy as np

def main(p):
    # Stochastic matrix A
    A = np.array([[0.3, 0.7], 
                  [0.2, 0.8]])
    
    # Identity matrix
    B = np.array([[1, 0], 
                  [0, 1]])
    
    # Check if powers of A are stochastic
    for i in range(p):
        B = B @ A
        print(B)
        if (sum(B[0,]) == 1) and (sum(B[1,]) == 1):
            print("Good")
        else:
            print("Bad")
    
    return(0)
    
main(5)

For this particular A, the program works just fine; it prints out "Good" every time, and I can check by hand the all the powers of A are stochastic.

However, if we instead use

A = np.array([[0.3, 0.7], 
              [0.5, 0.5]])

then something breaks; when I check manually, the powers of A are still stochastic. But the program prints out one "Good" and four "Bads". Could you please help me figure out what the issue is?

linear-algebra numpy python
1 Answers

I believe it's floating point precision. Lesson learned: avoid comparing two floats with ==.

Numpy has allclose method for this:

    if np.allclose(B.sum(axis=1),1):
        print("Good")
    else:
        print("Bad")
user_42384080
Related