L1正则化降噪，对偶函数的构造，求解含L1正则项的优化问题，梯度投影法

本文主要实现L1正则化降噪，L2 正则化降噪的文章在：
https://blog.csdn.net/IYXUAN/article/details/121565229

原创文章！转载需注明来源：©️ Sylvan Ding’s Blog ❤️

实验目的

1. 学会使用L1正则化项求解降噪问题；
2. 构造原问题的对偶问题，以求解含L1正则项的优化问题；
3. 使用梯度投影法，求解带约束的优化问题；
4. 在降噪过程中，注意对比L1和L2正则化的区别，并分析使用L1正则项在此降噪问题中的优势。

实验环境

MATLAB R2021a

实验内容

L2正则化

本文主要实现L1正则化降噪，L2 正则化降噪的文章在：
https://blog.csdn.net/IYXUAN/article/details/121565229

假设有一段被噪声污染的信号，即 $y=x+w,x,w,y\in {\mathbb{R}}^n$ . 其中， $x$ 是不含噪声的源信号，$w$ 是一未知噪声，$y$ 是观察到的含噪信号。

目标是从观察到的信号 $y$ 还原出源信号 $x$ . 在实验一中，我们假设最初的信号是一段光滑的曲线，选用二次惩罚（Quadratic Penalty）作为正则化项，则有

$$x^{\ast }=\arg \underset{x}{\min }\Vert x-y{\Vert }^2+\lambda \Vert Lx{\Vert }^2,\lambda >0$$

其中，$x^{\ast }$ 是所求得的最优解，$L\in {\mathbb{R}}^{(n-1)\times n}$ ，

$$L=\left[\begin{array}{ccccccc}1& -1& 0& 0& \cdots & 0& 0\\ 0& 1& -1& 0& \cdots & 0& 0\\ 0& 0& 1& -1& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& & ⋮& ⋮\\ 0& 0& 0& 0& \cdots & 1& -1\end{array}\right]$$

利用最小二乘法，可以解得 $x^{\ast }(\lambda )={\left(I+\lambda L^{T}L\right)}^{-1}y$ .

L1正则化

然而，这种使用二次惩罚作为正则项，接着使用最小二乘求解的方法具有一定的局限性，在某些情况下，无论正则化参数 $\lambda$ 如何设置，都不能得到理想的解。事实上，对于存在断点的信号（比如脉冲信号）来说，由于惩罚项 $\Vert Lx{\Vert }^2$ 是二次的，断点会导致惩罚项的值急剧增大。因此，这种二次惩罚项会让有噪信号的间断点处更加平滑，但这种间断点处的平滑并不是我们希望得到的。

为了缓解这一问题，削弱信号在跳跃处的惩罚项值增大的幅度，我们使用 L1 范数形式的惩罚项，而非二次惩罚，那么优化问题变为（P）

$$\underset{x}{\min }\Vert x-y{\Vert }^2+\lambda \Vert Lx{\Vert }_1$$ .

算法描述

构造对偶问题

优化问题（P）等价于如下问题：

$$\begin{array}{lc}\underset{x,z}{\min }& \Vert x-y{\Vert }^2+\lambda \Vert z{\Vert }_1\\ \mathrm{s}.\mathrm{t}.& z=Lx\end{array}$$

那么，拉格朗日函数 $L$ 为

$$\begin{array}{lc}L(x,z,\mu )& =\Vert x-y{\Vert }^2+\lambda \Vert z{\Vert }_1+{\mu }^T(Lx-z)\\ & =\Vert x-y{\Vert }^2+(L^T\mu {)}^{T}x+\lambda \Vert z{\Vert }_1-{\mu }^{T}z\end{array}$$

对偶目标函数（Dual objective function）为 $q(\mu )=\underset{x,z}{\min }L(x,z,\mu )$.

观察 $L$ 可以发现，可分别求 $x$ 和 $z$ 的最优解，即

$$q(\mu )=\underset{x}{\min }\left\{\Vert x-y{\Vert }^2+(L^T\mu {)}^{T}x\right\}+\underset{z}{\min }\left\{\lambda \Vert z{\Vert }_1-{\mu }^{T}z\right\}$$.

$q(\mu )$ 中关于 $x$ 部分的最小值

对于 $\Vert x-y{\Vert }^2+(L^T\mu {)}^{T}x$ 求关于 $x$ 的最小值，显然这是一个二次函数，在梯度消失时，取得最小值，即

$$2(x-y)+L^T\mu =0$$

解得 $x^{\ast }=y-\frac{1}{2}{L}^T\mu$ ，故

$$\underset{x}{\min }\Vert x-y{\Vert }^2+(L^T\mu {)}^{T}x=-\frac{1}{4}{\mu }^{T}L{L}^T\mu +{\mu }^{T}Ly.$$

$q(\mu )$ 中关于 $z$ 部分的最小值

