二次型规划实战：有效集法从入门到精通（附Python代码实现）

在优化算法的世界里，二次型规划（Quadratic Programming, QP）占据着一个独特而核心的位置。它不仅是许多机器学习模型（如支持向量机）的基石，更是金融投资组合优化、工程控制系统设计等领域的常见数学工具。对于已经掌握线性代数和基础微积分，却苦于如何将书本上的“KKT条件”、“有效集”等概念转化为实际代码的开发者来说，理论与实践之间的那道鸿沟常常令人望而却步。本文正是为你准备的桥梁。

我们将彻底抛开繁琐的纯理论推导，转而采用一种“手把手”的实战模式。想象一下，你面前有一个Jupyter Notebook，我们将一起用NumPy从零开始，一步步构建有效集法（Active Set Method）的求解器。你会看到约束条件如何被动态管理，步长如何精确计算，最优性条件又如何被程序化地检验。更重要的是，我们最终会用SciPy这样的工业级优化库来验证自己“手搓”的算法是否正确，这种从无到有、再与权威对比的过程，能带来最扎实的理解和信心。无论你是想深入优化算法内核的学生，还是需要在项目中自定义求解器的工程师，这篇文章都将提供一条清晰的路径。

1. 二次型规划与有效集法：核心思想重塑

在深入代码之前，我们有必要重新梳理一下二次型规划问题的本质及其求解思路。一个标准的二次型规划问题通常形如：

最小化： f(x) = 1/2 * x^T * G * x + c^T * x 满足： A * x = b (等式约束) 以及： C * x >= d (不等式约束)

其中，G 是一个对称矩阵（通常要求半正定以保证凸性），x 是我们的决策变量。有效集法的智慧，在于它巧妙地化繁为简。它不试图一次性处理所有约束，而是像一位经验丰富的侦探，在每次迭代中，只聚焦于当前最可能“起作用”的那些线索——即有效约束（Active Constraints）。

注意：所谓“有效约束”，是指在当前解 x_k 处，取等号的那些不等式约束（C_i * x_k = d_i）以及所有的等式约束。它们像墙壁一样，限制了搜索方向的自由。

有效集法的核心迭代思想可以概括为以下三步循环：

1. 识别战场：在当前点 x_k，确定所有有效约束的集合，记为 A_k。
2. 子问题求解：假设这些有效约束在下一步依然有效（即视为等式约束），在此前提下求解一个简化后的二次规划子问题（通常只有等式约束，易于求解）。这个子问题的解会给出一个候选的搜索方向 p_k。
3. 谨慎推进：
   • 如果沿 p_k 移动一个单位步长（α=1）后，新点 x_k + p_k 仍然满足所有原始约束，且目标函数下降，我们就接受这个移动。
   • 否则，我们需要计算一个更短的步长 α < 1，使得移动刚好触及一个新的约束边界（这个新约束随后会被加入有效集），或者发现当前点已经是最优点。

这个过程的精妙之处在于，有效集像一个动态的“工作清单”，随着迭代不断更新（增加新触及的约束或移除不再起作用的约束），直到找到满足所有原始问题最优性条件（KKT条件）的点。

2. 构建我们的求解器：Python核心模块拆解

现在，让我们进入实战环节。我们将不使用任何现成的优化求解器（除了最后的验证），而是完全依赖NumPy进行线性代数运算。我们将算法分解为几个核心函数模块，每个模块对应算法的一个逻辑步骤。

首先，定义问题。我们用一个类来封装二次型规划问题的所有参数，便于管理。

import numpy as np

