0 引言

1 正交投影非负矩阵分解

2 基于交替方向乘子的正交投影非负矩阵分解方法

3 仿真实验与应用

3.1 ADMM_OPNMF算法与乘性规则收敛性比较

为了比较ADMM算法和乘性规则在分解OPNMF时的性能，选取4种不同规模 (以行数×列数表示规模) 的随机矩阵 (分别为M1、M2、M3和M4) 进行模拟实验，不同矩阵的维数如表 1所示 (表 1中的Y表示待分解矩阵，W表示分解后的基矩阵)。对这4种维数的随机矩阵分别使用ADMM方法和乘性规则进行迭代更新，得到目标函数值的收敛曲线如图 1所示。

表 1 4种随机矩阵
Table 1 Four sizes of random matrixes

下载CSV

	M1	M2	M3	M4
Y	256×160	512×160	1 024×160	2 048×160
W	256×20	512×20	1 024×20	2 048×20

图 1 矩阵收敛曲线

Fig. 1 Matrix convergence curves ((a) convergence curve of M3; (b) convergence comparison of two methods)

收敛曲线图的纵坐标表示的是迭代过程中目标函数的值，横坐标为运行时间，实验过程中迭代次数固定为40 000次。图 1 (a) 是随机矩阵在不同参量ρ的情况下与乘性更新规则的比较 (以矩阵M3为例)，可以看出ADMM方法中的ρ越小，目标函数收敛越快，但当ρ较小时，收敛相对不稳定 (如图 1中红色的线所示)。与乘性规则的收敛曲线相比 (图 1(a)中黑色的虚线所示)，ADMM方法目标函数的收敛值更小，即达到了更高的收敛精度。图 1(b)为ρ=10时，ADMM方法与乘性规则分别在4种不同维数矩阵的收敛对比，实验结果表明，随着矩阵规模的增大，两种更新方法迭代40 000次所需的时间都在不断增加。其中，在矩阵维数较小时ADMM方法的收敛速度略慢乘性规则 (如图 1(b)中的红实线和红虚线所示)，这是因为在迭代过程中ADMM方法需要交替更新4个主变量X，H，W，W₊和3个对偶变量α_X，α_H，α_W，而乘性规则只需要更新一个变量W。但随着矩阵维数增加，乘性规则所需时间的增长越来越明显，收敛的速度越来越慢，ADMM方法收敛的速度以及迭代40 000次所需的时间要明显快于乘性规则 (如图 1(b)中的黑实线和黑虚线所示)。这是乘性规则在找每一个局部最优解时都是通过一定的步长慢慢趋近的，而ADMM方法每次迭代就能找到当前的一个局部最优解，因此找到最优解的时间较快。

3.2 ADMM_OPNMF算法在目标跟踪上的应用

人脸实验表明PNMF分解得到的基向量具有非常高的正交性，能够保证获取人脸各部分的正交表示^[10]。采用王栋等人^[20]提出的IOPNMF目标跟踪框架，利用基向量构建运动目标的观测模型，将每帧视频中的目标表示为基向量集的线性组合，并且以重构误差为基础设计了观测似然函数。本文目标跟踪的整体框架如图 2所示。

图 2 目标跟踪整体框图

Fig. 2 Target tracking framework

本文目标跟踪算法是基于贝叶斯粒子滤波框架的^[20]。粒子滤波是通过非参数化的蒙特卡罗模拟来实现贝叶斯滤波，即利用一组带有权重的随机样本近似描述系统状态的后验概率密度。目标在$ t-1 $时刻前的观测集$ {z_{1:t-1}} = \{ {z_1}, {z_2}, \cdots, {z_{t-1}}\} $，则$ t $时刻目标的最佳状态可由最大近似后验概率得出$ z_t^{^*} = {\rm{arg min}}\;p(x_{_t}^{^i}|{z_{1:t}}) $。其中$ x_{_t}^{^i} $表示时刻$ t $的第$ i $个采样粒子的系统状态，后验概率$ p(x_{_t}^{^i}|{z_{1:t}}) $可由贝叶斯理论递归求得，即

$ \begin{array}{l} p\left( {{x_t}|{z_{1:t}}} \right) \propto \\ p({z_t}|{x_t})\smallint p({x_t}|{x_{t-1}})p({x_{t-1}}|{z_{1:t-1}})d{x_{t - 1}} \end{array} $

(10)

式中，$ p({x_t}|{x_{t-1}})$是状态转移概率，用来描述动态模型，$ p({x_{t-1}}|{z_{1:t-1}}) $是上一时刻的后验概率，$ p({z_t}|{x_t}) $是根据当前观测模型计算的当前观测似然度。

