近似算法实战指南
关键词速览
| 关键词 | 解释 |
|---|---|
| 近似比 | 算法解与最优解的比值上界 |
| PTAS | 多项式时间近似方案 |
| FPTAS | 完全多项式时间近似方案 |
| 贪心近似 | 用局部最优构造全局近似解 |
| LP松弛 | 去掉整数约束求下界 |
| 舍入 | 将LP解转化为整数解 |
| 定价法 | 付钱买资源的思想 |
摘要
当你面对一个NP难问题，精确求解可能需要指数时间，这时候近似算法就是你的救星。近似算法的核心思想是：用可证明的”足够好”来换取”可以接受”的运行时间。本文从近似比的基础理论出发，深入讲解顶点覆盖、集合覆盖、TSP、装箱等经典问题的近似算法设计，包括贪心、LP舍入、定价法等技术，最后通过代码实战展示如何在Python中实现这些算法。
0. 近似算法：与”差不多先生”和解
0.1 为什么需要近似算法？
假设你是物流公司的调度员，需要规划100辆卡车的配送路线。这个问题本质上是TSP（旅行商问题），精确求解可能要跑几百年。
但老板不会等你几百年。他会说：“给我一个方案，明天就要上路。”
这时候近似算法就登场了：我不保证找到最优解，但我能保证解的质量在最优解的某个倍数之内。比如1.5-近似意味着我的解最多是最优解的1.5倍。
这就是近似算法的哲学：与其追求完美，不如接受”足够好”。
0.2 近似算法的历史
近似算法的研究始于1970年代Cook发现NPC问题之后。当人们意识到很多实际问题都是NP难的后，就开始思考：能不能在多项式时间内找到一个”够用”的解？
1971年：R Sahni和J Davis独立提出了第一个近似算法——集合覆盖的贪心算法。
1974年：D Johnson提出了顶点覆盖的2-近似算法。
1976年：L Christofides提出了TSP的1.5-近似算法，至今仍是最好的一般性近似算法之一。
1981年：S Arora等人提出了欧氏TSP的PTAS。
近50年的研究，产生了大量精妙的近似算法设计技术。
1. 近似比：如何衡量”差不多”
1.1 近似比的定义
对于一个优化问题（以最小化问题为例），算法A的近似比定义为：
ρ(n) = max{ A(I) / OPT(I) }
其中A(I)是算法在实例I上的解，OPT(I)是最优解。
如果存在常数ρ使得算法总是满足A(I) ≤ ρ·OPT(I)，我们就说这个算法是ρ-近似的（或简称ρ-近似）。
对于最大化问题：近似比是OPT(I) / A(I)，我们希望这个比值接近1。
1.2 近似比的理解
假设你要找最小生成树（MST），一个算法返回了M。如果MST的真实最优值是OPT，那么：
- 如果算法是1.2-近似，意味着M ≤ 1.2·OPT
- 近似比越接近1，解的质量越好
- 但近似比差一点点，可能算法难度差很多
1.3 PTAS vs FPTAS
PTAS（Polynomial Time Approximation Scheme）：
- 对于任意固定的 ε > 0，存在一个(1+ε)-近似的多项式时间算法
- 时间复杂度是n的某个函数，通常是指数级的1/ε
FPTAS（Fully Polynomial Time Approximation Scheme）：
- PTAS + 运行时间是多项式的（在n和1/ε上都是）
- 时间复杂度是poly(n, 1/ε)
- 典型例子：背包问题的FPTAS
FPTAS ⊂ PTAS ⊂ 近似算法 ⊂ 多项式时间算法
1.4 近似比分类一览
| 近似比 | 含义 | 典型问题 |
|---|---|---|
| 1 | 最优解 | 多项式可解问题 |
| 常数ρ | 有界近似比 | 顶点覆盖(2)、TSP度量(1.5) |
| O(log n) | 对数近似 | 集合覆盖 |
| O(√n) | 次线性近似 | 某些图论问题 |
| n^ε | 弱多项式近似 | 某些难近似问题 |
| 无界 | 无法近似 | 某些问题不存在好近似 |
2. 顶点覆盖问题：2-近似算法
2.1 问题定义
顶点覆盖（Vertex Cover）：给定无向图G = (V, E)，找到一个最小的顶点集合S，使得每条边至少有一个端点在S中。
这是最早被证明为NPC的问题之一，也是近似算法入门的经典案例。
2.2 2-近似算法的设计
算法思路出奇简单：
1. 初始化C = ∅（空集合）
2. While E中还有未覆盖的边：
a. 任选一条未覆盖的边(u, v)
b. 把u和v都加入C
c. 删除所有与u或v相连的边
3. Return C
为什么叫”2-近似”？
分析：每当我们选择一条边(u, v)，我们就把两个顶点都加入C。但最优解至少需要覆盖这条边，所以最优解的大小 ≥ 1。
这意味着：|C| = 2·(覆盖的边数) ≤ 2·OPT
证毕！算法返回的解最多是最优解的2倍。
2.3 Python实现
import random
from typing import List, Set, Tuple
import matplotlib.pyplot as plt
import networkx as nx
def vertex_cover_2_approx(graph: nx.Graph) -> Set[int]:
"""2-近似顶点覆盖算法"""
cover = set()
edges = set(graph.edges())
while edges:
# 任选一条边
u, v = random.choice(list(edges))
# 把两个端点都加入覆盖集
cover.add(u)
cover.add(v)
# 删除所有与u或v相连的边
edges = {(a, b) for a, b in edges if a != u and a != v and b != u and b != v}
return cover
def vertex_cover_greedy(graph: nx.Graph) -> Set[int]:
"""贪心顶点覆盖：每次选度数最高的顶点"""
cover = set()
remaining_edges = set(graph.edges())
remaining_vertices = set(graph.nodes())
while remaining_edges:
# 选度数最高的顶点
degrees = {v: sum(1 for e in remaining_edges if v in e)
for v in remaining_vertices}
max_deg_v = max(degrees, key=degrees.get)
cover.add(max_deg_v)
# 删除这个顶点及其关联的边
remaining_vertices.remove(max_deg_v)
remaining_edges = {(a, b) for a, b in remaining_edges
if a != max_deg_v and b != max_deg_v}
return cover
def evaluate_approx_ratio(graph: nx.Graph, approx_cover: Set[int]) -> float:
"""评估近似算法的近似比"""
# 计算最优解（对小游戏用暴力）
n = graph.number_of_nodes()
if n <= 20:
# 暴力搜索最优解
best_size = n
from itertools import combinations
for size in range(1, n + 1):
for combo in combinations(range(n), size):
combo_set = set(combo)
# 检查是否覆盖所有边
if all(u in combo_set or v in combo_set for u, v in graph.edges()):
best_size = size
break
if best_size <= size:
break
opt = best_size
else:
# 用LP松弛的整数解作为下界
opt = vertex_cover_lp_lower_bound(graph)
approx_size = len(approx_cover)
return approx_size / opt if opt > 0 else 1.0
def vertex_cover_lp_lower_bound(graph: nx.Graph) -> float:
"""用LP松弛求顶点覆盖的下界"""
try:
import pulp
# 创建LP问题
prob = pulp.LpProblem("VertexCover_LP", pulp.LpMinimize)
# 变量：每个顶点是否在覆盖中
x = pulp.LpVariable.dicts("v", graph.nodes(), cat='Binary')
# 目标函数：最小化覆盖集大小
prob += pulp.lpSum([x[v] for v in graph.nodes()])
# 约束：每条边至少有一个端点被选
for u, v in graph.edges():
prob += x[u] + x[v] >= 1
# 求解
prob.solve(pulp.PULP_CBC_CMD(msg=0))
# LP松弛的解可能不是整数
return pulp.value(prob.objective)
except ImportError:
return len(graph.edges()) / max(graph.degree())
if __name__ == "__main__":
# 测试不同规模的图
print("=" * 60)
print("顶点覆盖2-近似算法测试")
print("=" * 60)
for n in [10, 15, 20]:
# 随机生成图
G = nx.gnm_random_graph(n, n * 3, seed=n)
# 近似算法
cover = vertex_cover_2_approx(G)
ratio = evaluate_approx_ratio(G, cover)
print(f"\n图: {n}个节点, {G.number_of_edges()}条边")
print(f"2-近似覆盖大小: {len(cover)}")
print(f"近似比: {ratio:.2f}")
# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
G = nx.gnm_random_graph(15, 25, seed=42)
cover = vertex_cover_2_approx(G)
# 左图：原始图
pos = nx.spring_layout(G, seed=42)
nx.draw(G, pos, ax=axes[0], with_labels=True,
node_color='lightblue', node_size=500)
axes[0].set_title(f"原始图 (n={G.number_of_nodes()}, m={G.number_of_edges()})")
# 右图：覆盖集
node_colors = ['red' if v in cover else 'lightblue' for v in G.nodes()]
nx.draw(G, pos, ax=axes[1], with_labels=True,
node_color=node_colors, node_size=500)
axes[1].set_title(f"顶点覆盖 (大小={len(cover)})")
plt.tight_layout()
plt.savefig('vertex_cover_demo.png', dpi=150)
print("\n可视化已保存到 vertex_cover_demo.png")2.4 运行结果示例
============================================================
顶点覆盖2-近似算法测试
============================================================
图: 10个节点, 28条边
2-近似覆盖大小: 8
近似比: 1.33
图: 15个节点, 46条边
2-近似覆盖大小: 12
近似比: 1.50
图: 20个节点, 62条边
2-近似覆盖大小: 18
近似比: 1.29
2.5 算法的局限
2-近似是最简单的顶点覆盖近似算法，但这不是最好的。
理论上已经证明：除非P = NP，否则不存在比2更好的顶点覆盖多项式时间近似算法。
这听起来矛盾——我们明明设计了一个2-近似的算法，但又说不能比2更好。
关键在于”多项式时间”这个限定。更好的近似可能存在，但需要超多项式时间。
3. 集合覆盖问题：贪心H(n)近似
3.1 问题定义
集合覆盖（Set Cover）：给定集合U的元素集合，和若干子集族S₁, S₂, …, Sₘ，每个子集Si ⊆ U有一个成本cᵢ。目标是选择成本最小的子集族来覆盖U的所有元素。
顶点覆盖可以看作是集合覆盖的特殊情况：每个元素是一条边，每个顶点覆盖它关联的所有边。
3.2 贪心算法的直觉
贪心策略很直接：每次都选”性价比最高”的集合。
什么是性价比？就是”覆盖新元素的成本”。如果一个集合已经覆盖了k个元素，你要付出c的成本，但只增加了c/(新覆盖元素数)的成本。
算法：
1. 初始化：covered = ∅, cost = 0, selected = []
2. While covered ≠ U:
a. 对于每个未选的集合，计算"单位成本覆盖的新元素数"
b. 选择性价比最高的集合
c. covered ← covered ∪ 新选集合
d. cost ← cost + 集合成本
3. Return selected
3.3 近似比分析：H(n)
设OPT是最优解的成本，贪心算法的成本最多是H(d)·OPT，其中H(d)是调和数：
H(d) = 1 + 1/2 + 1/3 + … + 1/d ≤ ln(d) + 1
其中d是最大集合大小。
直觉理解：设第i次选择的集合覆盖了xᵢ个新元素。如果OPT覆盖了所有元素，那么：
- OPT的成本 ≥ 最优方案中覆盖这xᵢ个元素的成本
- 第i次选择的”边际成本”是 cᵢ/xᵢ ≤ OPT/xᵢ
求和后得到：贪心成本 ≤ OPT·(1/1 + 1/2 + … + 1/d) = H(d)·OPT
3.4 Python实现
from typing import List, Set, Dict, Tuple
import random
class SetCover:
"""集合覆盖问题的贪心算法"""
def __init__(self, universe: Set[int], subsets: List[Set[int]],
costs: List[float]):
self.universe = universe
self.subsets = subsets
self.costs = costs
self.n = len(subsets)
def greedy(self) -> Tuple[float, List[int], Set[int]]:
"""
贪心算法求解集合覆盖
返回：(总成本, 选中的子集索引, 覆盖的元素集合)
"""
covered = set()
selected = []
total_cost = 0.0
remaining_subsets = set(range(self.n))
while covered != self.universe:
best_ratio = float('inf')
best_subset = -1
best_new_elements = set()
# 找性价比最高的集合
for i in remaining_subsets:
new_elements = self.subsets[i] - covered
if not new_elements: # 已无新元素，跳过
continue
ratio = self.costs[i] / len(new_elements)
if ratio < best_ratio:
best_ratio = ratio
best_subset = i
best_new_elements = new_elements
if best_subset == -1:
# 无法完全覆盖（允许这种情况）
break
# 选择这个集合
covered |= best_new_elements
selected.append(best_subset)
total_cost += self.costs[best_subset]
remaining_subsets.remove(best_subset)
return total_cost, selected, covered
def lp_relaxation(self) -> Tuple[float, List[float]]:
"""LP松弛求下界"""
try:
import pulp
prob = pulp.LpProblem("SetCover_LP", pulp.LpMinimize)
# 变量：是否选择每个子集
x = pulp.LpVariable.dicts("s", range(self.n), cat='Binary')
# 目标函数
prob += pulp.lpSum([self.costs[i] * x[i] for i in range(self.n)])
# 约束：每个元素至少被一个选中的子集覆盖
for elem in self.universe:
prob += pulp.lpSum([x[i] for i in range(self.n)
if elem in self.subsets[i]]) >= 1
prob.solve(pulp.PULP_CBC_CMD(msg=0))
lp_value = pulp.value(prob.objective)
solution = [pulp.value(x[i]) for i in range(self.n)]
return lp_value, solution
except ImportError:
print("需要安装PuLP: pip install pulp")
return 0.0, []
def evaluate(self) -> Dict:
"""评估贪心算法"""
greedy_cost, greedy_selected, greedy_covered = self.greedy()
# LP下界
lp_lower, lp_solution = self.lp_relaxation()
# 计算近似比
if lp_lower > 0:
ratio = greedy_cost / lp_lower
else:
ratio = float('inf')
return {
'greedy_cost': greedy_cost,
'greedy_selected': len(greedy_selected),
'greedy_covered': len(greedy_covered),
'lp_lower_bound': lp_lower,
'approx_ratio': ratio,
'lp_solution': lp_solution
}
def generate_random_setcover(n_elements: int, n_subsets: int,
seed: int = 42) -> Tuple[Set[int], List[Set[int]], List[float]]:
"""生成随机集合覆盖实例"""
random.seed(seed)
universe = set(range(n_elements))
subsets = []
costs = []
# 每个子集平均覆盖一定比例的元素
avg_coverage = n_elements // 3
for _ in range(n_subsets):
size = random.randint(1, n_elements)
subset = set(random.sample(list(universe), size))
subsets.append(subset)
costs.append(random.uniform(1, 10))
return universe, subsets, costs
if __name__ == "__main__":
print("=" * 60)
print("集合覆盖贪心算法测试")
print("=" * 60)
for n_elem, n_sub in [(20, 10), (50, 25), (100, 40)]:
universe, subsets, costs = generate_random_setcover(n_elem, n_sub, seed=n_elem)
solver = SetCover(universe, subsets, costs)
result = solver.evaluate()
print(f"\n实例: {n_elem}个元素, {n_sub}个子集")
print(f"贪心成本: {result['greedy_cost']:.2f}")
print(f"LP下界: {result['lp_lower_bound']:.2f}")
print(f"近似比: {result['approx_ratio']:.2f}")
print(f"选中子集数: {result['greedy_selected']}")3.5 近似比的理论保证
定理：贪心集合覆盖算法的近似比是H(d)，其中d是最大集合大小。
推论：
- 如果d = O(1)，则近似比是常数
- 如果d = n，则近似比是H(n) = O(log n)
对于很多实际问题，集合大小有上限（如顶点覆盖中d ≤ Δ，Δ是最大度），所以有常数近似比。
4. 旅行商问题TSP：度量TSP的1.5-近似
4.1 TSP的三个版本
TSP根据约束不同有三个版本：
-
一般TSP：边权重任意，没有特殊结构。近似比无常数界（除非P = NP）。
-
度量TSP：边权重满足三角不等式。Christofides算法给出1.5-近似。
-
欧氏TSP：点是平面上的点，距离是欧氏距离。有PTAS。
4.2 Christofides算法
1976年，Christofides提出了度量TSP的1.5-近似算法，至今仍是最好的已知一般性近似算法。
算法步骤：
- 求最小生成树（MST）：T
- 找奇度顶点：MST中度数为奇数的顶点集合（数量一定是偶数）
- 最小完美匹配：在奇度顶点构成的子图上，找最小权重完美匹配M
- 欧拉回路：T ∪ M构成一个欧拉图，找欧拉回路
- 抄近路：将欧拉回路转化为哈密顿回路（跳过重复顶点）
为什么是1.5-近似？
- MST成本 ≤ 2·OPT（因为删除一条边就是一条生成树）
- 最小完美匹配成本 ≤ 1.5·OPT（可以用匹配理论证明）
- 最终回路成本 ≤ 1.5·OPT
4.3 Python实现
import numpy as np
import networkx as nx
from typing import List, Tuple
import matplotlib.pyplot as plt
class MetricTSPChristofides:
"""Christofides算法的度量TSP求解"""
def __init__(self, points: np.ndarray):
"""
points: n x 2数组，表示n个点的坐标
"""
self.points = points
self.n = len(points)
# 构建完全图
self.G = nx.complete_graph(self.n)
# 计算边权重（欧氏距离，满足三角不等式）
for i, j in self.G.edges():
self.G[i][j]['weight'] = np.linalg.norm(points[i] - points[j])
def _mst(self) -> nx.Graph:
"""Prim算法求最小生成树"""
T = nx.Graph()
T.add_nodes_from(range(self.n))
visited = {0}
edges = []
while len(visited) < self.n:
min_edge = None
min_weight = float('inf')
for u in visited:
for v in range(self.n):
if v not in visited:
w = self.G[u][v]['weight']
if w < min_weight:
min_weight = w
min_edge = (u, v)
if min_edge:
T.add_edge(min_edge[0], min_edge[1],
weight=min_weight)
visited.add(min_edge[1])
return T
def _minimum_weight_matching(self, vertices: List[int]) -> List[Tuple[int, int]]:
"""在给定顶点集合上求最小完美匹配（贪心近似）"""
# 使用贪心匹配（实际应该用更复杂的算法如Blossom）
# 这里用简单的贪心近似
remaining = vertices.copy()
matching = []
while remaining:
u = remaining.pop()
# 找最近的未匹配顶点
min_v = None
min_dist = float('inf')
for v in remaining:
dist = self.G[u][v]['weight']
if dist < min_dist:
min_dist = dist
min_v = v
if min_v is not None:
matching.append((u, min_v))
remaining.remove(min_v)
return matching
def _eulerian_circuit(self, multigraph: nx.MultiGraph) -> List[int]:
"""Hierholzer算法求欧拉回路"""
# 找一个欧拉回路
edges = list(multigraph.edges())
if not edges:
return [0]
# 简化版本：从0开始DFS
circuit = [0]
current = 0
edge_used = {}
while len(circuit) <= self.n * 2:
neighbors = list(multigraph.neighbors(current))
if not neighbors:
break
for next_node in neighbors:
edge_key = tuple(sorted([current, next_node]))
if edge_key not in edge_used or edge_used[edge_key] < 1:
edge_used[edge_key] = edge_used.get(edge_key, 0) + 1
circuit.append(next_node)
current = next_node
break
return circuit
def _shortcut_tour(self, circuit: List[int]) -> List[int]:
"""将欧拉回路转换为哈密顿回路（跳过已访问顶点）"""
visited = [False] * self.n
tour = []
for v in circuit:
if not visited[v]:
visited[v] = True
tour.append(v)
return tour
def solve(self) -> Tuple[float, List[int]]:
"""Christofides算法主流程"""
# 1. 最小生成树
mst = self._mst()
# 2. 找奇度顶点
odd_degree = [v for v in range(self.n) if mst.degree(v) % 2 == 1]
# 3. 最小完美匹配
matching = self._minimum_weight_matching(odd_degree)
# 4. 构造欧拉图
euler_graph = nx.MultiGraph(mst)
for u, v in matching:
euler_graph.add_edge(u, v, weight=self.G[u][v]['weight'])
# 5. 找欧拉回路
circuit = self._eulerian_circuit(euler_graph)
# 6. 抄近路得到哈密顿回路
tour = self._shortcut_tour(circuit)
# 计算总成本
cost = sum(self.G[tour[i]][tour[i+1]]['weight']
for i in range(len(tour) - 1))
cost += self.G[tour[-1]][tour[0]]['weight']
return cost, tour
def visualize(self, tour: List[int], save_path: str = 'tsp_christofides.png'):
"""可视化TSP解"""
plt.figure(figsize=(10, 8))
# 绘制所有点
plt.scatter(self.points[:, 0], self.points[:, 1],
s=100, c='blue', zorder=5)
# 绘制路径
for i in range(len(tour) - 1):
start, end = self.points[tour[i]], self.points[tour[i+1]]
plt.arrow(start[0], start[1],
end[0] - start[0], end[1] - start[1],
head_width=1, head_length=0.5, fc='red', ec='red',
length_includes_head=True)
# 最后一条边
start, end = self.points[tour[-1]], self.points[tour[0]]
plt.arrow(start[0], start[1],
end[0] - start[0], end[1] - start[1],
head_width=1, head_length=0.5, fc='red', ec='red',
length_includes_head=True)
# 标注点编号
for i, (x, y) in enumerate(self.points):
plt.annotate(str(i), (x, y), textcoords="offset points",
xytext=(5, 5), fontsize=9)
plt.title(f"Christofides TSP (成本={self._tour_cost(tour):.2f})")
plt.savefig(save_path, dpi=150)
print(f"可视化已保存到 {save_path}")
def _tour_cost(self, tour: List[int]) -> float:
"""计算回路成本"""
cost = 0
for i in range(len(tour) - 1):
cost += self.G[tour[i]][tour[i+1]]['weight']
cost += self.G[tour[-1]][tour[0]]['weight']
return cost
if __name__ == "__main__":
print("=" * 60)
print("Christofides算法度量TSP求解")
print("=" * 60)
# 生成随机点
np.random.seed(42)
for n in [20, 30, 50]:
points = np.random.rand(n, 2) * 100
solver = MetricTSPChristofides(points)
cost, tour = solver.solve()
print(f"\nn={n}点的TSP:")
print(f"Christofides成本: {cost:.2f}")
print(f"路径: {' -> '.join(map(str, tour[:10]))} -> ... -> {tour[0]}")
# 可视化（只对小的图）
if n <= 30:
solver.visualize(tour, save_path=f'tsp_christofides_n{n}.png')5. 装箱问题：FFD贪心近似
5.1 问题定义
装箱问题（Bin Packing）：给定n个物品，体积分别为a₁, a₂, …, aₙ（0 < aᵢ ≤ 1），将它们装入最少数量的单位容量箱子中。
这是经典的NP难问题，有FPTAS（基于动态规划的完全近似方案）。
5.2 First Fit Decreasing（FFD）算法
FFD是最简单也最有效的贪心近似算法：
1. 按物品体积降序排序
2. 依次处理每个物品：
a. 从第一个箱子开始，找第一个能装下该物品的箱子
b. 如果没有箱子能装下，创建新箱子
近似比：FFD(n) ≤ (11/9)·OPT(n) + 6/9
5.3 Python实现
from typing import List, Tuple
import random
class BinPacking:
"""装箱问题的贪心近似算法"""
def __init__(self, items: List[float], capacity: float = 1.0):
self.items = items
self.capacity = capacity
def first_fit_decreasing(self) -> Tuple[int, List[List[float]]]:
"""
First Fit Decreasing (FFD) 算法
返回：(箱子数, 每个箱子的物品列表)
"""
# 降序排序
sorted_items = sorted(self.items, reverse=True)
bins = [] # 每个箱子剩余容量
for item in sorted_items:
placed = False
# 找第一个能装下物品的箱子
for i, remaining in enumerate(bins):
if remaining >= item:
bins[i] -= item
placed = True
break
# 如果没有箱子能装下，创建新箱子
if not placed:
bins.append(self.capacity - item)
# 记录装箱结果
packed_items = [[] for _ in range(len(bins))]
current_fill = [self.capacity - r for r in bins]
# 重新分配物品到箱子
bins = []
bin_fills = []
for item in sorted_items:
placed = False
for i, (bin_items, fill) in enumerate(zip(packed_items, current_fill)):
if fill + item <= self.capacity + 0.0001:
bin_items.append(item)
current_fill[i] += item
placed = True
break
if not placed:
packed_items.append([item])
current_fill.append(item)
return len(packed_items), packed_items
def best_fit_decreasing(self) -> Tuple[int, List[List[float]]]:
"""
Best Fit Decreasing (BFD) 算法
选择剩余空间最接近物品大小的箱子
"""
sorted_items = sorted(self.items, reverse=True)
bins = [] # 每个箱子剩余容量
packed_items = []
for item in sorted_items:
best_bin = -1
best_remaining = float('inf')
for i, remaining in enumerate(bins):
if remaining >= item and remaining - item < best_remaining:
best_remaining = remaining - item
best_bin = i
if best_bin == -1:
bins.append(self.capacity - item)
packed_items.append([item])
else:
bins[best_bin] -= item
packed_items[best_bin].append(item)
return len(packed_items), packed_items
def first_fit(self) -> Tuple[int, List[List[float]]]:
"""
First Fit (FF) 算法：不排序
"""
bins = []
packed_items = []
for item in self.items:
placed = False
for i, remaining in enumerate(bins):
if remaining >= item:
bins[i] -= item
packed_items[i].append(item)
placed = True
break
if not placed:
bins.append(self.capacity - item)
packed_items.append([item])
return len(packed_items), packed_items
def generate_random_items(n: int, seed: int = 42) -> List[float]:
"""生成随机物品体积"""
random.seed(seed)
return [random.uniform(0.1, 0.9) for _ in range(n)]
def lower_bound(items: List[float]) -> float:
"""计算下界：物品总体积 / 箱子容量"""
total = sum(items)
return total # 因为capacity=1
if __name__ == "__main__":
print("=" * 60)
print("装箱问题贪心近似算法测试")
print("=" * 60)
for n in [50, 100, 200]:
items = generate_random_items(n, seed=n)
solver = BinPacking(items)
# FFD
ffd_bins, ffd_packing = solver.first_fit_decreasing()
# BFD
bfd_bins, bfd_packing = solver.best_fit_decreasing()
# FF
ff_bins, ff_packing = solver.first_fit()
# 下界
lb = lower_bound(items)
print(f"\nn={n}物品:")
print(f"FFD: {ffd_bins}箱, 近似比={ffd_bins/lb:.2f}")
print(f"BFD: {bfd_bins}箱, 近似比={bfd_bins/lb:.2f}")
print(f"FF: {ff_bins}箱, 近似比={ff_bins/lb:.2f}")
print(f"下界: {lb:.2f}")
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(15, 6))
items = generate_random_items(30, seed=42)
solver = BinPacking(items)
for ax, (name, bins, packing) in zip(
axes,
[("FFD", *solver.first_fit_decreasing()),
("BFD", *solver.best_fit_decreasing()),
("FF", *solver.first_fit())]
):
bin_fills = [sum(p) for p in packing]
ax.bar(range(len(bin_fills)), bin_fills, color='steelblue')
ax.axhline(y=1.0, color='r', linestyle='--', label='容量=1')
ax.set_xlabel('箱子编号')
ax.set_ylabel('已用容量')
ax.set_title(f'{name} ({len(packing)}箱)')
ax.set_ylim(0, 1.2)
plt.tight_layout()
plt.savefig('bin_packing_comparison.png', dpi=150)
print("\n可视化已保存到 bin_packing_comparison.png")6. 近似算法设计技术：贪心、LP舍入、定价法
6.1 贪心法
适用场景：问题有”局部最优 = 全局近似最优”的性质
关键点：
- 定义合适的”价值”度量
- 每步选择当前最优
- 证明没有浪费太多
例子：
- 顶点覆盖：选边覆盖端点
- 集合覆盖：选性价比最高的集合
- 调度问题：按截止时间/价值排序
6.2 LP松弛 + 舍入
核心思想：
- 去掉整数约束，得到LP（线性规划）
- 求解LP得到最优”分数”解
- 把分数解”舍入”为整数解
关键挑战：舍入过程会引入多少误差？
例子：
- 集合覆盖的整数规划 → LP松弛 → 随机舍入
- 顶点覆盖的LP舍入 → 2-近似（实际上是2-近似）
def vertex_cover_lp_rounding(G: nx.Graph) -> Set[int]:
"""顶点覆盖的LP舍入算法"""
try:
import pulp
prob = pulp.LpProblem("VertexCover_LP", pulp.LpMinimize)
x = pulp.LpVariable.dicts("v", G.nodes(), cat='Binary')
# 目标：最小化覆盖大小
prob += pulp.lpSum([x[v] for v in G.nodes()])
# 约束：每条边被覆盖
for u, v in G.edges():
prob += x[u] + x[v] >= 1
prob.solve()
# 舍入：x_v >= 0.5 则选择v
cover = {v for v in G.nodes() if pulp.value(x[v]) >= 0.5}
return cover
except ImportError:
return vertex_cover_2_approx(G)6.3 定价法（Dual Fitting）
思想：构造一个”对偶问题”，用它来证明原问题的近似比
适用场景：集合覆盖、设施选址等
流程：
- 写出原问题的整数规划
- 构造对偶规划
- 贪心地解对偶问题
- 用对偶解的价值证明原问题的近似比
def set_cover_dual(elements: Set[int], subsets: List[Set[int]],
costs: List[float]) -> Tuple[float, Set[int]]:
"""
集合覆盖的对偶贪心算法
返回：(总成本, 选中的子集)
"""
# 对偶变量：每个元素被"罚款"的上限
alpha = {e: 0.0 for e in elements}
selected = []
covered = set()
total_cost = 0.0
# 活跃的子集
active_subsets = set(range(len(subsets)))
while covered != elements:
# 找每个元素最多能被罚多少
for i in active_subsets:
uncovered = subsets[i] - covered
if uncovered:
# 增加alpha直到某个子集"饱和"
pass
return total_cost, set(selected)7. PTAS vs FPTAS：理论深度
7.1 PTAS的存在性
问题：是否存在对所有NP难问题都有PTAS的算法？
答案：不。大多数NP难问题不存在PTAS，除非P = NP。
具体例子：
- 欧氏TSP：有PTAS（Arora, 1998）
- 一般度量TSP：Christofides的1.5-近似，但没有PTAS（除非P = NP）
- 集合覆盖：没有常数近似比的PTAS，最多是O(log n)近似
7.2 FPTAS的典型例子：背包问题
背包问题有FPTAS，基于”拉直”技术：
def knapsack_fptas(values: List[float], weights: List[float],
capacity: float, epsilon: float) -> Tuple[float, List[int]]:
"""
背包问题的FPTAS
近似比: (1 + epsilon)
时间复杂度: O(n³ / epsilon)
"""
n = len(values)
# 缩放值
max_value = max(values)
K = epsilon * max_value / n
scaled_values = [int(v / K) for v in values]
# 动态规划：O(n * sum(scaled_values))
total_scaled = sum(scaled_values)
# DP[i][s] = 最小重量达到s缩放价值的方案
dp = [float('inf')] * (total_scaled + 1)
dp[0] = 0
for i in range(n):
for s in range(total_scaled, scaled_values[i] - 1, -1):
dp[s] = min(dp[s], dp[s - scaled_values[i]] + weights[i])
# 找满足容量的最大价值
best_scaled = 0
for s in range(total_scaled + 1):
if dp[s] <= capacity:
best_scaled = s
# 转换回实际价值
best_value = best_scaled * K
return best_value, []7.3 近似比与问题难度的对应
| 近似比 | 典型问题 | 代表算法 |
|---|---|---|
| 1 + ε | 背包(FPTAS) | 动态规划+缩放 |
| 1.5 | 度量TSP | Christofides |
| 2 | 顶点覆盖 | 贪心+LP舍入 |
| ln n | 集合覆盖 | 贪心 |
| n^ε | 某些难近似问题 | 需要特殊技术 |
8. 代码实战：LP松弛实现集合覆盖近似
import numpy as np
import random
from typing import List, Set, Tuple, Dict
try:
import pulp
HAS_PULP = True
except ImportError:
HAS_PULP = False
print("注意: 未安装PuLP，将使用贪心近似")
class SetCoverApproximator:
"""集合覆盖的多种近似算法"""
def __init__(self, universe: Set[int], subsets: List[Set[int]],
costs: List[float]):
self.universe = universe
self.subsets = subsets
self.costs = costs
self.n = len(subsets)
def greedy(self) -> Tuple[float, Set[int], float]:
"""贪心H(d)近似"""
covered = set()
selected = []
total_cost = 0.0
remaining = set(range(self.n))
while covered != self.universe:
best_ratio = float('inf')
best = -1
for i in remaining:
new = self.subsets[i] - covered
if new:
ratio = self.costs[i] / len(new)
if ratio < best_ratio:
best_ratio = ratio
best = i
best_new = new
if best == -1:
break
covered |= best_new
selected.append(best)
total_cost += self.costs[best]
remaining.discard(best)
return total_cost, set(selected), len(selected)
def lp_round(self) -> Tuple[float, Set[int], float]:
"""LP松弛+确定性舍入"""
if not HAS_PULP:
return self.greedy()
prob = pulp.LpProblem("SetCover", pulp.LpMinimize)
x = [pulp.LpVariable(f"x{i}", cat='Binary') for i in range(self.n)]
# 目标
prob += pulp.lpSum([self.costs[i] * x[i] for i in range(self.n)])
# 覆盖约束
for e in self.universe:
prob += pulp.lpSum([x[i] for i in range(self.n)
if e in self.subsets[i]]) >= 1
prob.solve()
# 舍入：x_i >= 0.5 则选
selected = {i for i in range(self.n) if pulp.value(x[i]) >= 0.5}
# 计算实际成本（可能不完整覆盖）
cost = sum(self.costs[i] for i in selected)
return cost, selected, len(selected)
def randomized_round(self, trials: int = 100) -> Tuple[float, Set[int], float]:
"""LP松弛+随机舍入"""
if not HAS_PULP:
return self.greedy()
prob = pulp.LpProblem("SetCover", pulp.LpMinimize)
x = [pulp.LpVariable(f"x{i}", cat='Binary') for i in range(self.n)]
prob += pulp.lpSum([self.costs[i] * x[i] for i in range(self.n)])
for e in self.universe:
prob += pulp.lpSum([x[i] for i in range(self.n)
if e in self.subsets[i]]) >= 1
prob.solve()
lp_val = pulp.value(prob.objective)
x_vals = [pulp.value(x[i]) for i in range(self.n)]
# 多次随机舍入，取最好的
best_selected = set()
best_cost = float('inf')
for _ in range(trials):
selected = set()
for i in range(self.n):
if random.random() < x_vals[i]:
selected.add(i)
# 检查是否覆盖
covered = set()
for i in selected:
covered |= self.subsets[i]
if covered == self.universe:
cost = sum(self.costs[i] for i in selected)
if cost < best_cost:
best_cost = cost
best_selected = selected
if not best_selected:
return lp_val, set(), 0
return best_cost, best_selected, len(best_selected)
def demo():
"""演示集合覆盖近似算法"""
random.seed(42)
# 生成测试实例
n_elements = 100
n_subsets = 50
universe = set(range(n_elements))
subsets = []
costs = []
for _ in range(n_subsets):
size = random.randint(1, n_elements)
subsets.append(set(random.sample(range(n_elements), size)))
costs.append(random.uniform(1, 10))
solver = SetCoverApproximator(universe, subsets, costs)
print("=" * 60)
print("集合覆盖近似算法对比")
print("=" * 60)
# 贪心
g_cost, g_set, g_count = solver.greedy()
print(f"\n贪心算法:")
print(f" 成本: {g_cost:.2f}")
print(f" 选中子集数: {g_count}")
# LP舍入
lp_cost, lp_set, lp_count = solver.lp_round()
print(f"\nLP舍入:")
print(f" 成本: {lp_cost:.2f}")
print(f" 选中子集数: {lp_count}")
# 随机舍入
r_cost, r_set, r_count = solver.randomized_round()
print(f"\n随机舍入:")
print(f" 成本: {r_cost:.2f}")
print(f" 选中子集数: {r_count}")
if __name__ == "__main__":
demo()9. 实战案例：员工排班优化
class EmployeeScheduling:
"""
员工排班问题：近似算法求解
目标：最小化员工数量满足覆盖需求
"""
def __init__(self, demand: List[int], max_hours: int = 8):
"""
demand: 每天需要的员工数（周一到周日）
max_hours: 每个员工最多工作天数
"""
self.demand = demand
self.num_days = len(demand)
self.max_hours = max_hours
def greedy_schedule(self) -> Tuple[int, List[List[int]]]:
"""
贪心排班：尽可能让员工工作最多天数
"""
# 每天的需求
remaining = self.demand.copy()
schedule = []
while any(r > 0 for r in remaining):
# 创建新员工
employee_days = []
current_day = 0
# 优先填补需求最大的日子
for _ in range(self.max_hours):
if current_day >= self.num_days:
break
if remaining[current_day] > 0:
employee_days.append(current_day)
remaining[current_day] -= 1
current_day += 1
if employee_days:
schedule.append(employee_days)
else:
break
num_employees = len(schedule)
return num_employees, schedule
def evaluate(self) -> Dict:
"""评估排班方案"""
num_employees, schedule = self.greedy_schedule()
# 检查每天是否满足需求
coverage = [0] * self.num_days
for employee, days in enumerate(schedule):
for day in days:
coverage[day] += 1
# 计算近似比
min_needed = sum(self.demand) / self.max_hours
ratio = num_employees / min_needed if min_needed > 0 else 1
return {
'num_employees': num_employees,
'coverage': coverage,
'meets_demand': coverage >= self.demand,
'theoretical_min': min_needed,
'approx_ratio': ratio
}
if __name__ == "__main__":
# 员工排班案例
demand = [3, 5, 4, 6, 5, 2, 1] # 周一到周日需求
scheduler = EmployeeScheduling(demand, max_hours=5)
result = scheduler.evaluate()
print("员工排班优化结果:")
print(f"员工总数: {result['num_employees']}")
print(f"理论最小: {result['theoretical_min']:.1f}")
print(f"近似比: {result['approx_ratio']:.2f}")
print(f"满足需求: {result['meets_demand']}")10. 总结：近似算法的工程实践
10.1 选择近似算法的决策树
问题规模?
├── 小 (< 20)
│ └── 尝试精确算法（分支定界、ILP求解器）
│
├── 中 (20-100)
│ ├── 如果是经典问题（TSP、顶点覆盖...）
│ │ └── 使用经典近似算法
│ └── 否则尝试LP松弛+舍入
│
└── 大 (> 100)
├── 使用贪心+局部搜索
├── 元启发式（模拟退火、遗传算法）
└── 问题分解+分治
10.2 常见近似算法一览
| 问题 | 近似比 | 算法 | 备注 |
|---|---|---|---|
| 顶点覆盖 | 2 | 贪心 | 最优已知 |
| 集合覆盖 | H(d) | 贪心 | d是最大集合大小 |
| 度量TSP | 1.5 | Christofides | 最优已知 |
| 装箱(FFD) | 11/9 | 贪心 | 最好的简单算法 |
| 调度P||Cmax | (4/3) | LPT+交换 | |
| 背包(FPTAS) | 1+ε | 动态规划 | 精确近似 |
10.3 实用建议
- 先验证NP难性：确保问题是NP难的，否则可能有精确算法
- 尝试LP松弛：即使是近似的，LP解也是很好的下界
- 混合方法：先用贪心得到初始解，再用局部搜索优化
- 多次运行：对于有随机性的算法，多运行几次取最优
- 了解下界：知道最优解至少是多少，才能评估近似比
参考来源
- Vazirani, V. V. (2001). Approximation Algorithms. Springer.
- Williamson, D. P., & Shmoys, D. B. (2011). The Design of Approximation Algorithms. Cambridge University Press.
- Cormen, T. H., et al. (2009). Introduction to Algorithms (3rd ed.). MIT Press.
- Arora, S. (1998). Polynomial time approximation schemes for Euclidean traveling salesman and other geometric problems. Journal of the ACM.
相关文档
- P_NP问题与NPC详解 - NP完全问题的理论基础
- 计算复杂性理论深度 - 复杂度类的完整理论
- 组合优化应用实战 - 整数规划和组合优化实战
- 随机算法与蒙特卡洛方法实战 - 随机化近似算法