罚函数法

罚函数法把违反约束的程度变成目标函数中的附加成本，再求解一串较易处理的问题。它常用于把受约束优化转成无约束或约束更简单的优化，但“惩罚很大”与“已经严格满足原约束”是两回事。

一个一维约束的完整计算

先看原问题 $min x^{2}$ ，条件为 $x \geq 1$ 。可行区间上 $x^{2}$ 随 $x$ 增大，原最优点 $x^{*} = 1$ 、最优值为 1。把违反量记为 $(1 - x)_{+} = max {0, 1 - x}$ ，对惩罚系数 $ρ > 0$ 构造二次罚函数 $F_{ρ} (x) = x^{2} + ρ (1 - x)_{+}^{2} .$ 若 $x \geq 1$ ，罚项为零，此区间的最低值仍为 1。若 $x \leq 1$ ，把 $F_{ρ} = x^{2} + ρ (1 - x)^{2}$ 展开求导，得 ${F_{ρ}}^{'} (x) = 2 (1 + ρ) x - 2 ρ$ ；于是这半区间的唯一极小点为 $x_{ρ} = \frac{ρ}{1 + ρ}, F_{ρ} (x_{ρ}) = \frac{ρ}{1 + ρ}, (1 - x_{ρ})_{+} = \frac{1}{1 + ρ} .$ 因为 $ρ / (1 + ρ) < 1$ ，该点比可行区间里的任一点有更小的罚后目标，所以它就是整个实数轴的罚问题最优点。注意它对任何有限 $ρ$ 都不满足 $x \geq 1$ ；只有当 $ρ \to \infty$ 时才从左侧趋近原最优点。

误差怎样随惩罚增长

在这个例子， $ρ = 1$ 给 $x_{ρ} = 1 / 2$ 、违反量 $1 / 2$ ； $ρ = 9$ 给 $x_{ρ} = 9 / 10$ 、违反量 $1 / 10$ 。要把违反量压到 $0.001$ 以下，就需要 $ρ > 999$ 。惩罚越来越大时，目标曲面在不同方向上的曲率可能相差很悬殊，数值求解会变难；实际计算通常逐步增大惩罚，并同时监测原约束残差和原目标值。

更一般地，若原问题是 $min f (x)$ ，约束 $g_{i} (x) \leq 0$ 、 $h_{j} (x) = 0$ ，二次罚函数可取 $f (x) + ρ \sum_{i} max {0, g_{i} (x)}^{2} + ρ \sum_{j} h_{j} (x)^{2} .$ 只有正的 $g_{i}$ 才违反不等式；等式的正负偏差都要罚。若 $x_{ρ}$ 是罚问题的全局最小点，且原可行最优点 $x^{*}$ 存在，则与 $x^{*}$ 比较可得 $F_{ρ} (x_{ρ}) \leq f (x^{*})$ 。再假设 $f (x) \geq m$ 对所有考虑的点成立，就有 $\sum_{i} max {0, g_{i} (x_{ρ})}^{2} + \sum_{j} h_{j} (x_{ρ})^{2} \leq \frac{f (x^{*}) - m}{ρ} .$ 这给约束违反量一个粗略上界；它依赖全局求解和下界 $m$ ，不能直接替代停机后的残差核对。

在本例中，罚函数在 $x < 1$ 的二阶导数是 $2 (1 + ρ)$ ，在 $x > 1$ 则是 2。惩罚系数越大，边界两侧的曲率差越明显。高维模型会在某些违反约束的方向上变得很陡、另一些方向仍较平，数值优化更难选一个对所有方向都合适的步长。因此逐步提高 $ρ$ 和报告约束残差，比一次设置天文数字的惩罚更可检验。

二次罚与精确罚的区别

仍用原来的 $x \geq 1$ 问题，改成一次罚 $G_{λ} (x) = x^{2} + λ (1 - x)_{+}, λ > 0.$ 在 $x < 1$ ，导数为 $2 x - λ$ ；在 $x > 1$ ，导数为 $2 x$ 。若 $λ \geq 2$ ，左侧函数朝 $x = 1$ 递减，右侧朝远离 1 的方向递增，因此 $x = 1$ 已是有限惩罚系数下的全局最优。若 $0 < λ < 2$ ，最优点为 $x = λ / 2 < 1$ ，仍不满足原约束。这个阈值 2 恰与该问题的最优KKT 乘子相同。一般问题中“一次罚可精确”的结论需要适当的正则性与足够大的系数，不能只把平方换成绝对值就无条件宣布成功。

Python 3 核对

下列代码利用本例已经推导出的解析解，列出二次罚问题在不同系数下的精确最优点。它是本例的计算器，不声称会自动求解任意罚函数。

Python 3

def quadratic_penalty_path(penalties):
    rows = []
    for rho in penalties:
        if rho <= 0:
            raise ValueError("penalty must be positive")
        x = rho / (1.0 + rho)
        violation = max(0.0, 1.0 - x)
        penalized_value = x * x + rho * violation * violation
        rows.append((rho, x, violation, penalized_value))
    return rows


if __name__ == "__main__":
    rows = quadratic_penalty_path([1, 9])
    assert rows[0] == (1, 0.5, 0.5, 0.5)
    assert abs(rows[1][1] - 0.9) < 1e-12
    assert abs(rows[1][2] - 0.1) < 1e-12

对 $m$ 个系数，本程序用 $O (m)$ 时间与 $O (m)$ 输出空间，不改写输入。一般罚函数的内层最小化通常要由数值方法完成；代码中闭式公式只因一维平方模型简单而存在。

参考资料

MIT 15.084J《Penalty and Barrier Methods》，讨论二次罚、精确罚与障碍方法。
MIT 6.252J 第16讲，含二次罚与乘子方法。