基于Rayleigh 商的次子空间准则跟踪函数

徐中英，高迎彬，孔祥玉，杜伯阳

（火箭军工程大学控制工程系，陕西 西安 710025）

1 引言

在现代信号处理和数据分析领域，次成分分析是一种常用的分析方法。次成分代表了信号或数据具有最小偏差的方向，由多个次成分张成的空间称为次子空间。次成分分析是指利用次成分或次子空间进行信号处理的方法[1]。次成分分析方法可以应用在很多领域，如滤波器设计[2]、Pisarenko 频率估计[3]、阵列天线设计[4]、总体最小二乘估计[5]、波达方向（DoA,direction of arrival）估计[6]等。

早期的次成分分析是通过特征值分解（EVD,eigenvalue decomposition）或者奇异值分解（SVD,singular value decomposition）来实现的，这类算法本质上是采用批处理方式，存在计算量大和无法实时运算的问题[7]。因此，有研究者提出基于神经网络的算法，该类算法的优点有：1)避免了自相关矩阵的直接计算，2)能够实现非平稳信号的跟踪，3)能够处理高维数据[7]。基于神经网络的算法自提出以来迅速成为国内外研究的热点。

目前，学者们提出了许多次成分分析算法，如Möller 算法[8]、SDPM（stable data projection method）[9]、PAST（projection approximation subspace tracking）算法[10]、Douglas 算法[11]等，这些算法都是基于启发式推理得到的，缺乏对应的准则函数。信息准则规定了次成分分析方法的搜索方式，利用准则函数可以快速确定算法的全局收敛域，因此发展准则函数具有十分重要的意义。目前，已经提出的准则函数有AMEX（adaptive minor component extraction）函数[12]、Hasan 函数[13]、OJAm 函数[3]等。相比丰富的次成分分析算法而言，准则函数并不多见，需要进一步丰富和发展。

本文的主要工作是，通过对Rayleigh 商函数添加适当的惩罚项，提出了一个新型的次子空间准则函数；通过对所提准则函数的平稳点进行分析，证明了准则函数的最优解是次子空间的一组基，构建了准则函数与次子空间之间的联系；通过梯度上升法导出了一个次子空间跟踪算法；通过李雅普诺夫函数法确定了导出算法的全局收敛域。

2 新型准则函数的提出

2.1 预备知识

将矩阵R的特征值按照从小到大的顺序排列，即

特征值对应的特征向量也相应排列，可以得到R的另外一种特征值分解表示，如式(3)所示。

其中，U1是由最小的r个特征值对应的特征向量构成的矩阵，U2则是由剩余的n-r个特征值对应的特征向量构成的矩阵。根据信号处理理论[14]可知，自相关矩阵R的最小的r个特征值对应的特征向量张成的空间等于信号的次子空间。

2.2 新型准则函数

其中，W是神经网络的权矩阵，tr(·)是矩阵的迹。J(W)是一个无约束优化函数，它由两部分构成，第一部分是Rayleigh 商函数，作用是使算法能够收敛到需要的次子空间；第二部分是一个惩罚函数，作用是对权矩阵模值施加一个隐形约束，保证算法迭代过程中模值的收敛性。

3 准则函数的平稳点分析

对于准则函数式(4)而言，是否存在全局极大值、是否可以建立与次子空间的联系、是否存在其他局部极值是3个非常重要的问题，本文以定理1和定理2对上述问题进行回答。

定理1在域内，当且仅当时，W是准则函数J(W)的平稳点，其中Q∈ℝr×r是一个正交矩阵。

证明由于WTRW＞0和WTW≠0均是对称正定矩阵，因此都是可逆矩阵。由式(4)可得，J(W)对矩阵W的一阶微分存在，且有

反之，由平稳点定义可得，在J(W)的平稳点有∇J(W)=0，即

式(7)同时左乘WT并化简得

证毕。

定理2在域内，当且仅当时，准则函数J(W)达到全局极大值。在全局极大值处，其他所有的平稳点都是J(W)的鞍点，其中任意一个正交矩阵。

