通过前面的介绍，我们了解到求解惩罚性回归模型等价于求解一个优化问题，有大量的通用算法可以求解。这次介绍一个非常快速敏捷的算法——最小角回归（LARS）算法。LARS算法可以理解为一种改进的前向逐步回归算法。
前向逐步回归算法和最小角回归算法（LARS）：
前向逐步回归算法过程如下：
将回归系数的所有值初始化为0
在每一步使用已经选择的变量找到残差值
确定哪个未使用的变量能够最佳的解释残差，将该变量加入选择变量中。
LARS算法与前向逐步回归算法类似，LARS与前向逐步回归算法的主要差异是LARS在引入新属性时只是部分引入，引入属性过程并非不可逆。LARS算法过程如下：
将回归系数的所有值初始化为0
在每一步中，决定哪个属性与残差有最大关联
如果关联为正，小幅度增加关联系数；关联未负，小幅度减少关联系数；
下面以红酒数据集为例，介绍一下LARS算法：
import urllib.request, urllib.error, urllib.parseimport numpyfrom sklearn import datasets, linear_modelfrom math import sqrtimport matplotlib.pyplot as plot#读取数据target_url = "http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv"data = urllib.request.urlopen(target_url)xList = []labels = []names = []firstLine = Truefor line in data:if firstLine:names = line.decode().strip().split(";")firstLine = Falseelse:#分号分割row = line.decode().strip().split(";")#读取标签labels.append(float(row[-1]))#移除标签row.pop()#特征转为浮点型floatRow = [float(num) for num in row]xList.append(floatRow)#读取数据规模nrows = len(xList)ncols = len(xList[0])#计算均值和标准差xMeans = []xSD = []for i in range(ncols):col = [xList[j][i] for j in range(nrows)]mean = sum(col)/nrowsxMeans.append(mean)colDiff = [(xList[j][i] - mean) for j in range(nrows)]sumSq = sum([colDiff[i] * colDiff[i] for i in range(nrows)])stdDev = sqrt(sumSq/nrows)xSD.append(stdDev)#数据标准化xNormalized = []for i in range(nrows):rowNormalized = [(xList[i][j] - xMeans[j])/xSD[j] for j in range(ncols)]xNormalized.append(rowNormalized)#标签数据标准化meanLabel = sum(labels)/nrowssdLabel = sqrt(sum([(labels[i] - meanLabel) * (labels[i] - meanLabel) for i in range(nrows)])/nrows)labelNormalized = [(labels[i] - meanLabel)/sdLabel for i in range(nrows)]#系数初始化beta = [0.0] * ncols#初始化迭代系数矩阵betaMat = []betaMat.append(list(beta))#迭代参数nSteps = 350stepSize = 0.004nzList = []for i in range(nSteps):#计算残差residuals = [0.0] * nrowsfor j in range(nrows):labelsHat = sum([xNormalized[j][k] * beta[k] for k in range(ncols)])residuals[j] = labelNormalized[j] - labelsHat#计算特征与残差的相关性corr = [0.0] * ncolsfor j in range(ncols):corr[j] = sum([xNormalized[k][j] * residuals[k] for k in range(nrows)]) / nrowsiStar = 0corrStar = corr[0]for j in range(1, (ncols)):if abs(corrStar) < abs(corr[j]):iStar = j; corrStar = corr[j]beta[iStar] += stepSize * corrStar / abs(corrStar)betaMat.append(list(beta))nzBeta = [index for index in range(ncols) if beta[index] != 0.0]for q in nzBeta:if (q in nzList) == False:nzList.append(q)nameList = [names[nzList[i]] for i in range(len(nzList))]print(nameList)for i in range(ncols):#画迭代步数和残差相关性图coefCurve = [betaMat[k][i] for k in range(nSteps)]xaxis = range(nSteps)plot.plot(xaxis, coefCurve)plot.xlabel("Steps Taken")plot.ylabel(("Coefficient Values"))plot.show()
代码主要分为三部分：
第一部分读入数据以及列名称，将数据转换为属性列表以及标签。
第二部分对属性以及标签进行数据标准化（归一化）。
第三部分对问题进行求解获取结果。这里因为上面已经对数据进行了归一化，所以不需要再有截距项。截距项一般是用于表示标签值与加权属性值之间的差异。
LARS算法生成的结果如图所示，对图进行理解可以参考如下的方法：沿着迭代的方向想象有一个点，在该点上，一条垂直的线会穿过所有曲线，垂直线和系数曲线相交的值就是LARS算法在该步得到的系数。如果生成曲线需要350步，那么就生成350组系数集合。我们注意到对于前25步，只有一个系数值非零，第一个系数不为零的属性是酒精，在后此一段时间内，这是LARS回归唯一使用的变量。随着迭代步数的增多第二个变量开始出现，这样的过程一直持续到所有变量被包含到解中。
那么350个解的集合中，到底哪个解最优呢？正如前面章节中说的，模型的性能是指在测试集上的性能。前面也介绍过预留测试集的方法，本节中使用10折交叉验证法决定部署使用的最佳系数集合。代码如下：
import urllib.request, urllib.error, urllib.parseimport numpyfrom sklearn import datasets, linear_modelfrom math import sqrtimport matplotlib.pyplot as plot#读取数据target_url = "http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/winequality-red.csv"data = urllib.request.urlopen(target_url)xList = []labels = []names = []firstLine = Truefor line in data:if firstLine:names = line.decode().strip().split(";")firstLine = Falseelse:#分号分割row = line.decode().strip().split(";")#读取标签labels.append(float(row[-1]))#移除标签row.pop()#转为浮点型floatRow = [float(num) for num in row]xList.append(floatRow)#数据标准化nrows = len(xList)ncols = len(xList[0])#计算均值和标准差xMeans = []xSD = []for i in range(ncols):col = [xList[j][i] for j in range(nrows)]mean = sum(col)/nrowsxMeans.append(mean)colDiff = [(xList[j][i] - mean) for j in range(nrows)]sumSq = sum([colDiff[i] * colDiff[i] for i in range(nrows)])stdDev = sqrt(sumSq/nrows)xSD.append(stdDev)#特征数据标准化计算xNormalized = []for i in range(nrows):rowNormalized = [(xList[i][j] - xMeans[j])/xSD[j] for j in range(ncols)]xNormalized.append(rowNormalized)#标签数据标准化meanLabel = sum(labels)/nrowssdLabel = sqrt(sum([(labels[i] - meanLabel) * (labels[i] - meanLabel) for i in range(nrows)])/nrows)labelNormalized = [(labels[i] - meanLabel)/sdLabel for i in range(nrows)]#10折法nxval = 10#迭代参数设置nSteps = 350stepSize = 0.004#初始化errors = []for i in range(nSteps):b = []errors.append(b)for ixval in range(nxval):#设置训练集和测试集idxTest = [a for a in range(nrows) if a%nxval == ixval*nxval]idxTrain = [a for a in range(nrows) if a%nxval != ixval*nxval]#定义测试集和训练集xTrain = [xNormalized[r] for r in idxTrain]xTest = [xNormalized[r] for r in idxTest]labelTrain = [labelNormalized[r] for r in idxTrain]labelTest = [labelNormalized[r] for r in idxTest]#训练LARSnrowsTrain = len(idxTrain)nrowsTest = len(idxTest)#初始化beta系数beta = [0.0] * ncols#初始化系数矩阵betaMat = []betaMat.append(list(beta))for iStep in range(nSteps):#计算残差residuals = [0.0] * nrowsfor j in range(nrowsTrain):labelsHat = sum([xTrain[j][k] * beta[k] for k in range(ncols)])residuals[j] = labelTrain[j] - labelsHat#计算特征和残差相关性corr = [0.0] * ncolsfor j in range(ncols):corr[j] = sum([xTrain[k][j] * residuals[k] for k in range(nrowsTrain)]) / nrowsTrainiStar = 0corrStar = corr[0]for j in range(1, (ncols)):if abs(corrStar) < abs(corr[j]):iStar = j; corrStar = corr[j]beta[iStar] += stepSize * corrStar / abs(corrStar)betaMat.append(list(beta))#计算均方误差for j in range(nrowsTest):labelsHat = sum([xTest[j][k] * beta[k] for k in range(ncols)])err = labelTest[j] - labelsHaterrors[iStep].append(err)cvCurve = []for errVect in errors:mse = sum([x*x for x in errVect])/len(errVect)cvCurve.append(mse)minMse = min(cvCurve)minPt = [i for i in range(len(cvCurve)) if cvCurve[i] == minMse ][0]print("Minimum Mean Square Error", minMse)print("Index of Minimum Mean Square Error", minPt)xaxis = range(len(cvCurve))plot.plot(xaxis, cvCurve)plot.xlabel("Steps Taken")plot.ylabel(("Mean Square Error"))plot.show()
通过模型输出结果可以看出，最小值在低311步出现，所以选择311步的系数最为模型系数进行部署。
All things are difficult before they are easy.
如果觉得《【Python机器学习预测分析算法实战五】最小角回归LARS》对你有帮助，请点赞、收藏，并留下你的观点哦！
