正则化实现降噪，分别使用最小二乘、定步长梯度下降和回溯法的梯度下降求解最优解

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

实验目的

参考 INTRODUCTION TO NONELINEAR OPTIMIZATION. Amir Beck. 2014 的 3.4 Denoising 相关内容，分别使用 最小二乘法、定步长的梯度下降（constant）和回溯法的梯度下降（backtracking），实现对 Example 3.4 中带有噪声信号的降噪过程，对比分析采用不同方法的降噪效果。添加不同的噪声，或以不同的方式添加噪声，观察上述三种降噪方法的降噪效果，并简要分析结果。了解并掌握最小二乘法和惩罚问题（Penalized Problem），学会使用MATLAB求解并分析相关问题。

实验环境

MATLAB R2021a

实验内容

考虑一段信号 $x\in {\mathbb{R}}^{300}$ ，在该信号上添加一段高斯加性噪声（$\mu =0,\sigma =0.05$）得到 $b$，由下述MATLAB语句生成。要求使用正则化的最小二乘（RLS）、定步长的梯度下降、回溯法的梯度下降分别对加噪后对信号 $b$ 进行降噪处理，并分析各方法降噪效果的优劣。

t=linspace(0,4,300)';
x=sin(t)+t.*(cos(t).^2);
randn('seed',314);
b=x+0.05*randn(300,1);

真实信号与噪声信号的对比图

subplot(1,2,1);
plot(1:300,x,'LineWidth',2);
subplot(1,2,2);
plot(1:300,b,'LineWidth',2);

算法描述

降噪问题的一般表述为：给定含有噪声的信号 $b$，找到一个“足够好”的估计值 $x^{\ast }$，使得 $x^{\ast }\approx b$，即 $x^{\ast }=\arg \underset{x}{\min }\Vert x-b{\Vert }^2$. 显然，这个优化问题的解是 $x^{\ast }=b$。然而，这一最优解是没有意义的，即产生了过拟合的解。

因此，为了解决过拟合的问题，我们需要向目标函数中添加正则化项。不妨假设，最初的信号是一段光滑的曲线，那么，向量 $x$ 相邻两元素间的的差值也应当尽可能的小（接近于0但不等于0）. 这样，正则化项就应当选用二次惩罚（Quadratic Penalty），即 $R(x)={\sum }_{i=1}^{n-1}(x_i-x_{i+1}{)}^2,n=300$. 也可以用更加紧凑的向量形式表示：$R(x)=\Vert Lx{\Vert }^2,L\in {\mathbb{R}}^{(n-1)\times n},n=300$,

$$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 }=\arg \underset{x}{\min }\Vert x-b{\Vert }^2+\lambda \Vert Lx{\Vert }^2,\lambda >0$. 其中，$\lambda$ 为正则化参数（又称惩罚系数），$\lambda$ 越大，惩罚力度越强，$x^{\ast }$ 越“平滑”，反之亦然。当 $\lambda =0$ 时，问题转化为一般的最小二乘问题（LS），也就会产生上述无意义的最优解。

最小二乘解

$$x^{\ast }(\lambda )=\arg \underset{x}{\min }\left\{f(x)\equiv x^T\left(I+\lambda L^{T}L\right)x-2b^{T}x+\Vert b{\Vert }^2\right\}$$

显然，目标函数 $f(x)$ 是一个二次函数，而且必然有 ${\nabla }^{2}f(x)=2\left(I+\lambda L^{T}L\right)\succ 0$，因此，由定理 2.41 可知，该函数的任何驻点都是全局最小点。驻点满足 $\nabla f(x)=0$，即 $\left(I+\lambda L^{T}L\right)x=b$，所以最优解 $x^{\ast }(\lambda )={\left(I+\lambda L^{T}L\right)}^{-1}b$.

梯度下降

目标函数：$f(x)=x^T\left(I+\lambda L^{T}L\right)x-2b^{T}x+\Vert b{\Vert }^2$

目标函数的梯度：$g(x)=\nabla f(x)=\left(I+\lambda L^{T}L\right)x-b$

梯度下降算法描述：

设定 $ϵ$；

