使用临近点梯度法、交替方向乘子法、次梯度法优化一个具体问题

问题描述：

考虑线性测量b=Ax+e，其中b为50维的测量值，A为50×100维的测量矩阵，x为100维的未知稀疏向量且稀疏度为5，e为50维的测量噪声。从b与A中恢复x的一范数规范化最小二乘模型如下： $min(1/2)||Ax-b||_2^2+p||x||_1$ 其中p为非负的正则化参数。请设计下述算法求解该问题：

临近点梯度法
交替方向乘子法
次梯度法

在实验中，设x的真值中的非零元素服从均值为0方差为1的高斯分布，A中的元素均服从均值为0方差为1的高斯分布，e中的元素服从均值为0方差为0.1的高斯分布。对于每种算法，请给出每步计算结果与真值的距离以及每步计算结果与最优解的距离。此外，请讨论正则化参数p对计算结果的影响

准备工作：数据生成

数据的生成：生成相应的A和x_real以及e

用到的几个python函数的解析：

#生成服从正太分布的值
numpy.random.normal(loc=0.0, scale=1.0, size=None) 
#loc:均值    scale:标准差         size:形状，可以是int或者tuple
#从list中随机获取num个元素，作为一个片断返回
random.sample(list,num)

故生成数据的代码：

import numpy as np
import random
import matplotlib.pyplot as plt
#定义所用矩阵的维度
A_dim = (50, 100) #矩阵A维度为50×100
x_dim = 100    #X的维度为100
e_dim = 50     #e的维度为50
#定义x的稀疏度：即x中非0元素的数量
x_sparsity=5 
# 生成测量矩阵A，标准正态分布，50×100
A = np.random.normal(0, 1, A_dim)
# 生成测量噪声e，均值为0方差为0.1的正太分布，维度为50
e = np.random.normal(0, 0.1, e_dim)
#生成x的真值：
x_real= np.zeros(x_dim)  #初始化x_real
x_nonzero_index = random.sample(range(x_dim),x_sparsity)  #随机生成x中非0元素的下标
#将非0处的x的值更新为服从标准正态分布的随机值
for i in x_nonzero_index:
    x_real[i] = np.random.randn()
# 计算带噪声的b
b =np.dot(A,x_real) + e

生成相应的数据后，对上述数据分别进行临近点梯度法、交替方向乘子法、此梯度法和次梯度牛顿法进行求解。

临近点梯度法

数学原理

berkeley的一个关于临近点梯度下降的pdf：非常详细地讲述了什么是临界点梯度法（ Proximal Gradient Method ）

https://people.eecs.berkeley.edu/~elghaoui/Teaching/EE227A/lecture18.pdf

邻近点梯度法通常用于求解目标函数不可微的最优化问题，最常用于求解以下形式的优化问题：
$minf_0(x)=s(x)+r(x)$
其中s(x) 光滑可微，r(x) 是非光滑项。

邻近点梯度法:对光滑项做一次梯度下降，得到中间结果 $x_{k+1/2}$ ，然后将这个中间结果代入非光滑项中求其邻近点投影,这是一个迭代过程：
$x_{k+1/2}=x_k-\alpha\bigtriangledown_{}s(x_k)\\ x_{k+1}=argmin_x{}\{r(x)+\frac{1}{2\alpha}||x-x_{k+1/2}||^2\}$
分析可知在此题中：

$s(x)=\frac{1}{2}||Ax-b||_2^2$ ， $r(x)=p||x||_1$

带入上式，得：
$x_{k+1/2}=x_k-\alpha\bigtriangledown_{}s(x_k)=x_k-\alpha A^T(Ax_k-b)\\ x_{k+1}=argmin_x{}\{r(x)+\frac{1}{2\alpha}||x-x_{k+1/2}||^2\}=argmin_x{}\{p||x||_1+\frac{1}{2\alpha}||x-x_{k+1/2}||^2\}$

