Multiple Generalized Eigenvalue Decomposition Algorithm in Parallel Based on Weighted Matrix
-
摘要: 针对串行广义特征值分解算法实时性差的缺点, 提出基于加权矩阵的多维广义特征值分解算法. 与串行算法不同, 所提算法能够在一次迭代过程中并行地估计出多维广义特征向量. 平稳点分析表明: 当且仅当算法中状态矩阵等于所需的广义特征向量时, 算法达到收敛状态. 通过对比相邻时刻的状态矩阵模值证明了所提算法的自稳定特性. 所提算法参数选取简单, 实际实施较为容易. 数值仿真和实例应用进一步验证了算法的并行性、自稳定性和实用性.Abstract: In order to overcome the disadvantages of sequential algorithms, such as poor real time, a multiple generalized eigenvalue decomposition algorithm is proposed based on weighted matrix method. Unlike sequential algorithms, the proposed algorithm is able to estimate multiple generalized eigenvectors in parallel only through one iteration procedure. The stationary point analysis shows that the algorithm reaches convergence state if and only if the state matrix is equal to the desired generalized eigenvectors. The self-stabilization characteristics of the proposed algorithm is proved by comparing the state matrix module values of adjacent moments. The proposed algorithm parameters are simple to select and easy to implement in practice. Numerical simulation and example application further verify the parallelism, self-stability and practicality of the algorithm.
-
广义特征值分解是现代信号处理领域内重要的分析工具, 已经广泛应用于多个领域[1-6]. 在工程实际中, 不同的信号处理问题需要广义特征值分解的维数是不一样的. 根据提取广义特征向量的维数, 广义特征值分解可以分为3类[7]: 1)仅能估计出矩阵束最大(小)广义特征值对应的广义特征向量的单维分解算法[8]; 2)能够估计出若干个广义特征向量张成空间的子空间跟踪算法[9]; 3)能够准确估计任意个广义特征向量的多维分解算法[10].
假定存在矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$, 其中, ${{\boldsymbol{R}}_y}$和${{\boldsymbol{R}}_x}$分别是两个$n \times n$维对称正定矩阵, 则广义特征值分解就是要计算出一个$n \times 1$维向量${\boldsymbol{w}}$和一个标量$\lambda $, 使得
$${{\boldsymbol{R}}_y}{\boldsymbol{w}} = \lambda {{\boldsymbol{R}}_x}{\boldsymbol{w}}$$ (1) 满足上述方程的向量${\boldsymbol{w}}$和标量$\lambda $分别记为矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$的广义特征向量和广义特征值. 基于特征值分解的算法是进行广义特征值分解的传统方法, 而基于神经网络的算法是近些年来涌现出的一类新型算法. 相比传统算法, 神经网络类算法具有计算量低、实时性强、能够处理非平稳信号等优点[11], 已经成为广义特征值分解领域内一个主流的研究方向.
基于神经网络, 学者们提出了很多优秀的广义特征值分解算法. 例如RNNM (Recurrent neural network model)算法[12]、R-GEVE (Reduced-rank generalized eigen-decomposition extraction)算法[13]、基于NOCS (Nested orthogonal complement structure)架构的算法[9]、ANQ (Adaptive normalized quasi-Newton)算法[14]、GChen (Generalized Chen algorithm)算法[8]、GDM (Generalized Douglas minor component analysis)算法[10]等. 在上述算法中, RNNM算法、R-GEVE算法、ANQ算法和GCHen算法是单维广义特征值分解算法; NOCS算法属于子空间跟踪算法; 只有GDM算法是多维分解算法. 由于广义特征向量可以看作子空间的一组特殊基, 因此多维分解算法也适用于子空间算法的领域; 同时当多维分解算法提取维数为1时, 多维分解算法退化为单维分解算法. 综上可知, 多维分解算法具有最广泛的应用范围[10], 即研究多维分解算法更具有普适性.
在多维广义特征值分解算法领域主要存在两类算法: 串行算法和并行算法. 串行算法中多维广义特征向量是依次串行获取的, 即先用单维分解算法计算出最小广义特征值对应的广义特征向量, 然后利用膨胀技术再计算出次小广义特征值对应的广义特征向量, 依次类推[10]. 文献[8]提出了一种压缩技术, 文献[10]详细分析了该压缩技术并对其进行了改进, 实现了多维广义特征向量的串行提取. 相比并行算法, 串行算法有如下两点不足: 1)由于串行算法中广义特征向量是串行依次提取, 整个提取过程需要时间较长, 因此串行算法不适合实时性要求较高的场合; 2)串行算法在广义特征向量的估计过程中会用到上一次的估计结果, 而前一次的估计误差会累积到本次估计过程中, 因此会造成估计误差的累积传播[11]. 然而, 目前多维广义特征值分解算法还不多见, 因此发展并行多维算法仍是非常有意义的工作.
基于加权矩阵的思想, 本文提出一个多维广义特征值并行分解算法, 并详细分析了算法的收敛性和稳定性. 相比串行算法, 所提算法可以实现多维广义特征向量的并行计算, 实时性较强. 本文结构安排如下: 第1节提出新型广义特征值分解算法并对其进行收敛性和稳定性分析; 第2节是仿真实验和实例验证; 第3节是全文的总结.
1. 新型广义特征值分解算法
1.1 算法的提出
基于神经网络模型, Oja算法[15]首次提出了下述特征值分解算法:
$$\begin{split} &{\boldsymbol{w}}(k + 1) = {\boldsymbol{w}}(k) \;+ \\ &\qquad\eta \left[ {{\boldsymbol{Rw}}(k) - {{\boldsymbol{w}}^{\rm{T}}}(k){\boldsymbol{Rw}}(k){\boldsymbol{w}}(k)} \right] \end{split} $$ (2) 其中, ${\boldsymbol{R}} \in {{\bf{R}}^{n \times n}}$是信号的自相关矩阵, ${\boldsymbol{w}} \in {{\bf{R}}^{n \times 1}}$是神经网络的状态向量, $\eta $是算法的学习因子, $n$为输入信号的维数. 然而该算法只能估计单个矩阵${\boldsymbol{R}}$的特征向量, 并非矩阵束的广义特征向量; 此外该算法只能实现单维分解, 不能进行多维分解. Oja算法是特征值分解算法领域最为经典的算法, 很多现有算法在基础理论方面和Oja算法是相同的, 因此如果以Oja算法为代表进行研究并成功将其扩展为多维广义特征向量估计算法, 该研究结果对其他算法势必具有很强的借鉴意义.
假定${{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x} \in {{\bf{R}}^{n \times n}}$是两个信号的自相关矩阵, 其广义特征值和广义特征向量分别为${\lambda _i}$和${{\boldsymbol{v}}_i}$, 其中, $i = 1,2, \cdots , n$. 为后续方便使用, 这里将矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$的$n$个正广义特征值进行降序排列, 即
$${\lambda _1} > {\lambda _2} > \cdots > {\lambda _n} > 0$$ (3) 其对应的$n$个关于${{\boldsymbol{R}}_x}$正交的广义特征向量${{\boldsymbol{v}}_i}\;(i = 1, 2, \cdots ,n)$也相应进行排列. 根据矩阵理论[16]有
$${{\boldsymbol{R}}_y}{{\boldsymbol{v}}_i} = \lambda {{\boldsymbol{R}}_x}{{\boldsymbol{v}}_i}$$ (4) $${\boldsymbol{v}}_j^{\rm{T}}{{\boldsymbol{R}}_y}{{\boldsymbol{v}}_i} = {\delta _{ij}}, \;\; i,j \in \{ 1,2, \cdots ,n\} $$ (5) 其中, ${\delta _{ij}}$是Kronecker $\delta $函数.
为了计算矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$的多维广义特征向量, 通过对Oja算法添加加权矩阵, 提出如下的多维分解算法:
$$\begin{split} {\boldsymbol{W}}(k + 1) =\;& {\boldsymbol{W}}(k) + \eta [ {{{\boldsymbol{R}}_x}^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{W}}(k){\boldsymbol{D}}}\, -\\ & { {\boldsymbol{W}}(k){{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{\boldsymbol{W}}(k)} ] \end{split} $$ (6) 其中, ${\boldsymbol{W}} \in {{\bf{R}}^{n \times r}}$是神经网络的状态矩阵; ${\boldsymbol{D}} = {\rm{diag}} \{{d_1}, {d_2}, \cdots ,{d_r}\}$是一个对角矩阵, 其对角线元素满足${d_1} > {d_2} > \cdots > {d_r} > 0$, 这里将其称为加权矩阵; $r$是需要计算的广义特征向量的维数. 通过式(6)给定状态矩阵更新规则, 状态矩阵${\boldsymbol{W}}$进行更新迭代, 最终将收敛到矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$最大的$r$个广义特征值对应的广义特征向量.
通过对式(6)进行分析可以发现, 如果令式(6)中${{\boldsymbol{R}}_x} = {{\boldsymbol{I}}_n}$和${\boldsymbol{D}} = {{\boldsymbol{I}}_r}$, 且进行单维分解(即$r = 1$), 式(6)就会退化为Oja算法. 因此可以说式(6)是Oja算法的扩展. 式(2)是单维特征向量估计, 如果利用文献[10]中的压缩技术将其扩展为串行算法, 则其计算复杂度为 ${\rm{O}}({{{n^3}r} / 3} + 4{n^2}r \;+ 2nr)$, 而式(6)是多维广义特征向量估计, 其计算复杂度为${\rm{O}}({{{n^3}} / 3} + 2{n^2}r + 2n{r^2})$, 显然要低于串行算法的计算复杂度. 由于加权矩阵的加入, 使得状态矩阵在收敛过程中对广义特征向量与其构成子空间的夹角更为敏感[17], 进而通过在迭代过程中对神经网络状态矩阵施加的隐形Gram-Schmidt正交化(GS orthogonalization, GSO)操作, 使得状态矩阵能够收敛到广义特征向量, 而非其构成的子空间. 加权矩阵仅需满足约束条件${d_1} > {d_2} > \cdots > {d_r} > 0$, 在实际使用中该约束条件很容易达到(如可以取一等差数列), 因此, 式(6)是符合实际使用需要的.
1.2 算法的收敛性分析
算法是否能够准确地计算出矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$的广义特征向量, 是算法发展过程中首先要解决的问题, 本节将通过定理1证明算法的收敛性.
定理 1. 当且仅当神经网络的状态矩阵${\boldsymbol{W}} = {\boldsymbol{U}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}$时, 所提算法达到稳定的收敛状态, 其中, ${\boldsymbol{U}} = [{{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2}, \cdots , {{\boldsymbol{v}}_r}]$是由矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$最大的$r$个广义特征值对应的广义特征向量构成的矩阵.
证明. 定义矩阵函数
$$J({\boldsymbol{W}}) = {\boldsymbol{R}}_x^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{WD}} - {\boldsymbol{W}}{{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{W}}$$ (7) 则$J({\boldsymbol{W}})$代表着算法的学习步长. 当矩阵取值为${\boldsymbol{W}} = {\boldsymbol{U}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}$时
$$\begin{split} J({\boldsymbol{W}}) = \;&{\boldsymbol{R}}_x^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{U}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}{\boldsymbol{D}} - {\boldsymbol{U}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}{{\boldsymbol{U}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{U}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}= \\ \;& {\boldsymbol{U}}{{\boldsymbol{\varLambda }}_r}{{\boldsymbol{D}}^{\tfrac{3}{2}}} - {\boldsymbol{U}}{{\boldsymbol{\varLambda }}_r}{{\boldsymbol{D}}^{\tfrac{3}{2}}} = {\boldsymbol{0}} \\[-10pt] \end{split} $$ (8) 其中, ${{\boldsymbol{\varLambda }}_r} = {\rm{diag}}\{{\lambda _1},{\lambda _2}, \cdots ,{\lambda _r}\}$是一个对角矩阵, 其对角线元素由矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$最大的$r$个广义特征值构成.
当且仅当学习步长等于零时算法达到收敛状态, 即$J({\boldsymbol{W}}) = {\boldsymbol{0}}$, 进而有
$${{\boldsymbol{R}}_x}^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{WD}} = {\boldsymbol{W}}{{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{W}}$$ (9) 对式(9)同时左乘${{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_x}$, 并经过适当简化, 有
$${{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{W}}\left( {{\boldsymbol{D}} - {{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_x}{\boldsymbol{W}}} \right) = {\boldsymbol{0}}$$ (10) 由状态矩阵${\boldsymbol{W}}$的任意性, 可得
$${{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_x}{\boldsymbol{W}} = {\boldsymbol{D}}$$ (11) 令${\boldsymbol{P}} = {\boldsymbol{W}}{{\boldsymbol{D}}^{\tfrac{1}{2}}}$并将其代入式(11), 有: ${{\boldsymbol{P}}^{\rm{T}}}{{\boldsymbol{R}}_x}{\boldsymbol{P}} = {\boldsymbol{I}}$, 即矩阵${\boldsymbol{P}}$的列向量是关于矩阵${{\boldsymbol{R}}_x}$正交的. 矩阵${{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{W}}$的特征值分解为
$${{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{W}} = {{\boldsymbol{Q}}^{\rm{T}}}{\boldsymbol{\Sigma Q}}$$ (12) 其中, ${\boldsymbol{Q}} \in {{\bf{R}}^{r \times r}}$是一个正交矩阵, ${\boldsymbol{\Sigma }} \in {{\bf{R}}^{r \times r}}$是一个对角矩阵. 将式(12)代入式(9)并进行适当化简, 有
$${{\boldsymbol{R}}_y}{\boldsymbol{U}}_r = {{\boldsymbol{R}}_x}{\boldsymbol{U}_r{\boldsymbol{\Sigma}}}$$ (13) 其中, ${{\boldsymbol{U}}_r} = {\boldsymbol{P}}{{\boldsymbol{D}}^{ - \tfrac{1}{2}}}{{\boldsymbol{Q}}^{\rm{T}}} = {\boldsymbol{W}}{{\boldsymbol{Q}}^{\rm{T}}}$. 由于${{\boldsymbol{R}}_x}$是对称正定矩阵, 则必然可以分解为${{\boldsymbol{R}}_x} = {({{\boldsymbol{C}}^{ - 1}})^{\rm{T}}}{{\boldsymbol{C}}^{ - 1}}$, 其中${\boldsymbol{C}}$是一个可逆矩阵. 根据$\left( {{{\boldsymbol{C}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{C}}} \right){{\boldsymbol{C}}^{ - 1}}{{\boldsymbol{U}}_r} = {{\boldsymbol{C}}^{ - 1}}{{\boldsymbol{U}}_r}{\boldsymbol{\Sigma }}$和矩阵${{\boldsymbol{C}}^{ - 1}}{{\boldsymbol{U}}_r}$列满秩有: ${{\boldsymbol{C}}^{ - 1}}{{\boldsymbol{U}}_r}$必然是由${{\boldsymbol{C}}^{\rm{T}}}{{\boldsymbol{R}}_y}{\boldsymbol{C}}$的特征向量构成, 即${{\boldsymbol{U}}_r}$是由矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$最大的$r$个广义特征向量构成的矩阵.
□ 1.3 自稳定特性分析
自稳定特性是指不论初始化状态矩阵如何选择, 状态矩阵模值最终能够收敛到一个固定的常值[18]. 由于自稳定特性可以保证状态矩阵模值的稳定性, 因此已经成为算法发展中一个重要的研究内容. 接下来将对所提算法进行自稳定特性分析, 为后续方便使用, 首先定义如下矩阵范数.
定义 1. 假定${\boldsymbol{R}} \in {{\bf{R}}^{n \times n}}$是一个对称正定矩阵, 则令$\sqrt {{\rm{tr}}({{\boldsymbol{W}}^{\rm{T}}}{\boldsymbol{RW}})} $为矩阵${\boldsymbol{W}} \in {{\bf{R}}^{n \times r}}$关于矩阵${\boldsymbol{R}}$的范数, 记为${\left\| {\boldsymbol{W}} \right\|_{\boldsymbol{R}}} = \sqrt {{\rm{tr}}({{\boldsymbol{W}}^{\rm{T}}}{\boldsymbol{RW}})} $.
容易证明, 上述定义符合矩阵范数的规定. 基于定义1, 定理2给出所提算法的自稳定特性分析.
定理 2. 当输入信号是有界的且学习因子$\eta $足够小时, 所提算法(6)中状态矩阵${\boldsymbol{W}}$模值是自稳定的, 且状态矩阵模值${\left\| {\boldsymbol{W}} \right\|_{{{\boldsymbol{R}}_x}}}$的收敛值等于${\rm{tr}}({\boldsymbol{D}})$.
证明. 假定当前时刻状态矩阵模值为${\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}} = \sqrt {{\rm{tr}}\left[ {{{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_x}{\boldsymbol{W}}(k)} \right]}$, 则下一时刻的状态矩阵模值为
$$\begin{split} & \left\| {{\boldsymbol{W}}(k + 1)} \right\|_{{{\boldsymbol{R}}_x}}^2 = {\rm{tr}}\left[ {{{\boldsymbol{W}}^{\rm{T}}}(k + 1){{\boldsymbol{R}}_x}{\boldsymbol{W}}(k + 1)} \right]= \\ & \;\;\; {\rm{tr}}\{ {[ {{\boldsymbol{W}}(k) + \eta {{\boldsymbol{R}}_x}^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{W}}(k){\boldsymbol{D}}\; - } } \\ &\;\;\;{ {\eta {\boldsymbol{W}}(k){{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{\boldsymbol{W}}(k)} ]^{\rm{T}}}{{\boldsymbol{R}}_x} [ {\boldsymbol{W}}(k)\; + \\ &\;\;\; \eta {{\boldsymbol{R}}_x}^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{W}}(k){\boldsymbol{D}} - { {\eta {\boldsymbol{W}}(k){{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{\boldsymbol{W}}(k)} ]} \} =\\ & \;\;\; {\rm{tr}}\{ {{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_x}{\boldsymbol{W}}(k) -2\eta {{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{\boldsymbol{W}}(k)\;\times\\ &\;\;\;{\boldsymbol{W}}(k){{\boldsymbol{R}}_x}{\boldsymbol{W}}(k) + 2\eta {\boldsymbol{D}}{{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{\boldsymbol{W}}(k)\;+ \\ &\;\;\; {\eta ^2}{\boldsymbol{D}}{{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{{\boldsymbol{R}}_x}^{ - 1}{{\boldsymbol{R}}_y}{\boldsymbol{W}}(k){\boldsymbol{D}}+ {\eta ^2}{{\boldsymbol{W}}^{\rm{T}}}(k)\;\times \\ & \;\;\; { {{\boldsymbol{R}}_y}{\boldsymbol{W}}(k){\boldsymbol{W}}(k){{\boldsymbol{R}}_x}{\boldsymbol{W}}(k){{\boldsymbol{W}}^{\rm{T}}}(k){{\boldsymbol{R}}_y}{\boldsymbol{W}}(k)} \} \approx \\ & \;\;\;\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 - 2\eta \left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_y}}^2\left( {\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 - {\rm{tr}}({\boldsymbol{D}})} \right) \end{split} $$ (14) 由于$\eta $足够小, 因此省去式(14)中$\eta $的二阶项. 对比相邻时刻状态矩阵模值, 有
$$\begin{split} &\frac{{\left\| {{\boldsymbol{W}}(k + 1)} \right\|_{{{\boldsymbol{R}}_x}}^2}}{{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2}} = \\ &\frac{{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 - 2\eta \left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_y}}^2\left( {\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 - {\mathop{\rm tr}\nolimits} ({\boldsymbol{D}})} \right)}}{{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2}} = \\ &1 - 2\eta \frac{{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_y}}^2}}{{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2}}\left( {\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 - {\mathop{\rm tr}\nolimits} ({\boldsymbol{D}})} \right)\\[-20pt] \end{split} $$ (15) 由式(15)可得
$$\frac{{\left\| {{\boldsymbol{W}}(k + 1)} \right\|_{{{\boldsymbol{R}}_x}}^2}}{{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2}}\left\{ \begin{aligned} & > 1,&{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 < {\mathop{\rm tr}\nolimits} ({\boldsymbol{D}})} \\ & = 1,&{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 = {\mathop{\rm tr}\nolimits} ({\boldsymbol{D}})} \\ & <1,&{\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 > {\mathop{\rm tr}\nolimits} ({\boldsymbol{D}})} \end{aligned} \right. $$ 从式(15)可知, 如果当前时刻$\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 <{\rm{ tr}}({\boldsymbol{D}})$, 则在下一时刻状态矩阵模值将增加; 反之, 如果$\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 > {\rm{tr}}({\boldsymbol{D}})$, 则下一时刻状态矩阵模值减小, 即状态矩阵模值最终的收敛结果为$\left\| {{\boldsymbol{W}}(k)} \right\|_{{{\boldsymbol{R}}_x}}^2 = {\rm{tr}}({\boldsymbol{D}})$.
□ 2. 仿真实验和实例应用
本节将通过仿真实验和实例应用来验证所提算法的性能. 实验1考察算法多维广义特征向量的估计能力; 实验2考察算法的自稳定性; 实验3考察算法在信号处理中的实用性. 在仿真实验开始前, 首先随机生成两个对称正定矩阵, 该矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$最大的3个广义特征值对应的广义特征向量分别如式(16) ~ (18)所示(见本页下方).
2.1 算法能力验证
实验1分别利用所提算法和GDM算法[10]对矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$最大的3个广义特征值对应的广义特征向量进行估计, 即计算${{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2},{{\boldsymbol{v}}_3}$. 在算法的迭代过程中计算两个指标, 第1个指标是状态矩阵列向量与${{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2},{{\boldsymbol{v}}_3}$之间的方向余弦(Directional cosine, DC), 即
$$D{C_i}(k) = \frac{{\left| {{\boldsymbol{v}}_i^{\rm{T}}{{\boldsymbol{W}}_i}(k)} \right|}}{{\left\| {{{\boldsymbol{v}}_i}} \right\|\left\| {{{\boldsymbol{W}}_i}(k)} \right\|}}\tag{19}$$ (19) 其中, $ i = 1,2,3$, ${{\boldsymbol{W}}_i}(k)$表示第$k$次迭代时状态矩阵${\boldsymbol{W}}(k)$的第$i$列向量.
第2个指标是状态矩阵列向量之间关于矩阵${{\boldsymbol{R}}_x}$的范数, 即
$$Nor{m_{ij}}(k) = {\boldsymbol{W}}_i^{\rm{T}}(k){{\boldsymbol{R}}_x}{{\boldsymbol{W}}_j}(k) ,\;\; i,j = 1,2,3\tag{20}$$ (20) 两个算法的初始化状态矩阵是随机产生的, 学习因子$\eta = 0.25$. GDM算法中, 参数$\tau = 3$; 所提算法中, 加权矩阵取${\boldsymbol{D}} = {\rm{diag}}\{3,2,1\}$, 满足${d_1} > {d_2} > {d_3} > 0$的要求. 在仿真过程中, 为了更好地衡量算法性能, 这里进行了100次独立实验, 然后分别计算每次实验过程中数据的均值和上下界. 表1是两种算法完成3个广义特征向量所需的平均时间. 图1和图2分别是所提算法和GDM算法的方向余弦曲线, 其中实线是方向余弦的平均值, 虚线是数据上界, 点划线是数据下界. 为了清晰起见, 图3和图4分别只给出了所提算法的列向量模值曲线(平均值)和不同列向量之间内积关系曲线(平均值).
$${{\boldsymbol{R}}_y} = \left[ {\begin{array}{*{20}{r}} {0.7904 }&{-0.0920 }&{-0.0370 }& {0.0697 }&{-0.0563 }&{-0.0082 } \\ {-0.0920 }& {0.7087 }&{-0.0087 }&{-0.0555 }&{-0.1367 }& {0.2456 } \\ {-0.0370 }&{-0.0087 }& {0.5733 }&{-0.1282 }&{-0.2532 }&{-0.1875 } \\ {0.0697 }&{-0.0555 }&{-0.1282 }& {0.5665 }&{-0.2483 }& {0.0344 } \\ {-0.0563 }&{-0.1367 }&{-0.2532 }&{-0.2483 }& {0.5787 }& {0.0201 } \\ {-0.0082 }& {0.2456 }&{-0.1875 }& {0.0344 }& {0.0201 }& {0.2733} \end{array}} \right]\tag{16}$$ (16) $${{\boldsymbol{R}}_x} = \left[ {\begin{array}{*{20}{r}} {0.5539 }&{-0.0273 }& {0.0212 }& {0.1184 }&{-0.1045 }& {0.1350 } \\ {-0.0273 }& {0.5815 }& {0.0254 }&{-0.0481 }&{-0.1415 }&{-0.0006 } \\ {0.0212 }& {0.0254 }& {0.7190 }&{-0.0575 }&{-0.0469 }& {0.1114 } \\ {0.1184 }&{-0.0481 }&{-0.0575 }& {0.5436 }& {0.0689 }& {0.0281 } \\ {-0.1045 }&{-0.1415 }&{-0.0469 }& {0.0689 }& {0.4302 }& {0.1308 } \\ {0.1350 }&{-0.0006 }& {0.1114 }& {0.0281 }& {0.1308 }& {0.5010} \end{array}} \right]\tag{17}$$ (17) $$\left\{ \begin{array}{l} {{\boldsymbol{v}}_1} = {\left[ { - 1.0839,\; - 0.1148,\; - 0.0482,0.8737, - 1.4241,0.7618} \right]^{\rm{T}}} \\ {{\boldsymbol{v}}_2} = {\left[ {0.4091, - 0.8696,0.5280,0.2559, - 0.4955, - 0.6944} \right]^{\rm{T}}} \\ {{\boldsymbol{v}}_3} = {\left[ { - 0.9839, - 0.4039,0.3737, - 0.3484,0.2866, - 0.2464} \right]^{\rm{T}}} \\ \end{array} \right.\;\quad\qquad\tag{18}$$ (18) 表 1 两种算法的计算时间Table 1 The time cost of the two algorithms算法 时间 (ms) 所提算法 2.16 GDM算法 14.61 从图1中可以看出, 所提算法在经历了最快10次, 最慢140次, 平均90次左右的迭代运算后, 3条方向余弦曲线均能收敛到单位1. 表明不论收敛快慢, 神经网络的3个列向量均能准确地收敛到3个广义特征向量${{\boldsymbol{v}}_1},{{\boldsymbol{v}}_2},{{\boldsymbol{v}}_3}$的方向, 且其对应的广义特征值也是按照降序排列的. 从图3中可以看出, 在方向余弦收敛时, 状态矩阵的各个列向量模值也都收敛到一个常值, 即算法具有很好的收敛性. 综合图1和图3可以得知, 所提算法能够有效地估计矩阵束的广义特征向量.
从图2中可以看出, GDM算法是串行算法, 所有广义特征向量是一个接一个顺序估计的. 由于下次估计过程只有在上一次广义特征向量估计完成后才可以开始, 因此后两次迭代开始时刻分别是200和400, 而并非0. 显然这样会产生较长的等待时间, 且当估计广义特征向量维数增加时, 该等待时间会更长.
表1中, 本文所提算法的平均运行时间是2.16 ms, 远小于GDM算法的14.61 ms. 通过对比图1和图2以及综合表1可得, 由于所提算法仅需要一次迭代运算, 进而避免了等待时间, 因此计算时间较短. 需要注意的是, 由于串行算法需要信号采样完成后才可以进行运算, 而并行算法可以边采样边计算, 因此, 实际运行时并行算法在实时性方面会更有优势.
此外, 通过研究图3和图4中的收敛结果可以发现: 各个列向量模值的收敛值分别为$3,2,1$, 刚好等于加权矩阵对角线上的元素; 而两个不同的列向量关于${{\boldsymbol{R}}_x}$的内积最终收敛到了零, 即状态矩阵中不同列向量是关于矩阵${{\boldsymbol{R}}_x}$正交的. 综合图2和图3可知, ${{\boldsymbol{W}}^{\rm{T}}}{{\boldsymbol{R}}_x}{\boldsymbol{W}} = {\boldsymbol{D}} = {\rm{diag}}\{ 3,2, 1\}$, 显然该结果与定理1中的理论分析结果是一致的.
2.2 自稳定特性实验
实验2仍利用所提算法对矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$的广义特征向量进行计算. 为了验证算法的自稳定特性, 这里进行四次独立实验. 第1次取加权矩阵${{\boldsymbol{D}}_1} ={\rm{ diag}}\{10,7,5\}$, 则${\rm{tr}}({{\boldsymbol{D}}_1}) = \sqrt {22} $; 第2次取加权矩阵${{\boldsymbol{D}}_2} ={\rm{ diag}}\{3,2,1\}$, 则${\rm{tr}}({{\boldsymbol{D}}_2}) = \sqrt{ 6}$; 第3次取加权矩阵${{\boldsymbol{D}}_3} = {\rm{diag}}\{1,0.8,0.2\}$, 则${\rm{tr}}({{\boldsymbol{D}}_3}) = \sqrt 2 $; 第4次取加权矩阵${{\boldsymbol{D}}_4} ={\rm{ diag}}\{0.5,0.3, 0.2\}$, 则${\rm{tr}}({{\boldsymbol{D}}_4}) = 1$. 算法的初始化状态矩阵是随机产生的, 其模值${\left\| {{\boldsymbol{W}}(0)} \right\|_{{{\boldsymbol{R}}_x}}}$分别设置为大于、等于和小于${\rm{tr}}({{\boldsymbol{D}}_i}) (i = 1,2,3,4)$等三种情况. 四次实验中学习因子均设置为$\eta = 0.001$. 图5是四次实验中的状态矩阵模值曲线.
以图5(a)为例, 如果初始化状态矩阵${\left\| {{\boldsymbol{W}}(0)} \right\|_{{{\boldsymbol{R}}_x}}} < \sqrt {22}$时, 状态矩阵模值会不断增加, 最终收敛到$\sqrt {22} $; 反之, 当${\left\| {{\boldsymbol{W}}(0)} \right\|_{{{\boldsymbol{R}}_x}}} > \sqrt {22} $时, 模值就会减小, 最终收敛值也是$\sqrt {22} $. 从图5(b) ~ 5(d)中可以看到相类似的状态矩阵模值变化曲线, 不同之处在于: 状态矩阵模值收敛结果分别为对应的加权矩阵的迹. 所有曲线的变化规律是与定理2的分析结果相吻合的, 这也进一步证明了所提算法具有自稳定特性.
2.3 盲分离实验
实验3将所提算法应用于解决盲分离问题, 其中四个源信号取自ICALAB工具箱中的ABio7.mat文件[19]. 图6是四个源信号的波形曲线.
源信号经过混合矩阵${\boldsymbol{M}}$, 并添加加性白噪声以后得到观测信号${\boldsymbol{x}}$. 盲分离问题就是找出合适的分离矩阵${\boldsymbol{W}}$, 使得分离信号${\boldsymbol{s}} = {{\boldsymbol{W}}^{\rm{T}}}{\boldsymbol{x}}$彼此相互独立. 文献[20]提出了一个基于广义特征值分解的盲分离算法, 而分离矩阵${\boldsymbol{W}}$是由矩阵束$({{\boldsymbol{R}}_y},{{\boldsymbol{R}}_x})$所有广义特征向量构成的, 其中${{\boldsymbol{R}}_x} = {\rm{E}}[{\boldsymbol{x}}{{\boldsymbol{x}}^{\rm{T}}}]$和${{\boldsymbol{R}}_y} = {\rm{E}}[{\boldsymbol{y}}{{\boldsymbol{y}}^{\rm{T}}}]$, ${\boldsymbol{y}}$是一个滤波器的输出. 接下来利用所提算法求解该分离矩阵, 算法的初始化参数设置与实验1相同, 加权矩阵取${\boldsymbol{D}} = {\rm{diag}}\{4,3,2,1\}$, 混合矩阵为
$${\boldsymbol{M}} = \left[ {\begin{array}{*{20}{r}} {0.547}&{0.276}&{ - 1.071}&{0.615} \\ {0.911}&{ - 0.327}&{0.512}&{ - 1.158} \\ { - 0.273}&{ - 0.418}&{ - 0.571}&{ - 0.595} \\ {0.757}&{ - 1.623}&{ - 2.361}&{0.408} \end{array}} \right]$$ (21) 图7和图8分别是混合后的信号和所提算法获得的分离信号. 从图7中可以看出, 经过混合矩阵后, 源信号与观测信号具有非常大的差别. 而对比图8和图6可以发现, 分离信号与源信号具有相似的波形. 通过计算, 四个分离信号与源信号之间的相关系数分别为$0.9999,\;1.0000, \;0.9989,\;0.9665$, 即分离信号与源信号相似程度非常高, 表明所提算法能够有效地解决盲分离问题.
3. 结束语
多维广义特征值分解是应用范围最广的一类广义特征值分解算法. 针对串行算法的缺点, 本文利用加权矩阵法提出了一个并行广义特征值分解算法, 并分析了算法的收敛性和自稳定特性. 与以往算法相比, 所提算法可以在一次迭代过程中完成多维广义特征向量的同时估计, 因此具有很好的实时性和准确性. 仿真实验和应用实验表明所提算法能够有效地估计出所需的广义特征向量, 而且能够应用于解决盲分离问题. 后续研究将会考虑如何能够进一步降低算法的复杂度和加权矩阵法推广应用的问题, 这将有利于多维分解算法的发展.
-
表 1 两种算法的计算时间
Table 1 The time cost of the two algorithms
算法 时间 (ms) 所提算法 2.16 GDM算法 14.61 -
[1] Kong X Y, Du B Y, Feng X W, Luo J Y. Unified and self-stabilized parallel algorithm for multiple generalized eigenpairs extraction. IEEE Transactions on Signal Processing, 2020, 68: 3644-3659. doi: 10.1109/TSP.2020.2997803 [2] Rippl M, Lang B, Huckle T. Parallel eigenvalue computation for banded generalized eigenvalue problems. Parallel Computing, 2019, 88: 102542. doi: 10.1016/j.parco.2019.07.002 [3] Sun S L, Xie X J, Dong C. Multiview learning with generalized eigenvalue proximal support vector machines. IEEE Transactions on Cybernetics, 2019, 49(2): 688-697. doi: 10.1109/TCYB.2017.2786719 [4] Miyata T. A Riccati-type algorithm for solving generalized Hermitian eigenvalue problems. The Journal of Supercomputing, 2021, 77(2): 2091-2102. doi: 10.1007/s11227-020-03331-w [5] 郭亚宁, 林伟, 潘泉, 赵春晖, 胡劲文, 马娟娟. 基于推广流形学习的高分辨遥感影像目标分类. 自动化学报, 2019, 45(4): 720-729 doi: 10.16383/j.aas.2017.c170318Guo Ya-Ning, Lin Wei, Pan Quan, Zhao Chun-Hui, Hu Jin-Wen, Ma Juan-Juan. Generalized manifold learning for high resolution remote sensing image object classification. Acta Automatica Sinica, 2019, 45(4): 720-729 doi: 10.16383/j.aas.2017.c170318 [6] 陈晓云, 廖梦真. 基于稀疏和近邻保持的极限学习机降维. 自动化学报, 2019, 45(2): 325-333 doi: 10.16383/j.aas.2018.c170216Chen Xiao-Yun, Liao Meng-Zhen. Dimensionality reduction with extreme learning machine based on sparsity and neighborhood preserving. Acta Automatica Sinica, 2019, 45(2): 325-333 doi: 10.16383/j.aas.2018.c170216 [7] 孔祥玉, 冯晓伟, 胡昌华. 广义主成分分析算法及应用. 北京: 国防工业出版社, 2018.Kong Xiang-Yu, Feng Xiao-Wei, Hu Chang-Hua. General Principal Component Analysis and Application. Beijing: National Defense Industry Press, 2018. [8] Gao Y B, Kong X Y, Zhang Z X, Hou L A. An adaptive self-stabilizing algorithm for minor generalized eigenvector extraction and its convergence analysis. IEEE Transactions on Neural Networks and Learning Systems, 2018, 29(10): 4869-4881 doi: 10.1109/TNNLS.2017.2783360 [9] Nguyen T D, Takahashi N, Yamada I. An adaptive extraction of generalized eigensubspace by using exact nested orthogonal complement structure. Multidimensional Systems and Signal Processing, 2013, 24(3): 457-483 doi: 10.1007/s11045-012-0172-9 [10] Li H Z, Du B Y, Kong X Y, Gao Y B, Hu C H, Bian X H. A generalized minor component extraction algorithm and its analysis. IEEE Access, 2018, 6: 36771-36779 doi: 10.1109/ACCESS.2018.2852060 [11] Kong X Y, Hu C H, Duan Z S. Principal Component Analysis Networks and Algorithms. Singapore: Springer, 2017. [12] Liu L J, Shao H M, Nan D. Recurrent neural network model for computing largest and smallest generalized eigenvalue. Neurocomputing, 2008, 71(16-18): 3589-3594 doi: 10.1016/j.neucom.2008.05.005 [13] Attallah S, Abed-Meraim K. A fast adaptive algorithm for the generalized symmetric eigenvalue problem. IEEE Signal Processing Letters, 2008, 15: 797-800 doi: 10.1109/LSP.2008.2006346 [14] Nguyen T D, Yamada I. Adaptive normalized quasi-newton algorithms for extraction of generalized Eigen-pairs and their convergence analysis. IEEE Transactions on Signal Processing, 2013, 61(6): 1404-1418 doi: 10.1109/TSP.2012.2234744 [15] Qiu J L, Wang H, Lu J B, Zhang B B, Du K L. Neural network implementations for PCA and its extensions. International Scholarly Research Notices, 2012, 2012: 847305. [16] Lewis D W. Matrix Theory. Singapore: World Scientific, 1991. [17] 杜柏阳, 孔祥玉, 冯晓伟. 次成分提取信息准则的加权规则方向收敛分析. 通信学报, 2020, 41(3): 25-32 doi: 10.11959/j.issn.1000-436x.2020014Du Bo-Yang, Kong Xiang-Yu, Feng Xiao-Wei. Direction convergence analysis of weighted rule for minor component extraction information criteria. Journal on Communications, 2020, 41(3): 25-32 doi: 10.11959/j.issn.1000-436x.2020014 [18] Möller R. Derivation of Coupled PCA and SVD Learning Rules From a Newton Zero-finding Framework, Computer Engineering, Faculty of Technology, Bielefeld University, Berlin, 2017. [19] Gao Y B, Kong X Y, Hu C H, Li H Z, Hou L A. A generalized information criterion for generalized minor component extraction. IEEE Transactions on Signal Processing, 2017, 65(4): 947-959 doi: 10.1109/TSP.2016.2631444 [20] Zhang W T, Lou S T, Feng D Z. Adaptive quasi-newton algorithm for source extraction via CCA approach. IEEE Transactions on Neural Networks and Learning Systems, 2014, 25(4): 677-689 doi: 10.1109/TNNLS.2013.2280285 -