Extend Coordinate Descent for LASSO to Adaptive LASSO

62 views Asked by At

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?

1

There are 1 answers

1
Riccardo Bucco On BEST ANSWER

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.