求 $argmin_x{}\{p||x||_1+\frac{1}{2\alpha}||x-x_{k+1/2}||^2\}$ 时采用软门限法 (Soft Thresholding）

软门限法的一个简单介绍：https://blog.csdn.net/lanyanchenxi/article/details/50448206

利用软门限法我们可以得到最终的迭代策略：
$x_{k+1/2}=x_k-\alpha A^T(Ax_k-b)\\\\ x_{k+1}= \begin{cases} x_{k+\frac{1}{2}}+\alpha p,& x_{k+\frac{1}{2}}<-\alpha p \\ 0,& -\alpha p\leq x_{k+\frac{1}{2}}\leq \alpha p \\x_{k+\frac{1}{2}}-\alpha p,& x_{k+\frac{1}{2}}>\alpha p \end{cases}$
根据该迭代策略coding
代码求解

代码如下:（代码的很多细节在注释中解释）

#-----------计算迭代过程-----------------

#步长α：为固定步长
alpha = 0.005  
#p为正则化参数
p = 10
#迭代的停止条件，当|x[k+1]-x[k]|<stop时显式终止
stop = 1e-5  

xk = np.zeros(x_dim)#存储xk
xk_last = np.zeros(x_dim)#存储上一步的xk
 
iteration_res = []  # 每步计算结果

#迭代次数,当最优解不收敛时也可以作为停止条件，并且记录迭代次数用于比较收敛的效果
times = 1  
 
while times < 1e6:
    #首先对光滑部分做梯度下降，求出x[k+1/2],存到xk_temp中
    xk_temp = xk - alpha* np.dot(A.T,(np.dot(A,xk) - b)) 
    
    # 然后做临近点投影，soft threshold:
    for i in range(x_dim):
        if xk_temp[i] < -alpha * p:
            xk[i] = xk_temp[i] + alpha * p
        elif xk_temp[i] > alpha * p:
            xk[i] = xk_temp[i] - alpha * p
        else:
            xk[i] = 0
    # 记录每步计算结果    
    iteration_res.append(xk.copy())  
    #查看是否终止
    if np.linalg.norm(xk - xk_last) < stop:
        break
    else:
        xk_last = xk.copy()  
        times += 1
        
x_optimal = xk.copy()  # 最优解
dist_real = []  # 每步结果与真值的距离
dist_optm = []  # 每步结果与最优解的距离
#最优解的稀疏度
x_optimal_sparsity=100
for each in x_optimal:
    if each<0.00001 and each>-0.00001:
        x_optimal_sparsity-=1
for iter in iteration_res:
    dist_real.append(np.linalg.norm(iter - x_real))
    dist_optm.append(np.linalg.norm(iter - x_optimal))
#---------------作图------------------- 
plt.plot(dist_real, 'r', label='Distance from real value')
plt.plot(dist_optm, 'g', label='Distance from optimal solution')
plt.title('Proximal Gradient Descent')
plt.xlabel('Iteration times')
plt.ylabel('Distance')
plt.legend()
plt.show()

结果分析

迭代过程中修改正则化参数p的值，查看结果：理论上p值为正则化参数，随着p值的增大，在迭代过程中为了最小化结果，将使得 $|x|_1$ 尽可能地小，因此最终结果的稀疏度将会下降，(在学习神经网络时，解决过拟合时经常采用正则化的方法，并且使用一范数正则化将使得结果的稀疏度降低，实际上就是对神经网络进行了剪枝，和dropout有着相似的作用，当然操作方式不同)

用表格的形式说明一下不同的p和α下的每步计算结果与真值的距离以及每步计算结果与最优解的距离，以及每一个p值对应的最优解稀疏度:(注意的是当行列式值小于 $10^{-5}$ 时认为其已经为0)

步长α	正则化参数p	最优解稀疏度	收敛次数times
0.001	0.01	56	61155
0.001	0.1	46	5087
0.001	1	12	632
0.001	10	5	235
0.005	0.1	46	1814
0.005	1	13	182
0.005	10	3	51

根据结果分析，我们可以发现，在一定范围内，随着正则化参数p的增大，解的稀疏度会降低，并且迭代至收敛的次数明显下降，但若p过大时，可能会导致最终解与真实解偏差较大，。当然，调整步长α，适当地增大α会使得最优解不失太大精度的情况下收敛速度提升，总的来说， $\alpha=0.001,p=1$ 和 $\alpha=0.005,p=1$ 均为不错的结果，在较少的迭代次数内便有一个很好的收敛效果，并且收敛的最优解与真实解的差距较小，该情况下的稀疏度也是可以接收的。

交替方向乘子法

关于交替方向乘子法的详细介绍和推导：

https://www.cnblogs.com/kailugaji/archive/2019/02/25/10433774.html

https://www.zhihu.com/question/36566112

数学原理：

交替方向乘子法（Alternating Direction Method of Multipliers, ADMM),是一种求解具有可分结构的凸优化问题的重要方法，通常用于求解带有如下等式约束的凸优化问题：
$\min f(x)+g(y)\\ \text{s.t. }Ax+By=c$
函数 $f(x),g(y)$ 为关于 $x,y$ 的凸函数。

