携手创作,共同成长!这是我参与「掘金日新计划 · 8 月更文挑战」的第1天,点击查看活动详情
我们把lasso回归损失函数分为两部分,然后依次求导:
第一部分:
其中:
第二部分求导:
两部分合并得:
令导数等于0得:
注:Lagrange方程的求导等于0,等价于求出差方和在附加的指定区域内的极值。等式的右边是典型的坐标下降法。也就是通过坐标下降法来更新回归系数。
编程思想:
1. 通过上面的公式,每次选取一维进行优化,优化出一维,然后再循环到下一维,直到全部优化成功。
2. 然后优化出的回归系数矢量代入到损失函数计算,与初始损失函数值比较,如果差别没有达到阈值,继续循环计算。
3. 优化出新的回归系数,再次代入计算出损失函数值,与上一次的损失函数值作比较,如果绝对差别小于阈值,循环结束,返回矢量回归系数值。
import itertools
from math import exp
import numpy as np
import matplotlib.pyplot as plt
def lasso_regression(X_array, y, lambd, threshold=0.1):
#通过坐标下降(coordinate descent)法获取LASSO回归系数
# 计算残差平方和
X=np.column_stack((np.ones( X_array.shape[0]),X_array))
square_diff_sum = lambda X, y, w: np.dot((np.dot(X,w)-y).T,(np.dot(X,w)-y))
#数据的初始化,参数是0数组,差方的初始化通过w的0数组来实现的
m, n = X.shape
w = np.zeros(n)#创造承载器
r = square_diff_sum(X, y, w)
# 迭代器
itertoll = itertools.count(1)
for it in itertoll:
for k in range(n):
#print(k)
z_k = np.dot(X[:,k].T,X[:,k])
#print(z_k)
p_k = 0
for i in range(m):
p_k += X[i,k]*(y[i] - sum([X[i,j]*w[j] for j in range(n) if j != k]))
#上面三行可直接改成矢量化
#p_k=X[:,k]@(y-(np.delete(X,k,1)@np.delete(w,k)))
#print(p_k,lambd/2)
if p_k < -lambd/2:
w_k = (p_k + lambd/2)/z_k
elif p_k > lambd/2:
w_k = (p_k - lambd/2)/z_k
else:
w_k = 0
w[k] = w_k
square_diff_sum_indirect = square_diff_sum(X, y, w)
error_abs = abs(square_diff_sum_indirect- r)#绝对差别
r = square_diff_sum_indirect
if error_abs< threshold:
break
return w
Proximal Gradient descent