粒子滤波在跟踪问题中大致可分为两个部分：运动模型和目标观测模型。运动模型的作用是根据上一时刻的目标状态预测当前时刻的目标状态。本文所使用的运动模型是仿射参数采样模型，使用6个仿射参数描述目标的运动状态变化，即

$ \left( \begin{array}{l} {u'}\\ {v'} \end{array} \right) = \left[ {\begin{array}{*{20}{c}} {\cos (\theta )}&{\sin (\theta )}\\ { - \sin (\theta )}&{\cos (\theta )} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} s&\beta \\ 0&{\alpha s} \end{array}} \right]\left( \begin{array}{l} u + x\\ v + y \end{array} \right) $

(11)

式中，$ (u, v) $和$ ({u^{'}}, {v^{'}}) $分别表示仿射变换前后的坐标位置，$ x $和$ y $分别代表水平和垂直方向的位移，θ表示旋转角度，s表示水平尺度，α表示宽高比，β表示切变系数。利用这6个仿射参数表示目标在$t $时刻的运动状态$ {\mathit{\boldsymbol{x}}_t} = ({x_t}, {y_t}, {\theta _t}, {s_t}, {\alpha _t}, {\beta _t}) $，假定前一时刻的目标状态到当前时刻状态的变化符合独立高斯分布，系统运动模型可表示为

$ p({\mathit{\boldsymbol{x}}_t}{\mathit{\boldsymbol{x}}_{t-1}}) = N({\mathit{\boldsymbol{x}}_t};{\mathit{\boldsymbol{x}}_{t-1}}, \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}) $

(12)

式中，Ψ为仿射变换参数方差 (即$\sigma _{_x}^{^2}、\sigma _{_y}^{^2}、\sigma _{_\theta }^{^2}、\sigma _{_s}^{^2}、\sigma _{_\alpha }^{^2}、\sigma _{_\beta }^{^2} $) 的对角协方差矩阵。

目标观测模型的主要作用从每一帧的数百个粒子中评估每个采样粒子点是目标的可能性，即计算采样粒子$\mathit{\boldsymbol{x}}_{_t}^i $对应的图像向量服从目标观测模型的概率$ p(\mathit{\boldsymbol{y}}_{_t}^{^i}|\mathit{\boldsymbol{x}}_{_t}^{^i}) $，即

$ p(\mathit{\boldsymbol{y}}_{_t}^{^i}\mathit{\boldsymbol{x}}_{_t}^{^i}) \propto {\rm{exp(-||}}{\mathit{\boldsymbol{y}}^i}-\mathit{\boldsymbol{W}}\overline {{\mathit{\boldsymbol{x}}^i}} ||_2^2 $

(13)

式中，$ i $代表粒子号，$ t$代表帧号，∝表示正比例于。针对跟踪过程中目标图像受到遮挡或异常噪声的情况，引入稀疏噪声项^[24]，候选目标的观测似然函数表示为

${\rm{arg}}\;\mathop {{\rm{min}}}\limits_{{x^i}, {e^i}} \frac{1}{2}||y_{_t}^{^i}-\mathit{\boldsymbol{W}}{\mathit{\boldsymbol{h}}^i}-{\mathit{\boldsymbol{e}}^i}||_2^2 + \lambda ||{\mathit{\boldsymbol{e}}^i}|{|_1} $

(14)

式中，W是目标图像的基向量集，$ \overline {{\mathit{\boldsymbol{e}}^i}} $和$ \overline {{\mathit{\boldsymbol{x}}^i}} $分别表示第$ i $个粒子的稀疏噪声项和投影系数。

在跟踪过程中，目标的外观不断发生变化，需要及时更新观测模型，以防跟踪漂移。当目标受到遮挡时，直接使用新的观测样本更新模型将遮挡信息加入新的观测模型，会使得模型退化。

借鉴王栋等人的IOPNMF更新策略^[20]，结合OPNMF在图像基于部分表示的特点以及ADMM在更新速度上的优势，提出一种ADMM_OPNMF更新方法。该方法在每次更新样本时，将当前第$ i $个样本的前$ k $帧并入样本集，以避免在发生遮挡或异常噪声时，更新后的样本只含有当前的噪声信息，使得跟踪漂移。第$ i $个样本的基图像$ {\mathit{\boldsymbol{W}}_i} $可以利用第2节的方法得到更新为

$ \begin{array}{l} {\mathit{\boldsymbol{W}}_i} \leftarrow {({\mathit{\boldsymbol{H}}_{i- 1}}\mathit{\boldsymbol{H}}_{^{i- 1}}^{\rm{T}} + \mathit{\boldsymbol{I}})^{- 1}}({\mathit{\boldsymbol{H}}_{i - 1}}{({\mathit{\boldsymbol{X}}_{i - 1}} + \frac{1}{\rho }{\mathit{\boldsymbol{\alpha }}_X})^{\rm{T}}} + \\ ({\mathit{\boldsymbol{H}}_{i - 1}} + \frac{1}{\rho }{\mathit{\boldsymbol{\alpha }}_H}\mathit{\boldsymbol{Y}}_{_{[i-k, i]}}^{^{\rm{T}}} - {\mathit{\boldsymbol{H}}_{i - 1}}\mathit{\boldsymbol{Y}}_{_{[i-k, i]}}^{\rm{T}} + {({\mathit{\boldsymbol{W}}_{ +, i -1}} -\frac{1}{\rho }{\mathit{\boldsymbol{\alpha }}_W})^{\rm{T}}}) \end{array} $

(15)

式中，$ {\mathit{\boldsymbol{Y}}_{[i-k, i]}} $表示第$ i $个样本的前$ k $帧样本集合。

本文跟踪算法的主要步骤如下：

1) 输入：输入视频帧。

2) 输出：跟踪目标位置。

