算法设计技巧: Rounding

最新推荐文章于 2023-05-04 20:50:59 发布

胡拉哥

最新推荐文章于 2023-05-04 20:50:59 发布

阅读量3.2k

点赞数 4

分类专栏：优化算法文章标签：算法

本文链接：https://blog.csdn.net/qx3501332/article/details/105545915

版权

优化算法专栏收录该内容

21 篇文章 55 订阅

订阅专栏

本文介绍一种近似算法的设计技巧: Rounding. 具体来说, 它有两种实现思路:

Rounding solution. 以整数规划为例, 可以先求其线性规划松弛解, 然后对解进行四舍五入(同时保证可行性), 从而得到整数解.
Rounding instance. 修改原问题实例然后求解, 同时保证解的可行性(对原问题而言). 修改实例的目的是使得问题更容易求解.

背包问题

考虑 $n$ 个商品, 它们的大小为 $w_1, w_2, \ldots, w_n \in \mathbb{N}^+$ , 价值为 $p_1, p_2, \ldots, p_n \in \mathbb{N}^+$ . 给定一个大小为 $W$ 的背包, 我们需要找到一组商品装入背包中(商品总大小不超过 $W$ )使得装入背包的商品总价值达到最大.

设决策变量 $x_i\in \{0,1\}$ 代表商品 $i$ 是否被装入背包中. 背包问题可以写成如下整数规划:
$\begin{aligned} \max & \sum_i p_ix_i \\ \text{s.t. } & \sum_{i=1}^n w_i x_i \leq W \\ & x_i \in \{0,1\}. \end{aligned}$

虽然背包问题可以用动态规划求得最优解, 但从时间复杂度上来看它是指数时间的算法. 当问题规模(encoding length of input)非常大时, 我们可以考虑对它进行近似地求解. 下面从Rounding的角度介绍几种近似求解的思路.

LP-Rounding

考虑上述规划的松弛问题(Linear programming relaxation):
$\begin{aligned} \max & \sum_i p_ix_i \\ \text{s.t. } & \sum_{i=1}^n w_i x_i \leq W \\ & 0 \leq x_i \leq 1. \end{aligned}$

令 $(x_1, x_2, \ldots, x_n)^T$ 代表松弛问题的最优解. 如果 $x$ 是整数解则直接输出, 否则需要把它Rounding成整数解.

Lemma. 存在松弛问题的最优解 $x$ ,它最多包含一个非整数的分量.

证明. 假设存在两个变量 $x_i, x_j\in (0,1)$ . 不失一般性, 不妨假设 $p_i/w_i \geq p_j/w_j$ . 增大 $x_i$ 的同时减小 $x_j$ 直到 $x_i = 1$ 或 $x_j$ 等于0.

因此, 我们得到如下算法:

