简述
凸优化会很详细地讲解这三个算法,这个学期刚好有这门课。
这里以期末的大作业的项目中的一个题目作为讲解。
题目
考虑线性测量b=Ax+e,其中b为50维的测量值,A为50*100维的测量矩阵,x为100维的未知稀疏向量且稀疏度为5,e为50维的测量噪声。从b和A中恢复x的一范数规范化最小二乘法模型(任务!!)
- p为非负的正则化参数。
- 提示:
在本实验中,设x的真值中的非零元素 服从标准正态分布,A中的元素服从标准正态分布,e中的元素服从均值为0,方差为0.1的高斯分布。
要求使用的算法:
- 邻近点梯度下降法
- 交替方向乘子法
- 次梯度法
实验部分
生成数据
generate-data.py
- 保存数据到二进制文件中
import numpy as np
import random
ASize = (50, 100)
XSize = 100
A = np.random.normal(0, 1, ASize)
X = np.zeros(XSize)
e = np.random.normal(0, 0.1, 50)
XIndex = random.sample(list(range(XSize)), 5) # 5 稀疏度
for xi in XIndex:
X[xi] = np.random.randn()
b = np.dot(A, X) + e
np.save("A.npy", A)
np.save("X.npy", X)
np.save("b.npy", b)
邻近点梯度下降法
- s为光滑项
- r为非光滑项
算法过程:
解析:
这一问本质上就是Lasso的邻近点梯度下降问题。
代入之前的算法过程,得到下面的结果(相比于原来的Lasso简单的变下就好了~)
- 求解argmin的方法–软门限法
关于光滑的部分,直接求导,在不光滑的部分,就求次梯度。
然后关于每个分量部分不进行。可以参照代码中的片段,在中间就为0,在区间之外就是对应的一个正数~
实现领近点梯度法
import numpy as np
A = np.load('A.npy')
b = np.load('b.npy')
X = np.load('X.npy')
ASize = (50, 100)
BSize = 50
XSize = 100
alpha = 0.001
P_half = 0.01
Xk = np.zeros(XSize)
zero = np.zeros(XSize)
while True:
Xk_half = Xk - alpha * np.dot(A.T, np.dot(A, Xk) - b)
# 软门限算子
Xk_new = zero.copy()
for i in range(XSize):
if Xk_half[i] < - alpha * P_half:
Xk_new[i] = Xk_half[i] + alpha * P_half
elif Xk_half[i] > alpha * P_half:
Xk_new[i] = Xk_half[i] - alpha * P_half
if np.linalg.norm(Xk_new - Xk, ord=2) < 1e-5:
break
else:
Xk = Xk_new.copy()
print(Xk)
print(X)
附加上画图的代码
- 蓝色的是算法最优值的距离(最终的收敛点的距离)
- 黄色的是预测的真实值的距离(一开始生成的数据)
- 距离都用的是二范数
import matplotlib.pyplot as plt
import numpy as np
A = np.load('A.npy')
b = np.load('b.npy')
X = np.load('X.npy')
ASize = (50, 100)
BSize = 50
XSize = 100
alpha = 0.005
P_half = 0.01
Xk = np.zeros(XSize)
zero = np.zeros(XSize)
X_opt_dst_steps = []
X_dst_steps = []
while True:
Xk_half = Xk - alpha * np.dot(A.T, np.dot(A, Xk) - b)
# 软门限算子
Xk_new = zero.copy()
for i in range(XSize):
if Xk_half[i] < - alpha * P_half:
Xk_new[i] = Xk_half[i] + alpha * P_half
elif Xk_half[i] > alpha * P_half:
Xk_new[i] = Xk_half[i] - alpha * P_half
X_dst_steps.append(np.linalg.norm(Xk_new - X, ord=2))
X_opt_dst_steps.append(Xk_new)
if np.linalg.norm(Xk_new - Xk, ord=2) < 1e-5:
break
else:
Xk = Xk_new.copy()
print(Xk)
print(X)
X_opt = X_opt_dst_steps[-1]
for i, data in enumerate(X_opt_dst_steps):
X_opt_dst_steps[i] = np.linalg.norm(data - X_opt, ord=2)
plt.title("Distance")
plt.plot(X_opt_dst_steps, label='X-opt-distance')
plt.plot(X_dst_steps, label='X-real-distance')
plt.legend()
plt.show()
交替方向乘子法
算法过程:
第一步中有交叉项的话,可以类似的换成下面的方式
上面说了,这次的问题是一个Lasso问题。
为了使用ADMM算法,这里引入一个新的变元
所以一致性约束,就变成了
-
出自于 增广拉格朗日算法中所提出的增广拉格朗日函数(augmented lagrangian function)
也就是在拉格朗日函数的基础上,在加上一个二范数作为penalty。
- 其中 c是一个常数(c>0)
对于这道题目,具体为:
代入之后,有(做适当地化简),
当然,我们注意到,x的更新的话,可以直接给出结果(因为是光滑的),通过矩阵求导,我们可以得到。
而z的计算也可以再做一次变换,
这个就用软门限的方式来进行求解。
同样关于 的正负性做分类。然后把x和v作为一个整体看,就是跟之前的软门限一模一样的了~
import matplotlib.pyplot as plt
import numpy as np
A = np.load('A.npy')
b = np.load('b.npy')
X = np.load('X.npy')
ASize = (50, 100)
BSize = 50
XSize = 100
P_half = 0.01
c = 0.005
Xk = np.zeros(XSize)
Zk = np.zeros(XSize)
Vk = np.zeros(XSize)
X_opt_dst_steps = []
X_dst_steps = []
while True:
Xk_new = np.dot(
np.linalg.inv(np.dot(A.T, A) + c * np.eye(XSize, XSize)),
c*Zk + Vk + np.dot(A.T, b)
)
# 软门限算子
Zk_new = np.zeros(XSize)
for i in range(XSize):
if Xk_new[i] - Vk[i] / c < - P_half / c:
Zk_new[i] = Xk_new[i] - Vk[i] / c + P_half / c
elif Xk_new[i] - Vk[i] / c > P_half / c:
Zk_new[i] = Xk_new[i] - Vk[i] / c - P_half / c
Vk_new = Vk + c * (Zk_new - Xk_new)
# print(np.linalg.norm(Xk_new - Xk, ord=2))
X_dst_steps.append(np.linalg.norm(Xk_new - X, ord=2))
X_opt_dst_steps.append(Xk_new)
if np.linalg.norm(Xk_new - Xk, ord=2) < 1e-5:
break
else:
Xk = Xk_new.copy()
Zk = Zk_new.copy()
Vk = Vk_new.copy()
print(Xk)
print(X)
X_opt = X_opt_dst_steps[-1]
for i, data in enumerate(X_opt_dst_steps):
X_opt_dst_steps[i] = np.linalg.norm(data - X_opt, ord=2)
plt.title("Distance")
plt.plot(X_opt_dst_steps, label='X-opt-distance')
plt.plot(X_dst_steps, label='X-real-distance')
plt.legend()
plt.show()
生成的图
次梯度法
其中 为 的次梯度。
对于这个问题其实分成简单,对于一范数来说,其实次梯度就是我们之前说的软门限算法。
这里重新描述一下:
- x不为0的情况下,为符号函数
- x为0的情况下,为[-1, 1]上的任意数
然后前半部分,就是用正常的梯度了。
实现 + 效果
import matplotlib.pyplot as plt
import numpy as np
A = np.load('A.npy')
b = np.load('b.npy')
X = np.load('X.npy')
def g_right(x):
Xnew = x.copy()
for i, data in enumerate(x):
if data == 0:
Xnew[i] = 2 * np.random.random() - 1
else:
Xnew[i] = np.sign(x[i])
return Xnew
ASize = (50, 100)
BSize = 50
XSize = 100
alpha = 0.001
p_half = 0.001
alphak = alpha
i = 0
g = lambda x: 2 * np.dot(A.T, (np.dot(A, x) - b)) + p_half * g_right(x)
Xk = np.zeros(XSize)
X_opt_dst_steps = []
X_dst_steps = []
while True:
Xk_new = Xk - alphak * g(Xk)
alphak = alpha / (i + 1)
i += 1
X_dst_steps.append(np.linalg.norm(Xk_new - X, ord=2))
X_opt_dst_steps.append(Xk_new)
print(np.linalg.norm(Xk_new - Xk, ord=2))
if np.linalg.norm(Xk_new - Xk, ord=2) < 1e-5:
break
else:
Xk = Xk_new.copy()
print(Xk)
print(X)
X_opt = X_opt_dst_steps[-1]
for i, data in enumerate(X_opt_dst_steps):
X_opt_dst_steps[i] = np.linalg.norm(data - X_opt, ord=2)
plt.title("Distance")
plt.plot(X_opt_dst_steps, label='X-opt-distance')
plt.plot(X_dst_steps, label='X-real-distance')
plt.legend()
plt.show()