UA MATH567 高维统计专题1 稀疏信号及其恢复4 Basis Pursuit的算法 Projected Gradient Descent

前三讲完成了对sparse signal recovery的modelling（即L0L_0L0-minimization建模，但考虑到它很难用于实际计算，再用L1L_1L1-minimization作为L0L_0L0-minimization的convex relaxation，并且讨论了二者full recovery的性质），这一讲介绍能实际用于求解L1L_1L1-minimization
min⁡∥x∥1s.t.y=Ax\min \ \ \left\| x\right\|_1 \\ s.t. \ \ y = Axmin ∥x∥1s.t. y=Ax也就是basis pursuit问题的算法。

首先考虑凸优化中最常用的Gradient descent：要求解min⁡xf(x)\min_xf(x)minxf(x)，只需要按下面的规则update即可
xk+1=xk−tk∇f(xk)x_{k+1}=x_k-t_k\nabla f(x_k)xk+1=xk−tk∇f(xk)

其中tk↓0t_k \downarrow 0tk↓0，表示每一次update的步长；但对于basis pursuit问题，它的形式要更复杂一点：

有等式约束y=Axy=Axy=Ax
目标函数∥x∥1\left\| x\right\|_1∥x∥1在角点（也就是∃j,xj=0\exists j,x_j=0∃j,xj=0的地方）处不可微

于是，我们先处理第一个难点，等式约束y=Axy=Axy=Ax，定义一个线性流形
L={x∈Rn:y=Ax}L=\{x \in \mathbb{R}^n:y=Ax\}L={x∈Rn:y=Ax}

然后为了保证每一次update之后的xk+1x_{k+1}xk+1满足等式约束，我们可以取xk−tk∇f(xk)x_k-t_k\nabla f(x_k)xk−tk∇f(xk)在LLL上的投影作为xk+1x_{k+1}xk+1；

然后处理第二个难点，在角点处不可微，我们可以用次梯度(subgradient)替代梯度。次梯度的定义是

f:S→Rf:S \to \mathbb{R}f:S→R is concave where SSS is a nonempty convex set. Then ξ\xiξ is called subgradient of fff at xˉ\bar xxˉ if
f(x)≤f(xˉ)+ξT(x−xˉ),∀x∈Sf(x) \le f(\bar x)+\xi^T(x-\bar x),\forall x \in Sf(x)≤f(xˉ)+ξT(x−xˉ),∀x∈S

f:S→Rf:S \to \mathbb{R}f:S→R is convex where SSS is a nonempty convex set. Then ξ\xiξ is called subgradient of fff at xˉ\bar xxˉ if
f(x)≥f(xˉ)+ξT(x−xˉ),∀x∈Sf(x) \ge f(\bar x)+\xi^T(x-\bar x),\forall x \in Sf(x)≥f(xˉ)+ξT(x−xˉ),∀x∈S

不清楚次梯度的计算和意义的同学可以看一下这个例题。

与梯度相比，虽然次梯度无法提供一个最速的“下降”方向（只能提供一种正确的下降方向），但它对目标函数的可微性没有要求，所以在basis pursuit中用次梯度代替梯度是非常合理的。

根据以上两条推理，我们可以修正Gradient decent算法，使之适用于basis pursuit，修正后的这个算法一般被称为projected gradient decent。

下面是这个算法的伪代码以及用R写的一个代码示例：

在上面的伪代码中，有一些值得注意的细节。x~=A†y\tilde x=A^{\dag}yx~=A†y中，AAA本身是不可逆的，但是我们希望x~∈L\tilde x \in Lx~∈L，所以取了AAA的pseudo-inverse A†A^{\dag}A†，可以验证
Ax~=AA†y=AA∗(AA∗)−1y=yA \tilde x =AA^{\dag}y=AA^*(AA^*)^{-1}y=yAx~=AA†y=AA∗(AA∗)−1y=y

也就是x~∈L\tilde x \in Lx~∈L，其中A∗A^*A∗表示AAA的共轭转置；Γ\GammaΓ矩阵是一个投影矩阵，作用是把一个向量投影到AAA的核空间；在while循环中，
xt=x~+Γ(xt−1−sign(xt−1)t)x_t=\tilde x+\Gamma(x_{t-1}-\frac{sign(x_{t-1})}{t})xt=x~+Γ(xt−1−tsign(xt−1))

这里其实是把可行解分为了两部分。因为对于y=Axy=Axy=Ax的解集L={x:y=Ax}L=\{x:y=Ax\}L={x:y=Ax}而言，它可以表示为特解x~\tilde xx~和AAA的核空间Null(A)={x:0=Ax}Null(A)=\{x:0=Ax\}Null(A)={x:0=Ax}构成的线性流形
L=x~+Null(A)L=\tilde x+Null(A)L=x~+Null(A)