class QuadraticProgram:
    """
    二次型规划问题类
    最小化: 1/2 * x.T @ G @ x + c.T @ x
    约束: A_eq @ x = b_eq
         A_ineq @ x >= b_ineq
    """
    def __init__(self, G, c, A_eq=None, b_eq=None, A_ineq=None, b_ineq=None):
        self.G = np.array(G, dtype=float)  # 二次项系数矩阵 (n x n)
        self.c = np.array(c, dtype=float).flatten()  # 一次项系数向量 (n,)
        self.n = self.G.shape[0]  # 决策变量维度

        # 处理等式约束（可以为空）
        self.A_eq = np.array(A_eq, dtype=float) if A_eq is not None else np.empty((0, self.n))
        self.b_eq = np.array(b_eq, dtype=float).flatten() if b_eq is not None else np.empty(0)

        # 处理不等式约束（可以为空）
        self.A_ineq = np.array(A_ineq, dtype=float) if A_ineq is not None else np.empty((0, self.n))
        self.b_ineq = np.array(b_ineq, dtype=float).flatten() if b_ineq is not None else np.empty(0)

    def objective(self, x):
        """计算目标函数值"""
        return 0.5 * x.T @ self.G @ x + self.c.T @ x

接下来是最关键的部分：求解等式约束二次规划子问题。当给定一个有效集（即一组被视为等式约束的索引）时，我们需要求解对应的子问题。这可以通过求解其KKT系统来实现。

def solve_equality_qp(G, c, A, b):
    """
    求解等式约束二次规划子问题。
    最小化 1/2 x.T @ G @ x + c.T @ x
    满足 A @ x = b

    参数:
        G: 二次项矩阵 (n x n)
        c: 一次项向量 (n,)
        A: 等式约束矩阵 (m x n)
        b: 等式约束右端项 (m,)

    返回:
        x_opt: 最优解 (n,)
        lambda_opt: 对应的拉格朗日乘子 (m,)
    """
    n = G.shape[0]
    m = A.shape[0]

    # 构建KKT系统矩阵
    # [ G  A.T ] [ x ] = [ -c ]
    # [ A   0  ] [ λ ]   [  b ]
    KKT_top = np.hstack([G, A.T])
    KKT_bot = np.hstack([A, np.zeros((m, m))])
    KKT_matrix = np.vstack([KKT_top, KKT_bot])

    KKT_rhs = np.hstack([-c, b])

    # 求解线性方程组
    try:
        solution = np.linalg.solve(KKT_matrix, KKT_rhs)
    except np.linalg.LinAlgError:
        # 如果矩阵奇异，使用最小二乘解（更稳健）
        solution = np.linalg.lstsq(KKT_matrix, KKT_rhs, rcond=None)[0]

    x_opt = solution[:n]
    lambda_opt = solution[n:]
    return x_opt, lambda_opt

这个函数是整个有效集法的心脏。它接收当前问题的 G, c 和由有效集定义的 A, b，直接求解出子问题的最优点 x_opt 和对应的拉格朗日乘子 lambda_opt。拉格朗日乘子的符号至关重要，它将用于判断最优性。

3. 有效集法主循环：步长计算与集合管理

有了求解子问题的能力，我们就可以构建主算法循环了。算法的每一步都需要处理几个关键逻辑：计算搜索方向、确定最大可行步长、更新迭代点以及管理有效集。