初始化 $x_0=0,x_0\in {\mathbb{R}}^n,n=300$；

FOR k = 0, 1, 2, …

1. 选取下降方向 $d_k=g(x_k)$；
2. 选取步长 $t_k$，使得 $f(x_k+t_k{d}_k)<f(x_k)$；
3. 设置 $x_{k+1}=x_k+t_k{d}_k$；
4. IF $\Vert g(x_{k+1})\Vert \le ϵ$ THEN STOP；

在算法循环的 2. 步长选取中，对于定步长的梯度下降来说，$t_k$ 为一常量，即 $t_k=\bar{t}$. 而对于采用回溯法的非精确线搜索来说，令步长的初值 $t_k=s,s>0$，更新步长 $t_k\leftarrow \beta t_k,\beta \in (0,1)$，直到满足 $f(x_k)-f(x_k+t_k{d}_k)\ge -\alpha t_k\nabla f(x_k{)}^T{d}_k,\alpha \in (0,1)$，可以证明，当 $t\in [0,ϵ]$，上述不等式总成立。

实验步骤

% 生成L
L=zeros(299,300);
for i=1:299L(i,i)=1;L(i,i+1)=-1;
end

function draw_pic(x,b,optimum,lambda,pic_title)
% 绘图函数
% 输入：
% x是不含噪声的源信号
% b是含有噪声的原信号
% optimum的每列是对原信号b降噪后的一个解
% lambda是选用的正则化惩罚因子
% pic_title是图片的标题
color_prism=colormap(prism);
plot(1:300,b,'LineWidth',1,'DisplayName','signal with noise',...'Color',color_prism(1,:));
hold on
plot(1:300,x,'LineWidth',0.7,'DisplayName','signal without noise',...'Color','b');
[~,sz]=size(optimum);
for j=1:szp=plot(1:300,optimum(:,j),'LineWidth',0.8,'DisplayName',...['denoised signal with \lambda=',num2str(lambda(j)),...', SSE=',num2str(norm(optimum(:,j)-x)^2)],...'Color',color_prism(2*j,:));p.Color(4)=0.5; % set transparent
end
hold off
xlim([0 300])
ylim([-0.5 3.5])
legend('Location','northwest')
title(pic_title)

最小二乘

lambda=[1 10 100 1000];
x_rls=[];
for l=lambdax_rls=[x_rls,(eye(300)+l*L'*L)\b];
end
draw_pic(x,b,x_rls,lambda,'RLS solutions')

定步长梯度下降

lambda=[0.5 5 50];
x0=zeros(300,1);
t=0.005;
eps=5;
x_constgd=[];
for l=lambdaf=@(x) x'*(eye(300)+l*L'*L)*x-2*b'*x+norm(b)^2;g=@(x) (eye(300)+l*L'*L)*x-b;[x_constgd_next,~]=...gradient_method_constant(f,g,x0,t,eps);x_constgd=[x_constgd,x_constgd_next];
end
draw_pic(x,b,x_constgd,lambda,[...'GMC solutions, t=',num2str(t),...'x_0=(0,0,\dots)',...' \epsilon=',num2str(eps)])

function [x,fun_val]=gradient_method_constant(f,g,x0,t,epsilon)
% Gradient method with constant stepsize
%
% INPUT
%=======================================
% f ......... objective function
% g ......... gradient of the objective function
% x0......... initial point
% t ......... constant stepsize
% epsilon ... tolerance parameter
% OUTPUT
%=======================================
% x ......... optimal solution (up to a tolerance)
% of min f(x)
% fun_val ... optimal function value
x=x0;
grad=g(x);
iter=0;
while (norm(grad)>epsilon)iter=iter+1;x=x-t*grad;fun_val=f(x);grad=g(x);fprintf('iter_number = %3d norm_grad = %2.6f fun_val = %2.6f \n',...iter,norm(grad),fun_val);
end

回溯法梯度下降

