本文地址：https://www.jeddd.com/article/use-python-to-solve-lasso-by-pgd-admm-and-subgd.html
问题描述
考虑线性测量 $b=Ax+e$，其中 $b$ 为 50 维的测量值，$A$ 为 50×100 维的测量矩阵，$x$ 为 100 维的未知稀疏向量且稀疏度为 5，$e$ 为 50 维的测量噪声。从 $b$ 与 $A$ 中恢复 $x$ 的一范数规范化最小二乘模型如下：
$$ \min \dfrac{1}{2}\left\|Ax-b\right\|_2^2+p\left\|x\right\|_1 $$
其中 $p$ 为非负的正则化参数。请设计下述算法求解该问题：
- 邻近点梯度下降法；
- 交替方向乘子法；
- 次梯度法。
在实验中，设 $x$ 的真值中的非零元素服从均值为 0 方差为 1 的高斯分布，$A$ 中的元素服从均值为 0 方差为 1 的高斯分布，$e$ 中的元素服从均值为 0 方差为 0.1 的高斯分布。对于每种算法，请给出每步计算结果与真值的距离以及每步计算结果与最优解的距离。此外，请讨论正则化参数 $p$ 对计算结果的影响。
思路
该实验的思路是，首先随机生成数据，然后依次用邻近点梯度下降法、交替方向乘子法、次梯度法在生成的数据上进行求解，同时对于每个方法每次迭代的效果进行绘图。
在本文中，未经特殊说明的范数均为2-范数；同样地，如果范数符号没有下标，也代表2-范数。
实验环境为 Python 3.7 和 Jupyter Notebook，用到的外部包只有 numpy 和 matplotlib。
完整代码：https://github.com/Jed-Z/optimization-theory/blob/master/Opt1_LASSO.ipynb
准备工作
首先，引入所需要的包，并进行一些配置；然后定义一些常量，包括矩阵和向量的维度、稀疏度。
import numpy as np
import random
import matplotlib.pyplot as plt
# 使Matplotlib输出矢量图
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
# 常量
A_dim = (50, 100)
x_dim = 100
x_sparsity = 5 # 稀疏度
e_dim = 50本文地址：https://www.jeddd.com/article/use-python-to-solve-lasso-by-pgd-admm-and-subgd.html
生成数据
依次生成 $x$ 的真值、矩阵$A$、向量 $e$，然后计算出向量 $b$。“$x$ 的稀疏度为 5”的意思是，该向量的 100 个元素中仅有 5 个非零，其余 95 个元素都为 0。另外注意高斯分布就是正态分布，可以用np.random.normal来生成服从给定 $\mu$ 和 $\sigma ^2$ 的正态分布。
# 生成x的真值
x_nonzero_index = random.sample(range(x_dim), 5) # 非0元素的下标
x_nonzero_element = np.random.normal(0, 1, x_sparsity)
x_real = np.zeros(x_dim)
for i in range(len(x_nonzero_element)):
x_real[x_nonzero_index] = x_nonzero_element
# 生成测量矩阵A
A = np.random.normal(0, 1, A_dim)
# 生成测量噪声e
e = np.random.normal(0, 0.1, e_dim)
# 计算带噪声的b
b = A @ x_real + e上面生成的一组数据将被用于邻近点梯度下降法、交替方向乘子法、次梯度法这三种方法分别求解。
邻近点梯度下降法
算法通用公式
邻近点梯度下降法（Proximal Gradient Method），更常见的译名为“近端梯度法”，常缩写为 PGD（“D”为 descent，意为“下降”）。邻近点梯度法常用于求解以下形式的优化问题：
$$ \min f_0(x)=s(x)+r(x) $$
其中 $s(x)$ 光滑可微，$r(x)$ 非光滑不可微。
邻近点梯度法的过程是，首先对光滑项做一次梯度下降，得到中间结果 $x^{k+\frac{1}{2}}$，然后将这个中间结果代入非光滑项中求其邻近点投影（proximal map），即完成一次迭代。算法如下（上标代表迭代次数而非指数）：
$$ \begin{equation} \left\{ \begin{array}{lr} x^{k+\frac{1}{2}}=x^k-\alpha \times \nabla s(x^k) \\ x^{k+1}=\arg\min_x\{ r(x)+\dfrac{1}{2\alpha}\left\|x-x^{k+\frac{1}{2}}\right\|^2 \}\\ \end{array} \right. \end{equation} $$
将算法用于本题
在本问题中，$s(x)=\dfrac{1}{2}\left\|Ax-b\right\|_2^2$，$r(x)=p\left\|x\right\|_1$，将它们带入上面的公式中，得到：
$$ \begin{equation} \left\{ \begin{array}{lr} x^{k+\frac{1}{2}}=x^k-\alpha \times A^T(Ax^k-b) \\ x^{k+1}=\arg\min_x\{ p\left\|x\right\|_1+\dfrac{1}{2\alpha}\left\|x-x^{k+\frac{1}{2}}\right\|^2 \} \end{array} \right. \end{equation} $$
求解上述迭代方程组第二行的 $\arg\min_x$ 时用到了软阈值（Soft Thresholding，又译做“软门限”），软阈值函数的推导较为简单，只需做分类讨论即可，这里就不过多涉及了，我们直接将软阈值的结论应用在 $x^{k+1}$ 的更新中。
$$ \begin{equation} x^{k+1}= \begin{cases} x^{k+\frac{1}{2}}+p,& x^{k+\frac{1}{2}}<-\alpha \times p \\ 0,& |x^{k+\frac{1}{2}}|\le \alpha \times p \\ x^{k+\frac{1}{2}}-p,& x^{k+\frac{1}{2}}>\alpha \times p \end{cases} \end{equation} $$
编写代码
Python 变量与数学符号的对应关系如下表：
| Python 变量 | 数学符号 | 含义 |
|---|---|---|
| k | $k$ | 迭代次数 |
| x_k | $x^k$ | 第 k 次迭代的解 |
| x_k_temp | $x^{k+\frac{1}{2}}$ | 对光滑部分梯度下降后的中间结果 |
| x_k_old | $x^{k-1}$ | 上次迭代的解 |
| x_optm | $x^*$ | 最优解 |
| epsilon | $\epsilon$ | 容许的最大误差，用于判断是否收敛 |
| alpha | $\alpha$ | 步长，也称为学习率 |
| p | $p$ | L1 正则化参数 |
下面这段代码按以上公式进行迭代，直到最优解收敛。我使用$\left\|x^k-x^{k-1}\right\|<\epsilon$作为收敛条件，为了防止出现迭代不收敛的情况，我为迭代次数设置了一个上限 100 万次，保证能够跳出循环而不会产生死循环。
alpha = 0.001 # 固定步长
p = 0.01 # 正则化参数
epsilon = 1e-5 # 最大允许误差
x_k = np.zeros(x_dim)
x_k_old = np.zeros(x_dim)
plot_iteration = [] # 每步计算结果
k = 1 # 迭代次数
while k < 1e6:
x_k_temp = x_k - alpha * A.T @ (A@x_k - b) # 对光滑部分做梯度下降
# 临近点投影（软阈值）
for i in range(x_dim):
if x_k_temp[i] < -alpha * p:
x_k[i] = x_k_temp[i] + alpha * p
elif x_k_temp[i] > alpha * p:
x_k[i] = x_k_temp[i] - alpha * p
else:
x_k[i] = 0
plot_iteration.append(x_k.copy()) # 记录每步计算结果
if np.linalg.norm(x_k - x_k_old) < epsilon:
break
else:
x_k_old = x_k.copy() # 深拷贝
k += 1
x_optm = x_k[:] # 最优解用以下语句查看最优解的稀疏度（即最优解中非零元素的个数）：
np.count_nonzero(x_optm)上述代码中，将每次迭代的结果 $x^k$ 都加入了列表plot_iteration，这是为了方便作图观察 $x^k$ 与最优解 $x^*$ 的距离和 $x^k$ 与真值 x_real 的距离。使用以下代码作图：
plot_dist_real = [] # 每步结果与真值的距离
plot_dist_optm = [] # 每步结果与最优解的距离
for iter in plot_iteration:
plot_dist_real.append(np.linalg.norm(iter - x_real))
plot_dist_optm.append(np.linalg.norm(iter - x_optm))
# 作图
plt.title('Proximal Gradient Descent')
plt.xlabel('Iteration times')
plt.ylabel('Distance')
plt.plot(plot_dist_real, 'r', label='Distance from true value of x')
plt.plot(plot_dist_optm, 'b', label='Distance from the optimal solution')
plt.grid()
plt.legend()
plt.show()结果分析
改变正则化参数 $p$，多次重复以上实验，结果如下：
| 正则化参数 $p$ | 图像 | $x^*$ 的稀疏度 |
|---|---|---|
| 0.01 |  | 60 |
| 0.1 |  | 47 |
| 1 |  | 11 |
| 10 |  | 4 |
可见正则化参数 $p$ 的大小影响最优解的稀疏度，$p$ 越大，最优解越稀疏；随着 $p$ 的增大，收敛的迭代次数变少。在上面的实验中，当 $p=1$ 时，收敛速度足够快，回归效果也足够好。（由于数据是随机生成的，因此以上结果并不是固定的。）
交替方向乘子法
算法通用公式
交替方向乘子法（Alternating Direction Method of Multipliers, ADMM）常用于求解以下形式的优化问题：
$$ \min f_1(x)+f_2(y)\\ \text{s.t. }Ax+By=0 $$
算法如下（上标代表迭代次数而非指数）：
$$ \begin{equation} \left\{ \begin{array}{lr} x^{k+1}=\arg\min_x L_c(x,y^k,v^k) \\ y^{k+1}=\arg\min_y L_c(x^{k+1},y,v^k) \\ v^{k+1}=v^k+c(Ax^{k+1}+By^{k+1}) \end{array} \right. \end{equation} $$
其中$L_c(x,y,v)$为原目标函数的增广拉格朗日函数，$c$ 为惩罚项系数。注意，有些文献中，用字母 $c$ 表示原问题等式约束的常数项（上式中该常数项为 0），用字母 $\rho$ 表示增广拉格朗日函数的惩罚项系数，务必注意区分。
以上算法可以等价变形为如下形式（与 $\arg\min$ 无关的项可以直接删去，内积与范数可以合并）：
$$ \begin{equation} \left\{ \begin{array}{lr} x^{k+1}=\arg\min_x f_1(x)+\dfrac{c}{2}\left\|Ax+By^k+\frac{v^k}{c} \right\|^2 \\ y^{k+1}=\arg\min_y f_2(y)+\dfrac{c}{2}\left\|Ax^{k+1}+By+\frac{v^k}{c} \right\|^2 \\ v^{k+1}=v^k+c(Ax^{k+1}+By^{k+1}) \end{array} \right. \end{equation} $$
将算法用于本题
通过引入一个新的变元 $y$，将原问题变形乘适用于交替方向乘子法的形式：
$$ \begin{align*} \min& \dfrac{1}{2}\left\|Ax-b\right\|_2^2+p\left\|y\right\|_1 \\ \text{s.t. }& x-y=0 \end{align*} $$
然后按照以上算法得出迭代方程组：
$$ \begin{equation} \left\{ \begin{array}{lr} x^{k+1}=\arg\min_x \dfrac{1}{2}\left\|Ax-b\right\|_2^2+\dfrac{c}{2}\left\|x-y^k+\frac{v^k}{c} \right\|_2^2 \\ y^{k+1}=\arg\min_y p\left\|y\right\|_1+\dfrac{c}{2}\left\|Ax^{k+1}+By+\frac{v^k}{c} \right\|_2^2 \\ v^{k+1}=v^k+c(x^{k+1}-y^{k+1}) \end{array} \right. \end{equation} $$
上述迭代方程组中，$x^{k+1}$ 的 $\arg\min$ 是可微的，可以通过令其导数等于 0 的方式解出 $x$：
$$ x^{k+1}=(A^TA+cI)^{-1}(A^Tb+cy^k-v^k) $$
$y^{k+1}$ 中有1-范数，因此它不可微，但它可通过软阈值（Soft Thresholding）的方法求解：
$$ \begin{equation} y^{k+1}= \begin{cases} x^{k+1}+\frac{v^k}{c}+\frac{p}{c},& x^{k+1}+\frac{v^k}{c}<-\frac{p}{c} \\ 0,& \left| x^{k+1}+\frac{v^k}{c}+\frac{p}{c}\right|\le \frac{p}{c} \\ x^{k+1}+\frac{v^k}{c}-\frac{p}{c},& x^{k+1}+\frac{v^k}{c}>-\frac{p}{c} \end{cases} \end{equation} $$
$v^{k+1}$ 可以对 $v^k$、$x^{k+1}$ 和 $y^{k+1}$ 通过简单的线性运算得到。
编写代码
Python 变量与数学符号的对应关系与前面“邻近点梯度下降法”中的对应关系基本相同，不再赘述。
# 调参区
c = 0.005
alpha = 1/c
p = 0.01 # 正则化参数
epsilon = 1e-5 # 最大允许误差
作图：
plot_dist_real = [] # 每步结果与真值的距离
plot_dist_optm = [] # 每步结果与最优解的距离
for iter in plot_iteration:
plot_dist_real.append(np.linalg.norm(iter - x_real))
plot_dist_optm.append(np.linalg.norm(iter - x_optm))
# 作图
plt.title('Alternating Direction Method of Multipliers')
plt.xlabel('Iteration times')
plt.ylabel('Distance')
plt.plot(plot_dist_real, 'r', label='Distance from true value of x')
plt.plot(plot_dist_optm, 'b', label='Distance from the optimal solution')
plt.grid()
plt.legend()
plt.show()结果分析
改变正则化参数 $p$，多次重复以上实验，结果如下：
| 正则化参数 $p$ | 图像 |
|---|---|
| 0.01 |  |
| 0.1 |  |
| 1 |  |
| 10 |  |
本文地址：https://www.jeddd.com/article/use-python-to-solve-lasso-by-pgd-admm-and-subgd.html
次梯度法
算法通用公式
用 $g_0(x)\in \partial f_0(x)$ 表示 $f_0(x)$ 在 $x$ 处的次梯度。所谓次梯度法，就是在不可微点处用次梯度代替梯度，然后使用梯度下降法。当函数可微时，与梯度下降法无异。
$$ x^{k+1}=x^{k}-\alpha\times g_0(x^k) $$
次梯度法的步长选取：
- 固定步长：$\alpha^k=\alpha$
- 不可加，但平方可加：$\alpha^k=\dfrac{1}{k}$
- 不可加、平方不可加，但能收敛至 0：$\alpha^k=\dfrac{1}{\sqrt{k}}$
将算法用于本题
若令原问题的目标函数 $f(x)=\dfrac{1}{2}\left\|Ax-b\right\|_2^2+p\left\|x\right\|_1$，则 $\dfrac{\partial f(x)}{\partial x}=A^T(Ax+b)+ p\dfrac{\partial \left\|x\right\|_1}{\partial x}$。L1正则化项在其元素为 0 处不可微，它的次梯度为 $[-1, 1]$ 之间的任意值，即：
$$ \begin{equation} \begin{cases} \left(\dfrac{\partial \left\|x\right\|_1}{\partial x}\right)_i=-1,& x_i<0 \\ -1\le \left(\dfrac{\partial \left\|x\right\|_1}{\partial x}\right)_i\le 1,& x_i=0 \\ \left(\dfrac{\partial \left\|x\right\|_1}{\partial x}\right)_i=1,& x_i>0 \end{cases} \end{equation} $$
我选取的递减步长为 $\alpha^k=0.001\times \dfrac{1}{\sqrt{k}}$。
编写代码
alpha_0 = 0.001 # 步长初值
p = 0.01 # 正则化参数
epsilon = 1e-5 # 最大允许误差
x_k = np.zeros(x_dim)
x_k_old = np.zeros(x_dim)
plot_iteration = [] # 每步计算结果
k = 1 # 迭代次数
while k < 1e6:
alpha = alpha_0 / np.sqrt(k) # 递减步长
# 计算目标函数次梯度
l1_subgradient = np.zeros(x_dim) # L1范数的次梯度
for i in range(len(x_k)):
if x_k[i] != 0:
l1_subgradient[i] = np.sign(x_k[i])
else:
l1_subgradient[i] = np.random.uniform(-1, 1) # 随机取[-1, 1]内的值作为次梯度
subgradient = A.T @ (A@x_k - b) + p * l1_subgradient
x_k = x_k - alpha * subgradient
plot_iteration.append(x_k.copy()) # 记录每步计算结果
if np.linalg.norm(x_k - x_k_old) < epsilon:
break
else:
x_k_old = x_k.copy() # 深拷贝
k += 1
x_optm = x_k.copy() # 最优解作图：
plot_dist_real = [] # 每步结果与真值的距离
plot_dist_optm = [] # 每步结果与最优解的距离
for iter in plot_iteration:
plot_dist_real.append(np.linalg.norm(iter - x_real))
plot_dist_optm.append(np.linalg.norm(iter - x_optm))
# 作图
plt.title('Subgradient Method')
plt.xlabel('Iteration times')
plt.ylabel('Distance')
plt.plot(plot_dist_real, 'r', label='Distance from true value of x')
plt.plot(plot_dist_optm, 'b', label='Distance from the optimal solution')
plt.grid()
plt.legend()
plt.show()结果分析
| 正则化参数 $p$ | 图像 |
|---|---|
| 0.01 |  |
| 0.1 |  |
| 1 |  |
| 10 |  |
在次梯度法中，$p$ 对结果的影响似乎并不明显。将次梯度法的图像与邻近点梯度法、交替方向乘子法相比较，可以看出次梯度法的回归效果并不是很理想，尽管最优解能够收敛，然而最优解与真值的差距较大。
L1 正则化的作用
通过查阅资料以及简单的数学分析，发现 L1 正则化可以起到使参数更加系稀疏的作用。当正则化惩罚项 $p\left\|x\right\|_1$ 被加入到目标函数中后，要极小化原目标函数与惩罚项之和，会受到该惩罚项的约束。当 $p$ 较大时，为了使目标函数尽可能小，需要使惩罚项尽可能小，而 $\left\|x\right\|_1=\sum_i \left|x_i \right|$，所以 $x$ 越稀疏，惩罚项的值就可以越小。
References
本文地址：https://www.jeddd.com/article/use-python-to-solve-lasso-by-pgd-admm-and-subgd.html
博主你好，在ADMM算法里的软阈值那里，关于最后的argminx的值是否少乘了一项alpha？x(k+1/2)+alpha*p?
博主您好，subgradient算法那里，你好像打错了。应该是 if x_k[i] !=0 : l1_sub[i] = sgn(x_k[i])
感谢博主
请问次梯度代码那里的第17行是不是应该是 np.sign(x_k[i])?
感谢指正
感谢博主，不过链接好像挂了
已修复
博主，已下载，很好用，谢谢
博主，链接挂了啊 OωO
没有啊，我这边下载正常 你再试一下吧