def active_set_method(qp, x0, max_iter=100, tol=1e-8):
    """
    有效集法主函数。

    参数:
        qp: QuadraticProgram 问题实例
        x0: 初始可行点 (n,)
        max_iter: 最大迭代次数
        tol: 收敛容差

    返回:
        x_opt: 找到的最优解
        active_set_hist: 有效集历史（用于调试）
    """
    x_k = np.array(x0, dtype=float).flatten()
    n = qp.n
    # 合并所有约束，便于索引
    A_full = np.vstack([qp.A_eq, qp.A_ineq]) if qp.A_ineq.size > 0 else qp.A_eq
    b_full = np.hstack([qp.b_eq, qp.b_ineq]) if qp.A_ineq.size > 0 else qp.b_eq
    num_constraints = A_full.shape[0]
    num_eq = qp.A_eq.shape[0]

    # 初始有效集：所有等式约束 + 在x0处紧致的不等式约束
    active_set = set(range(num_eq))  # 等式约束永远有效
    for i in range(num_eq, num_constraints):
        if abs(A_full[i] @ x_k - b_full[i]) < tol:
            active_set.add(i)

    active_set_hist = [active_set.copy()]
    history = {'x': [x_k.copy()], 'obj': [qp.objective(x_k)]}

    for iter in range(max_iter):
        # 步骤1: 基于当前有效集构建子问题并求解
        A_active = A_full[list(active_set), :]
        b_active = b_full[list(active_set)]
        p_k, lambda_k = solve_equality_qp(qp.G, qp.c, A_active, b_active)

        # 搜索方向是子问题最优解与当前点的差： p = x* - x_k
        p_k = p_k - x_k

        # 检查收敛性：如果搜索方向非常小，则认为当前点可能是最优解
        if np.linalg.norm(p_k) < tol:
            # 步骤2: 检查不等式约束的拉格朗日乘子符号
            # 我们需要将lambda_k映射回原始约束索引
            # 这里简化处理：只检查不等式约束对应的乘子（索引在num_eq之后）
            is_optimal = True
            lambda_dict = {list(active_set)[j]: lambda_k[j] for j in range(len(active_set))}
            for idx in active_set:
                if idx >= num_eq:  # 这是一个不等式约束
                    if lambda_dict[idx] < -tol:  # 乘子为负，违反最优性条件
                        is_optimal = False
                        # 移除乘子最小的约束
                        min_lambda_idx = min(
                            [idx for idx in active_set if idx >= num_eq],
                            key=lambda i: lambda_dict[i]
                        )
                        active_set.remove(min_lambda_idx)
                        break
            if is_optimal:
                print(f"在 {iter+1} 次迭代后收敛。")
                break
        else:
            # 步骤3: 计算沿方向p_k的最大可行步长
            alpha_max = 1.0
            blocking_constraint = None
            # 只考虑当前不在有效集中的不等式约束
            for i in set(range(num_eq, num_constraints)) - active_set:
                a_i = A_full[i]
                b_i = b_full[i]
                denom = a_i @ p_k
                if denom < -tol:  # 分母为负，移动可能违反约束
                    alpha_i = (b_i - a_i @ x_k) / denom
                    if 0 < alpha_i < alpha_max:
                        alpha_max = alpha_i
                        blocking_constraint = i

            # 步骤4: 更新迭代点
            alpha_k = min(1.0, alpha_max)
            x_new = x_k + alpha_k * p_k
            history['x'].append(x_new.copy())
            history['obj'].append(qp.objective(x_new))

            # 步骤5: 更新有效集
            if alpha_k < 1.0 - tol:
                # 步长被一个约束阻挡，将该约束加入有效集
                if blocking_constraint is not None:
                    active_set.add(blocking_constraint)
            x_k = x_new

        active_set_hist.append(active_set.copy())
        if iter == max_iter - 1:
            print("达到最大迭代次数，可能未收敛。")

    return x_k, history, active_set_hist

这段代码实现了算法的完整逻辑流。其中几个细节值得关注：

• 拉格朗日乘子检验：当搜索方向 p_k 很小时，我们检查有效集中不等式约束对应的拉格朗日乘子。根据KKT条件，它们必须非负。如果发现负的乘子，说明对应的约束并不“活跃”在全局最优解中，我们将其从有效集中移除。
• 步长计算：公式 α_i = (b_i - a_i^T x_k) / (a_i^T p_k) 计算了对于第 i 个不在有效集中的约束，从当前点 x_k 沿方向 p_k 移动到该约束边界所需的步长。我们取所有正的 α_i 中的最小值作为 α_max，确保移动后不违反任何约束。
• 有效集更新：如果使用了 α_max（即 α_k < 1），说明我们碰到了一个“阻挡”约束，需要将其加入有效集。如果步长为1，则有效集保持不变。

