马尔可夫过程的估计问题是计算机科学、系统工程和数据科学等领域中的一个核心问题^[1]。计算个性化网页排序的状态转移矩阵问题^[2]、解决电子商务中的排序问题^[3]和分析城市出租车或公交车的运行轨迹问题^[4-5]等都可归结为马尔可夫过程的估计问题。源于上述问题的马尔可夫过程往往拥有很大的状态空间，但是它们的状态转移矩阵却被证明了是低秩或者近似低秩的矩阵^[1]。因此，学者们对低秩马尔可夫过程的状态转移矩阵的估计及其应用问题开展了研究^{[1, 6-10]}。

据我们所知，现有的估计方法都不能保证得到低秩的转移矩阵估计。譬如，ZHANG和WANG^[1]利用频率矩阵的经验估计的截断奇异值分解结合非负投影，提出了低秩马尔可夫过程的谱估计方法，并建立了估计的统计误差界，证明了估计误差与极小极大误差的下界相差一个马尔可夫链轨迹长度的对数因子。但是，由于谱估计方法利用了非负投影，导致该文献最后得到的估计矩阵不是低秩的。ZHU等^[8]提出了状态转移矩阵的核范数正则罚极大似然估计模型和秩约束极大似然估计模型，并建立了2种模型的统计误差界，证明了估计误差与极小极大误差的下界相差一个马尔可夫链状态空间维数的对数因子。然而，核范数正则优化问题的最优解不一定满足低秩条件，秩约束优化问题的最优解虽能满足秩约束条件，但其求解一般都是NP难。尽管文献[8]设计了一类DC (凸函数的差) 规划算法来近似求解秩约束极大似然估计模型，但不能保证算法的输出是一个低秩矩阵。特别地，该DC规划算法的每一步都需要进行奇异值分解，计算量非常大，因此不适用于大规模的马尔可夫过程估计问题。

另一方面，误差界研究一直以来都是最优化领域中的重点和难题^[11-15]。PANG^[11]证明了：凸多面体集合具有全局Lipschitz型误差界; 在一定的约束规范下，一般凸不等式系统具有Lipschitz型误差界; 对一般非凸不等式系统，全局误差界(即使是Hölder型的)都很难成立。目前已有的研究主要是对次解析不等式系统和多项式不等式系统建立了局部Hölder型误差界^[14]，对矩阵秩约束系统的误差界研究还很少。对于一个有界且其多值函数在原点满足calmness条件的矩阵秩约束系统，BI和PAN^[15]得到了该系统的局部和全局Lipschitz型误差界。但是，验证多值函数的calmness条件与建立误差界的难度基本一样。因此，我们需要寻求新的工具来研究矩阵秩约束系统的误差界。

受上述启发，本文试图寻求一个能够快速获得低秩转移矩阵的方法，以估计大规模的低秩马尔可夫过程。首先，建立秩约束状态转移矩阵集合的局部Lipschitz型误差界，并寻求该集合高质量的近似投影方法; 然后，提出一种低秩马尔可夫过程状态转移矩阵的低秩谱估计算法(LRSEA)，并进行数值实验。

1.   局部Lipschitz型误差界

本节为秩约束状态转移矩阵集合建立局部Lipschitz型误差界。首先，给出本文的记号。

(1) 记$\boldsymbol{I} \in \mathbb{R}^{d \times d}、\boldsymbol{E} \in \mathbb{R}^{d \times d}、\boldsymbol{e} \in \mathbb{R}^d$分别为单位矩阵、元素全为1的矩阵、元素全为1的向量。

(2) 对于$\boldsymbol{P} \in \mathbb{R}^{d \times d}$，定义$\|\boldsymbol{P}\|_{\infty}: =\max \limits_{i, j}\left|P_{i j}\right|$为无穷范数，$\|\boldsymbol{P}\|_F: =\sqrt{\sum\limits_{i=1}^d \sum\limits_{j=1}^d P_{i j}^2}$为Frobenius范数。

