企业官网建设流程全解析

• 科学计算

导读

本文以 CVXPY 官方示例 tv_inpainting.rst 为蓝本，系统讲解如何使用 CVXPY 对灰度图与彩色图进行总变分（Total Variation, TV）图像修复（in-painting）：给定一幅部分像素损坏/缺失的图像，通过最小化图像的 TV 范数并约束已知像素不变，重建出完整图像。读完本文，你将掌握 TV 修复的数学模型、CVXPY 中tv原子（atom）的底层实现原理、灰度与彩色两套完整可运行的求解代码，以及面对超大规模锥规划问题时如何选择求解器并读懂 SCS 的迭代日志。

一、问题背景：什么是图像修复（Inpainting）

图像修复指在已知部分像素值的前提下，推测并填补缺失像素，从而恢复整幅图像。典型应用场景包括：

• 文字/Logo 遮挡恢复：如示例中被白色文字“This is the Loki test image...”覆盖的图片；
• 随机噪声/划痕去除：如示例中约 70% 像素被随机置零的彩色图；
• 老旧照片修复、视频去字幕、医学影像补全等。

在数学上，图像修复是一个不适定（ill-posed）问题：缺失像素有无数种补法。为此需要引入先验正则——总变分，它偏好“分段平滑”（piecewise smooth）的图像，即在像素值变化剧烈处允许跳变（保留边缘），在平坦区域则抑制噪声。这正是 TV 修复在图像去噪、去遮挡中效果出色的原因。

二、灰度图 TV 修复：数学模型

2.1 图像与已知像素的表示

一幅灰度图表示为 $m \times n$ 的强度矩阵 $U^{\text{orig}}$（像素值通常落在 $[0, 255]$）。已知像素下标集合为 $\mathcal{K} \subset {1,\ldots,m} \times {1,\ldots,n}$，已知值记为 $U^{\text{orig}}_{ij}, (i,j)\in\mathcal{K}$。我们要重建的矩阵为 $U \in \mathbf{R}^{m\times n}$，且必须满足保真约束：

$$U_{ij} = U^{\text{orig}}_{ij}, \quad \forall (i,j) \in \mathcal{K}$$

2.2 $\ell_2$ 总变分定义

CVXPY 示例采用$\ell_2$ 总变分（对矩阵按离散梯度取 2-范数再求和）：

$$\mathop{\bf tv}(U) = \sum_{i=1}^{m-1}\sum_{j=1}^{n-1}\left|\begin{bmatrix} U_{i+1,j}-U_{ij} \ U_{i,j+1}-U_{ij}\end{bmatrix}\right|_2$$

即对每个像素计算“右邻像素差、下邻像素差”组成的梯度向量，取其欧几里得范数，再对全图求和。注意：范数不加平方——这是 TV 与 Tikhonov 正则（$|\nabla U|_2^2$）的本质区别：不加平方的 $\ell_2$ 范数使目标函数对边缘处的梯度惩罚是线性的，从而允许图像中存在锐利边缘，而不会像平方惩罚那样过度模糊。

2.3 优化问题

重建图像 $U$ 通过求解如下凸优化问题得到：

$$ \begin{aligned} & \underset{U}{\text{minimize}} && \mathop{\bf tv}(U) \ & \text{subject to} && U_{ij} = U^{\text{orig}}_{ij}, \quad (i,j)\in\mathcal{K} \end{aligned} $$

目标函数 $\mathop{\bf tv}(U)$ 是凸的（范数求和），约束是线性等式约束，因此整个问题是一个凸优化问题，可用 CVXPY 直接建模并交给锥求解器（conic solver）求解。

三、灰度图修复：CVXPY 实现与代码逐行解析

3.1 加载图像与构造 Known 矩阵

import matplotlib.pyplot as plt import numpy as np # 加载原始图与损坏图。 u_orig = plt.imread("data/loki512.png") u_corr = plt.imread("data/loki512_corrupted.png") rows, cols = u_orig.shape # known 为 1 表示像素已知，0 表示像素被损坏。 known = np.zeros((rows, cols)) for i in range(rows): for j in range(cols): if u_orig[i, j] == u_corr[i, j]: known[i, j] = 1 %matplotlib inline fig, ax = plt.subplots(1, 2, figsize=(10, 5)) ax[0].imshow(u_orig, cmap='gray') ax[0].set_title("Original Image") ax[0].axis('off') ax[1].imshow(u_corr, cmap='gray'); ax[1].set_title("Corrupted Image") ax[1].axis('off');

