一种新的相关系数矩阵稀疏估计的对偶方法

薛文娟, 刘 潇

(1.上海电力大学 数理学院, 上海 200090; 2.上海理工大学 理学院, 上海 200439)

近年来,利用半正定协方差矩阵的稀疏逼近在声纳数据处理[1]、细胞信号恢复[2]、基因聚类[3]、语音信号分类[4-5]等领域都有着广泛的应用。正如文献[6]提到的,协方差矩阵不是尺度不变的。稀疏半正定相关系数矩阵在基因聚类的实际应用中也发挥着重要作用[6-7]。因此,寻求一种新的、高效求解半正定相关系数矩阵稀疏逼近问题的算法是非常必要的。

CUI Y等人[6]首次提出采用对偶方法,通过加速近端梯度来求解稀疏逼近问题的某种对偶问题。该方法是对已被广泛研究的协方差矩阵稀疏逼近[1-5]的推广。本文针对稀疏逼近优化问题推导出一个新的对偶问题,并利用邻近梯度法求解该新对偶问题。

本文首先将引入稀疏逼近优化问题的一个新对偶,并给出一阶最优性条件;然后提出一种新的基于近邻梯度的对偶方法来求解所得到的新对偶问题,并证明其具有全局收敛性;最后通过数值对比验证了该算法的有效性。

1 对偶问题及其最优性条件

1.1 稀疏逼近优化问题

本文主要考虑如下形式的稀疏优化问题

s.t.〈Ai,X〉=bi,

〈Bj,X〉≤uj,

X-

(1)

若对于矩阵P中某些元素,Pij=0,那么对Xij不要求具有稀疏性。若要求矩阵X的对角线元素为1,即定义Ai的第s行第t列中的元素为

(2)

其中:s=1,2,3,…,n;t=1,2,3,…,n。b的第i个元素bi=1,i=1,2,3,…,n。则可以得到式(1)模型的一个特例,即具有以下形式的相关系数矩阵的稀疏逼近:

s.t.Xii=1,i∈[n],

X-

(3)

其中:[n]表示集合{1,2,3,…,n}。

(4)

s.t.A(X)=b,B(X)≤u,

X-

(5)

1.2 一个新的对偶问题

经过计算可得问题式(5)的对偶问题如下

minΨ(Λ,γ,λ)=

s.t.(Λ,λ)∈Ω

(6)

其中:Ω是可行集,即

Λ≤P,λ≥o}

(7)

∇ΛΨ(Λ,γ,λ)=(C+A*(γ)-B*(λ)+

Λ-I)++I

∇γΨ(Λ,γ,λ)=A(C+A*(γ)-B*(λ)+

Λ-I)++A(I)-b

∇λΨ(Λ,γ,λ)=B(C+A*(γ)-B*(λ)+

Λ-I)++B(I)+u

(8)

式(8)中的梯度都是Lipschitz连续的,且当Slater约束规范条件成立时,问题式(5)和对偶问题式(6)之间的对偶间隙为零。从而,只需求解对偶问题式(6)即可。

1.3 最优性条件

由于不等式约束-Pij≤Λij≤Pij,i,j∈[n],λ≥o的特殊结构,Mangasarian-Fromovitz约束规范条件在对偶问题式(6)的任意可行点都成立[12],于是可以得到对偶问题式(6)的一阶最优性条件。

(9)

(10)

(11)

(12)

此外,容易验证一阶最优性条件式(9)等价于

r(Λk,γk,λk)=0

(13)

因此,KKT误差rk=r(Λk,γk,λk)可以用来度量最优性。

2 一种新的对偶方法及其收敛性分析

本节将应用邻近梯度法求解对偶问题式(6)。

2.1 求解对偶问题的邻近梯度法

邻近梯度法[13]一般用于求解可分解的、非光滑的无约束优化问题,即

minf1(x)+f2(x)

(14)

其中:f1是凸的且可微的;f2是凸的邻近映射(但不一定是可微的)。在当前迭代点,邻近梯度法的迭代公式为

xk+1=proxαkf2(xk-αk∇f1(xk))

(15)

其中:αk>0是步长。

min〈∇ΛΨ(Λk,γk,λk),Λk+DΛ〉+

(16)

(17)

其中:δΩ(Λ,λ)是Ω的指示函数,即

(18)

注意到式(16)中目标函数的前4项是连续可微的,最后1项是凸函数,其邻近映射可由式(19)计算得到

proxδΩ(Λk+DΛ,λk+dλ)=

ΠΩ(Λk+DΛ,λk+dλ)

(19)

其中:ΠΩ(·)由式(11)计算可得。容易推出式(16)和式(17)闭合形式的解分别如下

λk-∇λΨ(Λk,γk,λk))-

(Λk,λk)

(20)

(21)

2.2 主要算法