(3) 秩约束状态转移矩阵集合定义为：$\varPi: =\{\boldsymbol{P} \in \left.\mathbb{R}^{d \times d}: \operatorname{rank}(\boldsymbol{P}) \leqslant r, \boldsymbol{P e}=\boldsymbol{e}, P_{i j} \geqslant 0, 1 \leqslant i, j \leqslant d\right\}$，其中$r \in [1, d-1]$是一个给定整数。对于$\boldsymbol{Z} \in \mathbb{R}^{d \times d}$，定义Z到集合Π的距离为:

(4) 给定$\boldsymbol{x} \in \mathbb{R}^d$，定义l₁范数为$\|\boldsymbol{x}\|_1:=\sum\limits_{i=1}^d\left|x_i\right|$，无穷范数为$\|\boldsymbol{x}\|_{\infty}: =\max \limits_i\left|x_i\right|$，记$\operatorname{diag}(\boldsymbol{x}) \in \mathbb{R}^{d \times d}$是第i个对角元为x_i(i=1, 2, …, d)的对角矩阵。

(5) 定义秩r约束矩阵集合为$\mathcal{R}: =\left\{\boldsymbol{Z} \in \mathbb{R}^{d \times d}\right.$ : $\operatorname{rank}(\boldsymbol{Z}) \leqslant r\}$。对任意$\boldsymbol{P} \in \mathbb{R}^{d \times d}$，定义

其中，$\sigma_i(\boldsymbol{P})(i=1, 2, \cdots, r)$是P的第i个最大奇异值，$u_i(\boldsymbol{P}) \in \mathbb{R}^d、v_i(\boldsymbol{P}) \in \mathbb{R}^d$分别是σ_i(P)对应的左、右奇异向量。众所周知，在Frobenius范数距离意义下，$\varUpsilon_{\mathcal{R}}(\boldsymbol{P})$是P在$\mathcal{R}$上的一个投影矩阵。

(6) 定义集合$\varOmega: =\left\{\boldsymbol{Z} \in \mathbb{R}^{d \times d}: Z \boldsymbol{e}=\boldsymbol{e}, \beta / d \geqslant Z_{i j} \geqslant\right.$ $\alpha / d, 1 \leqslant i, j \leqslant d\}$，其中$\alpha \in(0, 1)$和$\beta \in(1, d)$是给定常数。

下面给出建立矩阵$\boldsymbol{P} \in \varOmega$ 与投影矩阵$\varUpsilon_{\mathcal{R}}(\boldsymbol{P})$之间关系的引理。

引理1   令$\gamma \in(0, 1)$是给定常数。任取$\boldsymbol{P} \in \varOmega$，若$\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_i \geqslant \gamma(i=1, 2, \cdots, d)$，则有：

其中$\boldsymbol{D}: =\operatorname{diag}\left(\frac{1}{\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_1}, \cdots, \frac{1}{\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_d}\right)$。

证明   由于($\varUpsilon_{\mathcal{R}}(\boldsymbol{P})$e)_j≥γ>0 (j=1, 2, …, d)且Pe=e，利用$\langle\boldsymbol{x}, \boldsymbol{y}\rangle \leqslant\|\boldsymbol{x}\|_1\|\boldsymbol{y}\|_{\infty}\left(\boldsymbol{x}, \boldsymbol{y} \in \mathbb{R}^d\right)$，有

因此，不等式(1)成立。利用不等式(1)、I－D是对角矩阵以及$d\|\boldsymbol{P}\|_{\infty} \geqslant 1>\gamma$，可得

于是，不等式(2)成立。当$0 \leqslant \min \limits_{1 \leqslant j \leqslant d}\left(\boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_j$时，不等式(3)显然成立。当$\min \limits_{1 \leqslant j \leqslant d}\left(\boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})\right)_j < 0$时, 由$\|\boldsymbol{D}\|_{\infty} \leqslant 1 / \gamma, \min \limits_{1 \leqslant j \leqslant d}\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_j < 0$且$\boldsymbol{P}_{i j} \geqslant 0, 1 \leqslant i, j \leqslant d$，可得

不等式(3)成立。证毕。

令$\gamma>0, c \in(0, 1)$为给定常数。对任意矩阵$\boldsymbol{P} \in \varOmega$, 定义如下矩阵

