pythonlasso-regression

Extend Coordinate Descent for LASSO to Adaptive LASSO


I currently have a simple coordinate descent algorithm to solve LASSO in python:

def lasso_nb(X, y, alpha, tol=0.001, maxiter=10000):
    n, p = X.shape
    beta = np.zeros(p)
    R = y.copy()
    norm_cols_X = (X ** 2).sum(axis=0)
    resids = []
    prev_cost = 10e10
    for n_iter in range(maxiter):
        for ii in range(p):
            beta_ii = beta[ii]
            if beta_ii != 0.:
                R += X[:, ii] * beta_ii
            tmp = np.dot(X[:, ii], R)
            beta[ii] = fsign(tmp) * max(abs(tmp) - alpha, 0) / (.00001 + norm_cols_X[ii])
            if beta[ii] != 0.:
                R -= X[:, ii] * beta[ii]
        cost = (np.sum((y - X @ beta)**2) + alpha * np.sum(np.abs(beta))) / n
        resids.append(cost)
        if prev_cost - cost < tol:
            break
        else:
            prev_cost = cost
    return beta

how can this be extended to do Adaptive LASSO? Or is it not possible with coordinate descent?


Solution

  • Here is a possible solution (I commented the lines of code that I added/edited):

    def adaptive_lasso_cd(X, y, alpha, tol=0.001, maxiter=10000):
        n, p = X.shape
        beta = np.zeros(p)
        w = np.ones(p)  # Initialize weights for adaptive LASSO
        R = y.copy()
        norm_cols_X = (X ** 2).sum(axis=0)
        resids = []
        prev_cost = 10e10
        for n_iter in range(maxiter):
            for ii in range(p):
                beta_ii = beta[ii]
                if beta_ii != 0.:
                    R += X[:, ii] * beta_ii
                tmp = np.dot(X[:, ii], R)
                z = np.abs(tmp)
                w[ii] = 1.0 / (z + 1e-10)  # Update weights based on current estimate
                beta[ii] = np.sign(tmp) * max(z - alpha, 0) / (.00001 + norm_cols_X[ii]) * w[ii]  # Use updated weights
                if beta[ii] != 0.:
                    R -= X[:, ii] * beta[ii]
            cost = (np.sum((y - X @ beta)**2) + alpha * np.sum(np.abs(beta))) / n
            resids.append(cost)
            if prev_cost - cost < tol:
                break
            else:
                prev_cost = cost
        return beta
    

    Basically, all you need to do is updating the regularization parameter with adaptive weights based on the current estimate of the coefficients.