计算松弛问题的最优解 $x$ ;
保证 $x$ 最多包含一个非整数分量 $x_k$ (单纯形算法直接保证), 然后令 $x_k=0$ ;
令 $\sum_{i}w_i$ 代表背包剩余的空间, 从剩下的商品中挑选一个价值最大且大小不超过 $W^{'}$ 的商品加入背包.

Python实现

from ortools.linear_solver import pywraplp


class KnapsackLPRounding(object):
    """ 背包问题LP Rounding(近似)解法.
    """

    def __init__(self, w, p, W):
        """
        :param w: 物品大小, list
        :param p: 物品价值, list
        :param W: 背包大小, int
        """
        self._w = w
        self._p = p
        self._W = W
        self._n = len(self._w)
        self._result = None

    def _solve_lp(self):
        solver = pywraplp.Solver('MasterModel', pywraplp.Solver.GLOP_LINEAR_PROGRAMMING)
        n = len(self._w)
        # 决策变量
        x = [solver.NumVar(0, 1, 'x[%d]' % i) for i in range(n)]
        # 约束
        ct = solver.Constraint(0, self._W)
        for i in range(n):
            ct.SetCoefficient(x[i], self._w[i])
        # 目标
        obj = solver.Objective()
        for i in range(n):
            obj.SetCoefficient(x[i], self._p[i])
        obj.SetMaximization()
        # 求解
        solver.Solve()
        # 得到计算结果
        return [x[i].solution_value() for i in range(n)]

    def _rounding(self, x):
        """ 对分数解x取整, 然后从剩余商品中挑一个价值最大的商品装入背包(如果可行）
        """
        sol = [i for i in range(len(x)) if abs(x[i] - 1) < 1e-6]
        # 背包剩余的空间
        available_space = self._W - sum([self._w[i] for i in sol])
        # 剩下的物品
        left_over = set(range(self._n)) - set(sol)
        # 剩下物品的价值
        left_over_profits = [self._p[i] for i in left_over]
        # 按从大到小排序
        left_over_items = list(sorted(zip(left_over, left_over_profits), key=lambda item: item[1], reverse=True))
        left_over = [item[0] for item in left_over_items if self._w[item[0]] <= available_space]

        return sol + [left_over[0]] if left_over else sol

    def solve(self):
        # 求解背包问题的松弛解
        x = self._solve_lp()
        # Rounding solution
        self._result = self._rounding(x)
        return self

完整代码

PTAS

PTAS的全称是Polynomial Time Approximation Scheme. 给定参数 $\epsilon>0$ , 它是多项式时间的 $(1+\epsilon)$ 近似算法, 即算法对应的目标函数值与最优目标函数值之差不超过最优值的 $\epsilon$ 倍.

算法-k

枚举所有不超过 $k$ 个商品的集合 $\mathcal{S}$ , 共有 $O(kn^k)$ 种情况.
从 $\mathcal{S}$ 中选择价值最大的集合 $S^*\in \mathcal{S}$ .
用贪心算法填充剩下的空间.

评估

设OPT代表最优解, ALG代表算法解. 可以证明¹: $\text{OPT} / \text{ALG} \leq 1+1/k$ . 时间复杂度是 $O(kn^{k+1})$ . 由于 $k$ 是常数, 该算法是多项式时间. 当参数 $k$ 越大, 算法解与最优解越相近, 但时间复杂度越高. 时间复杂度的增长是关于 $k$ 的指数函数.

Python实现

枚举算法的实现. 分成两个函数考虑

choose_exact(n, k): 从 $n$ 个物品里选 $k$ 个, 枚举所有情况.
choose_at_most(n, k): 从 $n$ 个物品里选至少 $1$ 个至多 $k$ 个, 枚举所有情况.

def choose_exact(n, k):
    """
    从n个物品[0, 1, ... n-1]里选择k个, 枚举所有情况. 例如:
    >>> choose_exact(4, 2)
    [[0, 1], [0, 2], [0, 3], [1, 2], [1, 3], [2, 3]]
    >>> choose_exact(3, 3)
    [[0, 1, 2]]
    """
    if k == 1:
        return [[i] for i in range(n)]
    res_mid = choose_exact(n, k - 1)
    result = []
    for p in res_mid:
        for i in range(max(p)+1, n):
            result.append(p + [i])
    return result


def choose_at_most(n, k):
    """ 从n个物品[0, 1, ..., n-1]里选择至少1个至多k个, 枚举所有情况.
    """
    result = []
    for i in range(1, k+1):
        result += choose_exact(n, i)
    return result

背包问题的PTAS.

class KnapsackPTAS(object):
    """
    背包问题的PTAS.
    近似比 = (1+1/k), 计算时间复杂度 = O(kn^{k+1}).
    """

    def __init__(self, w, p, W):
        """
        :param w: 物品大小, list
        :param p: 物品价值, list
        :param W: 背包大小, int
        """
        self._w = w
        self._p = p
        self._W = W
        self._n = len(self._w)
        self._result = None

    def solve(self, k):
        # step1. 枚举所有不超过k个物品的解, 取价值最大的可行解
        solutions = choose_at_most(self._n, k)
        max_sol = None
        max_p = 0
        for sol in solutions:
            if sum([self._w[i] for i in sol]) > self._W:
                continue
            p = sum([self._p[i] for i in sol])
            if p > max_p:
                max_p = p
                max_sol = sol
        # step2. 背包的剩余空间用贪心算法填充
        self._result = self._greedy(max_sol)
        return self

    def _greedy(self, sol):
        """
        给定可行解sol, 把背包剩余的空间用贪心算法填充.
        """
        
        # 背包剩余的空间
        available_space = self._W - sum([self._w[i] for i in sol])
        new_items = []
        # 剩下的物品
        left_over = set(range(self._n)) - set(sol)
        # 计算剩下物品的性价比
        left_over_values = [self._p[i] / self._w[i] for i in left_over]
        # 按性价比从大到小排序
        left_over_items = sorted(zip(left_over, left_over_values), key=lambda x: x[1], reverse=True)
        left_over = [item[0] for item in left_over_items]
        # 把剩余的物品依次装入背包(如果能装下)
        for i in left_over:
            if self._w[i] <= available_space:
                new_items.append(i)
                available_space -= self._w[i]

        return sol + new_items

完整代码

FPTAS

FPTAS的全称是Fully Polynomial Time Approximation Scheme. 给定 $\epsilon>0$ , 算法解ALG与最优解OPT的关系也满足 $\text{OPT}/\text{ALG} \leq 1+\epsilon$ . 它与PTAS区别的地方在于算法的时间复杂度是关于 $1/\epsilon$ 的多项式函数(PTAS是关于 $1/\epsilon$ 的指数函数).

该算法基于动态规划: 设计原问题的一个动态规划算法，然后修改问题实例, 通过动态规划求解从而得到原问题的解. 修改之后的实例可以保证算法的计算复杂度是多项式时间的.

动态规划

设 $f (i, j)$ 代表背包中物品价值等于 $j$ 所需要的最小的总体积(total weight), 其中物品来自集合 $\{1, 2, \ldots, i\}$ . 令 $P=\max(p_i)$ , 因此 $n P$ 是 $j$ 的上限.

考虑 $\ldots, m$ , $\ldots, nP$ .

递归式

$\begin{cases} \min\{f(i,j), f(i, j - p_{i+1}) + w_{i+1}\}, & \text{ if } p_{i+1}\leq j \\ f(i, j) & \text{ otherwise} \end{cases}$

初始条件
$\begin{aligned} & f(1, p_1) = w_1 \\ & f(1, 0) = 0 \\ & f(i, j) = \infty, \quad \text{for the rest } i, j \end{aligned}$

Python实现

下面是上述动态规划的实现. 更多关于动态规划的标准实践可以参考算法设计技巧: 动态规划 (Dynamic Programming).

import math


class KnapsackDP(object):
    """ 背包问题的动态规划算法.
    """

    def __init__(self, w, p, W):
        """
        :param w: 物品大小, list
        :param p: 物品价值, list
        :param W: 背包大小, int
        """
        self._w = w
        self._p = p
        self._W = W
        self._n = len(self._w)
        self._f = self._init_recurrence_formula()
        self._result = None

    def _init_recurrence_formula(self):
        n = len(self._w)
        f = [[]] * n
        max_p = max(self._p)
        for i in range(n):
            f[i] = [math.inf] * n * max_p
            f[0][0] = 0  # !
            f[0][self._p[0]] = self._w[0]
        return f

    def solve(self):
        n = len(self._w)
        max_p = max(self._p)
        # result_items保存计算的中间结果
        # key = profit, value = 达到此profit所包含的一个item
        result_items = {self._p[0]: 0}  # 初始化
        for i in range(n-1):
            for j in range(n * max_p):
                if self._p[i+1] <= j:
                    self._f[i+1][j] = min(self._f[i][j],
                                          self._f[i][j-self._p[i+1]] + self._w[i+1])
                    if self._f[i][j-self._p[i+1]] + self._w[i+1] < self._f[i][j]:
                        result_items[j] = i+1
                else:
                    self._f[i+1][j] = self._f[i][j]

            self._result = self._get_result(result_items, self._get_profit())
        return self

FPTAS

给定 $\epsilon > 0$ , 令 $\epsilon P / n$ .
令 $p'_i = \lfloor p_i/K\rfloor$ , $\forall i=1,2,\ldots,n$ .
用上述动态规划求解新的实例(对应价值 $p'_i$ ), 然后输出结果.

可以证明¹: $\text{ALG} \geq (1-\epsilon) \text{OPT}$ , 且算法的时间复杂度为 $O(n^2\lfloor n/\epsilon\rfloor)$ .

Python实现

class KnapsackFPTAS(object):
    """ 动态规划FPTAS.
    近似比: ALG >= (1-epsilon)OPT, 时间复杂度 = O(n^2 * floor(n/epsilon))
    """
    def __init__(self, w, p, W):
        """
        :param w: 物品大小, list
        :param p: 物品价值, list
        :param W: 背包大小, int
        """
        self._w = w
        self._p = p
        self._W = W
        self._n = len(self._w)
        self._result = None

    def solve(self, epsilon):
        k = epsilon * max(self._p) / len(self._w)
        p1 = [int(x/k) for x in self._p]
        dp = KnapsackDP(self._w, p1, self._W).solve()
        self._result = dp.get_result()
        return self

完整代码

参考文献

K. Lai and M. X. Goemans. The Knapsack Problem and Fully Polynomial Time Approximation Schemes (FPTAS), lecture notes, 2006. ↩︎ ↩︎

胡拉哥

关注

4
点赞
踩
25

收藏

觉得还不错? 一键收藏
4
评论
算法设计技巧: Rounding

本文介绍一种近似算法的设计技巧: Rounding. 具体来说, 它有两种实现思路:Rounding solution. 以整数规划为例, 可以先求其线性规划松弛解, 然后对解进行四舍五入(同时保证可行性), 从而得到整数解.Rounding instance. 修改原问题实例然后求解, 同时保证解的可行性(对原问题而言). 修改实例的目的是使得问题更容易求解.背包问题考虑nnn个商品...
复制链接

扫一扫