其中，$\boldsymbol{D}: =\operatorname{diag}\left(\frac{1}{\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_1}, \cdots, \frac{1}{\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_d}\right), \tilde{\boldsymbol{v}}=$ $\left(\frac{1}{(\widetilde{\boldsymbol{P}} \boldsymbol{e})_1}, \frac{1}{(\widetilde{\boldsymbol{P}} \boldsymbol{e})_2}, \cdots, \frac{1}{(\widetilde{\boldsymbol{P}} \boldsymbol{e})_d}\right)^{\mathrm{T}}, \widetilde{\boldsymbol{P}}=\boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})-\boldsymbol{t} \boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})$, $\boldsymbol{t} \in \mathbb{R}^d$，定义为$t_i=\min \left(0, \min \limits_{1 \leqslant j \leqslant d} \frac{\left(\boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_{i j}}{\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_j}\right)(i=1, 2$, $\cdots, d)$。

下面给出建立集合Π的局部Lipschitz型误差界的命题。

命题1   令$\gamma \in(0, 1 / 2], c \in(0, (\sqrt{5}-1) / 2), \beta \geqslant 1 \geqslant \alpha>0$是给定常数。任取$\boldsymbol{P} \in \varOmega$，有$\widetilde{\boldsymbol{P}}_{\varPi} \in \varPi$，且

证明   设$\boldsymbol{P} \in \varOmega$, 下面分5种情况证明。

(1) 假设$\left\|\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty} \leqslant 10\|\boldsymbol{P}\|_{\infty}, \left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_j \geqslant \gamma$且$\beta / c \geqslant\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_j \geqslant c \alpha(1 \leqslant j \leqslant d)$。由$\widetilde{\boldsymbol{P}}$和$t_i$的定义，可得$\widetilde{P}_{i l}=\left(\boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})\right)_{i l}-t_i\left(\boldsymbol{e} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_l \geqslant 0 (1 \leqslant i, l \leqslant d)$。注意$-t_i\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_l \geqslant 0$，可得$\widetilde{P}_{i l} \geqslant\left(\boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})\right)_{i l}(1 \leqslant i, l \leqslant d)$，因此$(\widetilde{\boldsymbol{P}} \boldsymbol{e})_i \geqslant \left(\boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_i=1(i=1, 2, \cdots, d)$，进而得到$0 < \tilde{v}_i \leqslant 1(i=1, 2, \cdots, d)$。定义$\widetilde{\boldsymbol{P}}_{\varPi}=\operatorname{Diag}(\tilde{\boldsymbol{v}}) \widetilde{\boldsymbol{P}}$，由$\operatorname{rank}(\widetilde{\boldsymbol{P}}) \leqslant r$，可得$\operatorname{rank}\left(\widetilde{\boldsymbol{P}}_{\varPi}\right) \leqslant r, \widetilde{\boldsymbol{P}}_{\varPi} \geqslant 0, \widetilde{\boldsymbol{P}}_{\varPi} \boldsymbol{e}=\boldsymbol{e}$，即$\widetilde{\boldsymbol{P}}_{\varPi} \in \varPi$。根据假设$\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} Y_{\mathcal{R}}(\boldsymbol{P})\right)_j \geqslant c \alpha(1 \leqslant j \leqslant d)$，可得$\frac{\left|\min \left(0, \min _{1 \leqslant i \leqslant d}\left(\boldsymbol{D} Y_{\mathcal{R}}(\boldsymbol{P})\right)_i\right)\right|}{c \alpha} \geqslant\|\boldsymbol{t}\|_{\infty}$。利用假设$\beta / c \geqslant\left\|\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty}$，可得

由不等式(6)、$\|\boldsymbol{D}\|_{\infty} \leqslant 1 / \gamma, \left\|\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty} \leqslant 10\|\boldsymbol{P}\|$及$\boldsymbol{P}_I=\operatorname{Diag}(\widetilde{\boldsymbol{v}}) \widetilde{\boldsymbol{P}}$，可得

因此，利用$d\|\boldsymbol{P}\|_{\infty} \leqslant \beta$及不等式(2)、(3)、(7)，可得

