何希平,杨 劲,刘 波
(1.重庆工商大学 电子商务及供应链系统重庆市重点实验室,重庆400067;2.重庆工商大学 计算机科学与信息工程学院,重庆400067;3.重庆工商大学 重庆市检测检测控制集成系统工程实验室,重庆400067)
对大规模数据的统计分析、稀疏低秩近似[1,2]、降维表达[3]与压缩感知等处理,是目前深受研究人员青睐的技术和方法。线性运算的高时效性使得如何有效地刻画研究对象之间局部或整体的近似线性关系成为数据分析过程中倍受重视并起着决定性作用的理论与技术问题。矩阵低秩逼近[4-8]由于其简洁性,已成为数据处理的基本手段之一。
在高维数据处理中,高维数据的相关性不容忽视,正是由于奇异值向量直接反映矩阵秩的大小,即矩阵中行或列向量的相关度,使得基于高维数据矩阵的奇异值分解(singular value decomposition,SVD)成为了信号主成分分析、高维数据降维、线性判别分析、稀疏表示、矩阵低秩分解等应用的重要技术手段。用奇异值分解技术实现信号降噪[9-11]即为线性近似的代表应用之一。
本文重点研究信号SVD重建的通用技术,其预期的主要目标包括:①构建信号基于SVD分解的低秩线性表达与重建的通用模型;②研究并实验一维信号的强相关结构矩阵的构建与SVD线性化表达模型;③通过含噪声的一维数字信号和二维数字图像的低秩近似与去噪声实验展示模型的适应性。
设X= (x1,…,xn)∈Rm×n为m 个n维信号的抽样矩阵,M为X的列均值矩阵,X的中心化信号矩阵Y=X
对信号X的SVD分解可转化为对Y的分解。
定理 设Y∈Rm×n为X的中心化信号矩阵,秩为r,则必存在正交矩阵U= (u1,…,um)∈Rm×m和正交矩阵V=(v1,…,vn)∈Rn×n,使得
称式 (1)为矩阵Y的奇异值分解 (简记为SVD),σi称为Y的奇异值。由式 (2)和式 (3)可见,YTY为X的协方差矩阵的m-1倍,且YTY的非零特征值σ2i的对应特征向量为vi,YYT的非零特征值σ2i的对应特征向量为ui(i=1,…,r),从而Y的SVD分解可通过对YTY和YYT进行特征分解实现。u1,…,um构成Rm空间的一组规范正交基,v1,…,vn构成Rn空间的一组规范正交基。
从中心化信号矩阵Y的奇异值分解式 (1)可以得到精确重建X的一些非常有用的信息
其中
式 (7)为X的低秩重建近似信号,式 (6)给出的是一个非凸优化模型,求解困难。借助拉格朗日乘数法,可以将秩不超过X的秩r的约束转化为在目标函数中引入一个核范数的正则化项
则这一优化问题可以借助拟梯度分析获得解的解析表达式[7]。而信号的恢复率η则可按式 (9)给出的基础成分的系数向量Λ与满秩奇异值分解的系数向量∑的相似性度量进行近似估算
显然,1≥η≥0。当Λ=∑时,η=1,则可通过式 (7)获得信号X的精确重建;否则,通过式 (8)的解析解抑制基础成分Yi的系数λi的大小,控制对原始信号的近似表达程度;同时由于奇异值与信号方差的紧密关系,也可以通过抑制λi有效地抑制信号噪声。
最常见的源信号是一维和二维信号,这也是本文重点研究的对象。根据不同的目的,如去噪声、稀疏表示、数据压缩等,可以通过式 (8)选择不同的线性系数来组合原始信号的SVD基础成分,进而利用式 (7)完成对信号X的重建。可见,重建算法的第一步应该是构建原始信号X的中心化信号矩阵Y,为利用式 (7)的信号重建模型,首先还必须把一维信号向量变换成对应的结构矩阵。利用一维信号可以构造出很多种矩阵,如Hankel矩阵、Toeplitz矩阵、Circulant矩阵等。研究实践结果表明,矩阵构造方式不同,则SVD的信号处理效果就会不一样。
信号重建的重要依据之一就是相邻信号时序的结构相关性,因此,重建一维信号也可通过构建一维数字信号时序的邻域相关结构矩阵,再利用结构矩阵的SVD低秩逼近模型实现重建。
设一维离散数字信号的时序x= (x1,…,xm)T∈Rm,则可以利用此信号构造如下的Hankel矩阵X (各行时序间或者列时序间依次超前一个时间步长)
则X的第一列和最后一行 (或第一行和最后一列)元素构成原离散数字信号x。Hankel矩阵X的相邻元素较好地反应了信号时序的相邻相关的特点。
也可用一维离散数字信号x构造Toeplitz矩阵X (同一条对角线上的元素相同)
则X的第一列和第一行元素构成离散数字信号的时间序列x。
而式 (12)表示的信号x的Circulant矩阵完全由它的第一列即x(或第一行即x′)确定,其余各列均为其前一列元素顺时针 (或逆时针)循环移一位得到,即以上一列第二 (或最后)元素为首元素依序循环排列 (相当于上一时序循环提前或延迟一个时间步长)
对于一维离散数字信号x,有了结构矩阵X以后,就可通过式 (1)~式 (4)计算结构矩阵X的中心化信号矩阵Y的SVD分解和基础成分。由式 (4)~式 (8)不难发现,X的重建信号是Y的SVD基础成分Yi的线性组合。若分别记Yi中的、与所构造的结构矩阵X中构成信号时序x的列和行对应的列和行向量 (行向量不含第一元素)为pi和qi,并记则x的重建信号可表示为
式中:μ——由结构矩阵X的列均值构成的与x同维度的向量。实际应用中,μ可直接取x的均值构成。
直观观察可见,在式 (10)~式 (12)表示的3种结构矩阵中,Hankel矩阵元素的相邻相关性最强,Toeplitz矩阵次之,Circulant矩阵元素的相邻相关性相对较弱。使用Hankel矩阵与Toeplitz矩阵的计算复杂性与参数k的选取有密切关系,k越小,SVD计算复杂性越低;反之,复杂性越高。当信号时序较长时,Circulant矩阵的阶数较大,SVD的计算复杂度呈指数级增长,效率较低。
根据前面分析,重建算法的具体步骤如下:
(1)由信号矩阵X求得中心化矩阵Y;
(2)对中心化矩阵Y进行SVD分解,得到奇异值向量∑和Y的基础成分Y1,…,Yr;
(3)求解优化模型 (8),得到Λ= (λ1,…,λr)T;
(4)根据式 (7)计算重建近似信号。
在重建模型求解过程中,重要的是对优化模型 (8)中的Λ= (λ1,…,λr)T的选择。针对不同应用可选择不同的信号恢复率η0,并在恢复率η0的控制下,基于主成分表达的思想,使用贪心法首先选取一个非零的λi,使得仅由基础成分Yi重建的X的恢复率η最接近η0,且η0≥η,并令λk=0 (k<i);然后再选取下一个非零的λj(i<j),使得由基础成分Yi和Yj重建的X的恢复率η最接近η0,且η0≥η,并令λk=0 (i<k<j);如此继续选取下一Yj。在η还未接近η0的大小时,均可取λj=σj,以便算法快速收敛,同时使Λ取值表现出稀疏特性。
利用SVD分解模型 (4)和重建模型 (7),可以对信号进行主成分分析,或者利用基础成分对信号进行近似或者稀疏压缩表示 (选取较少的几个基础成分近似表示原信号),或者消除信号噪声等。
对一维数字信号去噪,一般考虑如下加性白噪声模型
其中,ni~N(0,σ2)(i=1,2,…,m)为独立同分布的高斯白噪声。
由式 (1)、式 (2)、式 (4)可见,信号SVD分解模型实质上是对信号矩阵左右两边各乘以一个正交矩阵而得到信号的奇异值矩阵,而作为其逆变换的信号重建过程正好是以奇异值为所求基础成分的权的线性表出过程,从而噪声在变换过程中保持可加性未变,并且噪声通常表现为高频尖峰数据,因此在去噪声应用中,可以利用软硬阈值模型[12]作用于奇异值向量∑,求得重建模型 (7)中的基础成分的系数向量Λ,进而恢复原始信号。
去噪声的过程实质上就是选择一种误差度量模型,通常是均方误差 (MSE)
寻求与含噪信号x以MSE最小为目标的原始信号x0的估计
在SVD去噪中,x是根据式 (13)表示的模型得到的结果。
由于图像噪声与一维信号相似,更常见的也是加性高斯白噪声,而在SVD重建模型中,二维图像即构成了信号矩阵X,所以,可以直接在类似于式 (16)的MSE最小化的优化目标模型控制下,利用SVD重建算法实现图像去噪声,也可以对图像进行分块,然后通过分块线性近似实现。
基于主成分分析和贪心算法的思想,可以按式 (9)表示的恢复率计算模型,根据稀疏度要求按从前至后的顺序选择基础成分及其对应奇异值,而舍弃后面重要性相对较低的成分,可以有效实现信号的稀疏k近似表示
其中
借助式 (17)的模型,可以简便地实现对原始信号进行压缩与稀疏表示。事实上,由于k个基础成分 (通常k<<r)仅需要由k对ui、vi生成,所以信号存储要求由o(mn)降低为o(k(m+n)),可以大幅压缩信号的存储量,实现信号的稀疏近似表示。
对一维信号去噪声,我们选择了长度为2048的Heavy sine离散数字信号作为原始信号,并按式 (14)添加了高斯噪声,且取σ=1.5。
在具体算法中,我们对Λ的选取采用的一种策略是如下主成分选取模型
为进一步分析结构矩阵的特性,我们对同一信号的3种结构矩阵SVD低秩重建结果做了统计计算,并与小波软阈值自动去噪结果进行对比,结果见表1,其信号处理直观对比如图1所示。
表1 一维信号去噪结果的统计对比
其中小波变换选取的是5层sym8,阈值选取方法为SURE。从表中数据可见:①Hankel矩阵、Toeplitz矩阵、Circulant矩阵分别表达信号矩阵时,前两种矩阵比第3种更有利于用SVD对信号进行阈值去噪声和近似表达;②借助Hankel矩阵或者Toeplitz矩阵表达信号时,通过对信号进行SVD阈值去噪的效果明显优于小波阈值去噪的效果,且两者通过SVD实现的效果非常相近;③由式 (9)表达的恢复率计算模型只能用于近似表达信号间的相似程度,可作为算法选择时简化的近似计算依据,但不能严格把它作为信号相似性好坏的判据;④Circulant矩阵由于误差偏大的原因,不适合于信号的SVD去噪与最佳近似表达应用场合。
图1 一维信号去噪效果对比
虽然去噪声并非本文研究的最终目标,但是图1所示的直观效果表明对信号进行SVD阈值去噪的效果明显优于小波阈值去噪的效果。另一方面,SVD阈值去噪的效果与阈值选择和线性近似时系数调节算法也有关。实验中按式(19)进行系数选择,是一种比较粗略的处理方法,实用中可以在优化模型 (16)的控制下动态调节各系数以达到更加理想的效果。
为检验SVD阈值去噪对二维图像信号处理的效果,我们选择了两幅512×512的灰度图像,均添加了均值为0、标准差为16的高期白噪声,然后分别用小波软阈值自适应去噪与SVD阈值去噪进行效果对比,直观效果如图2所示,而统计结果见表2。其中,与一维信号类似,对经加高斯白噪声的图像选用sym8小波做3层小波分解,对所获得的细节系数选用自适应软阈值小波收缩;而SVD阈值去噪声时阈值为奇异值的均值。
图2 图像阈值去噪效果对比
表2 噪声图像阈值去噪结果对比
由表2所示的3种统计图像特征不难发现,尽管简化了系数优化过程,使得SVD去噪算法中利用奇异值均值阈值对近似系数的选择只是一种粗略估计,然而SVD去噪效果也与自适应软阈值小波去噪声结果非常接近,如果进一步优化SVD近似的逼近算法,提高去噪效果是必然的。
同时,为了测试k近似压缩表示的效果和优化近似系数选择的意义,我们对图2中的两幅加噪图像在相似度模型 (17)的控制下,分别取η0=0.85、0.90和0.95的3种不同值时,自适应地选取了不同的k值,将对应的重建结果与原始图像、自适应软阈值小波去噪结果图像进行了统计与直观对比。实验所得到的统计特征见表3,直观对比效果如图3所示。
表3 噪声图像k近似结果的统计特征对比
图3 图像k近似效果与小波去噪对比
由表3所示的图像的3种统计特征不难发现:①SVD逼近中,即使要求恢复率η0≥0.85,相对于奇异值个数来说,信号的k近似中的线性组合基础成分数k也比较小,说明SVD线性组合模型对信号逼近与压缩表示非常有效;②选择恰当的恢复率η0,不仅可以通过SVD低代价压缩近似表示噪声信号,而且还可能达到比自适应软阈值小波去噪声更好的信号处理效果;③表3第5行数据也说明,当一味追求过高的恢复率时,不仅线性表达的代价会提高,而且消除噪声的效果也会降低。
本文提出的基于信号矩阵的SVD分解的线性近似的优化重建模型,就其本质是对信号进行相关性分析,寻求信号线性表达的低秩矩阵。文中所给出的数学模型不仅可以用于一维和二维信号去噪声,而且可以用于信号的稀疏与压缩表示,表达式简单易于计算;所给出的模型的SVD奇异值均值阈值去噪声和k近似逼近算法的模型都简便易行。实验结果表明,新优化重建模型对一维信号和二维图像的噪声抑制均能获得很好的统计特征,可靠与稳定性高,实用性强。
[1]Meng Deyu,Zhao Qian,Xu Zongben.Improve robustness of sparse PCA by L1-norm maximization [J].Pattern Recognition,2012,45 (1):487-497.
[2]Haipeng Shena,Jianhua Z Huang.Sparse principal component analysis via regularized low rank matrix approximation [J].Journal of Multivariate Analysis,2008,99 (6):1015-1034.
[3]David Nelson,Siamak Noorbaloochi.Information preserving sufficient summaries for dimension reduction [J].Journal of Multivariate Analysis,2013,115 (3):347-358.
[4]Ma S,Goldfarb D,Chen L.Fixed point and bregman iterative methods for matrix rank minimization [J]. Math Progra,2011,128 (1-2):321-353.
[5]Cai J F,Candès E J,Shen Z W.A singular value thresholding algorithm for matrix completion [J].SIAM J Optim,2010,20:1956-1982.
[6]RAO Guo,PENG Yi,XU Zongben.Robust sparse and lowrank matrix decomposition based on S1/2modeling [J].Science China Information Sciences,2013,43 (6):733-748 (in Chinese).[饶过,彭毅,徐宗本.基于S1/2建模的稳健稀疏-低秩矩阵 分 解 [J]. 中 国 科 学: 信 息 科 学,2013,43 (6):733-748.]
[7]ZHANG Zhenyue,ZHAO Keke.Regularization methods for low-rank factorization of matrices with missing data [J].Scientia Sinica Mathematica,2013,43 (3):249-271 (in Chinese).[张振跃,赵科科.数据缺损矩阵低秩分解的正则化方法 [J].中国科学:数学,2013,43 (3):249-271.]
[8]PENG Yigang,SUO Jinli,DAI Qionghai,et al.From compressed sensing to low-rank matrix recovery:Theory and applications [J].ACTA Automatica Sinica,2013,39 (7):981-994(in Chinese).[彭义刚,索津莉,戴琼海,等.从压缩传感到低秩矩阵恢复:理论与应用 [J].自动化学报,2013,39(7):981-994.]
[9]HUANG Jianzhao,XIE Jian,LI Feng,et al.Com pound denoising method based on discrete stationary wavelet transform and optimized SVD [J].Piezoelectrics & Acoustooptics,2013,35 (3):448-451 (in Chinese). [黄建招,谢建,李锋,等.基于优化SVD和平稳小波的复合降噪方法 [J].压电与声光,2013,35 (3):448-451.]
[10]CHEN Enli,ZHANG Xi,SHEN Yongjun,et al.Fault diagnosis of rolling bearings based on SVD denoising and blind signals separation [J].Journal of Vibration and Shock,2012,31 (23):185-190 (in Chinese). [陈恩利,张玺,申永军,等.基于SVD降噪和盲信号分离的滚动轴承故障诊断[J].振动与冲击,2012,31 (23):185-190.]
[11]ZHAO Xuezhi,YE Bangyan,CHEN Tongjian.Principle of singularity detection based on SVD and its application [J].Journal of Vibration and Shock,2008,27 (6):11-14 (in Chinese). [赵学智,叶邦彦,陈统坚.基于SVD的奇异性信号检测原理及其应用 [J].振动与冲击,2008,27 (6):11-14.]
[12]HE Xiping,YANG Jin.A novel type of threshold functions for wavelet shrinkage [J].Computer Engineering &Science,2013,35 (8):114-119 (in Chinese). [何希平,杨劲.一类新的小波收缩阈值函数 [J].计算机工程与科学,2013,35 (8):114-119.]