Generalized Dantzig Selector

Question

I intended to write code for the generalized Dantzig selector based on the paper (https://www.jstor.org/stable/27798827, on Page 328), but I have encountered a problem. I believe the issue is related to an error indicating that the dimension of V is null. Initially, I started with the Dantzig selector and planned to extend it to a generalized version for Bernoulli response. I assumed the link function to be the logistic function and used its inverse to calculate $\mu$. Subsequently, I determined the variance based on $\mu$ and the distribution type.

I would greatly appreciate it if you could assist me with this code. I may have made a critical mistake. I've tried as the following

import cvxpy as cp
import numpy as np

def dantzig_selector(X, y, lmbda,W):
    n, p = X.shape
    beta = cp.Variable(p)


    objective = cp.Minimize(cp.norm(beta, 1))


    constraints = [cp.norm(cp.matmul(X[:, j].T, W @ (X @ beta - y)), 'inf') <= lmbda for j in range(p)]


    problem = cp.Problem(objective, constraints)
    problem.solve()
    


    selected_features = np.where(np.abs(beta.value) > 1e-2)[0]
    coefficients = beta.value[selected_features]

    return selected_features, coefficients,beta.value



def generalized_dantzig_selector(X, y, lmbda, max_iter=4, tol=1e-7):
    n, p = X.shape
    beta = cp.Variable(p)
    t = cp.Variable()
    Beta=np.zeros((max_iter,p))
    




    for i in range(max_iter):
        print("iter: ",i)

        mu = 1 / (1 + cp.exp(-X @ beta))  
        V = mu * (1 - mu) 
        W = cp.diag(V) 

        Z = X @ beta + (y - mu) / V

        print("Dimension of mu:", mu.shape)
        print("Dimension of V:", V.shape)
        print("Dimension of Z:", Z.shape)


        Z_star = cp.multiply(cp.sqrt(V), Z)
        X_star = cp.multiply(cp.sqrt(V), X)

        print("Dimension of Z_star:", Z_star.shape)
        print("Dimension of X_star:", X_star.shape)
        selected_features, coefficients, Beta[i]=dantzig_selector(X_star,Z_star,lmbda,W)
        
        X=X_star
        Z=Z_star

        if (i>0):
          if np.linalg.norm(Beta[i] - Beta[i-1]) < tol:
            break
        beta.value=Beta[i]

    return beta.value

n, p = 800, 20
X = np.random.randn(n, p)
true_beta = np.zeros(p)
true_beta[[3,1,2]]=0.8
true_beta[:int(np.ceil(p/5))] = 0.9

lmbda=0.01



random_index = np.random.randint(0, n)

y = np.zeros(n)
y[random_index] = 1

Final_beta= generalized_dantzig_selector(X, y, lmbda, max_iter=4, tol=1e-7)