lambda=[0.1 5 50];
x0=zeros(300,1);
s=0.1;
alpha=0.5;
beta=0.3;
eps=0.7;
x_constgd=[];
for l=lambdaf=@(x) x'*(eye(300)+l*L'*L)*x-2*b'*x+norm(b)^2;g=@(x) (eye(300)+l*L'*L)*x-b;[x_constgd_next,~]=...gradient_method_backtracking(f,g,x0,s,alpha,beta,eps);x_constgd=[x_constgd,x_constgd_next];
end
draw_pic(x,b,x_constgd,lambda,[...'GMB solutions, s=',num2str(s),...' x_0=(0,0,...)',...' \alpha=',num2str(alpha),...' \beta=',num2str(beta),...' \epsilon=',num2str(eps)])

function [x,fun_val]=gradient_method_backtracking(f,g,x0,s,alpha,...
beta,epsilon)
% Gradient method with backtracking stepsize rule
%
% INPUT
%=======================================
% f ......... objective function
% g ......... gradient of the objective function
% x0......... initial point
% s ......... initial choice of stepsize
% alpha ..... tolerance parameter for the stepsize selection
% beta ...... the constant in which the stepsize is multiplied
% at each backtracking step (0<beta<1)
% epsilon ... tolerance parameter for stopping rule
% OUTPUT
%=======================================
% x ......... optimal solution (up to a tolerance)
% of min f(x)
% fun_val ... optimal function value
x=x0;
grad=g(x);
fun_val=f(x);
iter=0;
while (norm(grad)>epsilon)iter=iter+1;t=s;while (fun_val-f(x-t*grad)<alpha*t*norm(grad)^2)t=beta*t;endx=x-t*grad;fun_val=f(x);grad=g(x);fprintf('iter_number = %3d norm_grad = %2.6f fun_val = %2.6f \n',...iter,norm(grad),fun_val);
end

结果分析

对于降噪后结果的评估使用和方差（SSE），其计算公式为 $SSE=\Vert x^{\ast }-x{\Vert }^2$.

最小二乘（RLS）

结果表明，当 $\lambda$ 取值越大时，RLS的解变得越平滑。当 $\lambda =1$ 时，RLS的解非常接近于有噪声的信号 $b$，不够平滑。当 $\lambda =100$ 时，我们得到了一条更加平滑的RLS信号，但是，相较于 $\lambda =10$ 的RLS结果来说，精确度更低，特别是在靠近边界的地方。RLS的解在 $\lambda =1000$ 时变得非常平滑，但是也更加偏离真实信号。

定步长梯度下降（GMC）

$x_0=(0,0,\dots ),\bar{t}=0.005,ϵ=5$

$\lambda$	Epochs
0.5	370
5	370
50	368

定步长梯度下降的精度整体偏低，其解严重偏离真实信号，特别是在波峰处。在实验时，发现定步长 $\bar{t}$ 的选取不能过大，否则目标函数发散，函数值持续增长到正无穷，梯度下降偏离下降中心，相反，若 $\bar{t}$ 选取过小，将导致收敛速度极慢。$ϵ$ 选取过小，导致目标函数不收敛，$ϵ$ 选取过大，将造成上述结果不精确到情况。

回溯法梯度下降（GMB）

$x_0=(0,0,\dots ),s=0.1,\alpha =0.5,\beta =0.3,ϵ=0.7$

$\lambda$	Epochs
0.1	37
5	41
50	335

回溯法的梯度下降精确度较高。在 $\lambda$ 较小时，需要的迭代步数较少，随着 $\lambda$ 增大，迭代步数随之增加，但GMB曲线也更加平滑。

总结

本次试验对比了最小二乘、定步长梯度下降和回溯法的梯度下降法对添加高斯加性噪声的信号的降噪情况，通过添加正则项达到降噪的目的。最小二乘法的结果受惩罚系数 $\lambda$ 的影响较大，但当 $\lambda$ 选取合适时，将实现较好的降噪效果，得到比较光滑且接近原信号的RLS解；定步长梯度下降的降噪结果严重偏离真实信号，且对 $\bar{t}$ 和 $ϵ$ 的选取要求较高，目标函数易发散。回溯法的梯度下降所得降噪的结果不仅精度高，而且可通过增大 $\lambda$ ，提高降噪效果，且不会使降噪后的信号严重偏离原信号。

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