关键点：

• u_orig与u_corr是形状为(rows, cols)的 NumPy 数组（示例图片为 512×512）；
• known是 0/1 掩码矩阵：像素未被损坏（两图相等）记为 1，被损坏记为 0；
• 掩码矩阵known后续既作为约束的乘子，也用于从u_corr中“抠出”已知像素值。

3.2 建模并求解

# 使用总变分修复重建原始图像。 import cvxpy as cp U = cp.Variable(shape=(rows, cols)) obj = cp.Minimize(cp.tv(U)) constraints = [cp.multiply(known, U) == cp.multiply(known, u_corr)] prob = cp.Problem(obj, constraints) # 使用 SCS 求解。 prob.solve(verbose=True, solver=cp.SCS) print("optimal objective value: {}".format(obj.value))

代码中三个要素对应数学模型：

• U = cp.Variable(shape=(rows, cols))：重建图像变量，形状与图像一致；
• cp.tv(U)：TV 原子，即目标函数 $\mathop{\bf tv}(U)$；
• cp.multiply(known, U) == cp.multiply(known, u_corr)：等式约束。cp.multiply是逐元素 Hadamard 乘法（见 cvxpy/atoms/affine/binary_operators.py），乘上 0/1 掩码后，约束等价于“在已知像素处 $U$ 等于原图值，在未知像素处 $0=0$ 无约束”，巧妙地将已知像素约束与缺失像素自由度写入同一个等式。

prob.solve(verbose=True, solver=cp.SCS)使用SCS（Splitting Conic Solver）求解并打印详细日志；求解结果存入U.value，最优目标值由obj.value读取。

3.3 SCS 求解日志解读（灰度图）

示例在 512×512 图上运行 SCS v2.0.2 的输出摘要如下：

Lin-sys: sparse-indirect, nnz in A = 1554199, CG tol ~ 1/iter^(2.00) eps = 1.00e-05, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00 Variables n = 523265, constraints m = 1045507 Cones: primal zero / dual free vars: 262144 soc vars: 783363, soc blks: 261121 ... Status: Solved Solve time: 2.55e+02s ... c'x = 11044.2661, -b'y = 11044.2813 optimal objective value: 11044.28989542425

这些数字揭示出问题的真实规模与结构：

• 变量数 $n = 523265$：512×512 = 262144 个像素变量，加上锥分解引入的松弛变量（soc vars 783363 等），总计约 52 万；
• 约束数 $m = 1045507$：包含 262144 个 primal zero / dual free 变量（对应保真等式约束的松弛），以及 783363 个二阶锥（SOC）分量、261121 个 SOC 块——每个像素对应一个二维 SOC 块，正是 $\ell_2$ 范数 $|[U_{i+1,j}-U_{ij};,U_{i,j+1}-U_{ij}]|_2$ 被 SCS 转换为 SOC 约束的结果；
• 求解耗时约 255 秒（该示例运行时的硬件环境下的实测值），收敛到primal res ≈ 9.0e-06、dual res ≈ 8.2e-06、rel gap ≈ 6.9e-07，均在默认容差eps=1e-05之下，判定Status: Solved；
• 最优目标值 ≈ 11044.29，即重建图像的 TV 值。

文档明确提示：这里选择 SCS 是因为它能扩展到比 ECOS 更大的问题规模。ECOS 基于内点法，对 52 万变量、105 万约束的锥规划内存开销过大。

3.4 展示修复结果与差异图

fig, ax = plt.subplots(1, 2, figsize=(10, 5)) # 展示修复后的图像。 ax[0].imshow(U.value, cmap='gray'); ax[0].set_title("In-Painted Image") ax[0].axis('off') img_diff = 10*np.abs(u_orig - U.value) ax[1].imshow(img_diff, cmap='gray'); ax[1].set_title("Difference Image") ax[1].axis('off');

• 修复结果直接来自U.value（Variable求解后的数值）；
• 差异图为10 * |u_orig - U.value|，将差异放大 10 倍以便肉眼观察。修复图与原始图几乎一致，差异图仅在被文字遮挡的区域显示微弱残影——说明 TV 修复成功抹去了文字而保留了底层图像结构。

四、彩色图 TV 修复：三通道变量与随机掩码

4.1 数学模型扩展