首先写出问题的增广拉格朗日函数：
$L_p(x,y,\lambda)=f(x)+g(y)+\lambda^T(Ax+By-c)+\frac{\rho}{2}||Ax+By-c||_2^2$ 其中 $\lambda$ 是对偶变量， $\rho$ 是惩罚系数

则ADMM法的迭代策略为：( $x_k$ 表示的是x的第k次迭代结果，此处均表示迭代结果)
$\begin{array}{lr} x_{k+1}=\arg\min_x L_\rho(x,y_k,\lambda_k) \\ y_{k+1}=\arg\min_y L_\rho(x_{k+1},y,\lambda_k) \\ \lambda_{k+1}=\lambda_k+\rho(Ax_{k+1}+By_{k+1}-c) \end{array}$ 将增广拉格朗日函数带入(对结果无影响的项均省略)
$\begin{array}{lr} x_{k+1}=\arg\min_x f(x)+\dfrac{\rho}{2}\left\|Ax+By_k+\frac{\lambda_k}{\rho} \right\|^2_2 \\ y_{k+1}=\arg\min_y g(y)+\dfrac{\rho}{2}\left\|Ax_{k+1}+By+\frac{\lambda_k}{\rho} \right\|^2_2 \\ \lambda_{k+1}=\lambda_k+\rho(Ax_{k+1}+By_{k+1}-c) \end{array}$
思考：当前优化问题为 $min(1/2)||Ax-b||_2^2+p||x||_1$ ，如何使用ADMM法？

当前问题为无约束问题，且只有一个变量 $x$ ，为了使用ADMM法，我们引入变量 $y$ ,令 $y=x$ ,则原问题将转换为有约束问题：
$min\frac{1}{2}||Ax-b||_2^2+p||y||_1\\ s.t. x-y=0$
将该问题带入式子(9)，得：
$\begin{array}{lr} x_{k+1}=\arg\min_x \dfrac{1}{2}\left\|Ax-b\right\|_2^2+\dfrac{\rho}{2}\left\|x-y_k+\frac{\lambda_k}{\rho} \right\|_2^2 \\ y_{k+1}=\arg\min_y p\left \|y\right\|_1+\dfrac{\rho}{2}\left\|x_{k+1}-y+\frac{\lambda_k}{\rho} \right\|_2^2 \\ \lambda_{k+1}=\lambda_k+\rho(x_{k+1}-y_{k+1}) \end{array}$ 当求解 $x_{k+1}=\arg\min_x \dfrac{1}{2}\left\|Ax-b\right\|_2^2+\dfrac{\rho}{2}\left\|x-y_k+\frac{\lambda_k}{\rho} \right\|_2^2$ 时，由于等式右边是两个凸函数，故可以直接求导找导数等于0的点即可：

令
$H(x)=\dfrac{1}{2}\left\|Ax-b\right\|_2^2+\dfrac{\rho}{2}\left\|x-y_k+\frac{\lambda_k}{\rho} \right\|_2^2\\$ 故 $H'(x)=A^T(Ax-b)-\rho(Ix-y_k+\frac{\lambda_k}{\rho})$

令 $H'(x)=0$

则： $x=\frac{A^Tb+\rho y_{k}-\lambda_k}{A^TA+\rho I}$

因此 $x_{k+1}=\frac{A^Tb+\rho y_{k}-\lambda_k}{A^TA+\rho I}$ ，即该式是一个显式数值解

求解 $y_{k+1}=\arg\min_y p\left\|y\right\|_1+\dfrac{\rho}{2}\left\|x_{k+1}-y+\frac{\lambda_k}{\rho} \right\|_2^2$ 时，由于一阶范数不可导，故不可以在此处直接求导等于0来解决，因此在此处仍使用软门限(Soft Thresholding)的方法：
$y_{k+1}= \begin{cases} x_{k+1}+\frac{\lambda _k}{\rho}+\frac{p}{\rho},& x_{k+1}+\frac{\lambda _k}{\rho}<-\frac{p}{\rho} \\ 0,& \left| x_{k+1}+\frac{\lambda _k}{\rho}\right|\le \frac{p}{\rho} \\ x_{k+1}+\frac{\lambda_k}{\rho}-\frac{p}{\rho},& x_{k+1}+\frac{\lambda _k}{\rho}>\frac{p}{\rho} \end{cases}$
根据该迭代策略可以写出代码
代码求解

#------------生成数据的过程与上述一致，在此省略

