2022 11月24 Ridge/LASSO Regression学习笔记
创始人
2024-02-12 09:33:16
  • 解决过拟合现象,减少高次项的影响,使曲线更加平滑。利用正则化。岭回归和LASSO都是一种正则化。

  • 岭回归是将代价函数正则化

  • LASSO回归是将高价的项正则化,让他们的影响不那么大。

  • 岭回归和LASSO回归是最小二乘法的优化,解决了最小二乘法的局限性

岭回归

  • 局限性: w=(X(T)X)(-1)X(T)y,如果 X(T)X 是奇异矩阵则无法进行求解

    • w=(X(T)X)(-1)X(T)y :在高维 x 的情况下,满足这个条件的w使得残差最小(和真实值的距离最小)
    • 两种情况
      1. X本身存在相关性 行列式为0 奇异矩阵
      2. m
  • 中心化:所有数据平移到原点

  • 标准化:让所有数据不同特征有相同尺度,数据减方差除标准差(类似正态变量标准化)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9y0K4cqS-1669271227239)(Ridge:LASSO Regression.assets/image-20221121205756485.png)]

  • 岭回归的作用:减少最小二乘法可能得过拟合现象,或者,引入惩罚因子。
    • w=(X(T)X+λI)-1XTy
    • X(T)X不可逆,则它们的行列式为0,则它们的特征值相乘等0。通过给他们加上一个正的单位向量解决这个问题。

岭回归主要是对传统的最小二乘法拟合曲线矩阵奇异的处理,通过引入乘法因子,这个算法是我第一次对机器学习进行的实验,首先我对岭回归的来龙去脉进行了学习,之后网上寻找相关代码和数据,这个代码进行了三部分,第一部分是对数据的导入,第二部分是岭回归系数的计算,第三部分是画图。对numpy和matplotlib画图有了一点点认识。对于算法的每部分含义写了注释,对于岭矩阵然后画图还是有点不太理解再做什么。

import random
from math import expimport numpy as np
from matplotlib import pyplot as pltdef loadDataSet(fileName):"""加载数据:param fileName: 文件名:return:xArr - x数据集yArr - y数据集"""# open函数主要运用到两个参数,文件名和mode,文件名是添加该文件对象的变量,# 返回文件对象# mode是告诉编译器和开发者文件通过怎样的方式进行使用。因此在Python中打开文件的代码如下:# mode参数可以不写,默认mode参数是“r”。mode参数如下:# ‘r’ – 只读模式,当文件处在”只读“的模式时使用。# readline() 读取一行numFeat = len(open(fileName).readline().split('\t')) - 1xArr = []yArr = []f = open(fileName)for line in f.readlines():lineArr = []# strip()表示删除掉数据中的换行符curLine = line.strip().split('\t')for i in range(numFeat):lineArr.append(float(curLine[i]))xArr.append(lineArr)yArr.append(float(curLine[-1]))x = np.array(xArr)y = np.array(yArr)return x,ydef ridgeRegres(xMat, yMat, lam=0.2):"""岭回归:param xMat: x数据集:param yMat: y数据集:param lam: 缩减系数:return:ws - 回归系数"""xTx = xMat.T.dot(xMat)  # np.dot(xMat) 矩阵乘法# 构造单位矩阵 np.eye(n) n为列数# np.shape(xMat)[1] 得到xMat矩阵的列数denom = xTx + np.eye(np.shape(xMat)[1]) * lamif np.linalg.det(denom) == 0.0:print("矩阵为奇异矩阵,不能转置")returnws = np.linalg.inv(denom).dot(xMat.T.dot(yMat))  # np.linalg.inv(denom) 矩阵求逆return wsdef ridge_traj(X, y, ntest=30):''' 获取岭轨迹矩阵'''_, n = X.shapews = np.zeros((ntest, n))X = np.array(X)y = np.array(y)for i in range(ntest):w = ridgeRegres(X, y, exp(i-10))ws[i, :] = w.Treturn wsif '__main__' == __name__:ntest = 30# 加载数据X,y = loadDataSet("data.txt")# 绘制岭轨迹ws = ridge_traj(X, y, ntest)# 创建自定义子图fig = plt.figure()#大致意思就是输入 三位数字,abc,根据abc将原图划分成a行,b列。# 那么就会将原图划分成a*b个子图,然后c就是我们的下标。# 我们通过c来指定展示我们要的子图。ax = fig.add_subplot(111)lambdas = [i-10 for i in range(ntest)]ax.plot(lambdas, ws)    # 确定横纵坐标plt.show()

数据下载地址

LASSO回归