证明根据定理1可知，任意的均是准则函数J(W)的平稳点，其中Ur是由R的任意r个特征向量构成的矩阵。令Ur中特征向量的位置序号构成的集合为S1，即S1={i1,i2,…,ir}。同理，S2={1,2,…,r}为U1中特征向量的位置序号构成的集合。

对于任意的S1(S1≠S2)而言，该集合中必定存在某一个元素j，满足

将矩阵Ur中的第j列向量替换为vj+εvk（k∈S2且k∉S1），则形成新矩阵。显然，有λj＞λk。令M1=[0,…,0,vk,0,…,0]，即矩阵M1中只有第j列为vk，而其他列都是0，则

进一步有

将式(13)和式(14)代入式(15)有

由于λj＞λk，则根据式(15)可得，在平稳点处，沿着向量vk的方向，J(W)是递增的。此外，将矩阵Ur中的vj替换为vj+εvj，则可以获得新的矩阵。令，即矩阵M2中只有第j列为vj，而其他列都是0，则

其中，D1=diag(0,…,0,1,0,…,0)是一个对角矩阵。

进一步有

根据式(14)和式(19)可得

证毕。

定理1和定理2表明，J(W)存在全局极大值，在全局极大值时，W是次子空间的一组基，且WTRW=I可以被唯一确定。由于J(W)只有全局极大值而没有其他局部极值，因此通过迭代法（如梯度法）导出的算法可以保证收敛到需要的次子空间。

4 次子空间跟踪算法及其收敛域分析

4.1 子空间跟踪算法

通过定理1可得，式(5)是J(W)的一阶微分。假设在第k次迭代时，神经网络的权矩阵为Wk，则通过利用梯度上升法可以给出下一时刻的权矩阵更新式，如式(22)所示。

其中，μ是算法的学习因子，满足0＜μ＜1。经过多次迭代运算后，神经网络权矩阵将最终收敛到次子空间的一组基。

4.2 算法的收敛域

本节对导出算法式(22)的收敛性进行讨论。当学习因子μ足够小时，离散差分式(22)可以近似为与之相对应的连续微分方程[15]。

其中，t=kμ。由此，算法的全局收敛性可以通过对该连续微分方程对应的动态系统分析来完成。根据李雅普诺夫第二定律，动态系统的收敛性分析主要涉及以下两方面：1)动态系统是否能够全局收敛到次子空间；2)动态系统可以吸引的范围是多少，即什么样的初始化条件能够保证系统的全局收敛性。

定理3将给出上述2个问题的答案。

定理3给定常微分方程式(23)，假如初始化全矩阵满足W(0)∈Ω，则当t→∞时，W(t)必将全局渐进收敛到集合中的一个点，其中Q∈ℝr×r任意一个正交矩阵。

证明令L(W)=-J(W)。根据定理1和定理2，L(W)与J(W)有相同的平稳点，且在处，L(W)取得全局极小值。利用求导链式法，则

将式(5)和式(23)代入式(24)，则很容易得到，对于任意的W∈Ω，均有。因此L(W)构成了式(23)的李雅普诺夫函数。根据李雅普诺夫第二定律，从任意初始化权矩阵出发的动态系统都收敛到相同的不变集由于是不稳定的鞍点，因此W(t)必将全局渐进收敛到中的一个点。

证毕。

5 仿真实验

本节通过3个仿真实验来验证所提准则函数和导出算法的正确性。第一个实验提供了一种特殊情况下的准则函数曲线，第二个实验展示了导出算法在次成分分析方面的优越性，第三个实验是算法在DoA 估计中的应用。

5.1 准则函数曲线

令准则函数J(W)中R=diag(2,1)。当提取子空间维数为1，即r=1时。如果W=[W11,W12]T，则通过适当简化有