下面将给出完整的对偶邻近梯度算法(Dual Proximal Gradient Algorithm,DPGA)。令(Λk,γk,λk)为当前迭代点,通过式(16)和式(17)(或式(20)和式(21))可得邻近梯度步。定义Ψ(Λ,γ,λ)的预计线性下降量Δk为

(22)

利用线搜索回溯技术,在{1,τ,τ2,…}中寻找最大的αk,使得Ψ(Λ,γ,λ)满足充分下降条件

Ψ(Λk,γk,λk)+σαkΔk

(23)

其中:τ∈(0,1);σ∈(0,1)是给定的常数。从而,可得DPGA算法的迭代公式为

(24)

对于算法的终止条件,将类似于文献[6],用ηgap来表示原始-对偶的间隙,可以衡量原始目标与对偶目标函数值之间的差异,即

(25)

其中:jprob和jdob分别表示原始目标函数值和对偶目标函数值。

当式(26)成立时,DPGA算法停止。

(26)

定理1由DPGA算法生成的序列{(Λk,γk,λk)}的任何聚点都是对偶问题式(6)的最优解。

证明:应用邻近梯度方法收敛性结论立即可得到所提算法的全局收敛性。

3 数值实验

为了验证DPGA算法的有效性,本文将其与文献[6]中的相关系数矩阵稀疏估计(Sparse Estimation of the Correlation matrix,SEC)算法进行了比较。为了公平比较,本文还将SEC算法代码中的EIG(·)替换为Mexeig(·)。

3.1 设置问题

本文实验中的测试问题形式如下

s.t.Xii=1,i∈[n],

X-

(27)

取ε=10-6。当DPGA算法满足终止条件式(26),或达到最大迭代次数(5 000次)时停止。DPGA算法中的其他参数设置如下:Λ0=O,γ0=o,σ=10-4,τ=0.5。此外,SEC算法采用默认参数。

3.2 合成数据集

本文生成矩阵C有以下两种方法。

(1) 测试问题E1式(27)中的矩阵C由MATLAB命令实现生成。具体为:C=rand(n,n);C=(C+C′)/2;C(1:n+1:end)=1。

(2) 测试问题E2式(27)中的矩阵C是一个样本相关系数矩阵,是由一组N(0,Σ)样本计算生成的,均值为零。其中:Σ是一个的n×n的协方差矩阵。特别地,定义

i,j∈[n]

(28)

接下来,将从N(0,Σ)中抽取p个独立同分布的样本,用D∈Rp×n来表示。最后,应用zscore方法对D进行归一化,使得每个特征的均值为零,单位方差为1。在归一化数据上,样本协方差矩阵与样本相关系数矩阵相同。在本文的数值实验中,对所有给定的矩阵C,ρ=0.01,将通过求解式(27)来估计相关系数矩阵。

本文从式(27)对偶问题的残差、CPU时间、总迭代次数、原始目标值,以及估计的相关系数矩阵的稀疏性等5个方面与算法SEC进行比较。实验结果表明,SEC算法和DPGA算法都能在预设误差小于10-6和最大迭代次数内成功地解决所有测试问题。

图1显示了SEC和DPGA算法在E1上的数值结果。其中,φ和φmin分别为目标函数值和其最小值。由图1(a)可以看出,在n∈[1 000,2 800]之间所有10个测试问题上,尽管差距很小,但DPGA算法得到的原始目标值都优于SEC算法。对于残差,图1(b)显示DPGA算法在所有情况下都比SEC算法更精确。此外,由图1(c)和图1(d)可以看出,DPGA算法在E1上的迭代次数不仅小于SEC算法的1/4,而且CPU时间也少了一半,这说明本文所提出的DPGA算法更具优越性。

图1 SEC和DPGA算法在E1上的数值结果

此外,通过设置n∈[1 000,2 000]和样本数p∈[1 000,2 000],对E2进行了实验研究。图2和表1给出了E2的相关数值结果。由图2可以看出,与SEC算法相比,DPGA算法在E2上几乎都能以更少的迭代次数获得更高精度的解。 从表1可以观察到,SEC和DPGA算法在E2上都可以达到几乎相同的原始目标值及相同的稀疏程度。

图2 SEC和DPGA算法在E2上的数值结果

表1 E2上的数值结果

综上所述,针对本文所提出的新对偶问题式(6),DPGA算法可以较低的计算量获得相关系数矩阵稀疏逼近的最优解。

4 结 语

本文提出了一种Frobenius范数下稀疏矩阵逼近问题的对偶邻近梯度算法,用于求解由原始问题导出的一种新的对偶问题。该算法利用邻近梯度产生搜索方向,并使用线搜索技术来寻找合适的步长以确保算法具有全局收敛性。此外,本文还对相关模型随机生成数据进行了初步的数值实验,并与文献[6]中的SEC算法进行了比较。通过对数值结果比较,结果显示本文提出的求解对偶问题的DPGA算法优于SEC算法。