Python/cvxpy中基于快速矩阵向量乘法的凸优化问题求解优化

1. 先修正CVXPY代码的约束错误 #

你当前代码中的cvx.norm(vx, 2) <= R约束不符合原问题要求，原问题约束是||x||₂=1（其中x = vx[:n] - vx[n:]），正确的约束应为cvx.norm(vx[:n] - vx[n:], 2) == 1。错误的约束不仅会导致结果偏离，还可能额外增加求解开销。

2. 在CVXPY中引入快速傅里叶运算 #

CVXPY本身不直接支持FFT加速，但可以通过自定义线性算子封装傅里叶矩阵的乘法/转置乘法，替代大矩阵的直接存储与运算，具体步骤如下：

自定义傅里叶线性算子 #

利用Scipy的LinearOperator封装FFT/IFFT操作，避免显式存储大规模傅里叶矩阵：

import cvxpy as cp
import numpy as np
from scipy.sparse.linalg import LinearOperator

# 假设已知参数：m(A的行数), n(x的维度), R(追加值), y(符号向量)
m = ...
n = ...
R = ...
y = ...  # 元素为±1的向量

# 定义A[:,:n]的乘法：等价于m点FFT
def matvec(x):
    return np.fft.fft(x, m)

# 定义A[:,:n]的转置乘法：等价于m*IFFT(v)[:n]的实部（x为实数向量）
def rmatvec(v):
    return m * np.fft.ifft(v)[:n].real

# 封装为Scipy线性算子
A_op = LinearOperator((m, n), matvec=matvec, rmatvec=rmatvec)
# 转换为CVXPY兼容的线性算子
A_cvx = cp.LinearOperator(A_op.shape, A_op.matvec, A_op.rmatvec)

# 预先计算A[:,n]*R（等价于FFT([0,...,0,R])）
AR = np.fft.fft(np.concatenate([np.zeros(n), [R]]), m).real

修改CVXPY求解代码 #

用自定义算子替换原矩阵乘法，修正约束：

# 定义变量与目标函数
vx = cp.Variable(2*n)
x = vx[:n] - vx[n:]
objective = cp.Minimize(cp.norm(vx, 1))

# 修正后的约束
constraints = [
    vx >= 0,
    cp.multiply(A_cvx @ x + AR, y) >= 0,
    cp.norm(x, 2) == 1
]

# 选择适合大规模问题的求解器（如ECOS/SCS）
prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.SCS, verbose=True)

solution = x.value

3. 更高效的方案：自定义一阶算法（ADMM） #

对于超大规模傅里叶矩阵问题，CVXPY的内点法效率有限，推荐实现ADMM（交替方向乘子法），每次迭代仅需几次FFT/IFFT操作，复杂度为O(m log m)，远快于内点法：

import numpy as np

# ADMM参数设置
rho = 1.0
max_iter = 1000
tol = 1e-6

# 初始化变量
x = np.random.randn(n)
x = x / np.linalg.norm(x)  # 初始化为||x||₂=1
z = np.fft.fft(np.concatenate([x, [R]]), m).real
u = np.zeros(m)
z_prev = z.copy()

for k in range(max_iter):
    # 1. x-update：L1正则化求解+投影到单位球
    c = z - u
    # 计算A^T(c - AR)
    AT_c_minus_AR = m * np.fft.ifft(c - AR)[:n].real
    # 软阈值算子求解L1正则化问题
    x_unproj = np.sign(AT_c_minus_AR) * np.maximum(np.abs(AT_c_minus_AR) - 1/(rho * m), 0)
    # 投影到||x||₂=1
    x = x_unproj / np.linalg.norm(x_unproj)
    
    # 2. z-update：投影到z⊙y ≥0的约束集
    AxR = np.fft.fft(np.concatenate([x, [R]]), m).real
    z = np.where(y == 1, np.maximum(AxR + u, 0), np.minimum(AxR + u, 0))
    
    # 3. u-update：对偶变量更新
    u = u + AxR - z
    
    # 收敛判断
    primal_res = np.linalg.norm(AxR - z)
    dual_res = np.linalg.norm(-rho * (z - z_prev)) if k > 0 else np.inf
    z_prev = z.copy()
    
    if primal_res < tol and dual_res < tol:
        break

solution = x

4. 其他可选库 #

• PyTorch/TensorFlow：适合需要自动求导的场景，可利用内置FFT操作和优化器（如Adam）实现一阶算法，结合投影操作处理约束。
• pyproximal：专门针对近端算法的工具库，支持自定义线性算子，适合大规模稀疏优化问题。

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