​ LASSO回归解决如果需要拟合的数据特征值远大于样本值,或者数据特征之间的相关性太大,利用LASSO回归正则化特征值,给特征加上权重,使得多重共线性问题不会那么明显。但是LASSO回归相比于岭回归,是利用特征值的绝对值乘以λ,所以不可以利用求导的方式来求代价函数的最小值,所以这里介绍两种方式一种是坐标轴下降法,另一种是最小角回归法。坐标轴下降法是固定一个参数系数对其余的参数进行迭代然后求出这个参数情况下代价函数的最小值,然后对每个参数都算一次就能求出总体的最小值。最小角回归法结合了向前选择算法和前向梯度算法,算法思路是找到和目标y最接近的自变量x然后走一大步直到和倒数第二接近的自变量x有相同的残差,再沿着残差走,循环以上思路直到最后获得最小残差。

下面是网上的代码,没有细看。同时使用这个手写的代码和库里的代码,结果不同。手写的代码还是会受到异常值的影响。

def lassoUseCd(X, y, lambdas=0.1, max_iter=1000, tol=1e-4):"""Lasso回归,使用坐标下降法(coordinate descent)args:X - 训练数据集y - 目标标签值lambdas - 惩罚项系数max_iter - 最大迭代次数tol - 变化量容忍值return:w - 权重系数"""# 初始化 w 为零向量w = np.zeros(X.shape[1])for it in range(max_iter):done = True# 遍历所有自变量for i in range(0, len(w)):# 记录上一轮系数weight = w[i]# 求出当前条件下的最佳系数w[i] = down(X, y, w, i, lambdas)# 当其中一个系数变化量未到达其容忍值,继续循环if (np.abs(weight - w[i]) > tol):done = False# 所有系数都变化不大时,结束循环if (done):breakreturn wdef down(X, y, w, index, lambdas=0.1):"""cost(w) = (x1 * w1 + x2 * w2 + ... - y)^2 + ... + λ(|w1| + |w2| + ...)假设 w1 是变量,这时其他的值均为常数,带入上式后,其代价函数是关于 w1 的一元二次函数,可以写成下式:cost(w1) = (a * w1 + b)^2 + ... + λ|w1| + c (a,b,c,λ 均为常数)=> 展开后cost(w1) = aa * w1^2 + 2ab * w1 + λ|w1| + c (aa,ab,c,λ 均为常数)"""# 展开后的二次项的系数之和aa = 0# 展开后的一次项的系数之和ab = 0for i in range(X.shape[0]):# 括号内一次项的系数a = X[i][index]# 括号内常数项的系数b = X[i][:].dot(w) - a * w[index] - y[i]# 可以很容易的得到展开后的二次项的系数为括号内一次项的系数平方的和aa = aa + a * a# 可以很容易的得到展开后的一次项的系数为括号内一次项的系数乘以括号内常数项的和ab = ab + a * b# 由于是一元二次函数,当导数为零时,函数值最小值,只需要关注二次项系数、一次项系数和 λreturn det(aa, ab, lambdas)def det(aa, ab, lambdas=0.1):"""通过代价函数的导数求 w,当 w = 0 时,不可导det(w) = 2aa * w + 2ab + λ = 0 (w > 0)=> w = - (2 * ab + λ) / (2 * aa)det(w) = 2aa * w + 2ab - λ = 0 (w < 0)=> w = - (2 * ab - λ) / (2 * aa)det(w) = NaN (w = 0)=> w = 0"""w = - (2 * ab + lambdas) / (2 * aa)if w < 0:w = - (2 * ab - lambdas) / (2 * aa)if w > 0:w = 0return w
def lassoUseLars(X, y, lambdas=0.1, max_iter=1000):"""Lasso回归,使用最小角回归法(Least Angle Regression)论文:https://web.stanford.edu/~hastie/Papers/LARS/LeastAngle_2002.pdfargs:X - 训练数据集y - 目标标签值lambdas - 惩罚项系数max_iter - 最大迭代次数return:w - 权重系数"""n, m = X.shape# 已被选择的特征下标active_set = set()# 当前预测向量cur_pred = np.zeros((n,), dtype=np.float32)# 残差向量residual = y - cur_pred# 特征向量与残差向量的点积,即相关性cur_corr = X.T.dot(residual)# 选取相关性最大的下标j = np.argmax(np.abs(cur_corr), 0)# 将下标添加至已被选择的特征下标集合active_set.add(j)# 初始化权重系数w = np.zeros((m,), dtype=np.float32)# 记录上一次的权重系数prev_w = np.zeros((m,), dtype=np.float32)# 记录特征更新方向sign = np.zeros((m,), dtype=np.int32)sign[j] = 1# 平均相关性lambda_hat = None# 记录上一次平均相关性prev_lambda_hat = Nonefor it in range(max_iter):# 计算残差向量residual = y - cur_pred# 特征向量与残差向量的点积cur_corr = X.T.dot(residual)# 当前相关性最大值largest_abs_correlation = np.abs(cur_corr).max()# 计算当前平均相关性lambda_hat = largest_abs_correlation / n# 当平均相关性小于λ时,提前结束迭代# https://github.com/scikit-learn/scikit-learn/blob/2beed55847ee70d363bdbfe14ee4401438fba057/sklearn/linear_model/_least_angle.py#L542if lambda_hat <= lambdas:if (it > 0 and lambda_hat != lambdas):ss = ((prev_lambda_hat - lambdas) / (prev_lambda_hat - lambda_hat))# 重新计算权重系数w[:] = prev_w + ss * (w - prev_w)break# 更新上一次平均相关性prev_lambda_hat = lambda_hat# 当全部特征都被选择,结束迭代if len(active_set) > m:break# 选中的特征向量X_a = X[:, list(active_set)]# 论文中 X_a 的计算公式 - (2.4)X_a *= sign[list(active_set)]# 论文中 G_a 的计算公式 - (2.5)G_a = X_a.T.dot(X_a)G_a_inv = np.linalg.inv(G_a)G_a_inv_red_cols = np.sum(G_a_inv, 1)     # 论文中 A_a 的计算公式 - (2.5)A_a = 1 / np.sqrt(np.sum(G_a_inv_red_cols))# 论文中 ω 的计算公式 - (2.6)omega = A_a * G_a_inv_red_cols# 论文中角平分向量的计算公式 - (2.6)equiangular = X_a.dot(omega)# 论文中 a 的计算公式 - (2.11)cos_angle = X.T.dot(equiangular)# 论文中的 γgamma = None# 下一个选择的特征下标next_j = None# 下一个特征的方向next_sign = 0for j in range(m):if j in active_set:continue# 论文中 γ 的计算方法 - (2.13)v0 = (largest_abs_correlation - cur_corr[j]) / (A_a - cos_angle[j]).item()v1 = (largest_abs_correlation + cur_corr[j]) / (A_a + cos_angle[j]).item()if v0 > 0 and (gamma is None or v0 < gamma):gamma = v0next_j = jnext_sign = 1if v1 > 0 and (gamma is None or v1 < gamma):gamma = v1next_j = jnext_sign = -1if gamma is None:# 论文中 γ 的计算方法 - (2.21)gamma = largest_abs_correlation / A_a# 选中的特征向量sa = X_a# 角平分向量sb = equiangular * gamma# 解线性方程(sa * sx = sb)sx = np.linalg.lstsq(sa, sb)# 记录上一次的权重系数prev_w = w.copy()d_hat = np.zeros((m,), dtype=np.int32)for i, j in enumerate(active_set):# 更新当前的权重系数w[j] += sx[0][i] * sign[j]# 论文中 d_hat 的计算方法 - (3.3)d_hat[j] = omega[i] * sign[j]# 论文中 γ_j 的计算方法 - (3.4)gamma_hat = -w / d_hat# 论文中 γ_hat 的计算方法 - (3.5)gamma_hat_min = float("+inf")# 论文中 γ_hat 的下标gamma_hat_min_idx = Nonefor i, j in enumerate(gamma_hat):if j <= 0:continueif gamma_hat_min > j:gamma_hat_min = jgamma_hat_min_idx = iif gamma_hat_min < gamma:# 更新当前预测向量 - (3.6)cur_pred = cur_pred + gamma_hat_min * equiangular# 将下标移除至已被选择的特征下标集合active_set.remove(gamma_hat_min_idx)# 更新特征更新方向集合sign[gamma_hat_min_idx] = 0else:# 更新当前预测向量cur_pred = X.dot(w)# 将下标添加至已被选择的特征下标集合active_set.add(next_j)# 更新特征更新方向集合sign[next_j] = next_signreturn w
  • 利用scikit-learn 实现
