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!
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 :-).