CVXPY求解TV最小化压缩感知问题遇DCPError报错求解决方案

问题背景 #

我要解决一个压缩感知问题，目标是最小化目标向量的Total Variation(TV)，变量向量是图像的反射率，将图像行像素堆叠成向量来计算TV。原代码实现如下：

import numpy as np
import matplotlib.pyplot as plt
import os
import numpy as np
from scipy.optimize import minimize
import scipy.io
import pandas as pd
import math

# A and Y are (2323,6561) and (2323,1) complex array respectively.
A = scipy.io.loadmat(r'C:\\Users\\mohammad\\Desktop\\feko1\\feko\\A.mat')
Y = scipy.io.loadmat(r'C:\\Users\\mohammad\\Desktop\\feko1\\feko\\Y.mat')

A = np.array(A["A"])
Y = np.array(Y["Y"])

import cvxpy as cp
delta = 10

# Define the variables
s_L1 = cp.Variable(A.shape[1])

# Define the objective function

xydim = int(math.sqrt(A.shape[1]))
            
print(xydim)

s_L1_2d = s_L1.reshape((xydim,xydim))

dim1diff = s_L1_2d[1:,:] - s_L1_2d[:-1,:]
dim2diff = s_L1_2d[:,1:]- s_L1_2d[:,:-1]
dim1diff = cp.square(dim1diff)
dim2diff = cp.square(dim2diff)
dim2diff2 = dim2diff.T

tot_var = cp.sum(cp.sqrt(dim1diff + dim2diff2))

s_L1 = s_L1.reshape((A.shape[1],1))
objective = tot_var

constraints = [cp.norm(A@s_L1 - Y,2) <= delta]

# Define the optimization problem
problem = cp.Problem(cp.Minimize(objective), constraints)

# Solve the optimization problem
problem.solve()

运行错误 #

在Jupyter中运行时出现如下DCPError：

---------------------------------------------------------------------------
DCPError                                  Traceback (most recent call last)
~\\AppData\\Local\\Temp\\ipykernel_15028\\3075150186.py in <module>
      5 
      6 # Solve the optimization problem
----> 7 problem.solve()

~\\AppData\\Roaming\\Python\\Python39\\site-packages\\cvxpy\\problems\\problem.py in solve(self, *args, **kwargs)
    501         else:
    502             solve_func = Problem._solve
--> 503         return solve_func(self, *args, **kwargs)
    504 
    505     @classmethod

~\\AppData\\Roaming\\Python\\Python39\\site-packages\\cvxpy\\problems\\problem.py in _solve(self, solver, warm_start, verbose, gp, qcp, requires_grad, enforce_dpp, ignore_dpp, canon_backend, **kwargs)
   1070                 return self.value
   1071 
--> 1072         data, solving_chain, inverse_data = self.get_problem_data(
   1073             solver, gp, enforce_dpp, ignore_dpp, verbose, canon_backend, kwargs
   1074         )

~\\AppData\\Roaming\\Python\\Python39\\site-packages\\cvxpy\\problems\\problem.py in get_problem_data(self, solver, gp, enforce_dpp, ignore_dpp, verbose, canon_backend, solver_opts)
    644         if key != self._cache.key:
    645             self._cache.invalidate()