$$\underset{z}{\min }\lambda \Vert z{\Vert }_1-{\mu }^{T}z=\{\begin{array}{cc}0,& \Vert \mu {\Vert }_{\infty }\le \lambda \\ -\infty & else.\end{array}$$

综上，对偶目标函数为

$$q(\mu )=\underset{x,z}{\min }L(x,z,\mu )=\{\begin{array}{cc}-\frac{1}{4}{\mu }^{T}L{L}^T\mu +{\mu }^{T}Ly,& \Vert \mu {\Vert }_{\infty }\le \lambda \\ -\infty & else.\end{array}$$

那么，对偶问题就为

$$\begin{array}{lc}\max & -\frac{1}{4}{\mu }^{T}L{L}^T\mu +{\mu }^{T}Ly\\ \mathrm{s}.\mathrm{t}.& \Vert \mu {\Vert }_{\infty }\le \lambda .\end{array}$$

使用梯度投影法求解对偶问题

投影的定义

约束域 $C:=\left\{\mu \in {\mathbb{R}}^{n-1}:\Vert \mu {\Vert }_{\infty }\le \lambda \right\}$ 是“盒状的”，即

$$-\lambda \le {\mu }_i\le \lambda ,i=1,2,\dots ,n-1.$$

那么，$\mu$ 在 $C$ 上的投影 $P_C(\mu )\in {\mathbb{R}}^{n-1}$ 为

$${\left[P_C(\mu )\right]}_i:=\{\begin{array}{cc}-\lambda ,& {\mu }_i\le -\lambda \\ {\mu }_i,& -\lambda \le {\mu }_i\le \lambda \\ \lambda ,& \lambda \le {\mu }_i\end{array}$$

其中，$i=1,2,\dots ,n-1$. 特别的，当 $\lambda =1$ 时，即 lambda=1 ，有PcMu=lambda*mu./max(abs(mu),lambda); ，PcMu 是 $\mu$ 在 $C$ 上的投影，即 $P_C(\mu )$.

求解梯度的Lipschitz常数

对偶目标函数 $q(\mu )$ 中，$L{L}^T\in {\mathbb{R}}^{(n-1)\times (n-1)}$ 是一个对称矩阵，对 $\forall \mu \ne 0$ ，有瑞利商（Rayleigh Quotient）

$$\begin{array}{lc}{R}_{L{L}^T}(\mu )& =\frac{{\mu }^{T}L{L}^T\mu }{\Vert \mu {\Vert }^2}=\frac{\Vert L\mu {\Vert }^2}{\Vert \mu {\Vert }^2}=\frac{{\sum }_{i=1}^{n-1}({\mu }_i-{\mu }_{i+1}{)}^2}{\Vert \mu {\Vert }^2}\\ & \le \frac{2\left({\sum }_{i=1}^{n-1}{\mu }_i^2+{\sum }_{i=1}^{n-1}{\mu }_{i+1}^2\right)}{\Vert \mu {\Vert }^2}\le \frac{4\Vert \mu {\Vert }^2}{\Vert \mu {\Vert }^2}=4\end{array}$$

由 定理1.11 知，$R_{L{L}^T}(\mu )\le {\lambda }_{\max }(L{L}^T)$ ，即 ${\lambda }_{\max }(L{L}^T)\le 4$.

${\lambda }_{\max }(L{L}^T)$ 是 $L{L}^T$ 的最大特征值

$L{L}^T$ 的诱导范数（Induced norm）是一个谱范数（Spectral norm），令 $A=L{L}^T$ ，则有

$$\Vert A{\Vert }_{2,2}=\sqrt{{\lambda }_{\max }\left(A^{T}A\right)}={\lambda }_{\max }(A).$$

因此，二次函数 $-\frac{1}{4}{\mu }^{T}L{L}^T\mu +{\mu }^{T}Ly$ 的梯度Lipschitz常数为

$$L^{\prime }=2\times \Vert \frac{1}{4}A{\Vert }_{2,2}=\frac{1}{2}{\lambda }_{\max }(A)=2.$$

综上，可用梯度投影法求解 ${\mu }^{\ast }$ ，令 $t_k=\bar{t}=\frac{1}{L^{\prime }}$，则有

$${\mu }_{k+1}=P_C\left({\mu }_k-\frac{1}{4}L{L}^T{\mu }_k+\frac{1}{2}Ly\right).$$

定理 9.16 和 9.18 保证了 ${\mu }^{\ast }$ 的收敛性.

假设 ${\mu }^{\ast }$ 为梯度投影法所得结果，那么原问题的解为

$$x^{\ast }=y-\frac{1}{2}{L}^T{\mu }^{\ast }.$$

实验步骤

在L1正则化中，设定 $\lambda =1$ ，梯度投影的迭代次数为 $1000$ 次.

%% L2正则化
L=sparse(n-1,n);
for i=1:n-1
L(i,i)=1;
L(i,i+1)=-1;
end
lambda=100;
xde=(speye(n)+lambda*L'*L)\y;
figure(3)
plot(t,xde);
randn('seed',314);
x=zeros(1000,1);
x(1:250)=1;
x(251:500)=3;
x(501:750)=0;
x(751:1000)=2;
figure(1)
plot(1:1000,x,'.')
axis([0,1000,-1,4]);
y=x+0.05*randn(size(x));
figure(2)
plot(1:1000,y,'.')

%% L1正则化
lambda=1;
mu=zeros(n-1,1);
for i=1:1000
mu=mu-0.25*L*(L'*mu)+0.5*(L*y);
mu=lambda*mu./max(abs(mu),lambda);
xde=y-0.5*L'*mu;
end
figure(5)
plot(t,xde,'.');
axis([0,1,-1,4])

结果分析

Figure 12.4. True signal (top image) and noisy signal (bottom image).

使用 L2 惩罚项时，不能正确处理函数图像的三个间断点。但是，使用 L1 惩罚项时，在间断点处，比任何 $\lambda$​ 时的 L2 惩罚项的表现情况都要好。（This result is much better than any of the quadratic regularization reconstructions, and it captures the breakpoints very well.）

参考文献

1. INTRODUCTION TO NONELINEAR OPTIMIZATION. Amir Beck. 2014 : 12.3.10 Denoising
2. THE GRADIENT PROJECTION ALGORITHM, https://sites.math.washington.edu/~burke/crs/408/notes/nlp/gpa.pdf

原创文章！转载需注明来源：©️ Sylvan Ding’s Blog ❤️