论文将AST(原子范数软阈值化)推广到非连续采样数据(不完整数据)的情况。提出了一种无网格的SPICE(Gridless SPICE,简称GLS),它既适用于完整数据,也适用于未知噪声水平的不完整数据。
NLS(nonlinear least squares) 非线性最小二乘
IAA(iterative adaptive approach) 迭代自适应
SORTE(second order statistic of eigenvalues) 特征值二阶统计量
SSR(sparse signal representation) 稀疏信号表示
CS(compressed sensing) 压缩感知
SAGE(generalized expectation-maximization) 广义期望最大化
ADMM(alternating direction method ofmultipliers) 交换方向乘子法
WCF(weighted covariance fitting) 加权协方差拟合
AST(atomic norm soft thresholding ) 原子范数软阈值
AST(atomic norm soft thresholding)是一种基于原子范数的信号处理方法,主要用于稀疏信号的恢复和去噪。它是通过对信号的原子范数进行软阈值处理来实现的。
软阈值处理是一种经典的信号处理技术,它可以抑制信号中的噪声并保留信号中的重要信息。在AST中,软阈值处理被应用于信号的原子范数,以实现信号的稀疏性。
Soft ( x , λ ) \operatorname{Soft}(\boldsymbol{x},\lambda) Soft(x,λ)表示对向量 x \boldsymbol{x} x进行软阈值处理, λ \lambda λ是一个阈值参数。
软阈值处理的具体形式为: Soft ( x , λ ) = s i g n ( x ) max { ∣ x ∣ − λ , 0 } \operatorname{Soft}(\boldsymbol{x},\lambda)=sign(\boldsymbol x)\max\{|\boldsymbol x|-\lambda,0\} Soft(x,λ)=sign(x)max{∣x∣−λ,0}
在AST中,软阈值处理被应用于信号的原子范数,即: x ^ = Soft ( A T ( y − A x ) , λ ∣ ∣ Φ ∣ ∣ ∗ ) \hat{\boldsymbol x}= \operatorname{Soft}(\boldsymbol{A^T(\boldsymbol y-Ax)},\lambda|| \boldsymbol{\Phi}||_*) x^=Soft(AT(y−Ax),λ∣∣Φ∣∣∗).
其中, ∣ ∣ Φ ∣ ∣ ∗ ||\boldsymbol{\Phi}||_{*} ∣∣Φ∣∣∗是矩阵 Φ \boldsymbol{\Phi} Φ的核范数, λ \lambda λ是软阈值参数。AST方法通过使用原子范数和软阈值处理相结合的方式,实现了稀疏信号的高效恢复。
SPICE(sparse iterative covariance-based estimation) 稀疏迭代协方差估计
SPICE(Sparse Iterative Covariance-based Estimation)是一种用于稀疏信号处理和估计的算法,主要用于估计协方差矩阵。
\space \space 基本思想是利用协方差矩阵的稀疏性来降低算法的复杂度,并通过迭代方法逐步逼近协方差矩阵的真实值。在迭代过程中,SPICE算法利用了一种基于协方差矩阵估计的方法,即通过最小化原始数据和协方差矩阵之间的误差来估计协方差矩阵。在估计协方差矩阵的同时,SPICE算法还实现了稀疏性,即只有少数非零元素的矩阵表示。
GLS:gridless version of SPICE 无网格SPICE
假设我们观察到一个有噪声的正弦信号(索引为
j
j
j)
y
i
=
∑
k
=
1
K
s
k
e
i
2
π
(
j
−
1
)
f
k
+
e
j
y_i=\displaystyle\sum_{k=1}^K s_ke^{i2 \pi (j-1)f_k} +e_j
yi=k=1∑Kskei2π(j−1)fk+ej
在索引集
[
M
]
=
△
{
1
,
2
,
.
.
.
,
M
}
[M]\xlongequal{\triangle}\{1,2,...,M\}
[M]△{1,2,...,M} 或者子集
Ω
⊂
[
M
]
\mathbf{\varOmega} \subset[M]
Ω⊂[M] 上,
Y
∈
C
M
\mathcal{Y}\in\cnums^M
Y∈CM。
f
k
f_ k
fk 和
s
k
s_k
sk 分别表示第
k
k
k 个正弦信号的归一化频率,(复)振幅。
e
j
∈
C
e_j \in \cnums
ej∈C 是测量干扰。
正弦曲线数
K
<
M
K<M
K<M 通常被称为模型阶数,在实践中通常是未知的。
当信号能在
[
M
]
[M]
[M] 上被观测到时的情况被称为完整数据情况,而当只有样本在
Ω
⊂
[
M
]
\mathbf{\varOmega} \subset[M]
Ω⊂[M] 上可用时的其它情况被称为不完整数据情况(或丢失数据情况),其中在
Ω
\mathbf{\varOmega}
Ω 的互补集合上
Ω
‾
=
[
M
]
/
Ω
\overline{\mathbf{\varOmega}} =[M]/\mathbf{\varOmega}
Ω=[M]/Ω 的样本被称为丢失数据。
缺失数据的情况是重要的,因为缺失样本在实践中是常见的,其可以由传感器故障、异常值、天气条件或其他物理约束引起。
频率估计和模型阶数选择是线谱估计中的两个重要课题。
压缩感知(Compressed Sensing )
采样网格不应该太密集,否则相邻原子(或导向矢量)之间几乎完全相关可能会降低稀疏恢复性能。
但是,密集网格可以使频率估计更精确,因为密集网格可以减少网格失配(在网格点和真实频率之间)和近似误差(观测模型的。
因此,现有的稀疏方法是否可以用无限密集的网格来实际实现,或者等效地,直接在连续区间 [0,1) 上而不进行网格化?
基的失配: 信号的实际参数值很可能不会恰好落在离散的网格上,在这种情况下,假设的变换基则无法表示稀疏信号。
ℓ 1 \ell1 ℓ1范数去噪
在离散字典
A
∈
C
M
×
N
\mathcal{A} \in \cnums^{M\times N}
A∈CM×N 下从有噪声的测量信号
Y
∈
C
M
\mathcal{Y}\in\cnums^M
Y∈CM 中恢复信号
z
∈
C
M
\mathcal{z}\in\cnums^M
z∈CM 的问题,即存在稀疏向量
s
∈
C
M
\mathcal{s}\in\cnums^M
s∈CM 使得
z
=
A
s
\mathcal{z}=\mathcal{A}\mathcal{s}
z=As。通过以下优化问题来求解:
m
i
n
s
μ
∥
s
∥
1
+
g
(
Y
−
A
s
)
\underset{\mathcal{s}}{min}\mu\lVert\mathcal{s}\rVert_1+g(\mathcal{Y}-\mathcal{A}\mathcal{s})
sminμ∥s∥1+g(Y−As)
其中,
g
(
⋅
)
g(\centerdot)
g(⋅) 表示拟合数据,以及正则项参数
μ
>
0
\mu>0
μ>0 平衡存在噪声的测量值
Y
\mathcal{Y}
Y 的保真度和稀疏性
s
\mathcal s
s。上述问题被称为
ℓ
1
\ell 1
ℓ1 范数去噪(L1ND),也叫做 Lasso 问题。此外,还有平方根-(SR-)Lasso和最小绝对偏差-(LAD-)Lasso问题。
从统计学角度看,Lasso算法适合于高斯噪声,而LAD-Lasso算法对离群点具有较强的鲁棒性。与Lasso相比,SR-Lasso需要对噪声分布进行宽松的假设,并可轻松选择 μ \mu μ。
Lasso(Least Absolute Shrinkage and Selection Operator)是一种用于回归分析的统计方法。它通过对目标变量和自变量之间的线性关系进行建模来预测目标变量的值。与普通最小二乘回归不同,Lasso回归通过在代价函数中添加一个正则化项来限制模型中的参数数量,从而避免过度拟合。
Lasso问题通常是指使用Lasso回归求解线性回归问题时所遇到的最小化目标函数的优化问题。其目标函数包括两个部分:第一个部分是普通最小二乘回归的平方误差损失函数,第二个部分是L1正则化项,这个项对应了Lasso回归中的缩减惩罚项。L1正则化项是指将每个参数的绝对值之和乘以一个常数λ,λ是一个可以通过交叉验证等技术来选择的超参数。通过最小化这个目标函数,Lasso回归可以得到一个稀疏的参数向量,即一些参数会被缩减为零,从而达到特征选择的目的。
平方根LASSO(Square-root Lasso)是一种对普通Lasso方法的改进,它使用平方根惩罚项而不是L1惩罚项。与Lasso的惩罚项不同,平方根Lasso的惩罚项是将每个参数的平方根之和乘以一个常数λ,λ是一个可以通过交叉验证等技术来选择的超参数。
与Lasso相比,平方根Lasso的主要优点是能够更好地保留小的参数值,并且在一些数据集中具有更好的预测性能。另外,由于平方根Lasso使用平方根惩罚项,因此它的目标函数在零处的导数为零,这意味着平方根Lasso可以更容易地处理共线性问题,避免了在Lasso中可能出现的随机性。
需要注意的是,平方根Lasso与Lasso一样,也是一种线性回归方法,因此可以使用类似于最小二乘回归的方法来求解。同时,与Lasso相同,平方根Lasso也可以用于特征选择和稀疏建模,从而提高模型的解释性和泛化性能。
最小绝对偏差LASSO(Least Absolute Deviations Lasso,简称LAD Lasso)是一种同时使用L1和L2惩罚项的线性回归方法。它的目标是最小化残差的绝对值之和,同时通过L1和L2惩罚项对模型参数进行约束。
LAD Lasso的惩罚项包括一个L1惩罚项和一个L2惩罚项,其中L1惩罚项用于稀疏建模和特征选择,而L2惩罚项用于控制参数的平滑性和泛化性能。与普通的Lasso不同,LAD Lasso的目标函数不是平方误差损失函数,而是残差的绝对值之和。 与其他线性回归方法相比,LAD Lasso的主要优点是可以在存在异常值或噪声的情况下提高模型的鲁棒性。由于LAD Lasso使用绝对值损失函数,它对于极端值的敏感度比平方损失函数较小。
此外,LAD Lasso还可以用于稀疏建模和特征选择,从而提高模型的解释性和泛化性能。 需要注意的是,LAD Lasso的求解问题不是凸优化问题,因此需要使用迭代算法来求解。常用的求解算法包括坐标下降算法、梯度下降算法、LARS算法等。
原子范数
原子范数的定义如下:
对于一个矩阵或向量
x
x
x,设 A 是一个包含向量或矩阵的原子集合,原子范数定义为:
∣
∣
x
∣
∣
A
=
inf
{
t
>
0
:
x
∈
t
c
o
n
v
(
A
)
}
||x||_A = \inf \{t > 0: x ∈ tconv(A)\}
∣∣x∣∣A=inf{t>0:x∈tconv(A)},其中,
t
c
o
n
v
(
A
)
tconv(A)
tconv(A) 是原子集合A的凸包,
t
c
o
n
v
(
A
)
tconv(A)
tconv(A) 是原子集合 A 的图像,即线性组合的集合。
简单来说,原子范数是所有包含
x
x
x 的原子集合 A 的最小凸包,也就是可以最小化
x
x
x 的线性组合。是一种用于描述矩阵或向量中的稀疏性的范数。与L1范数类似,原子范数也可以用于特征选择和稀疏建模,但它对于某些特定的稀疏结构更具有优势。
设
A
\mathcal{A}
A 为原子的集合,满足其凸包
c
o
n
v
(
A
)
conv(\mathcal A)
conv(A) 是紧致的,中心对称的,并且包含原点作为内点。然后,
c
o
n
v
(
A
)
conv(\mathcal A)
conv(A) 的规范函数定义了一个范数,称为原子范数,并表示为:
∣
∣
Y
∣
∣
A
=
{
t
>
0
:
y
∈
t
c
o
n
v
(
A
)
}
=
inf
{
∑
k
c
k
α
k
,
c
k
≥
0
,
α
∈
A
}
.
||\mathcal Y||_A= \{ t>0:y\in tconv(\mathcal A)\}=\inf\set{ \sum_{\substack{k}}c_k\alpha_k,c_k\geq0,\alpha\in\mathcal A}.
∣∣Y∣∣A={t>0:y∈tconv(A)}=inf{k∑ckαk,ck≥0,α∈A}.
原子范数的对偶范数为:
∣
∣
z
∣
∣
A
∗
=
sup
{
⟨
z
,
α
⟩
R
:
∣
∣
α
∣
∣
A
≤
1
}
||\mathcal z ||_{\mathcal A}^*=\sup \{\lang \mathcal z,\mathcal \alpha \rang_{\R}:||\alpha||_{\mathcal A} \le 1 \}
∣∣z∣∣A∗=sup{⟨z,α⟩R:∣∣α∣∣A≤1}
⟨
z
,
α
⟩
R
=
ℜ
⟨
z
,
α
⟩
=
ℜ
{
α
H
z
}
\lang\mathcal z,\mathcal \alpha \rang_{\R}=\Re\lang \mathcal z,\mathcal \alpha \rang=\Re\{ \alpha^H \mathcal z\}
⟨z,α⟩R=ℜ⟨z,α⟩=ℜ{αHz},
ℜ
\Re
ℜ表示取复数的实部。此外,
c
o
n
v
(
A
)
=
{
α
:
∣
∣
α
∣
∣
A
≤
1
}
conv(\mathcal A)=\{ \alpha :||\alpha||_{\mathcal A}\le1\}
conv(A)={α:∣∣α∣∣A≤1},并且
A
\mathcal A
A 包含
{
α
:
∣
∣
α
∣
∣
A
≤
1
}
\{\alpha :||\alpha||_{\mathcal A}\le1\}
{α:∣∣α∣∣A≤1} 的所有极值点,所以:
∣
∣
z
∣
∣
A
∗
=
sup
α
∈
A
⟨
z
,
α
⟩
R
||z||_{\mathcal A}^*=\sup_{\substack{\mathcal \alpha \in \mathcal A}}\lang\mathcal z,\mathcal \alpha \rang_{\R}
∣∣z∣∣A∗=α∈Asup⟨z,α⟩R
完整数据线谱估计的原子范数软阈值AST (atomic-norm soft thresholding) :
观测信号模型可以写成如下形式:
Y
=
∑
k
=
1
K
α
(
f
k
)
s
k
+
e
=
A
(
f
)
+
e
\mathcal Y=\displaystyle\sum_{k=1}^K \mathcal \alpha(f_k)\mathcal s_k+\mathcal e=\mathcal A(f)+\mathcal e
Y=k=1∑Kα(fk)sk+e=A(f)+e其中
α
(
f
k
)
=
[
1
,
e
i
2
π
f
k
,
.
.
.
,
e
i
2
π
(
M
−
1
)
f
k
]
∈
C
M
,
A
(
f
)
=
[
α
(
f
1
)
,
α
(
f
2
)
,
.
.
.
,
α
(
f
k
)
]
∈
C
M
×
K
,
Y
∈
C
M
\alpha(f_k)=[1,e^{i2\pi f_k},...,e^{i2\pi(M-1) f_k}]\in\cnums^M,\mathcal A(f)=[\alpha(f_1),\alpha(f_2),...,\alpha(f_k)]\in\cnums^{M\times K},\mathcal Y \in \cnums^M
α(fk)=[1,ei2πfk,...,ei2π(M−1)fk]∈CM,A(f)=[α(f1),α(f2),...,α(fk)]∈CM×K,Y∈CM表示所以
y
i
y_i
yi的向量集。
α
(
f
,
ϕ
)
=
α
(
f
)
ϕ
\alpha(f,\phi)=\alpha(f)\phi
α(f,ϕ)=α(f)ϕ,其中
ϕ
∈
S
1
=
△
{
ϕ
∈
C
:
∣
ϕ
∣
=
1
}
\phi \in \mathbb{S}^1\xlongequal{\triangle}\{\phi\in\mathbb{C}:|\phi|=1\}
ϕ∈S1△{ϕ∈C:∣ϕ∣=1}。因此原子集
A
\mathcal{A}
A 可以表示为:
A
=
△
{
α
(
f
,
ϕ
)
:
f
∈
[
0
,
1
)
,
ϕ
∈
S
1
}
\mathcal A\xlongequal{\triangle}\{\alpha(f,\phi):f\in[0,1),\phi\in\mathbb{S}^1\}
A△{α(f,ϕ):f∈[0,1),ϕ∈S1}这个原子范数可以通过半定规划(SDP)来计算:
∣
∣
Y
∣
∣
A
=
m
i
n
x
,
u
1
2
(
x
+
u
1
)
,
s
u
b
j
e
c
t
t
o
[
x
Y
H
Y
T
(
u
)
]
≥
0.
||\mathcal Y ||_\mathcal A=\underset{\mathcal{x,\mathcal u}}{min}\frac 1 2(x+u_1),subject\space to\space \begin{bmatrix} x & \mathcal Y ^H \\ \mathcal Y & T(u) \end{bmatrix}\geq 0.
∣∣Y∣∣A=x,umin21(x+u1),subject to [xYYHT(u)]≥0.其中
u
∈
C
M
u \in \cnums^M
u∈CM ,
T
(
u
)
∈
C
M
×
M
T(u)\in\cnums^{M\times M}
T(u)∈CM×M 表示一个 (Hermitian) Toeplitz 矩阵,
T
(
U
)
=
[
u
1
u
2
.
.
.
u
M
u
2
H
u
1
.
.
.
u
M
−
1
⋮
⋮
⋱
⋮
u
M
H
u
M
−
1
H
.
.
.
u
1
]
T(U)=\begin{bmatrix}u_1& u_2&...&u_M\\u_2^H & u_1&...&u_{M-1}\\ \vdots &\vdots & \ddots &\vdots \\u_M^H &u_{M-1}^H&...&u_1 \end{bmatrix}
T(U)=
u1u2H⋮uMHu2u1⋮uM−1H......⋱...uMuM−1⋮u1
估计无噪声正弦信号的原子软阈值(AST),令
z
=
△
A
(
f
)
s
z\xlongequal{\triangle}\mathcal A(f)s
z△A(f)s :
m
i
n
z
=
μ
∣
∣
z
∣
∣
A
+
1
2
∣
∣
Y
−
z
∣
∣
2
2
\underset{\mathcal{\mathcal z}}{min}=\mu||z||_\mathcal A+\frac 1 2||\mathcal Y-z||_2^2
zmin=μ∣∣z∣∣A+21∣∣Y−z∣∣22其中当
M
M
M 足够大时
μ
≈
M
ln
N
σ
0
1
2
\mu \approx\sqrt{\smash[b]{M\ln N }}\sigma_0^{\frac 1 2}
μ≈MlnNσ021,上面的问题可以转化为半定规划问题:
m
i
n
x
,
u
,
z
u
2
(
x
+
u
1
)
+
1
2
∣
∣
Y
−
z
∣
∣
2
2
,
s
u
b
j
e
c
t
t
o
[
x
z
H
z
T
(
u
)
]
≥
0.
\underset{\mathcal{x,u,z}}{min} \frac u2(x+u_1)+\frac 1 2||\mathcal Y-z||_2^2,subject \space to \begin{bmatrix} x & \mathcal z ^H \\ z & T(u) \end{bmatrix}\geq 0.
x,u,zmin2u(x+u1)+21∣∣Y−z∣∣22,subject to[xzzHT(u)]≥0.求得最优解
(
x
∗
,
u
∗
,
z
∗
)
(x^*,u^*,z^*)
(x∗,u∗,z∗)后,频率
f
^
\hat {f}
f^ 和幅度
s
^
\hat s
s^ 的估计值可以通过
T
(
u
∗
)
T(u*)
T(u∗) 的Vandermonde分解获得。
不完整数据的AST
不完整数据的原子范数
假设观察样本在子集
Ω
∈
[
M
]
\varOmega \in[M]
Ω∈[M] 上,假定
Ω
\varOmega
Ω 是按升序排序的。样本量
L
=
∣
Ω
∣
≤
M
L=|\varOmega| \le M
L=∣Ω∣≤M,采样周期的范围为:
M
‾
=
Ω
L
−
Ω
1
+
1
≤
M
\overline M=\varOmega_L-\varOmega_1+1\le M
M=ΩL−Ω1+1≤M。
缺失数据的原子集定义为:
A
(
Ω
)
≜
{
a
Ω
:
a
∈
A
}
=
a
Ω
(
f
,
ϕ
)
:
f
∈
[
0
,
1
)
,
ϕ
∈
S
1
.
\mathcal A(\varOmega)\triangleq{ \{ \boldsymbol a_{\varOmega}:\boldsymbol a\in\mathcal A \}}=\boldsymbol a_{\varOmega}(f,\phi):f\in[0,1),\phi\in\mathbb S^1.
A(Ω)≜{aΩ:a∈A}=aΩ(f,ϕ):f∈[0,1),ϕ∈S1.这里的
a
Ω
(
f
,
ϕ
)
=
a
Ω
(
f
)
ϕ
\boldsymbol a_{\varOmega}(f,\phi)=\boldsymbol a_{\varOmega}(f)\phi
aΩ(f,ϕ)=aΩ(f)ϕ和
a
Ω
(
f
)
\boldsymbol a_{\varOmega}(f)
aΩ(f)都是
a
(
f
)
\boldsymbol a(f)
a(f)的子向量。
lemma1:对于原子范数
∣
∣
⋅
∣
∣
A
(
Ω
)
||\centerdot||_{\mathcal A(\varOmega)}
∣∣⋅∣∣A(Ω),
∣
∣
Y
Ω
∣
∣
A
(
Ω
)
=
m
i
n
Y
Ω
‾
∣
∣
Y
∣
∣
A
=
m
i
n
x
,
u
,
Y
Ω
‾
1
2
∣
∣
x
+
u
1
∣
∣
,
s
u
b
j
e
c
t
t
o
[
x
Y
H
Y
T
(
u
)
]
≥
0.
||\mathcal Y_{\varOmega}||_{\mathcal A(\varOmega)}=\underset{\mathcal{\mathcal {\mathcal Y_{\overline\varOmega}}}}{min}||\mathcal Y||_\mathcal A=\underset{\mathcal{\mathcal x,u,{\mathcal Y_{\overline\varOmega}}}}{min}\frac 1 2||\mathcal x+u_1||,subject \space to \begin{bmatrix} x & {\mathcal Y} ^H \\ \mathcal Y& T(u) \end{bmatrix}\geq 0.
∣∣YΩ∣∣A(Ω)=YΩmin∣∣Y∣∣A=x,u,YΩmin21∣∣x+u1∣∣,subject to[xYYHT(u)]≥0.
对偶原子范数:
∣
∣
z
Ω
∣
∣
A
(
Ω
)
∗
=
sup
f
,
ϕ
∈
S
1
⟨
z
Ω
,
α
Ω
(
f
)
⟩
R
.
||z_\varOmega||_{\mathcal A(\varOmega)}^*=\sup_{\substack{f,\mathcal \phi\in \mathbb S^1}}\lang\mathcal z_\varOmega,\mathcal \alpha_\varOmega(f) \rang_{\R}.
∣∣zΩ∣∣A(Ω)∗=f,ϕ∈S1sup⟨zΩ,αΩ(f)⟩R.
不完整数据的AST
假设抽样数据被独立同分布的噪声污染,通过解决下面的AST问题来估计无干扰信号: min z Ω μ ∣ ∣ z Ω ∣ ∣ A ( Ω ) + 1 2 ∣ ∣ Y Ω − z Ω ∣ ∣ 2 2 , \min\limits_{z_\varOmega}\mu||z_\varOmega||_{\mathcal A(\varOmega)}+\frac 1 2||\mathcal Y_{\varOmega}-z_\varOmega||_2^2, zΩminμ∣∣zΩ∣∣A(Ω)+21∣∣YΩ−zΩ∣∣22,其中 μ > 0 \mu >0 μ>0为一个定值,可以将上式写为如下的半定规划问题: min x , u , z μ 2 ( x + u 1 ) + 1 2 ∣ ∣ Y Ω − z Ω ∣ ∣ 2 2 , s u b j e c t t o [ x z H z T ( u ) ] ≥ 0. \min\limits_{x,u,z}\frac \mu 2(x+u_1)+\frac 1 2||\mathcal Y_{\varOmega}-z_\varOmega||_2^2,subject\space to\begin{bmatrix} x & {\mathcal z} ^H \\ \mathcal z& T(u) \end{bmatrix}\geq 0. x,u,zmin2μ(x+u1)+21∣∣YΩ−zΩ∣∣22,subject to[xzzHT(u)]≥0.将无噪声信号表示为: z o = A ( f ) s z^o=\mathcal A(f)s zo=A(f)s。在 μ ≥ E ∣ ∣ e Ω ∣ ∣ A ( Ω ) ∗ \mu \geq E||e_\varOmega||_{\mathcal A(\varOmega)}^* μ≥E∣∣eΩ∣∣A(Ω)∗ 时, z o z^o zo 的估计值 z ^ \hat z z^ 具有期望(每个元素)均方误差(MSE) : 1 L E ∣ ∣ z ^ Ω − z Ω o ∣ ∣ 2 2 ≤ μ L ∑ k = 1 K ∣ s k ∣ . \frac 1 LE||\hat z_\varOmega-z_\varOmega^o||_2^2\le \frac \mu L \sum_{k=1}^K |s_k|. L1E∣∣z^Ω−zΩo∣∣22≤Lμk=1∑K∣sk∣.
无网络的SPICE(GLS)
A.SPICE
SPICE 是一种基于加权协方差拟合(WCF)的基于网格的线谱估计稀疏方法。在LASSO使用中,它自动估计通常预先不可用的噪声方差。假设在 Y = ∑ k = 1 K α ( f k ) s k + e = A ( f ) + e \mathcal Y=\displaystyle\sum_{k=1}^K \mathcal \alpha(f_k)\mathcal s_k+\mathcal e=\mathcal A(f)+\mathcal e Y=k=1∑Kα(fk)sk+e=A(f)+e 中的 s k s_k sk 的相位相互独立且均匀分布。因此, E [ s s H ] = d i a g ( ∣ s k ∣ 2 ) = △ d i a g ( p ) E[ss^H]=diag(|s_k|^2)\xlongequal {\triangle} diag(p) E[ssH]=diag(∣sk∣2)△diag(p),其中 p p p 叫做功率参数。进一步假设噪声 e e e 独立于信号 s k s_k sk 是独立的,并且满足 E [ e e H ] = d i a g ( σ ) E[ee^H]=diag(\sigma) E[eeH]=diag(σ),其中的 σ \sigma σ 表示噪声方差参数。然后 Y \mathcal Y Y 的协方差矩阵可以表示为: R = E [ Y Y H ] = A ( f ) d i a g ( p ) A H ( f ) + d i a g ( σ ) R=E[\mathcal Y \mathcal Y^H]=\mathcal A(f) diag(p)\mathcal A^H(f)+diag(\sigma) R=E[YYH]=A(f)diag(p)AH(f)+diag(σ)SPICE使用最小化WCF标准,如下所示: h ( f , p , σ ) = ∣ ∣ R − 1 2 ( Y Y H − R ) ∣ ∣ F 2 = t r ( R ) + ∣ ∣ Y ∣ ∣ 2 2 Y H R − 1 Y − 2 ∣ ∣ Y ∣ ∣ 2 2 . h(f,p,\sigma)=||R^{- \frac 1 2}(\mathcal Y \mathcal Y^H-R)||_F^2=tr(R)+||\mathcal Y||_2^2\mathcal Y ^HR^{-1}\mathcal Y-2||\mathcal Y||_2^2. h(f,p,σ)=∣∣R−21(YYH−R)∣∣F2=tr(R)+∣∣Y∣∣22YHR−1Y−2∣∣Y∣∣22.这个问题产生的优化问题 min f , p ⪰ 0 , σ ⪰ 0 t r ( R ) + ∣ ∣ Y ∣ ∣ 2 2 Y H R − 1 Y \min \limits_{f,p\succeq 0,\sigma \succeq 0} tr(R)+||\mathcal Y||_2^2\mathcal Y ^HR^{-1}\mathcal Y f,p⪰0,σ⪰0mintr(R)+∣∣Y∣∣22YHR−1Y是非凸问题,因为上式中关于 f f f 的数据协方差 R R R 是非线性的。与其他现有的稀疏方法一样,将离散化应用于连续频域,以消除 R R R 对 f f f 的依赖性。 可以通过 SPICE 来求解。
B.完整数据情况下的GLS
无网格版本的SPICE,即GLS。
GLS的一个关键技术是通过引入一个半正定Toeplitz矩阵
T
(
u
)
=
A
(
f
)
d
i
a
g
(
p
)
A
H
(
f
)
+
d
i
a
g
(
σ
)
T(u)=\mathcal A(f) diag(p)\mathcal A^H(f)+diag(\sigma)
T(u)=A(f)diag(p)AH(f)+diag(σ) 来重新参数化
R
R
R ,而不是线性化协方差矩阵
R
R
R 的离散化。因此
R
R
R 可以表示为:
R
=
T
(
u
)
+
d
i
a
g
(
σ
)
R=T(u)+diag(\sigma)
R=T(u)+diag(σ),其中
T
(
u
)
≥
0
,
σ
⪰
0
T(u)\geq 0,\sigma \succeq 0
T(u)≥0,σ⪰0。这与之前的
R
R
R 的表达式是等价的。因此,GLS 的优化问题可以等效为 SDP 问题:
min
u
,
σ
⪰
0
t
r
(
R
)
+
∣
∣
y
∣
∣
2
2
y
H
R
−
1
y
,
s
u
b
j
e
c
t
t
o
T
(
u
)
≥
0
=
min
x
,
u
,
σ
⪰
0
t
r
(
R
)
+
∣
∣
y
∣
∣
2
2
x
,
s
u
b
j
e
c
t
t
o
[
x
Y
H
Y
R
]
≥
0
T
(
u
)
≥
0.
\min\limits_{u,\sigma \succeq 0} tr(R)+||y||_2^2y^HR^{-1}y,subject\space to\space T(u)\geq 0 \\ =\min\limits_{x,u,\sigma \succeq0}tr(R)+||y||_2^2x,\\ subject \space to \begin{bmatrix} x & \mathcal Y^H \\ \mathcal Y & R \end{bmatrix} \geq 0\space T(u)\geq 0.
u,σ⪰0mintr(R)+∣∣y∣∣22yHR−1y,subject to T(u)≥0=x,u,σ⪰0mintr(R)+∣∣y∣∣22x,subject to[xYYHR]≥0 T(u)≥0.半定规划问题可以直接通过现有的求解器求解。求得解
(
u
∗
,
σ
∗
)
(u^*,\sigma^*)
(u∗,σ∗),进一步估计
(
f
^
,
p
^
,
σ
^
)
(\hat f,\hat p,\hat \sigma)
(f^,p^,σ^)。如果
T
(
u
∗
)
T(u^*)
T(u∗) 不是满秩的,那么
f
^
,
p
^
\hat f,\hat p
f^,p^ 可以通过范德蒙德分解唯一的确定,并且
σ
^
=
σ
∗
\hat \sigma=\sigma^*
σ^=σ∗;如果
T
(
u
∗
)
T(u^*)
T(u∗) 是满秩的,那么
f
^
,
p
^
\hat f,\hat p
f^,p^ 并不是唯一的。这就意味这优化问题具有多个最优解。 使
δ
=
λ
min
(
T
(
u
∗
)
)
\delta= \lambda_{\min}(T(u^*))
δ=λmin(T(u∗)) 为
T
(
u
∗
)
T(u^*)
T(u∗) 的最小特征值。然后通过 Vandermonde 分解
T
(
u
∗
)
−
δ
I
T(u^*)-\delta I
T(u∗)−δI 和
σ
^
=
σ
∗
+
δ
1
\hat \sigma = \sigma^* + \delta 1
σ^=σ∗+δ1 唯一地获得
f
^
,
p
^
\hat f , \hat p
f^,p^ 满足
∣
f
^
∣
=
∣
p
^
∣
≤
M
−
1
|\hat f|=|\hat p|\leq M-1
∣f^∣=∣p^∣≤M−1,
I
I
I 和
1
1
1 分别是单位矩阵和单位矩阵的一个向量。
在同方差噪声的假设下,数据协方差矩阵
R
R
R 本身是一个 Toeplitz 矩阵,因此可以表示为
R
=
T
(
u
˜
)
R=T (\~u)
R=T(u˜) ,
u
˜
∈
C
M
,
T
(
u
˜
)
≥
0.
\~u \in \mathbb C^M,T(\~u)\geq 0.
u˜∈CM,T(u˜)≥0.
C.缺失数据情况下的GLS
在缺失数据的情况下,仅观察上的样本
Ω
⊂
[
M
]
\varOmega \subset [M]
Ω⊂[M] 。采用相同的WCF标准,但仅应用于可用数据
y
Ω
y_\varOmega
yΩ 。在与SPICE相同的假设下,
R
Ω
R_\varOmega
RΩ 表示
y
Ω
y_\varOmega
yΩ 的协方差矩阵:
R
Ω
=
E
[
Y
Ω
Y
Ω
H
]
=
A
Ω
(
f
)
d
i
a
g
(
p
)
A
Ω
H
(
f
)
+
d
i
a
g
(
σ
Ω
)
.
R_\varOmega = E[\mathcal Y_\varOmega \mathcal Y_\varOmega^H] = A_\varOmega(f)diag(p)A_\varOmega^H(f)+diag(\sigma_\varOmega).
RΩ=E[YΩYΩH]=AΩ(f)diag(p)AΩH(f)+diag(σΩ).其中,
A
Ω
(
f
)
=
Γ
Ω
A
(
f
)
,
Γ
Ω
∈
{
0
,
1
}
L
×
M
A_\varOmega(f)=\Gamma_\varOmega A(f),\Gamma_\varOmega \in \{ 0,1 \}^{L \times M}
AΩ(f)=ΓΩA(f),ΓΩ∈{0,1}L×M,它的元素只有在
(
j
,
Ω
j
)
,
j
∈
[
L
]
(j,\varOmega_j),j\in [L]
(j,Ωj),j∈[L] 时为 1 。因此,
R
Ω
R_\varOmega
RΩ 在约束:
T
(
u
)
≥
0
a
n
d
σ
Ω
⪰
0
T(u)\geq 0 \space and \space \sigma _\varOmega \succeq 0
T(u)≥0 and σΩ⪰0 可以等效地重新参数化为:
R
Ω
=
Γ
Ω
T
(
u
)
Γ
Ω
T
+
d
i
a
g
(
σ
Ω
)
.
R_\varOmega=\Gamma_\varOmega T(u) \Gamma_\varOmega^T+diag(\sigma_\varOmega).
RΩ=ΓΩT(u)ΓΩT+diag(σΩ).其中
T
(
u
)
T(u)
T(u) 可以解释为“干净”完整数据的协方差。
因此,GLS解决了以下凸优化问题:
min
u
,
σ
Ω
⪰
0
t
r
(
R
Ω
)
+
∣
∣
Y
Ω
∣
∣
2
2
Y
Ω
H
R
Ω
−
1
Y
Ω
,
s
u
b
j
e
c
t
t
o
T
(
u
)
≥
0.
\min\limits_{u,\sigma_\varOmega \succeq0}tr(R_\varOmega)+||\mathcal Y_\varOmega||_2^2\mathcal Y_\varOmega^H R_\varOmega^{-1}\mathcal Y_\varOmega,subject \space to\space T(u)\geq 0.
u,σΩ⪰0mintr(RΩ)+∣∣YΩ∣∣22YΩHRΩ−1YΩ,subject to T(u)≥0.得到解
(
u
∗
,
σ
Ω
∗
)
(u^*,\sigma_\varOmega^*)
(u∗,σΩ∗) 后,可以通过范德蒙德分解
T
(
u
∗
)
T(u^*)
T(u∗) 或者
T
(
u
∗
)
−
λ
min
(
T
(
u
∗
)
)
I
.
T(u^*)-\lambda_{\min}(T(u^*)) I.
T(u∗)−λmin(T(u∗))I. 可以得到估计参数
(
f
^
,
p
^
,
σ
^
Ω
)
(\hat f,\hat p,\hat \sigma_\varOmega)
(f^,p^,σ^Ω)。
GLS是SPICE的无网格版本,并且可以应用于存在同方差或异方差噪声的完整数据和缺失数据的情况,其中噪声方差是未知的。
GLS的扩展:线谱估计的系统框架
A.GLS的两个局限性:模型阶次不准确和频率分裂
GLS 是一种无需模型阶数或噪声方差的先验知识即可估计线性模型参数的强大方法。
该方法仍然存在一些限制:
1.一种限制与高斯噪声假设有关。如果系统中的噪声是非高斯分布的,那么 GLS 可能不是参数估计的最佳方法。
2.在系统具有大量参数的情况下,GLS 可能表现不佳,因为估计协方差矩阵的计算成本会很高。
GLS通常产生比真实模型阶数大得多的长度的频率估计。数据协方差 R ≥ 0 R\geq 0 R≥0 可以表示为Toeplitz矩阵本身。因为 GLS 的 SDP 问题 Y \mathcal Y Y 位于 R R R 的值域中。则 R R R 的解 R ∗ R^* R∗ 满秩的概率为 1。因为具有随机噪声的 Y \mathcal Y Y 可以被分解为 M-1 个正弦波的叠加。因此,GLS 的估计正弦分量的数量等于 R ∗ − λ min ( R ∗ ) I R^*-\lambda_{\min}(R^*)I R∗−λmin(R∗)I,几乎肯定是 M-1。(最小特征值是多重特征值的概率为零)
GLS会产生频率分裂问题。与基于网格的方法中遇到的不同(一个正弦分量被分成支持在附近网格点上的几个)。是由网格和某些算法的收敛问题引起的。相比之下,GLS 的频率分裂部分地由于不存在噪声水平而引起(AST很少遭受频率分裂)
无网格方法不需要将频率离散化成网格上的点,而是直接在连续的频率域上进行估计,因此不会受到网格大小的限制和频率分裂的影响;相比之下,网格方法需要将频率离散化成网格上的点进行估计,容易受到频率分裂的影响。
实际上在一定条件下,Gridless方法和基于网格的方法可能会遇到相似的频率分裂现象,但两者的表现和原因可能有所不同。在一些情况下,Gridless方法可能比基于网格的方法更容易受到频率分裂的影响,因为它们更依赖于数据的分布和噪声模型。因此,需要根据具体问题来选择合适的方法,并注意它们可能存在的限制和局限性。
B.论文提出线谱估计框架
因为GLS的两个局限性:
1)粗糙的网格导致频率分辨率较差,只能检测到由网格引起的频率分裂;
2)高度密集的网格意味着相邻原子之间几乎完全相关,可能会产生数值问题,给检测带来困难。
为了克服GLS的上述两个局限性,提出了线谱估计的框架,该框架由三个步骤组成:
a. 基于GLS进行协方差估计 ;
b. 基于协方差估计进行模型阶数选择;
c. 基于协方差和模型阶估计进行频率估计 。
论文框架创新点
1.在存在缺失数据的情况下,如何估计数据的协方差通常是不清楚的,而GLS通过协方差拟合和利用其 Toeplitz 结构提供了一种解决方案 。
2.通过适当地选择时间窗,可以在完整数据情况下给出数据协方差估计,然而,时间窗缩短了数据长度并潜在地降低了分辨率极限。在论文中,选择 Sorte 算法用于第二步的模型阶数选择,选择MUSIC算法用于第三步的频率估计。 这些方法的选择会产生最好的性能。
在此框架内可以修正频率分裂,正确估计模型阶数
GLS与原子范数去噪的关系
将
l
1
l1
l1 范数去噪的优化问题中的
l
1
l1
l1 范数用
∣
∣
⋅
∣
∣
A
||\cdot ||_{\mathcal A}
∣∣⋅∣∣A 或
∣
∣
⋅
∣
∣
A
(
Ω
)
||\cdot||_{\mathcal A(\varOmega)}
∣∣⋅∣∣A(Ω) 代替,将会得到优化问题:
min
z
μ
∣
∣
z
∣
∣
A
+
g
(
Y
−
z
)
\min\limits_{z}\mu ||z||_\mathcal A+g(\mathcal Y-z)
zminμ∣∣z∣∣A+g(Y−z)称为原子范数去噪(AND)。
GLS优化问题相当于下列AND问题之一:
1)在异方差噪声假设下,
min
z
Ω
L
∣
∣
z
Ω
∣
∣
A
(
Ω
)
+
∣
∣
Y
Ω
−
z
Ω
∣
∣
1
;
\min\limits_{z_\varOmega}\sqrt{\smash[b]{L}}||z_\varOmega||_{\mathcal A(\varOmega)}+||\mathcal Y_\varOmega - z_\varOmega||_1;
zΩminL∣∣zΩ∣∣A(Ω)+∣∣YΩ−zΩ∣∣1;
2)在同方差噪声下,
min
z
Ω
∣
∣
z
Ω
∣
∣
A
(
Ω
)
+
∣
∣
Y
Ω
−
z
Ω
∣
∣
2
\min\limits_{z_\varOmega}||z_\varOmega||_{\mathcal A(\varOmega)}+||\mathcal Y_\varOmega - z_\varOmega||_2
zΩmin∣∣zΩ∣∣A(Ω)+∣∣YΩ−zΩ∣∣2
3)在方差已知为
σ
0
\sigma_0
σ0 的同方差噪声下 ,
min
z
Ω
L
∣
∣
Y
Ω
∣
∣
2
σ
0
∣
∣
z
Ω
∣
∣
A
(
Ω
)
+
1
2
∣
∣
Y
Ω
−
z
Ω
∣
∣
2
2
\min\limits_{z_\varOmega}\frac {\sqrt{\smash[b]{L}}}{||\mathcal Y_\varOmega||_2}\sigma_0 ||z_\varOmega||_{\mathcal A(\varOmega)}+\frac 1 2||\mathcal Y_\varOmega - z_\varOmega||_2^2
zΩmin∣∣YΩ∣∣2Lσ0∣∣zΩ∣∣A(Ω)+21∣∣YΩ−zΩ∣∣22此外,在不同假设条件下的GLS优化问题在给定上述相应 AND 问题的解
u
∗
u^*
u∗ 的情况下都有最优解
∣
∣
Y
Ω
∣
∣
2
L
\frac {||\mathcal Y_\varOmega||_2}{\sqrt{\smash[b]{L}}}
L∣∣YΩ∣∣2。这些问题都是凸的,可以转化为SDP问题。
在不同的噪声假设下,GLS可以解释为三种不同的方法。 在异方差噪声下,GLS采用GL-LAD-LASSO的形式,通过在数据拟合中使用范数来抑制重要的噪声项,并对离群点具有鲁棒性;在同方差噪声下,GLS采用GL-SR-LASSO的形式,将所有噪声项综合考虑,只反映噪声能量; 最后,在方差固定的同方差噪声下,GLS表现为GL-Lasso或AST形式,噪声方差反映在正则化参数中。 此外,值得注意的是,由于SPICE和L1ND都是基于网格的GLS和AND,它们之间存在类似的等价性。
S P I C E N ∗ SPICE_N^* SPICEN∗ 和 G L S ∗ GLS^* GLS∗ 分别表示 SPICE(具有均匀大小为 N 的网格)和 GLS 的最优目标函数值。得到: ( 1 − π M ‾ N ) S P I C E N ∗ ≤ G L S ∗ ≤ S P I C E N ∗ . (1-\frac{\pi \overline M}{N})SPICE_N^* \leq GLS^* \leq SPICE_N^*. (1−NπM)SPICEN∗≤GLS∗≤SPICEN∗.因此,SPICE和L1ND是GLS和AND的基于网格的版本,它们之间存在类似的等价性。 lim N → + ∞ S P I C E N ∗ = G L S ∗ . \lim\limits_{N\to +\infty} SPICE_N^*=GLS^*. N→+∞limSPICEN∗=GLS∗.
计算可行解
A.通过对偶性求解
GL-LAD-Lasso 可以写成下面的SDP问题: min x , u , z τ ( x + u 1 ) + ∣ ∣ Y Ω − z Ω ∣ ∣ 1 , s u b j e c t t o [ x z H z T ( u ) ] ≥ 0. \min\limits_{x,u,z}\tau (x+u_1)+|| \mathcal Y_\varOmega-z_\varOmega ||_1,subject \space to \space \begin{bmatrix} x & z^H \\ z & T(u) \end{bmatrix} \geq 0. x,u,zminτ(x+u1)+∣∣YΩ−zΩ∣∣1,subject to [xzzHT(u)]≥0.其中 τ = L 2 \tau=\frac {\sqrt L} {2} τ=2L。它的对偶问题由下面的SDP按照标准的拉格朗日分析给出: min u , W 2 ℜ { Y Ω H v Ω } , s u b j e c t t o { [ τ v H v W ] ≥ 0 , ∣ ∣ v Ω ∣ ∣ ∞ ≤ 1 2 , T ∗ ( W ) = τ e 1 , \min\limits_{u,W}2\Re \{ \mathcal Y_\varOmega^H v_\varOmega \},subject \space to \space \begin{cases} {\begin{bmatrix} \tau & v^H \\ v & W \end{bmatrix}}\geq0,\\ {|| v_\varOmega ||_{\infty}\leq {\frac 1 2}}, \\{T^*(W)=\tau e_1,} \end{cases} u,Wmin2ℜ{YΩHvΩ},subject to ⎩ ⎨ ⎧[τvvHW]≥0,∣∣vΩ∣∣∞≤21,T∗(W)=τe1,其中 e 1 = [ 1 , 0 , . . . , 0 ] T ∈ R M . e_1=[1,0,...,0]^T\in \R^M. e1=[1,0,...,0]T∈RM.
B.通过ADMM求解
基于 ADMM 提出了本文涉及的 SDP 的一阶算法,这是解决大规模问题的成熟方法。GL-LAD-LASSO的算法可以写为: min x , u , z , Q > 0 τ ( x + u 1 ) + ∣ ∣ Y Ω − z Ω ∣ ∣ 1 , s u b j e c t t o Q = [ x z H z T ( u ) ] . \min\limits_{x,u,z,Q>0}\tau(x+u_1)+||\mathcal Y_\varOmega-z_\varOmega||_1,\\ subject \space to \space Q=\begin{bmatrix} x & z^H \\ z & T(u) \end{bmatrix}. x,u,z,Q>0minτ(x+u1)+∣∣YΩ−zΩ∣∣1,subject to Q=[xzzHT(u)].引入拉格朗日乘子 Λ \varLambda Λ: L ( x , y , z , Q , Λ ) = τ ( x − u 1 ) + ∣ ∣ Y Ω − z Ω ∣ ∣ 1 + t r [ ( Q − [ x z H z T ( u ) ] ) Λ ] + β 2 ∣ ∣ Q − [ x z H z T ( u ) ] ∣ ∣ F 2 \mathcal L(x,y,z,Q,\varLambda)=\tau(x-u_1)+||\mathcal Y_\varOmega-z_\varOmega||_1+tr[(Q-\begin{bmatrix} x & z^H \\ z & T(u) \end{bmatrix}) \varLambda ]+\frac \beta 2 ||Q-\begin{bmatrix} x & z^H \\ z & T(u) \end{bmatrix}||_F^2 L(x,y,z,Q,Λ)=τ(x−u1)+∣∣YΩ−zΩ∣∣1+tr[(Q−[xzzHT(u)])Λ]+2β∣∣Q−[xzzHT(u)]∣∣F2 β > 0 \beta>0 β>0 是惩罚参数集 ,遵循ADMM的规则,变量 ( x , u , z , Q ) (x,u,z,Q) (x,u,z,Q) 和 Λ \varLambda Λ可以以封闭形式迭代更新。