即不等式(5)成立。

假设$\left\|\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty}>10\|\boldsymbol{P}\|_{\infty}$或$\min \limits_{1 \leqslant j \leqslant d}\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right) < \gamma$或$\min \limits_{1 \leqslant j \leqslant d}\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right) < c \alpha$或$\max \limits_{1 \leqslant j \leqslant d}\left(e^1 \boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})\right)>\beta / c$, 由式(4)可知$\widetilde{\boldsymbol{P}}_{\varPi}=\boldsymbol{E} / d$。显然1$\widetilde{\boldsymbol{P}}_{\varPi} \in \varPi$且

(2) 假设$\left\|\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty}>10\|\boldsymbol{P}\|_{\infty}$。显然

不等式(5)成立。

(3) 假设$\left\|\varUpsilon_{\mathcal{R}}(P)\right\|_{\infty} \leqslant 10\|\boldsymbol{P}\|_{\infty}$且$\min \limits_{1 \leqslant j \leqslant d}\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_j < \gamma。$由于$\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}$ 中存在小于$\gamma$ 的元素，不妨设$\left(\boldsymbol{Y}_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_i < \gamma_{\circ}$由$(\boldsymbol{P e})_i=1$和d‖P－$\varUpsilon_{\mathcal{R}}(\boldsymbol{P})$$d\|\boldsymbol{P}\|_{\infty} \leqslant \beta$，可得

则

由此可得不等式(5)成立。

(4) 假设$\left\|\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty} \leqslant 10\|\boldsymbol{P}\|_{\infty}, \min \limits_{1 \leqslant j \leqslant d}\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_j \geqslant \gamma$ 且$\min \limits_{1 \leqslant j \leqslant d}\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right) < c \alpha$。由于$\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})$中存在小于$c \alpha$的元素，不防设$\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})\right)_i < c \alpha < \alpha$。由$\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{P}\right)_i \geqslant \alpha$，可得

由不等式(8)、$d\|\boldsymbol{P}\|_{\infty} \leqslant \beta$ 及引理1，可得

由$c \in(0, (\sqrt{5}-1) / 2)$，可得不等式(5)成立。

(5) 假设$\left\|\varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right\|_{\infty} \leqslant 10\|\boldsymbol{P}\|_{\infty}, \min \limits_{1 \leqslant j \leqslant d}\left(\varUpsilon_{\mathcal{R}}(\boldsymbol{P}) \boldsymbol{e}\right)_j \geqslant$ $\gamma$。由于$\max \limits_{1 \leqslant j \leqslant d}\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \varUpsilon_{\mathcal{R}}(\boldsymbol{P})\right)_j>\beta / c_{\text {。}}$ 由于$\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})$中存在大于β/c的元素，不妨设$\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{D} \boldsymbol{\varUpsilon}_{\mathcal{R}}(\boldsymbol{P})\right)_i>\beta / c$。由$\left(\boldsymbol{e}^{\mathrm{T}} \boldsymbol{P}\right)_i \leqslant$ $\beta$，可得

由不等式(9)、$d\|\boldsymbol{P}\|_{\infty} \leqslant \beta$及引理1, 可得

由$c \in(0, (\sqrt{5}-1) / 2)$，可得不等式(5)成立。证毕。

由于$\widetilde{\boldsymbol{P}}_{\varPi} \in \varPi$，命题1给出了任意$\boldsymbol{P} \in \varOmega$到集合Π的距离估计，即

同时，在集合Π中找到了低秩状态转移矩阵$\widetilde{\boldsymbol{P}}_{\varPi I}$，使得P到$\widetilde{\boldsymbol{P}}_{\varPi}$的距离满足误差界不等式。所以，$\widetilde{\boldsymbol{P}}_I$可以看成是$\boldsymbol{P} \in \varOmega$ 在集合Π上的近似投影矩阵。

2.   LRSEA算法

