Lasso原理与Python实现

Lasso原理与Python实现

• lasso原理

• Python实现

Lasso 原理
    在一般的回归情形中，对于数据集$\left({\mathrm{x}}^i,y^i\right)$,最小二乘法估计通过计算预测值与真实值之间的最小平方和对数据进行预测，但是数据分析不常满足最小二乘估计，而提升最小二乘估计的两个标准技术：子集选择与岭回归都存在缺陷。子集选择通过从特征中选择部分特征的方法建立模型，但数据中很小的变化就能导致不一样的模型，因此子集选择降低了预测的准确度。而岭回归不把任何系数设置为0，并不提供易于解释的模型。
    针对以上问题，作者提出了Lasso方法，即“最小绝对值收敛和选择算子”，从本质上说，Lasso方法能够收缩系数并将一些系数设置为0，保留了子集选择和岭回归的优点。
    Lasso的公式为：
        ${\min }_{\beta }\frac{1}{2}{\Vert \mathbf{y}-{\sum }_{i=1}^n{\mathbf{x}}_i{\beta }_i\Vert }_2^2+\lambda \Vert \beta {\Vert }_1$
    其中，前一项为最小二乘估计公式，后一项实际为正则化项L1。为计算使得误差最小的beta值，对公式做如下推导：
    $f=RSS(w)+\lambda {\Vert w_1\Vert }_1$
    $={\sum }_{i=0}^N{\left(y_i-{\sum }_{j=0}^D{w}_j{x}_{ij}\right)}^2+\lambda {\sum }_{j=0}^D\left|w_j\right|$

    $\frac{\partial RSS}{\partial W_j}$ = $-2{\sum }_{i=0}^N{x}_{ij}\left(y_i-{\sum }_{j=0}^D{w}_j{x}_{ij}\right)$

    $=-2{\sum }_{i=0}^N{x}_{ij}$$\left(y_i-{\sum }_{k\ne j}{w}_k{x}_{ik}-W_j{x}_{ij}\right)$

    $=-2{\sum }_{i=0}^N{x}_{ij}$$\left(y_i-{\sum }_{k\ne j}{w}_k{x}_{ik}\right)$+2${\sum }_{i=0}^N{w}_j{x}_{ij}^2$
    令$p_j=$${\sum }_{i=0}^N{x}_{ij}$$\left(y_i-{\sum }_{k\ne j}{w}_k{x}_{ik}\right)$

    $z_j=$${\sum }_{i=0}^N{x}_{ij}^2$

    则有$\frac{\partial RSS}{\partial W_j}=-2p_j+2w_j{z}_j$
    此外，对于正则项的求导本质上即为y=|x|的形式，求导结果为：
    $$\{\begin{array}{c}\lambda ,w_j>0\\ [-\lambda ,\lambda ],w_j=0\\ -\lambda ,w_j<0\end{array}$$
    因此当$w_j>0$时，$\frac{\partial f}{\partial w}=2z_j{w}_j-2p_j-\lambda$ ，令其为0解得：$w_j=\left(p_j-\frac{\lambda }{2}\right)/z_j$，且$p_j>\frac{\lambda }{2}$。同理，当$w_j<0$可解得$w_j=\left(p_j+\frac{\lambda }{2}\right)/z_j$，且$p_j<-\frac{\lambda }{2}$。
    根据解得的$w_j$公式，再计算对应的$p_j$、$z_j$，即可得到lasso回归的权重矩阵。
    具体实现时，$z_j$即为向量的内积，只需取出矩阵中的该列，并以内积的形式相乘即可得到$z_j$。计算$p_j$时，需设置双重循环。
Python实现

import numpy as np
# 计算残差
def rss(X,y,w):
	return (y-X*w).T*(y-X*w)
	
# lasso 回归	
def lassoRegression(X,y,a):
	threshold = 0.1
	# init para
	m,n = X.shape

	w = np.matrix(np.zeros((n,1)))
	r = rss(X,y,w)
	z = [0 for i in range(n)]

	tempw = [0 for i in range(n)]
	# 记录迭代次数
	iteration = 0
	# 迭代以计算参数
	while True:
		iteration += 1
		for k in range(n):
			temp_X = X[:, k].reshape(-1,1)
			z[k] = (temp_X.T*temp_X)[0, 0]
			p = [0 for i in range(n)]
			p_part = [0 for i in range(n)]
			for i in range(m):
				for j in range(n):
					if k != j:
						# print w[j,0]
						p_part[j] = X[i, j]*w[j, 0]
				p[k] += X[i, k] * (y[i, 0] - sum(p_part))
			if p[k] < -a/2:
				tempw[k] = (p[k] + a/2)/z[k]
			elif p[k] > a/2:
				tempw[k] = (p[k] - a/2)/z[k]
			else:
				tempw[k] = 0
			w[k, 0] = tempw[k]
		tempr = rss(X, y, w)
		delta = abs(tempr - r)[0, 0]
		r = tempr

		if delta < threshold:
			print 'iteration', iteration
			print 'delta', delta
			break
	return w
	
def load_data(filename):
	X, Y = [], []
	with open(filename, 'r') as f:
		for line in f:
			splited_line = [float(i) for i in line.split()]
			x, y = splited_line[: -1], splited_line[-1]
			X.append(x)
			Y.append(y)
	X, Y = np.matrix(X), np.matrix(Y).T
	return X, Y
def std(X):
	std_deviation = np.std(X, 0)
	mean = np.mean(X, 0)
	return (X - mean)/std_deviation
if '__main__'==__name__:
	# 使用鲍鱼数据集
	X, y = load_data('abalone.txt')
	X, y = std(X), std(y)
	enum = 30
    w_all = lassoTrace(X, y, enum)

    lambdas = [i-10 for i in range(enum)]

    data = pd.DataFrame({
        'x1': w_all[:, 0],
        'x2': w_all[:, 1],
        'x3': w_all[:, 2],
        'x4': w_all[:, 3],
        'x5': w_all[:, 4],
        'x6': w_all[:, 5],
        'x7': w_all[:, 6],
        'x8': w_all[:, 7],
    })
    plt.figure(12)

    plt.subplot(221)
    data.boxplot()

    plt.subplot(222)
    plt.plot(lambdas, w_all)
    plt.show()

输出各系数随惩罚力度变化的曲线图以及各系数的箱形图。