4. 实战演练：从简单例子到复杂案例

理论说得再多，不如运行一段代码看得真切。让我们用一个经典的例子来测试我们的算法。考虑以下凸二次规划问题：

最小化： f(x) = (x1 - 1)^2 + (x2 - 2.5)^2 满足： -x1 + 2x2 <= 2 x1 + 2x2 <= 6 x1 - 2x2 <= 2 x1 >= 0 x2 >= 0

首先，我们需要将其转化为标准形式 1/2 * x^T G x + c^T x。展开目标函数： f(x) = x1^2 - 2x1 + 1 + x2^2 - 5x2 + 6.25 = x1^2 + x2^2 - 2x1 -5x2 + 7.25 常数项不影响优化，可以忽略。因此： G = [[2, 0], [0, 2]], c = [-2, -5]

不等式约束（已包含 x>=0）可以写成 A_ineq * x >= b_ineq 的形式。注意我们的函数要求是 >=，所以需要调整符号： -x1 + 2x2 <= 2 -> x1 - 2x2 >= -2 x1 + 2x2 <= 6 -> -x1 - 2x2 >= -6 x1 - 2x2 <= 2 -> -x1 + 2x2 >= -2 x1 >= 0 -> x1 >= 0 x2 >= 0 -> x2 >= 0

现在，我们用代码定义并求解它。

# 定义问题参数
G = np.array([[2, 0],
              [0, 2]], dtype=float)
c = np.array([-2, -5], dtype=float)

# 不等式约束 A_ineq * x >= b_ineq
A_ineq = np.array([[1, -2],   # 对应 -x1 + 2x2 <= 2
                   [-1, -2],  # 对应 x1 + 2x2 <= 6
                   [-1, 2],   # 对应 x1 - 2x2 <= 2
                   [1, 0],    # 对应 x1 >= 0
                   [0, 1]])   # 对应 x2 >= 0
b_ineq = np.array([-2, -6, -2, 0, 0], dtype=float)

# 无等式约束
A_eq = None
b_eq = None

# 创建问题实例
qp_example = QuadraticProgram(G, c, A_eq, b_eq, A_ineq, b_ineq)

# 选择一个初始可行点，例如原点 (0, 0) 通常是可行的
x0 = np.array([0.0, 0.0])

# 调用有效集法求解
x_opt, history, active_hist = active_set_method(qp_example, x0, max_iter=50, tol=1e-10)

print("自定义有效集法求得的最优解：", x_opt)
print("最优目标函数值：", qp_example.objective(x_opt))
print("迭代历史中的目标值：", history['obj'])

运行这段代码，你会看到算法在几次迭代后收敛。为了更直观地理解算法的搜索路径，我们可以用Matplotlib绘制等高线图和迭代轨迹。

import matplotlib.pyplot as plt

# 绘制目标函数等高线及约束区域
x1 = np.linspace(-0.5, 3, 400)
x2 = np.linspace(-0.5, 3, 400)
X1, X2 = np.meshgrid(x1, x2)
F = 0.5 * (2*X1**2 + 2*X2**2) + (-2)*X1 + (-5)*X2  # 忽略常数项

fig, ax = plt.subplots(figsize=(10, 8))
# 绘制等高线
contour = ax.contour(X1, X2, F, levels=20, colors='gray', alpha=0.5)
ax.clabel(contour, inline=True, fontsize=8)

# 绘制不等式约束线
# 约束1: -x1 + 2x2 <= 2 -> x2 <= (2 + x1)/2
x1_line = np.array([-0.5, 3])
ax.plot(x1_line, (2 + x1_line)/2, 'b-', label='-x1+2x2=2', linewidth=2)
# 约束2: x1 + 2x2 <= 6 -> x2 <= (6 - x1)/2
ax.plot(x1_line, (6 - x1_line)/2, 'r-', label='x1+2x2=6', linewidth=2)
# 约束3: x1 - 2x2 <= 2 -> x2 >= (x1 - 2)/2
ax.plot(x1_line, (x1_line - 2)/2, 'g-', label='x1-2x2=2', linewidth=2)
# 非负约束
ax.axvline(x=0, color='k', linestyle='--', alpha=0.5)
ax.axhline(y=0, color='k', linestyle='--', alpha=0.5)