彩色图表示为 $m\times n\times 3$ 的 RGB 矩阵 $U^{\text{orig}}$，每个像素 $U^{\text{orig}}_{ij} \in \mathbf{R}^3$ 是一个 RGB 向量。TV 定义与灰度版形式相同，但每个梯度分量变为三维向量：

$$\mathop{\bf tv}(U) = \sum_{i=1}^{m-1}\sum_{j=1}^{n-1}\left|\begin{bmatrix} U_{i+1,j}-U_{ij} \ U_{i,j+1}-U_{ij}\end{bmatrix}\right|_2$$

这里向量的每个“分量”本身是三维 RGB 差向量，范数取在拼接后的六维向量上（对应到 CVXPY 实现，则是对 R/G/B 三个通道的梯度统一堆叠取 2-范数，见下文源码解析）。

4.2 随机损坏掩码的构造

与灰度版“文字遮挡”不同，彩色版通过随机丢弃 70% 像素制造损坏：

import matplotlib.pyplot as plt import numpy as np np.random.seed(1) # 加载图像。 u_orig = plt.imread("data/loki512color.png") rows, cols, colors = u_orig.shape # known 为 1 表示像素已知，0 表示像素被损坏。 # known 矩阵随机初始化。 known = np.zeros((rows, cols, colors)) for i in range(rows): for j in range(cols): if np.random.random() > 0.7: for k in range(colors): known[i, j, k] = 1 u_corr = known * u_orig %matplotlib inline fig, ax = plt.subplots(1, 2, figsize=(10, 5)) ax[0].imshow(u_orig, cmap='gray'); ax[0].set_title("Original Image") ax[0].axis('off') ax[1].imshow(u_corr); ax[1].set_title("Corrupted Image") ax[1].axis('off');

• np.random.seed(1)固定随机种子，保证可复现；
• 每个像素以 30% 概率保留（np.random.random() > 0.7），三个通道共享同一保留/丢弃决策；
• u_corr = known * u_orig：被丢弃的像素值置 0，损坏图呈黑色斑点状。

4.3 三变量建模与求解

# 使用总变分修复重建原始图像。 import cvxpy as cp variables = [] constraints = [] for i in range(colors): U = cp.Variable(shape=(rows, cols)) variables.append(U) constraints.append(cp.multiply(known[:, :, i], U) == cp.multiply(known[:, :, i], u_corr[:, :, i])) prob = cp.Problem(cp.Minimize(cp.tv(*variables)), constraints) prob.solve(verbose=True, solver=cp.SCS) print("optimal objective value: {}".format(prob.value))

• 三个矩阵变量：R、G、B 各用一个(rows, cols)的cp.Variable，放入variables列表；
• 逐通道约束：对每个通道i，用该通道的掩码known[:, :, i]构造保真等式约束；
• cp.tv(*variables)：tv原子接受多个矩阵参数，将三个通道作为“第三维”统一处理（这正是彩色 TV 的实现方式，详见下一节源码）。

求解器同样选择 SCS，文档明确指出：ECOS 和 CVXOPT 无法扩展到如此大规模的问题（彩色版问题规模比灰度版更大）。

4.4 彩色版求解规模与日志要点

WARN: A->p (column pointers) not strictly increasing, column 523264 empty WARN: A->p (column pointers) not strictly increasing, column 785408 empty WARN: A->p (column pointers) not strictly increasing, column 1047552 empty ... Lin-sys: sparse-indirect, nnz in A = 3630814, CG tol ~ 1/iter^(2.00) Variables n = 1047553, constraints m = 2614279 Cones: primal zero / dual free vars: 786432 soc vars: 1827847, soc blks: 261121 ... Status: Solved Solve time: 6.99e+02s ... optimal objective value: 11465.652787130613

• 变量数升至约 105 万（262144×3 个通道像素变量加上锥松弛），约束数约 261 万；
• 日志开头的三条WARN: A->p ... column ... empty是 SCS 对稀疏矩阵中空列的提示，属于无害告警；
• 求解约 700 秒收敛，primal res ≈ 9.0e-06、dual res ≈ 9.7e-06、rel gap ≈ 2.4e-06，判定Status: Solved；
• 最优目标值约 11465.65。

4.5 彩色修复结果的可视化