3) 初始化：每一帧采样粒子数设为600，取前$ k $($ k $=20) 帧作为初始观测矩阵。

4) 从$ k $+1帧至总帧数

(1) 根据仿射采样运动模型式 (12)，获取600个粒子样本及其对应的观测图像矩阵$ \mathit{\boldsymbol{y}}_{_t}^{^i} $；

(2) 利用目标函数每个采样粒子矩阵，进行PNMF分解式 (4) 得到对应的基图像矩阵W；

(3) 通过观测模型式 (14) 求解观测目标似然函数的最大值，得到最优的候选目标，即得到当前帧所跟踪的目标；

(4) 利用ADMM_OPNMF方法式 (15) 进行观测模型更新，得到基图像矩阵${\mathit{\boldsymbol{W}}_i} $；

(5) 转到步骤3)，继续取下一帧图像。

由于OPNMF学习到的基向量是基于部分的表示，具有潜在的遮挡处理能力，在目标短时间内被遮挡时可以取得相对较好的跟踪效果。因此，本文选取Dudek，FaceOcc1，FaceOcc2，Signer1 4个视频集^[25]，这些视频集具有短时或部分遮挡、短时光照变化以及尺度变化等特点。同时，为了验证本文ADMM_OPNMF算法的性能，不仅选取乘性规则的IOPNMF算法 (Multi_IOPNMF)，还选取目前较为常用的3种目标跟踪算法IVT (incremental visual tracking)^[26]、LOT (locally orderless tracking)^[27]以及L1_APG (L1 tracker using accelerated proximal gradient approach)^[28]。图 3是不同算法在4个视频集上的跟踪过程图，可以直观的反应跟踪过程。而为了客观地对本文算法与其他4种算法性能进行评估，本文利用重叠率准则来评估跟踪算法的精度, 即$score = \frac{{area\left( {{\mathit{\boldsymbol{R}}_r} \cap {\mathit{\boldsymbol{R}}_g}} \right)}}{{area({\mathit{\boldsymbol{R}}_r} \cup {\mathit{\boldsymbol{R}}_g})}} $，其中${{\mathit{\boldsymbol{R}}_r}} $是每帧目标的跟踪框覆盖的区域，$ {{\mathit{\boldsymbol{R}}_g}} $是测试序列提供的真实位置所在区域。重叠率取值范围在0和1之间，数值越大代表结果越精确^[20]。

图 3和表 2所示为5种算法在视频集上的具体跟踪过程图和重叠率比较，首先可以看出对于短时或者部分遮挡的情况，IOPNMF算法的效果明显好于LOT、IVT以及L1_APG 3个算法，这是因为在对目标区域进行分解时得到的基图像是基于部分的表示，可以有效地避免少量遮挡对于模板的影响，获取更好的跟踪效果。而对于使用ADMM方法更新OPNMF也获得了和乘性规则相当的效果，验证了此方法的正确性和实用性。

图 3 跟踪结果图

Fig. 3 Tracking results ((a) Dudek; (b) FaceOcc1; (c) FaceOcc2; (d) Singer1)

