使用Cvxpy在Python实现HMLasso时遭遇DCP规则合规性问题

我尝试在Python中使用Cvxpy实现HMLasso（高缺失率Lasso），但遇到了DCP规则不满足的异常。

我的代码如下：

# 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)

运行后报错：

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

我尝试用cp.quad_form(beta, Sigma_opt)替换问题表达式，但问题依旧，请问该如何解决？

问题核心是第一次优化得到的Sigma_opt因数值误差，不再严格满足半正定（PSD）属性，导致Cvxpy无法确认二次型的凸性，违反DCP规则。以下是具体解决步骤：

1. 修正Sigma_opt的半正定性
   对Sigma_opt做半正定投影，通过特征值分解将微小负特征值替换为0，重构严格半正定矩阵：

import numpy as np

# 特征值分解修正半正定性
eigenvals, eigenvecs = np.linalg.eigh(Sigma_opt)
# 把极小的负特征值截断为0（阈值可根据精度调整）
eigenvals[eigenvals < 1e-8] = 0
Sigma_psd = eigenvecs @ np.diag(eigenvals) @ eigenvecs.T

2. 规范二次型表达式
   使用Cvxpy原生的cp.quad_form定义二次项，确保符合DCP规范：

beta = cp.Variable(n)
# 使用修正后的半正定矩阵构建目标函数
obj2 = cp.Minimize(0.5 * cp.quad_form(beta, Sigma_psd) - rho_pair.T @ beta + mu * cp.norm1(beta))
constraints2 = []
prob2 = cp.Problem(obj2, constraints2)
prob2.solve()
beta_opt = beta.value

3. 完整修正代码

import numpy as np
import cvxpy as cp

# 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
################

n = R.shape[0]
Sigma = cp.Variable((n, n), PSD=True)
obj = cp.Minimize(cp.sum_squares(cp.multiply(R, Sigma-S_pair)))
constraints = []
prob = cp.Problem(obj, constraints)
prob.solve()

Sigma_opt = Sigma.value

# 修正Sigma_opt的半正定性
eigenvals, eigenvecs = np.linalg.eigh(Sigma_opt)
eigenvals[eigenvals < 1e-8] = 0
Sigma_psd = eigenvecs @ np.diag(eigenvals) @ eigenvecs.T

print("修正后的Sigma:\n", Sigma_psd)

################
## Second problem
################
beta = cp.Variable(n)
obj2 = cp.Minimize(0.5 * cp.quad_form(beta, Sigma_psd) - rho_pair.T @ beta + mu * cp.norm1(beta))
constraints2 = []
prob2 = cp.Problem(obj2, constraints2)
prob2.solve()

beta_opt = beta.value

print("最优beta:\n", beta_opt)

内容的提问来源于stack exchange，提问作者Noomkwah