import matplotlib.pyplot as plt import matplotlib.cm as cm %matplotlib inline rec_arr = np.zeros((rows, cols, colors)) for i in range(colors): rec_arr[:, :, i] = variables[i].value rec_arr = np.clip(rec_arr, 0, 1) fig, ax = plt.subplots(1, 2, figsize=(10, 5)) ax[0].imshow(rec_arr) ax[0].set_title("In-Painted Image") ax[0].axis('off') img_diff = np.clip(10 * np.abs(u_orig - rec_arr), 0, 1) ax[1].imshow(img_diff) ax[1].set_title("Difference Image") ax[1].axis('off')

• 三个变量的.value分别对应 R/G/B 通道，重新拼装成(rows, cols, 3)数组rec_arr；
• 用np.clip(rec_arr, 0, 1)将浮点像素值裁剪回合法显示范围（plt.imread读入的 PNG 像素归一化到 $[0,1]$）；
• 差异图同样放大 10 倍并裁剪到 $[0,1]$。视觉上修复图与原始图几乎一致，但 RGB 差异图显示许多像素的通道值仍有可见差异——这是随机丢失 70% 像素后重建的必然结果。

五、源码纵深：CVXPY 中tv原子是如何实现的

TV 修复的核心依赖是 CVXPY 的tv原子，其完整实现位于 cvxpy/atoms/total_variation.py，并通过 cvxpy/atoms/init.py 的from cvxpy.atoms.total_variation import tv导出为cp.tv。

5.1 向量与矩阵的分支处理

value = Expression.cast(value) if value.ndim == 0: raise ValueError("tv cannot take a scalar argument.") # 向量使用 L1 范数。 elif value.ndim == 1: return norm(value[1:] - value[0:value.shape[0]-1], 1) # 矩阵使用 L2 范数。 elif value.ndim == 2: ... else: raise ValueError("tv cannot have input arrays with more than 2 dimensions.")

• 一维向量：相邻元素差的一范数 $\text{tv}(x) = \sum_i |x_{i+1}-x_i|$，即 L1 总变分；
• 二维矩阵：离散梯度的 $\ell_2$ 范数求和，即本文 2.2 节的公式；
• 标量（ndim == 0）与三维以上输入会抛出ValueError——彩色图正是通过“多个矩阵参数”而非三维数组来绕过该限制。

5.2 彩色 TV 的核心：多矩阵参数堆叠

rows, cols = value.shape args = map(Expression.cast, args) values = [value] + list(args) diffs = [] for mat in values: diffs += [ mat[0:rows-1, 1:cols] - mat[0:rows-1, 0:cols-1], mat[1:rows, 0:cols-1] - mat[0:rows-1, 0:cols-1], ] length = diffs[0].shape[0]*diffs[1].shape[1] stacked = vstack([reshape(diff, (1, length), order='F') for diff in diffs]) return sum(norm(stacked, p=2, axis=0))

实现分四步：

1. 取离散梯度：对每个矩阵（第一个参数value加上*args中的其余通道矩阵），计算两个差分切片——mat[0:rows-1, 1:cols] - mat[0:rows-1, 0:cols-1]是右邻差（对应公式中 $U_{i,j+1}-U_{ij}$），mat[1:rows, 0:cols-1] - mat[0:rows-1, 0:cols-1]是下邻差（对应 $U_{i+1,j}-U_{ij}$）；
2. 按列优先展平：每个差分矩阵用reshape(..., order='F')展平成一行，length = (rows-1)*(cols-1)是差分元素总数；
3. 堆叠：vstack将所有通道、所有方向的差分行堆成一个大矩阵，每一“列”恰好包含 $U_{i+1,j}-U_{ij}$ 与 $U_{i,j+1}-U_{ij}$ 在各通道上的分量；
4. 逐列取 2-范数并求和：norm(stacked, p=2, axis=0)对每列求 $\ell_2$ 范数（灰度图列长为 2，彩色图列长为 $2\times3=6$，正是 4.1 节公式的向量化实现），最后sum求和。

这也解释了求解日志中的锥结构：norm(..., p=2, axis=0)中的每个 $\ell_2$ 范数项都会被 CVXPY 的 DCP 到锥（dcp2cone）约简转换成一个二阶锥约束，最终由 SCS 处理（相关约简代码见 cvxpy/reductions/dcp2cone/canonicalizers/pnorm_canon.py）。灰度图soc blks: 261121（$261121 = (512-1)\times(512-1)$）与像素差分位置数完全一致，可作为实现正确性的交叉验证。