表 2 重叠率比较
Table 2 Overlap rate comparison

下载CSV

视频集	LOT	IVT	L1_APG	Multi_IOPNMF	ADMM_OPNMF
Dudek	0.536	0.753	0.692	0.767	0.761
FaceOcc1	0.439	0.746	0.702	0.793	0.790
FaceOcc2	0.458	0.727	0.564	0.612	0.623
Singer1	0.494	0.574	0.504	0.744	0.749
平均值	0.482	0.701	0.616	0.729	0.731
注：粗体为最优结果。

表 3是ADMM方法与乘性规则在4种视频集下每帧的更新时间比较，可以看出在跟踪过程中使用ADMM方法更新IOPNMF时的速度明显快于乘性规则。在视频集中IOPNMF的模板Y的维数大致为1 024×20，尽管王栋等人为解决乘性规则的实时性问题提出增量式更新方法，但是与本文提出的ADMM方法相比跟踪速度仍旧较慢。在本文实验条件下，ADMM_OPNMF方法的帧处理速度约是Multi_IOPNMF方法的3.8倍。

表 3 跟踪时间比较
Table 3 Time tracking comparison

下载CSV

/(帧/s)
视频集	Multi_IOPNMF	ADMM_OPNMF
Dudek	1.438	0.469
FaceOcc1	1.395	0.370
FaceOcc2	1.421	0.343
Singer1	1.598	0.351
平均值	1.463	0.383
注：粗体为最优结果。

4 结论

参考文献