#-----------计算迭代过程-----------------
#parameters:
rho = 0.001
p = 0.1       # 正则化参数

#迭代的停止条件，当|x[k+1]-x[k]|<stop时显式终止
stop = 1e-5  

#初始化
xk = np.zeros(x_dim)
yk = np.zeros(x_dim)
lambda_k = np.zeros(x_dim)
xk_last = np.zeros(x_dim)#存上一步的xk，用于作减法

iteration_res = []  # 每步计算结果
k = 1  # 记录迭代次数

while k < 1e5:  #避免结果不收敛导致死循环
    # 更新xk，直接带入公式即可
    xk = np.dot(np.linalg.inv(np.dot(A.T,A) + rho*np.eye(x_dim, x_dim)),(np.dot(A.T,b) + rho*yk - lambda_k))
    
    # 更新yk,用软门限法
    for i in range(x_dim):
        if xk[i] + lambda_k[i]/rho < -p/rho:
            yk[i] = xk[i] + lambda_k[i]/rho + p/rho
        elif xk[i] + lambda_k[i]/rho > p/rho:
            yk[i] = xk[i] + lambda_k[i]/rho - p/rho
        else:
            yk[i] = 0
            
    # 更新lambda
    lambda_k = lambda_k + rho * (xk - yk)
    
    # 记录每步计算结果
    iteration_res.append(xk.copy()) 
    #已经收敛
    if np.linalg.norm(xk - xk_last) < stop:
        break
    else:
        xk_last = xk.copy()  
        k += 1
        
x_optimal = xk[:]  # 最优解

#最优解的稀疏度
x_optimal_sparsity=100
for each in x_optimal:
    if each<0.0001 and each>-0.0001:
        x_optimal_sparsity-=1

dist_real = []  # 对应每步结果与真值的距离
dist_optm = []  # 对应每步结果与最优解的距离
for iter in iteration_res:
    dist_real.append(np.linalg.norm(iter - x_real))
    dist_optm.append(np.linalg.norm(iter - x_optimal))

#---------------作图------------------- 
plt.plot(dist_real, 'r', label='Distance from real value')
plt.plot(dist_optm, 'g', label='Distance from optimal solution')
plt.title('ADMM')
plt.xlabel('Iteration times')
plt.ylabel('Distance')
plt.legend()
plt.grid()
plt.show()

结果分析

用表格的形式说明一下不同的p和α下的每步计算结果与真值的距离以及每步计算结果与最优解的距离，以及每一个p值对应的最优解稀疏度:(注意的是当矩阵中元素值小于 $10^{-5}$ 时认为其已经为0)

正则化参数p	解的稀疏度	收敛次数times
0.01	54	1804
0.1	99	1073
0.5	100	1763
1	99	2630
5	99	10590

在此处我发现调整正则化参数p并没有使得稀疏度减小，这是令我困惑的。当然，改变正则化参数仍将对解的收敛速度和收敛效果产生较大影响，对上述进行对比，易得当p=0.5时效果最好，此时收敛速度非常快且最优解与真实值的差距较小，当然此时的x稀疏度仍是较高的。

次梯度法

次梯度法的介绍：(非常详细的介绍)

https://blog.csdn.net/qq_32742009/article/details/81704139

https://blog.csdn.net/quiet_girl/article/details/79648124

次梯度法通常用来处理不可导的凸函数，解决不可导梯度问题，该方法可以使得不可导点处仍可以用次梯度表示从而继续进行迭代，缺点是算法的收敛速度会较慢。

数学原理

次梯度：(当该点处可导时次梯度即为梯度，不可导时次梯度为一个范围内的任意数)

在点 $x_0$ 的次导数的集合是一个非空闭区间[a, b]，其中a和b是单侧极限:

$a=lim_{x−>x_0^-}\frac{f(x)-f(x_0)}{x-x_0}$ ， $b=lim_{x−>x_0^+}\frac{f(x)-f(x_0)}{x-x_0}$

所有次导数的集合[a, b]称为函数 $f(x)$ 在$ x_0 $的次导数

迭代策略：
$x^{k+1}=x^{k}-\alpha\times g_0(x^k)$
其中， $g_0(x)\in \partial f_0(x)$ ， $\partial f_0(x)$ 表示 $f_0(x)$ 在 $x$ 处的次梯度。

