UA MATH567 高維統計專題1 稀疏信號及其恢復4 Basis Pursuit的算法 Projected Gradient Descent

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

---

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

其中$t_k\downarrow 0$，表示每一次update的步長；但對于basis pursuit問題，它的形式要更復雜一點：

有等式約束$y=Ax$

目標函數${\Vert x\Vert }_1$在角點（也就是$?j,x_j=0$的地方）處不可微

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

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

然后處理第二個難點，在角點處不可微，我們可以用次梯度(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$$

不清楚次梯度的計算和意義的同學可以看一下這個例題。

與梯度相比，雖然次梯度無法提供一個最速的“下降”方向（只能提供一種正確的下降方向），但它對目標函數的可微性沒有要求，所以在basis pursuit中用次梯度代替梯度是非常合理的。

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

---

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

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

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

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

因此固定了一個特解$\tilde{x}$之后，要得到使得目標函數${\Vert x\Vert }_1$最小的解，我們只需要在$Null(A)$中進行下降即可，這就是$\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$是步長，$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) }

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。