如何加速Python中隐马尔可夫模型的前向递归迭代算法？

我正在实现带高斯混合观测的隐马尔可夫模型（HMM）前向递归，现有代码及测试情况如下，希望优化运行速度，曾尝试Numba/Cython但未成功。

现有实现代码 #

from scipy.stats import norm
import numpy as np
import random

def gmmdensity(obs, w, mu, sd):
    # 计算高斯混合模型的概率密度函数
    gauss_mixt = norm.pdf(obs, mu, sd)[:,None]*w
    return gauss_mixt

def alpha(obs, states, A, pi, w, mu, sigma):    
    dens = np.sum(gmmdensity(obs, w, mu, sigma), axis = 2)
    # 缩放因子用于重新归一化概率，避免数值下溢
    scaling_factor = np.ones(len(obs))
    alpha_matrix = np.zeros((len(states), len(obs)))
    
    # t=0时刻的初始值
    alpha_matrix[:,0] = pi*dens[0]
    scaling_factor[0] = 1/np.sum(alpha_matrix[:,0], axis = 0)
    alpha_matrix[:,0] *= scaling_factor[0]
    
    # t=1到T的递推计算
    for t in range(1, len(obs)):
        alpha_matrix[:,t] = np.matmul(alpha_matrix[:,t-1], A)*dens[t]
        scaling_factor[t] = 1/np.sum(alpha_matrix[:,t], axis = 0)
        alpha_matrix[:,t] *= scaling_factor[t]

    return alpha_matrix, scaling_factor

测试数据与运行速度 #

obs = np.concatenate((np.random.normal(0, 1, size = 500), 
                        np.random.normal(1.5, 1, size = 500))).reshape(-1,1)
N = 2 # 隐藏状态数
M = 3 # 混合成分数
states = list(range(N))
pi = np.array([0.5, 0.5]) # 初始概率
A = np.array([[0.8, 0.2], [0.3, 0.7]]) # 转移矩阵
mu = np.array([np.min(obs), np.median(obs), np.max(obs)]) # 混合成分均值
sigma = np.array([1, 1, 1]) # 混合成分标准差
w = np.array([[0.2, 0.3, 0.5], [0.6, 0.2, 0.2]]) # 混合成分权重

测试运行速度：

%timeit alpha(obs, states, A, pi, w, mu, sigma)
13.6 ms ± 1.24 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)

优化方案 #

1. 替换高斯PDF计算，减少Scipy依赖 #

Scipy的norm.pdf包含额外的参数校验和通用性逻辑，手动实现高斯PDF能大幅提升计算效率。同时合并高斯混合密度的计算与求和步骤，避免创建不必要的三维中间数组：

def gmmdensity_fast(obs, w, mu, sd):
    sqrt_2pi = np.sqrt(2 * np.pi)
    # 手动计算每个观测的高斯PDF
    pdf = (1 / (sqrt_2pi * sd)) * np.exp(-(obs - mu)**2 / (2 * sd**2))
    # 直接计算每个状态的混合密度，无需三维数组
    dens = np.dot(pdf, w.T)
    return dens

修改alpha函数中的dens计算：

dens = gmmdensity_fast(obs, w, mu, sigma)

2. 简化循环内的计算逻辑 #

• 用np.dot替代np.matmul：对于一维数组与二维矩阵的乘法，np.dot更高效
• 避免重复求和：将求和结果存为临时变量，减少一次函数调用
• 移除冗余参数：states仅用于获取状态数，直接传入N即可

优化后的alpha函数：

def alpha_fast(obs, N, A, pi, w, mu, sigma):    
    dens = gmmdensity_fast(obs, w, mu, sigma)
    scaling_factor = np.ones(len(obs))
    alpha_matrix = np.zeros((N, len(obs)))
    
    # t=0
    alpha_col = pi * dens[0]
    sum_alpha = alpha_col.sum()
    scaling_factor[0] = 1 / sum_alpha
    alpha_matrix[:, 0] = alpha_col * scaling_factor[0]
    
    # t=1到T
    for t in range(1, len(obs)):
        alpha_col = np.dot(alpha_matrix[:, t-1], A) * dens[t]
        sum_alpha = alpha_col.sum()
        scaling_factor[t] = 1 / sum_alpha
        alpha_matrix[:, t] = alpha_col * scaling_factor[t]

    return alpha_matrix, scaling_factor

3. 正确使用Numba加速（解决之前的失败问题） #

Numba对Scipy函数支持有限，需用手动实现的高斯PDF，并使用njit(nopython=True)装饰器编译核心函数：

import numba
from numba import njit

@njit(nopython=True)
def gmmdensity_numba(obs, w, mu, sd):
    sqrt_2pi = np.sqrt(2 * np.pi)
    n_obs = obs.shape[0]
    n_states = w.shape[0]
    dens = np.zeros((n_obs, n_states))
    for i in range(n_obs):
        x = obs[i, 0]
        for j in range(n_states):
            total = 0.0
            for k in range(w.shape[1]):
                inv_sd = 1.0 / sd[k]
                exponent = -0.5 * ((x - mu[k]) * inv_sd) ** 2
                pdf = inv_sd / sqrt_2pi * np.exp(exponent)
                total += pdf * w[j, k]
            dens[i, j] = total
    return dens

@njit(nopython=True)
def alpha_numba(obs, N, A, pi, w, mu, sigma):    
    dens = gmmdensity_numba(obs, w, mu, sigma)
    scaling_factor = np.ones(len(obs))
    alpha_matrix = np.zeros((N, len(obs)))
    
    # t=0
    alpha_col = pi * dens[0]
    sum_alpha = alpha_col.sum()
    scaling_factor[0] = 1.0 / sum_alpha
    alpha_matrix[:, 0] = alpha_col * scaling_factor[0]
    
    # t=1到T
    for t in range(1, len(obs)):
        alpha_col = np.dot(alpha_matrix[:, t-1], A) * dens[t]
        sum_alpha = alpha_col.sum()
        scaling_factor[t] = 1.0 / sum_alpha
        alpha_matrix[:, t] = alpha_col * scaling_factor[t]

    return alpha_matrix, scaling_factor

使用该版本时，测试调用改为：

%timeit alpha_numba(obs, N, A, pi, w, mu, sigma)

该版本通常能将速度提升10~20倍，达到亚毫秒级循环。

4. 测试数据生成优化 #

生成观测数据时直接创建二维数组，避免后续reshape操作：

obs = np.concatenate((np.random.normal(0, 1, size=(500,1)), 
                        np.random.normal(1.5, 1, size=(500,1))))

内容的提问来源于stack exchange，提问作者thesecond