因此将该问题带入：
$f(x)=min(1/2)||Ax-b||_2^2+p||x||_1$
因此对 $f(x)$ 求导得：
$\dfrac{\partial f(x)}{\partial x}=A^T(Ax+b)+ p\dfrac{\partial \left\|x\right\|_1}{\partial x}$
由于 $||x||_1$ 在 $x=0$ 处不可导，故对于 $||x||_1$ 的梯度需要分段考虑
$\begin{cases} \dfrac{\partial \left\|x\right\|_1}{\partial x}=-1,& x<0 \\ -1\le \dfrac{\partial \left\|x\right\|_1}{\partial x}\le 1,& x=0 \\ \dfrac{\partial \left\|x\right\|_1}{\partial x}=1,& x>0 \end{cases}$
故迭代策略：
$g_0(x^k)= \begin{cases} A^T(Ax+b)-p,& x^k<0 \\ A^T(Ax+b)+ps,s\in[-1,1],& x^k=0 \\ A^T(Ax+b)+p=1,& x^k>0 \end{cases}$ $x^{k+1}=x^{k}-\alpha\times g_0(x^k)$ 需要注意的是，由于次梯度的选择并不一定是下降的方向，故可能导致固定步长不收敛。因此通常会选择递减步长进行次梯度。
根据该迭代策略很容易完成编程。

代码求解

#生成数据与之前一致，在此省略
#---------------
#-----------计算迭代过程-----------------
#parameters:
alpha = 0.005
p = 0.1       # 正则化参数

#迭代的停止条件，当|x[k+1]-x[k]|<stop时显式终止
stop = 1e-5  

#初始化
xk = np.zeros(x_dim)
xk_last = np.zeros(x_dim)#存上一步的xk，用于作减法

iteration_res = []  # 每步计算结果
k = 1  # 记录迭代次数

while k < 1e5:  #避免结果不收敛导致死循环

    l1_subgradient = np.zeros(x_dim)
    # 更新l1的次梯度，分情况讨论
    for i in range(x_dim):
        if xk[i]!=0:  #可导则直接求导
            l1_subgradient[i] = np.sign(xk[i])
        else:      #否则随机取次梯度
            l1_subgradient[i]=np.random.uniform(-1, 1)
    #更新梯度
    subgradient=np.dot(A.T,(np.dot(A,xk))-b)+p * l1_subgradient
    
    xk=xk-alpha*subgradient
    # 记录每步计算结果
    iteration_res.append(xk.copy()) 
    #已经收敛
    if np.linalg.norm(xk - xk_last) < stop:
        break
    else:
        xk_last = xk.copy()  
        k += 1
#最优解        
x_optimal = xk[:] 

#最优解的稀疏度
x_optimal_sparsity=100
for each in x_optimal:
    if each<0.00001 and each>-0.00001:
        x_optimal_sparsity-=1

dist_real = []  # 对应每步结果与真值的距离
dist_optm = []  # 对应每步结果与最优解的距离
for iter in iteration_res:
    dist_real.append(np.linalg.norm(iter - x_real))
    dist_optm.append(np.linalg.norm(iter - x_optimal))

#---------------作图------------------- 
plt.plot(dist_real, 'r', label='Distance from real value')
plt.plot(dist_optm, 'g', label='Distance from optimal solution')
plt.title('SubGradient')
plt.xlabel('Iteration times')
plt.ylabel('Distance')
plt.legend()
plt.grid()
plt.show()

结果分析

固定步长/递减步长	正则化参数p的值	解的稀疏度	收敛次数times
固定 $\alpha$ =0.001	0.1	95	不收敛
固定 $\alpha$ =0.001	1	98	不收敛
递减步长 $\alpha_{k}=0.001/\sqrt{k}$	0.1	97	7822
递减步长 $\alpha_{k}=0.001/\sqrt{k}$	1	15	不收敛
递减步长 $\alpha_{k}=0.001/\sqrt{k}$	10	63	不收敛

分析上述结果，发现效果都不是特别好，当为固定步长时，解是不收敛的。当为递减步长时，可以收敛但收敛效果并不特别好，当p较大时又会导致不收敛。因此分析思考如果限制前n步为固定步长，之后为递减步长，效果将如何？尝试前n次迭代为固定步长 $\alpha=0.001$ ,之后为 $1/\sqrt{k-n}$ 程度的递减:

递减步长	正则化参数p的值	每步计算结果与真值的距离以及每步计算结果与最优解的距离	解的稀疏度	收敛次数times
前5000步固定步长，5000步后以 $0.001/\sqrt{k-5000}$ 递减步长	0.1		46	7659
前10000步固定步长，10000步后以 $0.001/\sqrt{k-10000}$ 递减步长	0.1		43	12346