CODEWITHSUNDEEP
pythonnumpymathscipyhidden-markov-models
Halsey
Halsey

Reputation: 129

Compute null space of a matrix in Python

I'm trying to understand how to implement equation (6) in this paper (in page 5). The equation is this:

equation6

With C, P are the emission matrix and state transition matrix in a HMM. The goal is to compute N.

Let's take the example 4 in this paper.
example input

The result in the paper is given as:
enter image description here

Here's the code I used:

import numpy as np
from scipy.linalg import null_space

n = 4
P = np.array([[1/2, 0, 1/3, 1/4], [0, 1/3, 1/3, 1/4], [1/2, 0, 1/3, 0], [0, 2/3, 0, 1/2]])
C = np.array([[1/4, 1/4, 1/2, 7/16], [3/4, 3/4, 1/2, 9/16]])

N = C
for i in range(1, n):
    N = np.vstack((N, np.dot(C, np.linalg.matrix_power(P, i))))

ns = null_space(N)
print(ns)

The result from the code is this:

[[ 0.71390642 -0.01189748]
 [-0.68971636 -0.18464902]
 [ 0.07257016 -0.58963951]
 [-0.09676022  0.78618601]]

Both spans are equal but I want to know if there is anyway to get result in the paper using python.

Upvotes: 1

Views: 172

Answers (1)

Serge de Gosson de Varennes
Serge de Gosson de Varennes

Reputation: 11532

Following what the paper does, your code should be

from scipy.linalg import svd

def compute_null_space(C, P, n):
    N = C
    for i in range(1, n):
        N = np.vstack((N, np.dot(C, np.linalg.matrix_power(P, i))))
    U, s, Vh = svd(N)
    tol = max(N.shape) * np.max(s) * np.finfo(s.dtype).eps
    null_mask = (s <= tol)
    null_space = Vh[null_mask, :]
    return null_space.T

P = np.array([[1/2, 0, 1/3, 1/4], [0, 1/3, 1/3, 1/4], [1/2, 0, 1/3, 0], [0, 2/3, 0, 1/2]])
C = np.array([[1/4, 1/4, 1/2, 7/16], [3/4, 3/4, 1/2, 9/16]])


ns_paper_method = compute_null_space(C, P, n=4)
ns_paper_method

which return

array([[ 0.16383871,  0.69495381],
       [-0.34842041, -0.62322319],
       [-0.5537451 ,  0.21519187],
       [ 0.7383268 , -0.2869225 ]])

Upvotes: 0

Related Questions