[1] Guo K L, Xu X M, Qiu F H, et al. A novel incremental weighted PCA algorithm for visual tracking[C]//Proceedings of 2013 IEEE International Conference on Image Processing. Melbourne, VIC, Australia:IEEE, 2013:3914-3918.[DOI:10.1109/ICIP.2013.6738806]
[2] Bartlett M S, Movellan J R, Sejnowski T J. Face recognition by independent component analysis[J]. IEEE Transactions on Neural Networks, 2002, 13(6): 1450–1464. [DOI:10.1109/TNN.2002.804287]
[3] Gao Q X, Liang Y, Pan Q, et al. The problem existed in face recognition using SVD and its solution[J]. Journal of Image and Graphics, 2006, 11(12): 1784–1791. [高全学, 梁彦, 潘泉, 等. SVD用于人脸识别存在的问题及解决方法[J]. 中国图象图形学报, 2006, 11(12): 1784–1791. ] [DOI:10.11834/jig.2006012312]
[4] Lee D D, Seung H S. Learning the parts of objects by non-negative matrix factorization[J]. Nature, 1999, 401(6755): 788–791. [DOI:10.1038/44565]
[5] Chen X R, Gu L, Li S Z, et al. Learning representative local features for face detection[C]//Proceedings of 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Kauai, HI, USA:IEEE, 2001, 1:I-1126-I-1131.[DOI:10.1109/CVPR.2001.990657]
[6] Barman P C, Lee S Y. Tree cluster of text data by NMF based neural network[C]//Proceedings of 2006 International Conference on Electrical and Computer Engineering. Dhaka, Bangladesh:IEEE, 2006:312-315.[DOI:10.1109/ICECE.2006.355634]
[7] Du H S, Zhang X D. Graph embedding regularized projective non-negative matrix factorization for face image feature extraction[J]. Journal of Image and Graphics, 2014, 19(8): 1176–1184. [杜海顺, 张旭东. 图嵌入正则化投影非负矩阵分解人脸图像特征提取[J]. 中国图象图形学报, 2014, 19(8): 1176–1184. ] [DOI:10.11834/jig.20140809]
[8] Guo X, Uhlich S, Mitsufuji Y. NMF-based blind source separation using a linear predictive coding error clustering criterion[C]//Proceedings of 2015 IEEE International Conference on Acoustics, Speech and Signal Processing. South Brisbane, QLD, Australia:IEEE, 2015:261-265.[DOI:10.1109/ICASSP.2015.7177972]
[9] Wu Y, Shen B, Ling H. Visual tracking via online nonnegative matrix factorization[J]. IEEE Transactions on Circuits and Systems for Video Technology, 2014, 24(3): 374–383. [DOI:10.1109/TCSVT.2013.2278199]
[10] Yuan Z J, Oja E. Projective nonnegative matrix factorization for image compression and feature extraction[C]//Proceedings of the 14th Scandinavian Conference on Image Analysis. Berlin Heidelberg:Springer, 2005:333-342.[DOI:10.1007/11499145_35]
[11] Yuan Z J, Yang Z R, Oja E. Projective nonnegative matrix factorization:sparseness, orthogonality and clustering[J/OL]. Neural Processing Letters, 2009.[2016-10-28]https://pdfs.semanticscholar.org/c9b6/b62ce526b6d5a12d324ec40309e0a63c4582.pdf
[12] Feng T, Li S Z, Shum H Y, et al. Local non-negative matrix factorization as a visual representation[C]//Proceedings of the 2nd International Conference on Development and Learning. Cambridge, MA, USA:IEEE, 2002:178-183.[DOI:10.1109/DEVLRN.2002.1011835]
[13] Hoyer P O. Non-negative matrix factorization with sparseness constraints[J]. The Journal of Machine Learning Research, 2004, 5: 1457–1469.
[14] Lin C J. Projected gradient methods for nonnegative matrix factorization[J]. Neural Computation, 2007, 19(10): 2756–2779. [DOI:10.1162/neco.2007.19.10.2756]
[15] Yang Z R, Oja E. Linear and nonlinear projective nonnegative matrix factorization[J]. IEEE Transactions on Neural Networks, 2010, 21(5): 734–749. [DOI:10.1109/TNN.2010.2041361]
[16] Wang D, Lu H C. On-line learning parts-based representation via incremental orthogonal projective non-negative matrix factorization[J]. Signal Processing, 2013, 93(6): 1608–1623. [DOI:10.1016/j.sigpro.2012.07.015]
[17] Zhang Z Y. Nonnegative matrix factorization:models, algorithms and applications[M]//Holmes D E, Jain L C. Data Mining:Foundations and Intelligent Paradigms:Volume 2:Statistical, Bayesian, Time Series and other Theoretical Aspects. Berlin Heidelberg:Springer, 2012:99-134.[DOI:10.1007/978-3-642-23241-1_6]
[18] Sun D L, Févotte C. Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[C]//Proceedings of 2014 IEEE International Conference on Acoustics, Speech and Signal Processing. Florence, Italy:IEEE, 2014:6201-6205.[DOI:10.1109/ICASSP.2014.6854796]
[19] Li Y N, Zhang X W, Sun M, et al. Speech enhancement based on robust NMF solved by alternating direction method of multipliers[C]//Proceedings of the 17th International Workshop on Multimedia Signal Processing. Xiamen, China:IEEE, 2015:1-5.[DOI:10.1109/MMSP.2015.7340815]
[20] Wang D, Lu H C. Incremental orthogonal projective non-negative matrix factorization and its applications[C]//Proceedings of the 18th International Conference on Image Processing. Brussels:IEEE, 2011:2077-2080.[DOI:10.1109/ICIP.2011.6115890]
[21] Gabay D, Mercier B. A dual algorithm for the solution of nonlinear variational problems via finite element approximation[J]. Computers & Mathematics with Applications, 1976, 2(1): 17–40. [DOI:10.1016/0898-1221(76)90003-1]
[22] Boyd S, Parikh N, Chu E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends® in Machine Learning, 2011, 3(1): 1–122. [DOI:10.1561/2200000016]
[23] Boyd S, Vandenberghe L, Faybusovich L. Convex optimization[J]. IEEE Transactions on Automatic Control, 2006, 51(11): 1859. [DOI:10.1109/TAC.2006.884922]
[24] Mei X, Ling H B. Robust visual tracking using ℓ1 minimization[C]//Proceedings of the 12th International Conference on Computer Vision. Kyoto, Japan:IEEE, 2009:1436-1443.[DOI:10.1109/ICCV.2009.5459292]
[25] Wu Y, Lim J, Yang M H. Object tracking benchmark[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2015, 37(9): 1834–1848. [DOI:10.1109/TPAMI.2014.2388226]
[26] Ross D A, Lim J, Lin R S, et al. Incremental learning for robust visual tracking[J]. International Journal of Computer Vision, 2008, 77(1-3): 125–141. [DOI:10.1007/s11263-007-0075-7]
[27] Bao C L, Wu Y, Ling H B, et al. Real time robust L1 tracker using accelerated proximal gradient approach[C]//Proceedings of 2012 IEEE Conference on Computer Vision and Pattern Recognition. Providence, RI:IEEE, 2012:1830-1837.[DOI:10.1109/CVPR.2012.6247881]
[28] Oron S, Bar-Hillel A, Levi D, et al. Locally orderless tracking[J]. International Journal of Computer Vision, 2015, 111(2): 213–228. [DOI:10.1007/s11263-014-0740-6]