企业官网建设流程全解析

• 科学计算

导读

本文以 CVXPY 官方示例 doc/source/examples/dqcp/minimum_length_least_squares.rst 为骨架，讲解如何在 CVXPY 中建模并求解"最小长度最小二乘"问题：在保证均方误差（MSE）不超过给定阈值的前提下，寻找长度（最后一个非零分量的下标）最短的向量解。该问题是典型的拟凸规划（Quasiconvex Program, QCP），可通过**有纪律的拟凸规划（DQCP）**框架声明式建模，并借助二分法（bisection）求解。读完本文，你将掌握cp.length拟凸原子、DQCP 建模语法、problem.solve(qcp=True)求解入口，以及 CVXPY 内部将 QCP 归约为参数化 DCP 问题并二分求解的完整原理，并可在自己的稀疏解恢复类任务中直接复用。

一、问题定义：什么是最小长度最小二乘

最小长度最小二乘问题的目标是：找到一个长度最短的向量 x ∈ Rⁿ，使其对给定最小二乘系统的均方误差足够小。其标准形式为：

minimize len(x) subject to (1/n)·‖Ax - b‖₂² ≤ ε

其中：

• 变量：x ∈ Rⁿ；
• 问题数据：矩阵 A、观测向量 b、样本数 n、误差容忍度 ε；
• 目标函数：len(x)，即向量 x 的"长度"——最后一个非零分量的下标（按 1 计数）；若 x 全为零向量则长度为 0；
• 约束：归一化均方误差(1/n)·‖Ax - b‖₂² ≤ ε。

直观理解：给定存在噪声或冗余的线性方程组 Ax ≈ b，我们希望用尽量"稀疏"（有效分量集中在开头、尾部尽快归零）的解逼近它，这与稀疏信号恢复、特征选择等场景目标一致。

从 CVXPY 源码看，length原子定义于 cvxpy/atoms/length.py，其数值实现使用容差ATOM_EVAL_TOL判定非零分量，再取最大非零下标加 1：

• 仅接受向量输入，对非向量会抛出ValueError("length can only be applied to vectors.")（cvxpy/atoms/length.py）；
• 数值计算中，np.abs(x) > ATOM_EVAL_TOL判定非零，若无非零元素返回 0，否则返回np.max(nz) + 1（cvxpy/atoms/length.py）。

二、为什么是 QCP/DQCP：曲率分析

关键问题是：len(x)与 MSE 约束组合后，问题本身不是凸问题，但属于拟凸问题。

• len(x)既不凸也不凹（源码中is_atom_convex与is_atom_concave均返回False），但它满足**拟凸（quasiconvex）**性质：is_atom_quasiconvex返回True（cvxpy/atoms/length.py）。其子水平集{x : len(x) ≤ t}是凸集——当 t ≥ n 时退化为全空间，否则等价于x[t:] == 0，显然为凸集；
• MSE 约束(1/n)·‖Ax - b‖₂² ≤ ε是凸约束（平方范数凸）；
• "拟凸目标 + 凸约束"构成的极小化问题即为拟凸规划（QCP）。

CVXPY 通过**有纪律的拟凸规划（DQCP）**语法支持此类问题：在length上还规定其单调性——对非负参数单调不减、对非正参数单调不增（is_incr/is_decr，见 cvxpy/atoms/length.py），这是 DQCP 组合规则（单调函数复合拟凸原子）成立的前提。

注意：length的单调性依赖参数符号。示例 cvxpy/tests/test_dqcp.py 验证了length(abs(x))单调不减故为 DQCP，而length(abs(x) - 1)不再单调，is_dqcp()返回False。

三、完整可运行示例：建模与求解

以下代码完整复刻官方示例（问题数据构造与示例一致），可直接在安装了 CVXPY 与 NumPy 的 Python 3 环境中运行。

# 安装最新版 CVXPY（如尚未安装） # !pip install --upgrade cvxpy import cvxpy as cp import numpy as np # 1. 构造问题数据 n = 10 np.random.seed(1) A = np.random.randn(n, n) x_star = np.random.randn(n) # 真实解（作为对照） b = A @ x_star # 无噪声观测 epsilon = 1e-2 # 均方误差容忍度 # 2. 声明变量与目标/约束表达式 x = cp.Variable(n) mse = cp.sum_squares(A @ x - b) / n problem = cp.Problem(cp.Minimize(cp.length(x)), [mse <= epsilon]) # 3. 验证 DQCP 属性并求解 print("Is problem DQCP?: ", problem.is_dqcp()) problem.solve(qcp=True) print("Found a solution, with length: ", problem.value)

求解输出（示例文档中的运行结果）：

Is problem DQCP?: True Found a solution, with length: 8.0

进一步查看均方误差与最优解：

print("MSE: ", mse.value) print("x: ", x.value) print("x_star: ", x_star)

示例运行输出：

MSE: 0.00926009328813662 x: [-2.58366030e-01 1.38434327e+00 2.10714108e-01 9.44811159e-01 -1.14622208e+00 1.51283929e-01 6.62931941e-01 -1.16358584e+00 2.78132907e-13 -1.76314786e-13] x_star: [-0.44712856 1.2245077 0.40349164 0.59357852 -1.09491185 0.16938243 0.74055645 -0.9537006 -0.26621851 0.03261455]

结果解读：

• problem.value为 8.0，即最优解 x 的最后一个非零分量下标为 8（第 9、10 个分量在数值容差内为 0），与官方示例结果一致；
• mse.value ≈ 0.00926 < ε = 1e-2，约束满足；
• 对照x_star可见解的前 8 个分量与真实解接近，后两个分量被压缩至 1e-13 量级（数值零），实现了"尾部截断"式的稀疏化。

