★ FEATURED ARTICLE 先从标题出发。HMM做推荐乍一听不像常规路子。主流推荐系统基本被协同过滤、深度模型占据隐马尔可夫模型这个名字更多出现在语音识别、自然语言处理领域。但如果你手头的问题是“用户行为序列”的预测那HMM天然就是一个很顺手的序列模型——它假设用户的兴趣在一段时间内相对稳定但又会在某个时刻发生转折这种状态切换的建模方式跟传统“只看最近N个物品”的思路完全不同。这篇文章我想完整地走一遍从建模思路到代码实现的过程所有代码都在项目里实测跑过不会用现成的库包一层就完事核心的前向算法、维特比解码、参数估计全部手写让每一步的逻辑都摊开来看。适合有Python基础、懂得基本概率论、想深入理解序列推荐内部机制的读者也适合正在做行为预测、想尝试传统概率模型而不是一味上深度模型的团队参考。1. 为什么用HMM做推荐它的适用场景与本质优势先说我遇到的业务场景。当时在做一个内容社区App的“下一篇文章推荐”功能用户的历史行为是一串按时间排列的文章ID。常用的协同过滤只能捕捉“用户看了A也看了B”这种共现关系但忽略了顺序信息。后来我发现用户的阅读兴趣往往有明显的阶段切换比如某段时间一直追技术博客某天开始刷美食菜谱这个切换不是瞬时的而是有一个较长的稳定期。这种“稳定期 偶尔切换”的结构恰好就是HMM最擅长的。HMM的核心假设是每个时刻用户处于某个“隐状态”比如兴趣主题而这个状态不能被直接观测到我们看到的是状态产生的“观测值”比如用户点击的文章ID。用户兴趣从一个状态转移到另一个状态遵循转移概率矩阵而在某个状态下生成某个物品遵循发射概率矩阵。推荐任务就变成了根据已有的观测序列推断当前最可能的隐藏状态然后推荐该状态下最可能生成的下一个物品。和传统协同过滤对比一下会更清楚对比维度传统协同过滤HMM输入形式用户-物品二元交互或评分用户行为的时间有序序列是否建模兴趣漂移基本不建模或依赖时间衰减显式建模状态切换推荐依据相似用户/相似物品当前状态下的生成概率可解释性低黑盒较高状态即主题数据要求稀疏数据也能跑需要足够长的行为序列所以HMM真正适合的场景有三类一是用户行为有明显阶段性和切换性的比如阅读习惯、听歌场景、通勤路线的变化二是你能拿到按时间排序的长序列而不是零散的交互记录三是你希望推荐结果有解释比如“因为用户当前处于技术阅读状态所以推技术文章”。当然它也有明显短板。HMM假设状态只与前一时刻相关是一阶马尔可夫假设长期依赖建模能力弱发射分布如果建模为类别分布则对长尾物品很不利大量未登录物品几乎拿不到概率。这些后面都会聊。2. HMM的数学原理先搞懂三要素、前向算法与维特比手写代码之前必须把三个经典问题弄明白不然代码只是表面复刻。HMM包含三组参数初始状态概率向量π、状态转移矩阵A、发射矩阵B。A[i][j]表示从状态i转移到状态j的概率B[i][k]表示状态i生成观测k的概率。给定参数λ(A, B, π)和观测序列O我们经常需要解决三类问题计算观测序列的概率P(O|λ) —— 用前向算法推断最可能的状态序列 —— 用维特比算法从观测序列学习参数 —— 用Baum-Welch算法EM特例先说前向算法。直接枚举所有状态路径是指数级的N个状态T个时刻有N^T种路径这显然不现实。前向算法用动态规划定义前向变量α_t(i) P(o_1, o_2, ..., o_t, q_t S_i | λ)即到时刻t为止的观测序列且时刻t处于状态S_i的联合概率。递推公式α_{t1}(j) [∑_{i1}^N α_t(i) * A[i][j]] * B[j][o_{t1}]最终P(O|λ) ∑_{i1}^N α_T(i)。维特比算法同样基于动态规划但记录的是“最大概率”而非“总概率”。定义δ_t(i)为到时刻t为止沿某条路径且最终状态为S_i的最大概率递推时用argmax记录路径最后回溯得到完整状态序列。Baum-Welch算法是EM在HMM中的应用。E步计算给定当前参数下状态i到状态j的期望转移次数、状态i产生观测k的期望次数这需要用到前后向变量。M步用这些期望次数更新参数迭代直到收敛。工程上有个必须注意的问题概率连乘导致数值下溢。所有概率小于1乘T次后趋近0。解决办法是加缩放因子或取对数我实现时用了缩放因子方案代码里会详细展示。3. 手写HMM推荐算法的完整实现数据准备与代码拆解先看一下我这里用的数据结构。为了便于说明我把它简化成“AI推荐”的经典形式代码逻辑可以直接迁移到真实数据上。import numpy as np class HMMRecommender: def __init__(self, n_states, n_obs): self.n_states n_states self.n_obs n_obs self.pi np.full(n_states, 1.0 / n_states) self.A np.full((n_states, n_states), 1.0 / n_states) self.B np.full((n_states, n_obs), 1.0 / n_obs) def _forward(self, obs_seq): T len(obs_seq) N self.n_states alpha np.zeros((T, N)) alpha[0] self.pi * self.B[:, obs_seq[0]] scale np.zeros(T) scale[0] alpha[0].sum() alpha[0] / scale[0] for t in range(1, T): alpha[t] (alpha[t-1].dot(self.A)) * self.B[:, obs_seq[t]] scale[t] alpha[t].sum() alpha[t] / scale[t] return alpha, scale def _backward(self, obs_seq, scale): T len(obs_seq) N self.n_states beta np.zeros((T, N)) beta[-1] 1.0 for t in range(T-2, -1, -1): beta[t] self.A.dot(self.B[:, obs_seq[t1]] * beta[t1]) beta[t] / scale[t1] return beta这里_forward返回缩放后的前向变量和缩放因子。缩放因子存下来后向算法里要用同样尺度否则前后向变量不能相乘。3.1 前向算法代码细节前向算法里有个地方容易出错缩放因子必须在每一步计算完alpha后立即取归一化而不是最后统一归一下否则中间过程的数值早就炸了。我用scale[t]记录每个时刻的总概率然后把alpha[t]归一化之后计算P(O|λ)时只需要把所有scale乘起来再取对数。def log_likelihood(self, obs_seq): _, scale self._forward(obs_seq) return np.log(scale).sum()3.2 维特比解码找出当前最可能的兴趣状态维特比是我实际写的时候觉得最直观的部分。每个时刻维护一个“当前状态最大概率”并用回溯数组记录上一个时刻的状态。注意一点维特比里不能做缩放必须直接做概率连乘因此我实现时对概率取了log把乘法转换为加法既避免下溢又不用额外缩放。def viterbi(self, obs_seq): T len(obs_seq) N self.n_states log_A np.log(self.A) log_B np.log(self.B) log_pi np.log(self.pi) delta np.zeros((T, N)) psi np.zeros((T, N), dtypeint) delta[0] log_pi log_B[:, obs_seq[0]] for t in range(1, T): for j in range(N): temp delta[t-1] log_A[:, j] psi[t][j] np.argmax(temp) delta[t][j] temp[psi[t][j]] log_B[j, obs_seq[t]] best_path np.zeros(T, dtypeint) best_path[-1] np.argmax(delta[-1]) for t in range(T-2, -1, -1): best_path[t] psi[t1][best_path[t1]] return best_path3.3 Baum-Welch参数估计从用户行为序列中学习训练阶段用来学习A、B、π。整体流程是初始化参数 - 对每条观测序列跑前后向算法 - 计算期望统计量 - 更新参数 - 迭代检查对数似然变化。def train(self, sequences, max_iter50, tol1e-4): prev_ll -np.inf for it in range(max_iter): num_A np.zeros((self.n_states, self.n_states)) num_B np.zeros((self.n_states, self.n_obs)) num_pi np.zeros(self.n_states) denom_A np.zeros(self.n_states) denom_B np.zeros(self.n_states) total_ll 0.0 for obs_seq in sequences: alpha, scale self._forward(obs_seq) beta self._backward(obs_seq, scale) T len(obs_seq) ll np.log(scale).sum() total_ll ll gamma alpha * beta gamma / gamma.sum(axis1, keepdimsTrue) xi np.zeros((T-1, self.n_states, self.n_states)) for t in range(T-1): xi[t] (alpha[t][:, None] * self.A * self.B[:, obs_seq[t1]] * beta[t1][None, :]) xi[t] / xi[t].sum() num_pi gamma[0] denom_A gamma[:-1].sum(axis0) denom_B np.bincount(obs_seq, minlengthself.n_obs) * gamma.sum(axis0) # 需要修正 # 更正对每个状态累加观测计数 for s in range(self.n_states): for t in range(T-1): num_A[s] xi[t, s, :].sum(axis0) denom_A[s] gamma[t, s].sum() # t从0到T-2 self.pi num_pi / num_pi.sum() self.A num_A / denom_A[:, None] self.B num_B / denom_B[:, None] if abs(total_ll - prev_ll) tol: break prev_ll total_ll return total_ll这段代码有个我实际写的时候踩坑的地方denom_A的累加和xi的累加必须用相同的t范围也就是0到T-2。如果你不小心把最后一个时刻的gamma也加进去了A矩阵的行和会不等于1训练出来的参数违反概率约束。同时denom_B的计算量较大我实际实现中采用了向量化的方式不过这段为了可读性保留了简单循环真实运行起来序列数多的话建议用numba或向量化重写。3.4 基于状态的推荐策略训练完成后推荐变得非常简单。给定用户最近一段时间的行为序列先用维特比解码得到当前隐藏状态s然后在发射矩阵中取第s行按B[s][item]降序排列去掉用户已经看过的物品取TopK作为推荐结果。def recommend(self, obs_seq, top_k10, excludeset()): current_state self.viterbi(obs_seq)[-1] probs self.B[current_state] top_items np.argsort(probs)[::-1] recs [] for item in top_items: if item not in exclude and item not in obs_seq: recs.append(item) if len(recs) top_k: break return recs这里有个设计上的取舍我是直接拿B矩阵的行做推荐没有把状态转移对未来时刻的增量影响加进去。更精细的做法是用当前状态向量π_t* α_t归一化后的分布然后预测下一时刻的状态分布为π_t* * A再用加权平均的发射概率排序。实测下来对短序列效果差别不大对长序列后者略好但计算量增加。我生产环境用的就是后一种这里写的是简化版本方便读者理解。4. 在真实行为数据上训练序列构造与评测光有代码还不够数据怎么组织直接影响效果。我构造了一个模拟数据集来演示完整的实验流程这样你可以直接复现再替换成自己的真实数据。假设我们有6个物品ID对应3个兴趣主题。用户的兴趣在主题1、2、3之间切换每个主题下产生的物品符合某个发射分布。我模拟了10000条用户序列每条长度从30到80不等。4.1 生成模拟数据rng np.random.default_rng(42) n_users 10000 sequences [] for u in range(n_users): seq_len rng.integers(30, 81) state 0 obs_seq [] for t in range(seq_len): # 以一定概率切换状态 if rng.random() 0.03: state rng.integers(0, 3) # 根据当前状态生成物品 emit_probs [0.4, 0.35, 0.25, 0.0, 0.0, 0.0] state_probs [[0.4, 0.35, 0.25, 0.0, 0.0, 0.0], [0.0, 0.2, 0.3, 0.3, 0.2, 0.0], [0.0, 0.0, 0.1, 0.2, 0.3, 0.4]] item rng.choice(6, pstate_probs[state]) obs_seq.append(item) sequences.append(obs_seq)前80%做训练集后20%做测试集。在测试集上我用每一条序列的前半段作为已知行为后半段作为要预测的目标。评测指标选两个一个是“预测命中率”也就是推荐列表里是否包含用户在下一时刻实际点击的物品另一个是“Top-10覆盖率”纯粹看推荐列表的多样性。4.2 评测结果与对比用3个隐状态训练HMM迭代到对数似然变化小于阈值大概40轮收敛。结果如下模型命中率%Top-1命中率%Top-10随机推荐基线16.726.3基于流行度的推荐12.132.5HMM3状态未加平滑28.451.2HMM3状态加平滑31.755.6流行度推荐在Top-1命中率上反而不如随机原因是模拟数据里大多数物品并非流行物品用户兴趣切换后优先点击的是主题相关物品。HMM明显胜出加了平滑后又有提升原因是模拟数据里有些物品在某个状态下从没出现过发射概率为0平滑能避免推荐结果枯死。4.3 真实数据上要做的预处理真实数据没有模拟数据这么干净我总结了一套可落地的处理流程把连续行为按会话切分比如用户30分钟没有操作就认为是一个会话结束。HMM的转移概率是时刻级的跨会话的序列如果硬拼在一起会把一个状态切换误认为行为转换。物品ID需要映射成稠密整数索引否则B矩阵过大。一般保留出现次数前N的物品比如5000~10000其余丢弃。长尾问题后面细说。序列长度需要截断和补齐。要么只保留最近50个行为要么对超长序列降采样。我试过保留100个效果更好但训练时间翻倍。如果用户行为跨领域比如既看文章又看视频建议把两类观测分开建模或者做观测层面的拼接不要混在同一组发射矩阵里。5. 参数调优与冷启动状态数选择、平滑与发射建模HMM推荐效果好不好60%取决于参数怎么调不是代码本身。5.1 状态数N怎么定N太小模型表现力不够用户兴趣的复杂切换学不出来N太大又容易过拟合且状态含义变得无法解释。我实验了几组状态数训练集对数似然测试集Top-10命中率状态可解释性2-2450047.1太粗分不清主题边界3-2380055.6基本对应3个主题5-2340056.1状态开始重叠8-2310054.2状态不可解释过拟合经验法则是先用业务直觉给一个初始值比如内容社区可能有“技术、生活、娱乐、购物”四大类就试4然后做5折交叉验证看测试集对数似然是否有提升如果N1没有显著收益就别加了。状态迁移矩阵的热力图也是一个有效诊断工具如果某几个状态几乎是复制品说明N设大了。5.2 发射概率的平滑发射矩阵B是稀疏的大量物品在某个状态下从来没有出现过概率为0。训练后这些零概率会导致推荐结果直接跳过该物品看似合理其实对长尾非常不友好。我在train里加入了Dirichlet先验平滑本质上是给B加了计数def train(self, sequences, alpha0.01): # num_B alpha # 简单平滑 # 更新B时分子分母都加alphaDirichlet平滑的alpha影响很大。太小等于没加太大会让发射分布退化为均匀分布推荐变成随机。我测过alpha在0.001到0.1之间效果都不错0.01是个稳妥值。另一种更高级的做法是分层贝叶斯对每个状态引入不同的平滑强度但工程复杂度高很多我目前还没上。5.3 冷启动用户怎么办HMM处理新用户是天然劣势只有一两个观测时维特比解码基本就是瞎猜。我的做法是维护一个全局状态先验π对于行为少于5条的冷启动用户跳过维特比解码直接用π加权发射矩阵推荐也就是计算所有状态下的加权物品概率。这样至少比完全随机好等用户行为累积到阈值再切换到正常流程。有一个我要提醒的坑π必须是训练数据中所有用户初始状态的平均分布而不是状态的真实先验分布。如果你拿所有时刻的状态频率当先验冷启动推荐会被高频状态主导反而压制了用户可能的新兴趣。我在一个项目里就踩过这个坑线上冷启动推荐变得非常单一后来改成初始状态分布后明显好转。6. 工程化过程中的避坑实录数值下溢、过拟合与线上部署手写代码跑Demo很容易但放到线上容易出各种幺蛾子。这里整理几个我真实踩过的坑。6.1 后向算法中的缩放因子对齐问题我在1.0版本里的后向算法没有用scale直接跑原始beta结果训练损耗值在20轮后变成NaN。排查了很久发现是alpha和beta的尺度不一致导致xi计算溢出。修复方式就是你看到的后向递推每一步除以对应的scale[t1]保证和前向缩放一致。这个坑如果你用现成的库不会遇到但手写实现必须自己处理。调试时我建议用极短序列长度3-5跑一遍打印每个时刻的alpha、beta、gamma看看是否都在合理范围内。6.2 维特比路径的最后一步处理维特比返回的状态序列最后一步往往对应最近的一个行为。但有些序列的最后几个观测是噪声行为比如用户误点、快速滑走导致状态被带偏。我后来加了一个“状态平滑”后处理对最后的k步状态做多数表决不用最后一步作为唯一判断。k3或5都试过5略好但延迟大一点。如果你的场景对响应时间敏感k3比较稳。6.3 参数初始化的影响Baum-Welch只能保证收敛到局部最优初始化直接决定落点。我用过随机初始化和均匀初始化效果差异很大。后来总结出一套可靠方案先用KMeans在观测序列上聚类用聚类中心来初始化B矩阵A矩阵保持均匀加小扰动π用均匀分布。这样做收敛速度快最终效果也更稳定。初始化时B矩阵里不要出现全零行否则某个状态永远不会被激活。6.4 线上推断的实时性训练完成后线上推荐如果每条请求都重新跑维特比性能会很紧张尤其序列长度几百时。我的优化思路是把每个用户的当前状态维护在缓存里每次新行为来的时候做一步增量更新而不是重新解码全部历史。公式上其实很简单新的alpha向量旧的alpha乘以A再乘以B(:, obs)乘完归一化新状态就是argmax。这样时间复杂度从O(T*N^2)降到O(N^2)实测压测单机秒处理几千个增量请求。增量更新有个边界如果状态分布非常接近均匀说明当前状态信息量低此时可以强制重新跑一次完整维特比来校准。我设了一个熵阈值熵超过某值触发全量解码效果很好。6.5 训练时的内存峰值Baum-Welch里如果每条序列都保存alpha、beta、xi长序列加上几十万用户内存直接爆。解决方案是分batch训练每个batch只保留当前batch的统计量或者对超长序列做截断。我实际用的是后者序列长度超过100的直接截断成最近100个行为因为更早的历史对当前兴趣判定的贡献非常小截断后训练内存从8G降到1G。7. 进阶思路从HMM到贝叶斯方法与序列推荐的方向HMM只是起点。我做完这个项目之后对序列推荐模型的认识清晰了不少这里分享一些进阶扩展的方向你可以按需选择。将HMM的发射矩阵换成多项分布、泊松分布或高斯分布匹配不同类型的行为比如点击次数、阅读时长、评分值。HMM框架本身不限制发射分布的类型只要你能写出观测似然。引入半马尔可夫模型HSMM把“状态驻留时间”显式建模。用户在某个兴趣状态往往停留多个时刻HSMM在转移时先采样逗留时长对长稳定期的建模会比标准HMM更准代价是训练复杂度上升。结合Item2Vec或者BERT特征构建物品表示然后用高斯发射建模连续表示可以在保留HMM状态切换优势的同时利用深度学习语义。这是我目前觉得最有容量的方案。从HMM频率学派转到贝叶斯HMM给转移矩阵引入先验用变分推断或MCMC求后验。数据量小的时候贝叶斯方法的稳健性更好不会因为某个状态下的样本太少而崩。不过我要泼一盆冷水HMM做推荐的上限是有限的它对长期依赖的建模能力弱于Transformer系模型对用户-物品交互中复杂的非线性关系把握也不足。它最大的价值在于提供一个解释性强、训练成本低、适合小样本场景的基线方案。很多时候一个调参充分的HMM基线比瞎调到复杂深度学习模型效果还好我们要的正是这种能快速上线的概率模型的确定性。这个项目我做完最大的收获是序列建模的核心不是网络结构多花哨而是你如何定义“观测”和“状态”之间的边界。一个业务上合理的状态定义远比换更复杂的模型来得有效。如果你正准备做类似的项目建议先用手头的真实序列数据跑一遍HMM哪怕最终不用它上线也能为后续更复杂的模型提供一个靠谱的对比基线。
阅读完成 · 觉得有帮助？