因此固定了一个特解x~\tilde xx~之后，要得到使得目标函数∥x∥1\left\| x \right\|_1∥x∥1最小的解，我们只需要在Null(A)Null(A)Null(A)中进行下降即可，这就是Γ(xt−1−sign(xt−1)t)\Gamma(x_{t-1}-\frac{sign(x_{t-1})}{t})Γ(xt−1−tsign(xt−1))这部分的作用，其中Γ\GammaΓ的作用是将(xt−1−sign(xt−1)t)(x_{t-1}-\frac{sign(x_{t-1})}{t})(xt−1−tsign(xt−1))投影到Null(A)Null(A)Null(A)中，如果不做投影，那么第ttt次update就是xt−1−sign(xt−1)tx_{t-1}-\frac{sign(x_{t-1})}{t}xt−1−tsign(xt−1)，其中1/t1/t1/t是步长，sign(xt−1)sign(x_{t-1})sign(xt−1)就是次梯度。

PS <- function(A,y,max_iter = 1000){m = nrow(A); n = ncol(A)A1 = t(A)%*%solve(A%*%t(A))Gamma = diag(n)-A1%*%Ax_tilde = A1%*%yx0 = rep(0,n); t = 0while(t<max_iter){t = t+1x1 = x_tilde + Gamma%*%(x0 - (1/t)*sign(x0))x0 = x1}return(x0)
}

需要注意的是这个R代码示例只是简单翻译了伪代码，还有很多可以优化的地方，比如矩阵求逆的计算、初始值的选取、while循环的stopping criteria等；另外，想更详细的了解这个算法的同学可以去看Wright and Ma (2020)那本高维数据分析的section 2.3.3；另外，在这本书上的算法都可以找到对应的matlab code，想学习一下怎么把一个用伪代码描述的算法变成一个完整的优化过的算法的同学可以自己去找一下。

UA MATH567 高维统计专题1 稀疏信号及其恢复4 Basis Pursuit的算法 Projected Gradient Descent相关推荐

UA MATH567 高维统计专题1 稀疏信号及其恢复7 LASSO的预测误差与变量选择一致性
UA MATH567 高维统计专题1 稀疏信号及其恢复7 LASSO的预测误差与变量选择一致性 Prediction Error Variable Selection Consistency Pred ...
UA MATH567 高维统计专题1 稀疏信号及其恢复6 随机设计矩阵下LASSO的估计误差
UA MATH567 高维统计专题1 稀疏信号及其恢复6 随机设计矩阵下LASSO的估计误差上一讲我们推导了noisy setting下LASSO估计误差的阶O(slog⁡d/n)O(\sqrt{s ...
UA MATH567 高维统计专题1 稀疏信号及其恢复5 LASSO的估计误差
UA MATH567 高维统计专题1 稀疏信号及其恢复5 LASSO的估计误差 Signal Recovery Noisy Setting LASSO的估计误差 Signal Recovery Noi ...
UA MATH567 高维统计专题1 稀疏信号及其恢复3 Coherence与RIP简介
UA MATH567 高维统计专题1 稀疏信号及其恢复3 Coherence与RIP简介 Pairwise inc oherence Mutual Coherence RIP 前两讲介绍了L0-min ...
UA MATH567 高维统计专题1 稀疏信号及其恢复2 用L1-norm作为L0-norm的convex relexation
UA MATH567 高维统计专题1 稀疏信号及其恢复2 用L1-norm作为L0-norm的convex relexation L1L_1L1-norm minimization L1L_1L1 ...
UA MATH567 高维统计专题1 稀疏信号及其恢复1 L0-norm minimization
UA MATH567 高维统计专题1 稀疏信号及其恢复1 L0-norm minimization L0L^0L0-norm L0L_0L0-norm minimization Exhaustive ...
UA MATH567 高维统计专题0 为什么需要高维统计理论？——理解稀疏向量与hard-threshold
UA MATH567 高维统计专题0 为什么需要高维统计理论?--理解稀疏向量与hard-threshold 稀疏向量的soft-threshold与hard-threshold近似引入hard- ...
UA MATH567 高维统计专题0 为什么需要高维统计理论？——高维统计理论的常用假设
UA MATH567 高维统计专题0 为什么需要高维统计理论?--高维统计理论的常用假设延续前三讲对线性判别分析的讨论,在高维时,根据中心极限定理 n(Xˉ−μ)→dN(0,Id)\sqrt{n} ...
UA MATH567 高维统计专题3 含L1-norm的凸优化6 Stochastic Gradient Descent简介
UA MATH567 高维统计专题3 含L1-norm的凸优化6 Stochastic Gradient Descent简介 Stochastic Gradient Descent的思想 Varian ...

UA MATH567 高维统计专题1 稀疏信号及其恢复4 Basis Pursuit的算法 Projected Gradient Descent

UA MATH567 高维统计专题1 稀疏信号及其恢复4 Basis Pursuit的算法 Projected Gradient Descent

UA MATH567 高维统计专题1 稀疏信号及其恢复4 Basis Pursuit的算法 Projected Gradient Descent相关推荐

最新文章

热门文章