# 填充可行域（近似）
# 这里需要一点技巧来定义多边形顶点，本例中可行域是一个凸多边形
# 顶点可以通过求解约束交点得到，此处简化为已知
feasible_vertices = np.array([[0,0], [0,1], [2,2], [4,1], [2,0]])
feasible_poly = plt.Polygon(feasible_vertices, alpha=0.2, color='y')
ax.add_patch(feasible_poly)

# 绘制算法迭代路径
x_history = np.array(history['x'])
ax.plot(x_history[:, 0], x_history[:, 1], 'ko-', markerfacecolor='red', markersize=8, linewidth=2, label='迭代路径')
ax.plot(x0[0], x0[1], 's', markersize=12, color='blue', label='初始点')
ax.plot(x_opt[0], x_opt[1], '*', markersize=15, color='gold', label='最优解')

ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_title('有效集法求解二次型规划迭代轨迹')
ax.set_xlim([-0.5, 3])
ax.set_ylim([-0.5, 3])
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.show()

这张图会清晰地展示算法如何从初始点出发，沿着目标函数下降的方向，在约束边界上“拐弯”，最终到达最优解。有效集的变化（在哪些边界上移动）也一目了然。

5. 验证与对比：用SciPy确保算法正确性

自己实现的算法，其正确性必须经过严格验证。最直接的方法是与一个被广泛认可的优化库进行结果对比。Python的SciPy库提供了强大的 minimize 函数，可以求解二次规划问题。我们将使用它来验证我们的解。

首先，需要将我们的问题转化为SciPy接受的格式。SciPy的 minimize 函数处理的是 <= 形式的不等式约束，所以我们需要再次调整符号。

from scipy.optimize import minimize, LinearConstraint, Bounds
import warnings
warnings.filterwarnings('ignore')  # 忽略一些可能的警告

# 定义SciPy格式的目标函数（注意：SciPy最小化 f(x)，而我们形式是 1/2 x^T G x + c^T x）
def objective_scipy(x):
    return 0.5 * x @ qp_example.G @ x + qp_example.c @ x

# 定义梯度（对于二次函数，梯度是 Gx + c）
def gradient_scipy(x):
    return qp_example.G @ x + qp_example.c

# 定义约束： SciPy使用 A_ub @ x <= b_ub 和 A_eq @ x = b_eq
# 我们的 A_ineq 是 >=，所以需要取负号
A_ub_scipy = -qp_example.A_ineq  # 因为我们的 A_ineq @ x >= b_ineq 等价于 -A_ineq @ x <= -b_ineq
b_ub_scipy = -qp_example.b_ineq

# 变量边界 (x>=0)
bounds_scipy = Bounds([0, 0], [np.inf, np.inf])  # 下界为0，上界无穷

# 调用SciPy的SLSQP或trust-constr求解器
result = minimize(objective_scipy, x0, method='trust-constr',
                  jac=gradient_scipy,
                  constraints=[LinearConstraint(A_ub_scipy, -np.inf, b_ub_scipy)],
                  bounds=bounds_scipy,
                  options={'verbose': 0, 'gtol': 1e-10, 'xtol': 1e-10})

print("\n--- SciPy 验证结果 ---")
print("SciPy 求得的最优解：", result.x)
print("SciPy 最优目标函数值：", result.fun)
print("自定义算法与SciPy解的差异（范数）：", np.linalg.norm(x_opt - result.x))
print("目标函数值差异：", abs(qp_example.objective(x_opt) - result.fun))