官方示例被 cvxpy/tests/test_dqcp.py 的test_length_example完整收录（注释Fix #1760，作为回归测试），断言problem.value与 8 数值接近，进一步印证了上述输出具有可复现性。

四、求解原理：DQCP → 参数化 DCP + 二分法

problem.solve(qcp=True)并非直接求解 QCP，而是走一条"归约 + 二分"的流水线：

1. DQCP 归约（Dqcp2Dcp）：将拟凸目标len(x)替换为带参数 t 的不等式len(x) ≤ t，原问题变成一族参数化的 DCP 可行性问题{x : mse ≤ ε, len(x) ≤ t}；若目标非负则 t 声明为非负 Parameter，否则为一般 Parameter（cvxpy/reductions/dqcp2dcp/dqcp2dcp.py）。其中"惰性约束（lazy constraints）"以可调用对象形式存在，仅在 t 取具体值时求值；
2. 子水平集改写：拟凸原子length的子水平集按 cvxpy/reductions/dqcp2dcp/sets.py 的length_sub展开：当 t 为 Parameter 时，若 t < 0 视为不可行、t ≥ n 视为无约束，否则生成约束x[floor(t):] == 0（对 t 向下取整）；当 t 为常数时直接生成该等式约束。这就是"len(x) ≤ t ⟺ 尾部 t 个分量全为 0"的 DCP 化表达；
3. 二分求解（bisect）：cvxpy/reductions/solvers/bisection.py 的bisect函数在 t 上做二分：先求解可行性问题判断原问题是否可行；若未给定区间，则通过_find_bisection_interval指数扩张寻找 [low, high]（cvxpy/reductions/solvers/bisection.py），随后在_bisect中反复取中点 t = (low+high)/2 求解参数化 DCP，直至区间宽度小于 eps（默认 1e-6）或达到最大迭代次数（默认 100）；
4. 整数目标加速：length属于整数取值原子（与ceil、floor并列），cvxpy/reductions/dqcp2dcp/tighten.py 为其提供收紧函数：不可行侧向上取整np.ceil、可行侧向下取整np.floor，使二分区间快速收敛到精确整数值，这正是示例能精确得到 8.0 而非 7.9 之类近似值的原因。

可选的verbose=True参数会逐次打印二分迭代的下界、上界与查询点（cvxpy/reductions/solvers/bisection.py）；cvxpy/tests/test_dqcp.py 还提供了verbose=True下length原子不崩溃的回归测试。

五、DQCP 建模规则速查与注意事项

从示例及 CVXPY 源码可归纳出以下使用要点：

• 入口开关：DQCP 问题的求解必须显式传qcp=True，即problem.solve(qcp=True)；否则 CVXPY 会按 DCP 路径处理并报错；
• 属性自检：建模后先用problem.is_dqcp()、expr.is_quasiconvex()/expr.is_quasiconcave()验证，避免进入求解阶段才发现语法错误。源码曲率测试见 cvxpy/tests/test_dqcp.py：length(x)曲率为 QUASICONVEX，-length(x)为 QUASICONCAVE；
• 单调性组合规则：DQCP 允许"单调函数 ∘ 拟凸原子 ∘ DCP 表达式"的复合，但单调性必须成立（见 cvxpy/reductions/dqcp2dcp/dqcp2dcp.py 的约束改写逻辑：单调函数通过求逆作用到约束两侧，拟凸原子替换为其子水平集/超水平集）。破坏单调性（如length(abs(x) - 1)）会导致is_dqcp()为False；
• length 的输入限制：length只能作用于向量（源码显式校验），矩阵输入会抛出异常；
• 容差语义：length的数值结果依赖ATOM_EVAL_TOL容差，非零判定是"按数值容差"而非"精确非零"，因此打印出的解尾部会出现 1e-13 量级的数值残差；
• 求解器选择：二分内部使用常规凸求解器求解每个子问题，默认选择可用求解器即可；若需手动指定，可传problem.solve(qcp=True, solver=cp.SCS)等。测试 cvxpy/tests/test_dqcp.py 展示了length加等式约束的小规模验证（期望值为 2），可作为最小复现样例。

六、延伸应用：从示例到实际场景

最小长度最小二乘示例虽小，其"拟凸目标 + 凸约束"的建模范式可推广到大量实际问题：

• 稀疏解恢复：将目标换成len(x)或拟凸的sign、dist_ratio等原子，配合 L2 拟合约束，可在不做 L1 松弛的情况下直接求"最短有效支撑"解；
• 特征/模型选择：在误差上限约束下最小化有效参数量，得到尾部截断的稀疏系数向量；
• 通用 DQCP 工具箱：CVXPY 在 cvxpy/reductions/dqcp2dcp/sets.py 中注册了multiply（乘积）、DivExpression（比率）、length、sign、dist_ratio、gen_lambda_max、condition_number等原子的子水平集/超水平集改写，同类拟凸目标（如条件数最小化、广义最大特征值）均可套用本文的建模与求解流程。

结语

本文从 doc/source/examples/dqcp/minimum_length_least_squares.rst 出发，完整给出了最小长度最小二乘问题的数学模型、CVXPY 可运行代码、DQCP 属性分析与结果解读，并深入源码揭示了其底层实现：length拟凸原子（cvxpy/atoms/length.py）→ DQCP 归约（cvxpy/reductions/dqcp2dcp/dqcp2dcp.py）→ 子水平集改写（cvxpy/reductions/dqcp2dcp/sets.py）→ 二分求解（cvxpy/reductions/solvers/bisection.py）与整数收紧（cvxpy/reductions/dqcp2dcp/tighten.py）。掌握这套"拟凸建模 + qcp=True + 二分归约"的方法论后，你可以在 CVXPY 中自信地处理更大规模的拟凸规划问题。

• 科学计算