5.3 仓库中的自动化测试印证

TV 修复流程在仓库测试中有最小化复现：cvxpy/tests/test_examples.py 的test_inpainting用 20×20 的随机图像、30% 保留率掩码，构造了与本文完全相同的模型：

U = cvx.Variable((rows, cols)) obj = cvx.Minimize(cvx.tv(U)) constraints = [cvx.multiply(Known, U) == cvx.multiply(Known, Ucorr)] prob = cvx.Problem(obj, constraints) prob.solve(solver=cvx.SCS)

该测试作为TestExamples套件的一部分随 CI 运行，保证cp.tv、cp.multiply与 SCS 链路长期可用。

六、求解器选择：为什么用 SCS 而不是 ECOS/CVXOPT

文档明确的两条结论：灰度版“SCS scales to larger problems than ECOS does”；彩色版“ECOS and CVXOPT don't scale to this large problem”。此外，从仓库求解器注册表 cvxpy/reductions/solvers/defines.py 可以看到，SCS、ECOS、CVXOPT 均为 CVXPY 内置支持的锥求解器，用户只需pip install scs等安装对应求解器包即可（详见 doc/source/install/index.rst）。

实用建议：

• 追求速度与规模：选 SCS（或同样一阶的 Clarabel、COSMO 等现代锥求解器，仓库中均有对应接口，见 cvxpy/reductions/solvers/conic_solvers）；
• 追求高精度小规模：选 ECOS/CVXOPT/MOSEK；
• 调试阶段可先用小图（如测试中的 20×20）快速验证建模正确性，再上全尺寸图像。

七、完整可运行代码（灰度版）

将以下代码保存为脚本（把两张示例图片放到当前目录下的data/子目录）即可端到端复现：

import matplotlib.pyplot as plt import numpy as np import cvxpy as cp # 1. 加载图像 u_orig = plt.imread("data/loki512.png") u_corr = plt.imread("data/loki512_corrupted.png") rows, cols = u_orig.shape # 2. 构造已知像素掩码 known = np.zeros((rows, cols)) for i in range(rows): for j in range(cols): if u_orig[i, j] == u_corr[i, j]: known[i, j] = 1 # 3. 建模 U = cp.Variable(shape=(rows, cols)) obj = cp.Minimize(cp.tv(U)) constraints = [cp.multiply(known, U) == cp.multiply(known, u_corr)] prob = cp.Problem(obj, constraints) # 4. 求解（SCS 可扩展至大规模） prob.solve(verbose=True, solver=cp.SCS) print("optimal objective value: {}".format(obj.value)) # 5. 可视化 fig, ax = plt.subplots(1, 2, figsize=(10, 5)) ax[0].imshow(U.value, cmap='gray') ax[0].set_title("In-Painted Image") ax[0].axis('off') img_diff = 10*np.abs(u_orig - U.value) ax[1].imshow(img_diff, cmap='gray') ax[1].set_title("Difference Image") ax[1].axis('off') plt.show()

彩色版仅需把第 3 步替换为 4.3 节的循环建模即可。更详细的完整示例与其余应用案例（水填充、信道容量、鲁棒卡尔曼滤波等）可在 doc/source/examples/applications 目录下找到。

八、总结与延伸阅读

本文完整复现了 CVXPY 官方文档的灰度与彩色 TV 图像修复示例：从 $\ell_2$ 总变分的数学定义出发，推导出以“TV 最小化 + 已知像素保真约束”为核心的凸优化模型；随后给出两套可运行的 CVXPY 代码，解读了 SCS 求解日志中的锥结构与规模信息；最后深入 cvxpy/atoms/total_variation.py 源码，说明tv原子如何通过离散梯度差分、按列堆叠与逐列 2-范数求和，将数学公式精确翻译为 DCP 表达式并最终转换为二阶锥规划。

延伸方向：

• TV 修复的目标函数可以扩展为加权 TV、各向异性 TV 等变体，CVXPY 的tv原子已支持向量 L1 与矩阵 L2 两种形态；
• 若需对更大图像更快求解，可考虑仓库中同样内置的 Clarabel/COSMO 等一阶求解器（cvxpy/reductions/solvers/conic_solvers）；
• 原始文档位置：doc/source/examples/applications/tv_inpainting.rst，完整示例索引见 doc/source/examples/applications。

• 科学计算