from sklearn.linear_model import Lasso# 初始化Lasso回归器,默认使用坐标下降法
reg = Lasso(alpha=0.1, fit_intercept=False)
# 拟合线性模型
reg.fit(X, y)
# 权重系数
w = reg.coef_
from sklearn.linear_model import LassoLars# 初始化Lasso回归器,使用最小角回归法
reg = LassoLars(alpha=0.1, fit_intercept=False)
# 拟合线性模型
reg.fit(X, y)
# 权重系数
w = reg.coef_

梯度下降法

#临时写的函数,要在引入一个copy包,进行深度拷贝
#大家写一份代码,把要引入的包全放在最前面
import copydef CoordinateDescent(x, y,epochs,learning_rate,Lambda):        m=x.shape[0]X = np.concatenate((np.ones((m,1)),x),axis=1)xMat=np.mat(X)yMat =np.mat(y.reshape(-1,1))w = np.ones(X.shape[1]).reshape(-1,1)for n in range(epochs):out_w = copy.copy(w)for i,item in enumerate(w):#在每一个W值上找到使损失函数收敛的点for j in range(epochs):h = xMat * w gradient = xMat[:,i].T * (h - yMat)/m + Lambda * np.sign(w[i])w[i] = w[i] - gradient* learning_rateif abs(gradient)<1e-3:breakout_w = np.array(list(map(lambda x:abs(x)<1e-3, out_w-w)))if out_w.all():breakreturn  w