考虑d (≥2)个状态{S₁, S₂, …, S_d}的离散时间马尔可夫过程，假设其状态转移矩阵为$\overline{\boldsymbol{P}} \in \mathbb{R}^{d \times d}$，频率矩阵为$\overline{\boldsymbol{Q}} \in \mathbb{R}^{d \times d}$，且满足$\operatorname{rank}(\overline{\boldsymbol{P}})=\operatorname{rank}(\overline{\boldsymbol{Q}}) \leqslant r \ll d$。源于现实生活当中的一些马尔可夫过程已经被证明是低秩或者近似低秩的，比如交通网络、网页及电子商务排序^[1-7]。本文旨在通过轨迹长为n+1的马尔可夫链{X₀, X₁, …, X_n}来估计其状态转移矩阵$\overline{\boldsymbol{P}}$。

首先给出本文的假设。

假设1   存在常数$\beta \in[1, d), \alpha \in(0, 1]$，使得:

(1) $\forall 1 \leqslant i, j \leqslant d, \alpha / d \leqslant \bar{P}_{i j} \leqslant \beta / d$;

(2) $\left\{X_0, X_1, \cdots, X_n\right\}$是遍历马尔可夫链，平稳分布为$\boldsymbol{\pi} \in \mathbb{R}^d$，满足$\pi_i \geqslant \alpha / d(i=1, 2, \cdots, d)$, 其中$\pi_i$是向量$\boldsymbol{\pi}$的第i个元素。

注意到文献[1]的谱估计矩阵虽然不是低秩的，但其是一个状态转移矩阵。因此，以谱估计矩阵为基础，利用式(4)计算秩-r状态转移矩阵作为$\overline{\boldsymbol{P}}$的估计，具体见算法1。

算法1   LRSEA算法

输入：谱估计矩阵$\hat{\boldsymbol{S}} \in\left\{\boldsymbol{Z} \in \mathbb{R}^{d \times d}: Z \boldsymbol{e}=\boldsymbol{e}, Z_{i j} \geqslant 0, 1 \leqslant i, j \leqslant\right.$ $d\}, \alpha \in(0, 1), \beta \in(1, d)$和r。

第1步：计算修正的谱估计

for i=1 to d

if $\min \limits_{1 \leqslant j \leqslant d} \hat{S}_{i j} <\frac{\alpha}{d}$

$n_{+}=\sum\limits_{j=1}^d \max \left(\frac{\alpha}{d}-\hat{S}_{i j}, 0\right); n_{-}=1-\alpha+n_{+}$;

for j=1 to d

$\widetilde{S}_{i j}=\max \left(\hat{S}_{i j}-\frac{n_{+}}{n_{-}}\left(\hat{S}_{i j}-\frac{\alpha}{d}\right), \frac{\alpha}{d}\right)$;

end for

else

$\tilde{\boldsymbol{S}}(i, : )=\hat{\boldsymbol{S}}(i, : )$;

end if

end for

for i=1 to d

if $\max \limits_{1 \leqslant j \leqslant d} \widetilde{S}_{i j}>\frac{\beta}{d}$

$m_{+}=\sum\limits_{j=1}^d \max \left(\widetilde{S}_{i j}-\frac{\beta}{d}, 0\right); m_{-}=m_{+}+\beta-1$;

for j=1 to d

$S_{i j}=\min \left(\widetilde{S}_{i j}+\frac{m_{+}}{m_{-}}\left(\frac{\beta}{d}-\widetilde{S}_{i j}\right), \frac{\beta}{d}\right)$

end for

else

$\boldsymbol{S}(i, : )=\widetilde{\boldsymbol{S}}(i, : );$

end if

end for

第2步：令$\boldsymbol{P}=\boldsymbol{S}, \gamma \in(0, 1 / 2], c \in(0, (\sqrt{5}-1) / 2)$，由式(4)得到$\overline{\boldsymbol{P}}$的低秩谱估计$\widetilde{\boldsymbol{S}}_{\varPi}$。

由于LRSEA算法产生的估计矩阵$\tilde{\boldsymbol{S}}_{\varPi} \in \varPi$是一个满足低秩约束条件的状态转移矩阵，因此称其为低秩谱估计算法。可以预见，由于LRSEA算法可以满足低秩性要求，其有效性将高于现有的谱估计方法。为了建立LRSEA算法的理论保证，定义马尔可夫链的$\frac{1}{4}$-混合时间^[16]为