如果我们的算法实现正确，那么两个解之间的差异应该在数值容差范围内（例如 1e-6 或更小）。这种对比不仅能验证结果的正确性，还能让我们对算法的数值稳定性有更深的了解。

为了更全面地评估，我们可以设计一个更复杂的测试案例，比如一个更高维度（例如10维）的随机凸二次规划问题，并比较两种方法的求解时间和精度。

import time

def test_random_qp(n=10, m_ineq=20):
    """生成并求解一个随机凸二次规划问题"""
    np.random.seed(42)  # 固定随机种子以便复现
    # 生成随机正定矩阵 G
    A = np.random.randn(n, n)
    G = A.T @ A + 0.1 * np.eye(n)  # 使其强凸
    c = np.random.randn(n)

    # 生成随机不等式约束 A_ineq @ x >= b_ineq
    A_ineq = np.random.randn(m_ineq, n)
    # 为了确保存在可行域，我们生成一个可行中心点x_center，并令 b_ineq = A_ineq @ x_center - abs(噪声)
    x_center = np.random.randn(n)
    b_ineq = A_ineq @ x_center - np.abs(np.random.randn(m_ineq))

    # 无等式约束
    qp_random = QuadraticProgram(G, c, A_ineq=A_ineq, b_ineq=b_ineq)

    # 初始点选择为 x_center（确保可行）
    x0 = x_center.copy()

    # 使用自定义有效集法求解
    print(f"\n=== 测试 {n} 维随机QP问题（{m_ineq}个不等式约束）===")
    start_time = time.time()
    x_opt_custom, hist_custom, _ = active_set_method(qp_random, x0, max_iter=200, tol=1e-12)
    custom_time = time.time() - start_time
    print(f"自定义算法耗时: {custom_time:.4f} 秒")
    print(f"最终目标值: {qp_random.objective(x_opt_custom):.6e}")

    # 使用SciPy验证
    # 转换约束格式
    A_ub_sp = -A_ineq
    b_ub_sp = -b_ineq
    bounds_sp = Bounds([-np.inf]*n, [np.inf]*n)  # 无显式边界

    def obj_sp(x): return 0.5 * x @ G @ x + c @ x
    def grad_sp(x): return G @ x + c

    start_time = time.time()
    result_sp = minimize(obj_sp, x0, method='trust-constr',
                         jac=grad_sp,
                         constraints=[LinearConstraint(A_ub_sp, -np.inf, b_ub_sp)],
                         bounds=bounds_sp,
                         options={'maxiter': 1000, 'verbose': 0, 'gtol': 1e-12})
    sp_time = time.time() - start_time
    print(f"SciPy算法耗时: {sp_time:.4f} 秒")
    print(f"SciPy最终目标值: {result_sp.fun:.6e}")

    diff_x = np.linalg.norm(x_opt_custom - result_sp.x)
    diff_f = abs(qp_random.objective(x_opt_custom) - result_sp.fun)
    print(f"解向量差异范数: {diff_x:.2e}")
    print(f"目标函数值差异: {diff_f:.2e}")
    return custom_time, sp_time, diff_f

# 运行测试
test_random_qp(n=5, m_ineq=10)

通过这样的对比测试，我们不仅能确认算法的正确性，还能对其计算效率有一个初步的认识。对于中小规模、稠密的凸二次规划问题，有效集法通常非常可靠和高效。当然，对于大规模或稀疏问题，可能需要考虑内点法等其他算法。

在实现和测试过程中，你可能会遇到一些典型的“坑”。例如，初始点的选择必须严格可行，否则算法可能无法启动。又比如，在判断约束是否“活跃”时，由于浮点数精度问题，需要使用一个容差（tol）而不是直接判断等于零。还有，当KKT系统矩阵奇异时（例如有效约束线性相关），我们的 solve_equality_qp 函数使用了 np.linalg.lstsq 来求最小二乘解，这是一种稳健的处理方式，但在严格的数学意义上，这种情况下子问题可能有无穷多解，需要更复杂的策略来处理。