UA MATH567 高維統(tǒng)計(jì)專題1 稀疏信號(hào)及其恢復(fù)4 Basis Pursuit的算法 Projected Gradient Descent

前三講完成了對(duì)sparse signal recovery的modelling（即$L_0$-minimization建模，但考慮到它很難用于實(shí)際計(jì)算，再用$L_1$-minimization作為$L_0$-minimization的convex relaxation，并且討論了二者full recovery的性質(zhì)），這一講介紹能實(shí)際用于求解$L_1$-minimization
$$\begin{gathered} \min {\Vert x\Vert }_1 \\ s.t.y=Ax \end{gathered}$$也就是basis pursuit問(wèn)題的算法。

---

首先考慮凸優(yōu)化中最常用的Gradient descent：要求解${\min }_{x}f(x)$，只需要按下面的規(guī)則update即可
$$x_{k+1}=x_k?t_k?f(x_k)$$

其中$t_k\downarrow 0$，表示每一次update的步長(zhǎng)；但對(duì)于basis pursuit問(wèn)題，它的形式要更復(fù)雜一點(diǎn)：

有等式約束$y=Ax$

目標(biāo)函數(shù)${\Vert x\Vert }_1$在角點(diǎn)（也就是$?j,x_j=0$的地方）處不可微

于是，我們先處理第一個(gè)難點(diǎn)，等式約束$y=Ax$，定義一個(gè)線性流形
$$L=\{x\in {\mathbb{R}}^n:y=Ax\}$$

然后為了保證每一次update之后的$x_{k+1}$滿足等式約束，我們可以取$x_k?t_k?f(x_k)$在$L$上的投影作為$x_{k+1}$；

然后處理第二個(gè)難點(diǎn)，在角點(diǎn)處不可微，我們可以用次梯度(subgradient)替代梯度。次梯度的定義是

$f:S\to \mathbb{R}$ is concave where $S$ is a nonempty convex set. Then $\xi$ is called subgradient of $f$ at $\bar{x}$ if
$$f(x)\le f(\bar{x})+{\xi }^T(x?\bar{x}),?x\in S$$

$f:S\to \mathbb{R}$ is convex where $S$ is a nonempty convex set. Then $\xi$ is called subgradient of $f$ at $\bar{x}$ if
$$f(x)\ge f(\bar{x})+{\xi }^T(x?\bar{x}),?x\in S$$

不清楚次梯度的計(jì)算和意義的同學(xué)可以看一下這個(gè)例題。

與梯度相比，雖然次梯度無(wú)法提供一個(gè)最速的“下降”方向（只能提供一種正確的下降方向），但它對(duì)目標(biāo)函數(shù)的可微性沒(méi)有要求，所以在basis pursuit中用次梯度代替梯度是非常合理的。

根據(jù)以上兩條推理，我們可以修正Gradient decent算法，使之適用于basis pursuit，修正后的這個(gè)算法一般被稱為projected gradient decent。

---

下面是這個(gè)算法的偽代碼以及用R寫的一個(gè)代碼示例：

在上面的偽代碼中，有一些值得注意的細(xì)節(jié)。$\tilde{x}=A^{?}y$中，$A$本身是不可逆的，但是我們希望$\tilde{x}\in L$，所以取了$A$的pseudo-inverse $A^{?}$，可以驗(yàn)證
$$A\tilde{x}=A{A}^{?}y=A{A}^{?}(A{A}^{?}{)}^{?1}y=y$$

也就是$\tilde{x}\in L$，其中$A^{?}$表示$A$的共軛轉(zhuǎn)置；$\Gamma$矩陣是一個(gè)投影矩陣，作用是把一個(gè)向量投影到$A$的核空間；在while循環(huán)中，
$$x_t=\tilde{x}+\Gamma (x_{t?1}?\frac{sign(x_{t?1})}{t})$$

這里其實(shí)是把可行解分為了兩部分。因?yàn)閷?duì)于$y=Ax$的解集$L=\{x:y=Ax\}$而言，它可以表示為特解$\tilde{x}$和$A$的核空間$Null(A)=\{x:0=Ax\}$構(gòu)成的線性流形
$$L=\tilde{x}+Null(A)$$

因此固定了一個(gè)特解$\tilde{x}$之后，要得到使得目標(biāo)函數(shù)${\Vert x\Vert }_1$最小的解，我們只需要在$Null(A)$中進(jìn)行下降即可，這就是$\Gamma (x_{t?1}?\frac{sign(x_{t?1})}{t})$這部分的作用，其中$\Gamma$的作用是將$(x_{t?1}?\frac{sign(x_{t?1})}{t})$投影到$Null(A)$中，如果不做投影，那么第$t$次update就是$x_{t?1}?\frac{sign(x_{t?1})}{t}$，其中$1/t$是步長(zhǎng)，$sign(x_{t?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) }

需要注意的是這個(gè)R代碼示例只是簡(jiǎn)單翻譯了偽代碼，還有很多可以優(yōu)化的地方，比如矩陣求逆的計(jì)算、初始值的選取、while循環(huán)的stopping criteria等；另外，想更詳細(xì)的了解這個(gè)算法的同學(xué)可以去看Wright and Ma (2020)那本高維數(shù)據(jù)分析的section 2.3.3；另外，在這本書上的算法都可以找到對(duì)應(yīng)的matlab code，想學(xué)習(xí)一下怎么把一個(gè)用偽代碼描述的算法變成一個(gè)完整的優(yōu)化過(guò)的算法的同學(xué)可以自己去找一下。

總結(jié)

以上是生活随笔為你收集整理的UA MATH567 高维统计专题1 稀疏信号及其恢复4 Basis Pursuit的算法 Projected Gradient Descent的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。