定理1   定义$\bar{r}: =\|\overline{\boldsymbol{Q}}\|_F^2 / \sigma_r^2(\overline{\boldsymbol{Q}})$，其中$\sigma_r(\overline{\boldsymbol{Q}})$是$\overline{\boldsymbol{Q}}$的第r个最大奇异值，且记$\kappa=d^2\|\overline{\boldsymbol{Q}}\|_{\infty}$。若假设1成立，$\tilde{S}, \boldsymbol{S}, \tilde{\boldsymbol{S}}_I \in \mathbb{R}^{d \times d}$由LRSEA算法产生，则有：

(a) $S \in \varOmega$且

(b) 令$C_1=\frac{11 \beta}{\gamma}\left(\frac{\beta}{c^2 \alpha \gamma}+1\right) \sqrt{d}+1$。存在常数C，使得

其中E表示数学期望。

证明   (a)由算法1，易证$n_{-}=\sum\limits_{j=1}^d \max \left(\hat{S}_{i j}-\frac{\alpha}{d}, 0\right)$, $m_{-}=\sum\limits_{j=1}^d \max \left(\frac{\beta}{d}-\tilde{S}_{i j}, 0\right)$。由于$\alpha < 1$，所以$n_{-}>n_{+} \geqslant 0$。注意到$\hat{\boldsymbol{S}} \in\left\{\boldsymbol{Z} \in \mathbb{R}^{d \times d}: Z \boldsymbol{e}=\boldsymbol{e}, Z_{i j} \geqslant 0, 1 \leqslant i, j \leqslant d\right\}$。若$\min \limits_{1 \leqslant j \leqslant d} S_{i j} \geqslant \alpha / d$，则$\tilde{\boldsymbol{S}}(i, j)=\hat{\boldsymbol{S}}(i, j) \geqslant \alpha / d(j=1, 2$, $\cdots, d)$ 且$\tilde{\boldsymbol{S}}(i, : ) \boldsymbol{e}=1$。若$\min \limits_{1 \leqslant j \leqslant d} \hat{S}_{i j} < \alpha / d$，利用$n_{-}>$ $n_{+} \geqslant 0$，$\hat{S}_{i j} \leqslant \alpha / d$ 时, 有$\tilde{S}_{i j}=\alpha / d$; 当$\hat{S}_{i j}>\alpha / d$ 时, 有

且

因此，可得$\alpha / d \leqslant \tilde{S}_{i j}$且$\tilde{\boldsymbol{S}}(i, : ) \boldsymbol{e}=1(i, j=1, \cdots, d)$。同理, 可以证明$\alpha / d \leqslant S_{i j} \leqslant \beta / d$ 且$\boldsymbol{S}(i, : ) \boldsymbol{e}=1(i, j=$ $1, \cdots, d)$。因此, $S \in \varOmega_{\text {。}}$

若$\min \limits_{1 \leqslant j \leqslant d} \hat{S}_{i j} \geqslant \alpha / d$，由于$\tilde{\boldsymbol{S}}(i, : )=\hat{\boldsymbol{S}}(i, : )$，有

现假设$\min \limits_{1 \leqslant j \leqslant d} \hat{S}_{i j} < \alpha / d$。利用$\alpha / d \leqslant \bar{P}_{i j} \leqslant \beta / d \quad(1 \leqslant i$, $j \leqslant d)$, 可得

由上述讨论，可得

同理可以证明

所以，

(b) 若假设1成立，则由文献[1]的定理1可知：存在一个常数C，使得

由(a)部分的结论，可知存在常数C，使得

由于$\boldsymbol{S} \in \varOmega$，由命题1可知

其中$C_2=\frac{11 \beta}{\gamma}\left(\frac{\beta}{c^2 \alpha \gamma}+1\right)$。利用$\varUpsilon_{\mathcal{R}}(S)$的定义可得$\left\|\boldsymbol{S}-\varUpsilon_{\mathcal{R}}(\boldsymbol{S})\right\|_F \leqslant\|\boldsymbol{S}-\overline{\boldsymbol{P}}\|_F$。于是

