Lasso回归算法

84 阅读2分钟

携手创作,共同成长!这是我参与「掘金日新计划 · 8 月更文挑战」的第1天,点击查看活动详情

我们把lasso回归损失函数分为两部分,然后依次求导:

第一部分:

图片.png

图片.png

其中:

图片.png
第二部分求导:

图片.png

两部分合并得:

图片.png

 

令导数等于0得:

图片.png

注: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