图1是利用Matlab 绘出上述情况下的J(W)曲线。从图1中可以看出，J(W)存在全局极大值，而没有局部极值，这与定理2的分析结果相一致。当权矩阵W=[0,±1]T时，准则函数J(W)取得全局极大值，此时W=[0,±1]T正好是R=diag(2,1)的次成分，从而证明了所提准则函数的正确性。

图1 准则函数的3D 曲线

5.2 次成分分析实验

5.1 节实验主要通过一个特例来证明所提信息准则的正确性，本节实验将考察准则函数导出算法进行次成分分析的能力。首先，利用文献[10]的方法随机生成一个5×5对称正定矩阵，如式(26)所示。

该矩阵的最小特征值为λ1=0.8348，其对应的次成分特征向量为

然后，分别利用导出算法、OJAm 算法[3]和Douglas 算法[11]对该矩阵的次成分进行提取。为了定量描述算法的性能，这里引入2个评价函数[10]。

1)方向余弦

方向余弦（DC,direction cosine）实际上是权矩阵与次成分方向之间的夹角，如式(28)所示。如果DC=1，则表示Wk与1v方向重合。

2)权矩阵过程模值

对于OJAm 算法和Douglas 算法，已经证明算法收敛时WTW=I，因此过程模值取；对于导出算法而言，根据式(8)有WTRW=I，因此取

在迭代过程中，为了保证算法公平比较，3个算法采用相同的初始化权矩阵和学习因子。本实验中，初始化权矩阵是随机产生的（矩阵的每个元素均是高斯白噪声），学习因子μ=0.01。3种算法的仿真结果如图2和图3所示，该结果曲线是100次独立仿真结果的平均值。

从图2中可以看出，所提算法的方向余弦曲线最终收敛到了单位1，即所提算法能够收敛到次成分的方向。从图3中可以看出，迭代过程中，权矩阵模值是有界的且最终也收敛到了1，与定理1的分析结果相吻合。因此，综合图2与图3结果可知，所提算法具备次成分分析的能力。对比图2和图3中3种算法的表现可知，相比OJAm 算法和Douglas算法，本文所提算法在方向余弦和权矩阵模值方面具有较快的收敛速度。

图2 方向余弦曲线

图3 权矩阵模值曲线

5.3 DoA 估计

假设一个具有8个阵子的线性等距天线阵列，阵子的间距为半个波长。3个远场信号的入射角（即DoA 角）分别为12°、19°和27°。信号中的噪声为加性高斯白噪声，且信噪比为20 dB。分别利用所提算法、OJAm 算法[3]和Douglas 算法[11]对信号的次子空间进行估计，通过MUSIC 法计算DoA 角。本节实验中，3种算法的初始化方法与实验2相同，初始化权矩阵是随机产生的，学习因子μ=0.1。为了评价算法对子空间的估计结果，算法迭代过程中分别计算3个算法估计子空间与单位矩阵之间的正交误差，如式(29)所示。

其中，对于 Douglas 算法和 OJAm 算法有A=WTW；对于所提算法有A=WTRW。图4是3个算法的正交误差曲线，图5是利用导出算法得到的谱密度曲线。

图4 正交误差曲线

图5 所提算法的DoA 估计曲线

从图4中可以看出，经过大约100次迭代运算后，所提算法的正交误差已经趋于0，即权矩阵已经收敛到了信号的次子空间。通过与其他2个算法对比可以发现，所提算法的收敛速度要快于其他2种算法。从图5中可以看出，谱密度曲线分别在3个DoA 角处取得极值。因此可以得出结论，所提算法能够有效地解决DoA 估计问题，而且收敛速度优于同类型其他算法。

6 结束语

次成分分析是信号处理和数据分析领域一门重要的工具，而寻找准则函数在次成分分析算法中具有十分重要的地位。本文首先提出了一个新型的次子空间跟踪准则函数，并通过分析准则函数平稳点建立起了准则函数与次子空间之间的关系；基于此信息准则构建了一个新的次子空间跟踪算法，并证明了算法的全局收敛性。仿真实验表明，与一些现有算法相比，本文提出的算法具有较快的收敛速度。