0

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!

Noomkwah
  • 133
  • 6
  • That's because ``Sigma_opt`` is not strictly PSD, it has an eigenvalue of order -1e-9. This is a sort of numerical issue that always can happen, the same way as a variable declared as ``x>=0`` can end up as -1e-9 in the solution. Life in floating point. You can fix it for instance by boosting ``Sigma_opt`` by ``np.eye(2)*1e-8`` or make some other tiny numerical fix to sanitize the data before the second problem. And of course use ``quad_form``. – Michal Adamaszek Mar 29 '23 at 19:54
  • Ok thank you! Why is using quad_form better than writing `beta.T @ Sigma_opt @ beta` ? – Noomkwah Mar 30 '23 at 12:15
  • The latter will not be recognized as convex. You have to follow the rules for informing CVXPY about the convex structure. quad_form does that, while ``x.T@S@x`` is just a case of ``y.T@S@x``, which is not convex in general. In other words, CVXPY does not check that both sides of ``Sigma_opt`` are multiplied with the same ``x``, you have to tell CVXPY that yourself. Maybe it could, but then where would be the limit of how deep it should go through the structure? – Michal Adamaszek Mar 30 '23 at 12:19
  • I see thank you! Actually, by adding Sigma >> 0 in the constraints, the code does converge. Is this a good practice or does that have any cons? – Noomkwah Mar 30 '23 at 12:22
  • It is a coincidence which you could also achieve or nivelate by changing the solver or other such small changes. – Michal Adamaszek Mar 30 '23 at 12:39
  • Hello @MichalAdamaszek! As you seem quite well knowledgeable on the subject, i was wondering if you would agree to help me on again on a related issue I am facing? I keep encountering the same error again :(. I explain it here: https://stackoverflow.com/questions/75965882/cvxpy-quadratic-programming-arpacknoconvergence-error-and-assertionerror – Noomkwah Apr 08 '23 at 17:42

1 Answers1

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

Noomkwah
  • 133
  • 6