--> 646             solving_chain = self._construct_chain(
    647                 solver=solver, gp=gp,
    648                 enforce_dpp=enforce_dpp,

~\\AppData\\Roaming\\Python\\Python39\\site-packages\\cvxpy\\problems\\problem.py in _construct_chain(self, solver, gp, enforce_dpp, ignore_dpp, canon_backend, solver_opts)
    896         candidate_solvers = self._find_candidate_solvers(solver=solver, gp=gp)
    897         self._sort_candidate_solvers(candidate_solvers)
--> 898         return construct_solving_chain(self, candidate_solvers, gp=gp,
    899                                        enforce_dpp=enforce_dpp,
    900                                        ignore_dpp=ignore_dpp,

~\\AppData\\Roaming\\Python\\Python39\\site-packages\\cvxpy\\reductions\\solvers\\solving_chain.py in construct_solving_chain(problem, candidates, gp, enforce_dpp, ignore_dpp, canon_backend, solver_opts, specified_solver)
    215     if len(problem.variables()) == 0:
    216         return SolvingChain(reductions=[ConstantSolver()])
--> 217     reductions = _reductions_for_problem_class(problem, candidates, gp, solver_opts)
    218 
    219     # Process DPP status of the problem.

~\\AppData\\Roaming\\Python\\Python39\\site-packages\\cvxpy\\reductions\\solvers\\solving_chain.py in _reductions_for_problem_class(problem, candidates, gp, solver_opts)
    130             append += ("\nHowever, the problem does follow DQCP rules. "
    131                        "Consider calling solve() with `qcp=True`.")
--> 132         raise DCPError(
    133             "Problem does not follow DCP rules. Specifically:\n" + append)
    134     elif gp and not problem.is_dgp():

DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP. Its following subexpressions are not:
power(power(reshape(var117, (81, 81), F)[1:81, 0:81] + -reshape(var117, (81, 81), F)[0:80, 0:81], 2.0) + power(reshape(var117, (81, 81), F)[0:81, 1:81] + -reshape(var117, (81, 81), F)[0:81, 0:80], 2.0).T, 0.5)

问题原因 #

CVXPY要求优化问题必须符合DCP（Disciplined Convex Programming）规则：

1. 原代码中手动计算sqrt(dim1diff + dim2diff2)的方式，本质是sqrt(a² + b²)，这种表达式不符合DCP的凸函数组合规则
2. 代码中对dim2diff做了不必要的转置dim2diff2 = dim2diff.T，导致维度不匹配，进一步破坏了表达式的合法性

解决方法 #

使用CVXPY内置的cp.norm函数计算梯度的L2范数，该函数经过封装，天然符合DCP规则；同时移除错误的转置操作，保证维度一致。修正后的代码如下：

import numpy as np
import scipy.io
import cvxpy as cp
import math

# 加载数据
A = scipy.io.loadmat(r'C:\\Users\\mohammad\\Desktop\\feko1\\feko\\A.mat')["A"]
Y = scipy.io.loadmat(r'C:\\Users\\mohammad\\Desktop\\feko1\\feko\\Y.mat')["Y"]

delta = 10

# 定义变量
s_L1 = cp.Variable(A.shape[1])
xydim = int(math.sqrt(A.shape[1]))
s_L1_2d = s_L1.reshape((xydim, xydim))

# 计算TV：用cp.norm计算每个梯度的L2范数，再求和
dim1diff = s_L1_2d[1:, :] - s_L1_2d[:-1, :]  # 行方向差分
dim2diff = s_L1_2d[:, 1:] - s_L1_2d[:, :-1]  # 列方向差分

# 对每个位置的梯度（行差分和列差分组成的向量）计算L2范数
# 注意：行差分的形状是(80,81)，列差分是(81,80)，需要分别处理后再求和
tv_rows = cp.sum(cp.norm(dim1diff, axis=0))  # 按列计算每行差分的L2范数并求和
tv_cols = cp.sum(cp.norm(dim2diff, axis=1))  # 按行计算每列差分的L2范数并求和
tot_var = tv_rows + tv_cols

# 构造约束和问题
constraints = [cp.norm(A @ s_L1.reshape((A.shape[1], 1)) - Y, 2) <= delta]
problem = cp.Problem(cp.Minimize(tot_var), constraints)

# 求解问题
problem.solve(verbose=True)  # 开启verbose可以查看求解过程

# 获取结果
reconstructed_image = s_L1.value.reshape((xydim, xydim))

关键改动说明 #

• 替换手动平方开根号的操作，改用cp.norm计算L2范数，确保表达式符合DCP规则
• 分别处理行方向和列方向的差分，通过axis参数指定范数计算维度，避免维度不匹配问题
• 简化变量定义和数据加载代码，移除不必要的库导入（如matplotlib、pandas等未用到的库）

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