利用这个代码调整Lambda来观察拟合曲线

#Lambda=100时;
w = CoordinateDescent(x_1,y_1,epochs=250,learning_rate=0.001,Lambda=100)
print(w)#计算新的拟合值
y_1_pred = x_1 * w[1] + w[0]ax1= plt.subplot()
ax1.scatter(x_1,y_1,label='样本分布')
ax1.plot(x,y_pred,c='y',label='原始样本拟合')
ax1.plot(x_1,y_1_pred,c='r',label='新样本拟合')
ax1.legend(prop = {'size':15}) #此参数改变标签字号的大小
plt.show()
正则项参数过大、过小都不好,
过小起不到惩罚效果,模型任然过拟合;
过大惩罚太大,会使得模型欠拟合,达不到要求;
我们选择参数的标准:模型在训练集、验证集、测试集上,评估效果接近时,这个正则项参数较好。
  1. 参考资料

  2. 参考资料

  3. 参考资料

  4. 参考资料

  5. 参考资料

相关内容

热门资讯

埃菲尔铁塔在哪 中国仿建埃菲尔... 2019年4月26日,广西南宁市,街头惊现一座巨型山寨版埃菲尔铁塔,高约20米,白色塔身,造型逼真,...
苗族的传统节日 贵州苗族节日有... 【岜沙苗族芦笙节】岜沙,苗语叫“分送”,距从江县城7.5公里,是世界上最崇拜树木并以树为神的枪手部落...
北京的名胜古迹 北京最著名的景... 北京从元代开始,逐渐走上帝国首都的道路,先是成为大辽朝五大首都之一的南京城,随着金灭辽,金代从海陵王...
长白山自助游攻略 吉林长白山游... 昨天介绍了西坡的景点详细请看链接:一个人的旅行,据说能看到长白山天池全凭运气,您的运气如何?今日介绍...
应用未安装解决办法 平板应用未... ---IT小技术,每天Get一个小技能!一、前言描述苹果IPad2居然不能安装怎么办?与此IPad不...
脚上的穴位图 脚面经络图对应的... 人体穴位作用图解大全更清晰直观的标注了各个人体穴位的作用,包括头部穴位图、胸部穴位图、背部穴位图、胳...
demo什么意思 demo版本... 618快到了,各位的小金库大概也在准备开闸放水了吧。没有小金库的,也该向老婆撒娇卖萌服个软了,一切只...
猫咪吃了塑料袋怎么办 猫咪误食... 你知道吗?塑料袋放久了会长猫哦!要说猫咪对塑料袋的喜爱程度完完全全可以媲美纸箱家里只要一有塑料袋的响...
世界上最漂亮的人 世界上最漂亮... 此前在某网上,选出了全球265万颜值姣好的女性。从这些数量庞大的女性群体中,人们投票选出了心目中最美...
埃菲尔铁塔在哪 中国仿建埃菲尔... 2019年4月26日,广西南宁市,街头惊现一座巨型山寨版埃菲尔铁塔,高约20米,白色塔身,造型逼真,...
苗族的传统节日 贵州苗族节日有... 【岜沙苗族芦笙节】岜沙,苗语叫“分送”,距从江县城7.5公里,是世界上最崇拜树木并以树为神的枪手部落...
北京的名胜古迹 北京最著名的景... 北京从元代开始,逐渐走上帝国首都的道路,先是成为大辽朝五大首都之一的南京城,随着金灭辽,金代从海陵王...
应用未安装解决办法 平板应用未... ---IT小技术,每天Get一个小技能!一、前言描述苹果IPad2居然不能安装怎么办?与此IPad不...