Implementing HMLasso in Python with Cvxpy leads to problem with DCP rules

Question

In an attempt to implement the HMLasso (Lasso with High missing rate, see https://www.ijcai.org/proceedings/2019/0491.pdf) in Python, I wanted to use Cvxpy. However, I've been facing an exception that I fail to solve.

This is the code:

# Known variables
rho_pair = np.array([1, 2])
R = np.array([[1, 2], [2, 4]])
S_pair = np.array([[3, 6], [7, 5]])
mu = 1

################
## First problem
################

# Variable to optimize
n = R.shape[0]
print(n)
Sigma = cp.Variable((n, n), PSD=True)

# Objective to minimize
obj = cp.Minimize(cp.sum_squares(cp.multiply(R, Sigma-S_pair)))

# Constraints
constraints = []

# Solve the optimization problem
prob = cp.Problem(obj, constraints)
prob.solve()

Sigma_opt = Sigma.value

print(Sigma_opt)


################
## Second problem
################
beta = cp.Variable(n)
obj2 = cp.Minimize(0.5 * beta.T @ Sigma_opt @ beta - rho_pair.T @ beta + mu * cp.norm1(beta))

# Constraints
constraints2 = []

# Solving
prob2 = cp.Problem(obj2, constraints2)
prob2.solve()

beta_opt = beta.value

print(Sigma_opt, beta_opt)

and here is the error I get:

DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP. Its following subexpressions are not:
Promote(0.5, (2,)) @ var331 @ [[6.11942834 5.65628775]
 [5.65628775 5.228199  ]] @ var331

I have tried to replace the problematic expression by cp.quad_form(beta, Sigma_opt) but it does not change anything. Do you know how to fix this issue?

Thank you!

Answer 1

I've been able to solve this error this way:

# It appears that, due to floating points exceptions, Sigma_opt is not always
# Positive semidefinite. Hence, we shall check it.
eigenvalues = np.linalg.eig(self.Sigma_opt)[0]
min_eigenvalue = min(eigenvalues)
if min_eigenvalue < 0:
    print(f"[Warning] Sigma_opt is not PSD, its minimum eigenvalue is 
    {min_eigenvalue}. Error handled by adding {-min_eigenvalue} to each 
    eigenvalue.")
    self.Sigma_opt = self.Sigma_opt - min_eigenvalue * np.eye(self.p, self.p)

In the process of implementing this HMLasso, I've solved many other issues. Now it works. My code is available on github for those interested :-).