由不等式(10)、(11)可得定理结论。证毕。

定理1建立了LRSEA算法产生的估计矩阵与真实状态转移矩阵$\overline{\boldsymbol{P}}$之间的距离估计。由定理1(b)可知，轨迹长度n越大，估计误差的期望越小。

3.   数值实验

在LRSEA算法中，取$\gamma=0.01, c=0.02, \alpha=\frac{1}{2} d \times \min \limits_{1 \leqslant i, j \leqslant d} \bar{P}_{i j}, \beta=\min \left(d, 2 d \max \limits_{1 \leqslant i, j \leqslant d} \bar{P}_{i j}\right)$。首先，通过人工合成数据实验来比较LRSEA算法、谱估计方法^[1]、经验估计方法^[17]的效果。然后，利用LRSEA算法与k-均值聚类方法结合来分析纽约曼哈顿岛出租车运行轨迹的公开数据集，揭露该城市的潜在交通模式。所有数值实验都是在配置为英特尔奔腾G4600处理器和CPU主频3.6 GHz的笔记本电脑上，通过运行MATLAB(2019年版本)完成的。

3.1   合成数据实验

首先，考虑具有平衡分布的低秩马尔可夫过程估计问题。假设$\boldsymbol{U}_0 \in \mathbb{R}^{d \times r}, V_0 \in \mathbb{R}^{d \times r}$，它们的元素由标准正态分布随机生成。定义矩阵$\widetilde{\boldsymbol{U}} \in \mathbb{R}^{d \times r}, \widetilde{\boldsymbol{V}} \in \mathbb{R}^{d \times r}$ : $\widetilde{\boldsymbol{U}}_{[i, i]}=\left(\boldsymbol{U}_0 \odot \boldsymbol{U}_0\right)_{[i, :]} /\left\|\left(\boldsymbol{U}_0\right)_{[i, i]}\right\|_2^2(i=1, \cdots, d)$, $\widetilde{\boldsymbol{V}}_{[:, j]}=\left(\boldsymbol{V}_0 \odot \boldsymbol{V}_0\right)_{[:, j]} /\left\|\left(\boldsymbol{V}_0\right)_{[:, j}\right\|_2^2(j=1, \cdots, r)$，这里$\widetilde{\boldsymbol{U}}_{[i, , ]}$表示$\widetilde{\boldsymbol{U}}$的第i行，$\widetilde{\boldsymbol{V}}_{[:, j]}$ 表示$\widetilde{\boldsymbol{V}}$ 的第$j$ 列, $\odot$ 表示矩阵Hadamard积。定义真实状态转移矩阵为$\overline{\boldsymbol{P}}=\widetilde{\boldsymbol{U}} \widetilde{\boldsymbol{V}}^{\mathrm{T}}$。本文利用状态转移矩阵$\overline{\boldsymbol{P}}$ 生成状态数为$d$、长度为$n=\operatorname{round}\left(q d r(\log d)^2\right)$ 的马尔可夫链$X_0$, $\cdots, X_n$, 这里$q$ 是常数。

下面比较文献[17]的经验估计方法、文献[1]的谱估计方法和LRSEA算法的估计误差。记P为相应的估计矩阵，本文利用以下2个数值来衡量其估计效果：

其中，$\boldsymbol{U} \in \mathbb{R}^{d \rtimes r}、\boldsymbol{V} \in \mathbb{R}^{d \times r}$分别是P的前r个最大奇异值对应的左奇异向量、右奇异向量，$\overline{\boldsymbol{U}} \in \mathbb{R}^{d \times r}、\overline{\boldsymbol{V}} \in \mathbb{R}^{d x_r}$ 分别是$\overline{\boldsymbol{P}}$的前r个最大奇异值对应的左奇异向量、右奇异向量。

取d=1 000、r=10、k$\in$^{[1, 10]}时，3种方法的估计效果(图 1)表明：LRSEA算法与谱估计方法的估计误差相差不大，小于经验估计方法的。

图  1  平衡分布下3种估计方法的比较

Figure  1.  The comparison of three estimators in balanced distribution

下载: 全尺寸图片 幻灯片

接下来，考虑具有非平衡分布的低秩马尔可夫过程的估计问题，即一部分状态的概率远小于其他状态的概率。此处构造矩阵$\boldsymbol{P}=\widetilde{\boldsymbol{U}} \widetilde{\boldsymbol{V}}^{\mathrm{T}} \boldsymbol{B}$，并将其标准化得到$\overline{\boldsymbol{P}}$，这里$\widetilde{\boldsymbol{U}}、\widetilde{\boldsymbol{V}}$与前面一致，B为对角矩阵，其对角元服从独立同分布的贝塔分布(Beta(α₁, α₂))。取α₁=α₂=0.1，d=1 000、r=10、轨迹长度n=round(qdr(log d)²)时, 3种方法的估计效果(图 2)表明：对于非平衡分布情况下的马尔可夫过程估计问题，LRSEA算法的估计误差小于谱估计方法和经验估计方法的。

图  2  非平衡分布下3种估计方法的比较

Figure  2.  The comparison of three estimators in imbalanced distribution

下载: 全尺寸图片 幻灯片

3.2   分析曼哈顿出租车运行轨迹

纽约市曼哈顿岛于2016年公开的黄色出租车运行轨迹数据集记录了r次乘客的行程(https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2016-01.csv)，本研究利用此数据集将曼哈顿岛分割成几个区域，满足同一区域中的乘客前往同一个目的地的概率相似。

类似文献[1]，本文将曼哈顿岛细分为大小一致的小正方形网格，将每个小网格近似地看成马尔可夫过程的一个状态，乘客从一个网格到另一个网格的行程视为该马尔可夫过程的一次状态转移。为了排除干扰项，本文只选取那些作为行程的起点或者终点的总次数超过2 000的网格作为有效状态。然后，利用LRSEA算法来估计该马尔可夫过程的状态转移矩阵，并利用k-均值聚类方法对估计矩阵的左奇异子空间进行聚类划分。由r=k分别为4, 5, 6, 7时的聚类结果(图 3)可知：当聚类数r增加时，LRSEA算法可以给曼哈顿岛的交通网络一个较好的分区，同一个区域的乘客前往同一个地点的概率相似。

图  3  LRSEA算法对曼哈顿岛交通网络的划分结果

注：1种颜色代表 1个区域。

Figure  3.  The result of a citywide partition of Manhattan traffic network with the low-rank spectral method

下载: 全尺寸图片 幻灯片

鉴于在不同时段，人们出行目的不同，将该出租车运行数据集化分为3个时段：早上(06:00~11:59)、中午(12:00~17:59)、晚上(18:00~23:59)，每个时段的有效状态数分别为769、1 029、1 147个。利用LRSEA算法来估计每个时段的马尔可夫过程的状态转移矩阵，并利用k-均值聚类方法对估计矩阵的左奇异子空间进行聚类划分，其中xietichuli=xietichuli=5。聚类结果所展示的在不同时段下LRSEA算法对曼哈顿岛的交通网络的分区结果(图 4)表明：同一时段下，同一个区域的乘客前往同一个目的地的概率相似。

图  4  LRSEA算法在不同时间段下对曼哈顿岛交通网络的划分结果

注：1种颜色代表 1个区域。

Figure  4.  The result of a per-time-segment citywide partition of Manhattan traffic network with the low-rank spectral method

下载: 全尺寸图片 幻灯片

4.   结论

针对马尔可夫过程的估计问题，利用秩约束状态转移矩阵集合的近似投影，本文对现有的谱方法进行低秩修正，提出了一个低秩谱估计算法(LRSEA)，以快速得到满足秩约束条件的状态转移矩阵。此外，通过建立秩约束状态转移矩阵集合的局部Lipschitz型误差界，给出该算法的统计误差界，建立了算法的理论保证。数值实验结果表明，对于具有非平衡分布的低秩马尔可夫过程的估计问题，LRSEA算法的估计误差小于谱估计方法和经验估计方法的。下一步，将把LRSEA算法应用到强化学习问题以及系统工程领域中的控制问题。