算法设计技巧: Primal-Dual

原始对偶(Primal-Dual)是一种求解优化问题的思想, 在凸优化和组合优化问题中有重要应用. 本文以线性规划问题为例解释原始对偶算法的设计思路, 进而介绍如何设计基于原始对偶的近似算法(一般用于求解NP-hard问题).

互补松弛条件

考虑线性规划的原始问题(P)和对偶问题(D)如下:

$$\begin{array}{rl}\min & c^{T}x\\ \text{s.t. }& Ax\ge b\\ & x\ge 0\end{array}\text{\hspace{1em}(P)}$$

$$\begin{array}{rl}\max & b^{T}y\\ \text{s.t. }& A^{T}y\le c\\ & y\ge 0\end{array}\text{\hspace{1em}(D)}$$

互补松弛条件(Complementary Slackness Condition)

• Primal complementary slackness conditions
  $$x_j>0\Rightarrow [A^{T}y{]}_j=c_j$$
• Dual complementary slackness conditions
  $$y_i>0\Rightarrow [Ax{]}_i=b_i$$
  其中$[\cdot {]}_i$是表达式的第$i$个分量.

当互补松弛条件成立且$x$和$y$分别是原问题(P)和对偶问题(D)的可行解, 那么它们也是对应问题的最优解. 回顾线性规划的单纯形算法(Simplex Algorithm), 它从原问题(P)的一个可行解出发, 在满足互补松弛条件的前提下, 使得其对偶变量朝着对偶可行解的方向迭代. Primal-Dual的思想是从对偶可行解出发, 在满足互补松弛条件的前提下, 使得原始变量朝着可行解的方向迭代. 在经典的组合优化问题例如Maximal Flow, Minimum Cost Flow, Maximum Weight Matching等都有Primal-Dual算法的身影.

下面我们介绍如何利用Primal-Dual的思想设计近似算法.

松弛的互补松弛条件

设$x$, $y$分别为(P)和(D)的可行解. 存在常数$\alpha \ge 1,\beta \ge 1$满足如下条件:

• Relaxed primal complementary slackness
  $$x_j>0\Rightarrow [A^{T}y{]}_j\ge c_j/\alpha$$
• Relaxed dual complementary slackness
  $$y_i>0\Rightarrow [Ax{]}_i\le \beta b_i$$

基于Primal-Dual的近似算法从对偶可行解出发, 始终满足松弛的(Relaxed)互补松弛条件, 且朝着原问题可行解的方向迭代.

由于$x$和$y$分别是(P)和(D)的可行解, OPT为原问题最优解的目标值. 此外, 根据上述 松弛的 互补松弛条件, 我们有
$$c^{T}x\le \alpha (A^{T}y{)}^{T}x=\alpha y^T(Ax)\le \alpha \beta y^{T}b=\alpha \beta b^{T}y\le \alpha \beta \cdot \text{OPT}.$$

换句话说, 上述算法得到的目标值不超过最优值的$\alpha \beta$倍. 当$\alpha \beta$越小, 算法的解与最优解越近.

基于Prima-Dual的近似算法

以最小化问题为例, 基于Primal-Dual的近似算法框架如下:

1. 把问题用整数规划建模, 然后考虑其松弛的线性规划问题(P)以及对应的对偶问题(D).
2. 初始化对偶可行解$y$以及原始问题的不可行解$x$. (一般来说, 初始化时$x=0,y=0$.)
3. 用某种方式增加$y_i$直到某个约束$[Ay{]}_j=c_j/\alpha$且始终保持$y$的可行性. 然后根据对偶约束增加$x_i$的值直到$x$成为原问题(整数规划问题)的可行解.
4. 对偶问题的解是OPT的下界. 算法的近似比是$\alpha \beta$.

(更多细节参考这里)

Set Cover

Set Cover是一个经典的组合优化问题. 我们已知:

• 元素的集合 E = { 1 , 2 , … , n } E=\{1, 2, \ldots, n\} E={1,2,…,n}
• $m$个集合$S_1,S_2,\dots ,S_m$, 其中$S_j\subseteq E$
• 每个集合的费用$c_j$, $j=1,2,\dots ,m$

目标是找到一些集合$S_{j_1},S_{j_2},\dots ,S_{j_k}$使得它们覆盖$E$, 即${\cup }_{i=1}^k{S}_{j_i}=E$, 且总成本最低.

例 如下图所示 E = { 1 , 2 , … , 9 } E=\{1, 2, \ldots, 9\} E={1,2,…,9}, S = { S 1 , S 2 , … , S 6 } \mathcal{S} = \{S_1, S_2, \ldots, S_6\} S={S1​,S2​,…,S6​}. 根据定义, C = { S 1 , S 2 , S 5 , S 6 } \mathcal{C} = \{S_1, S_2, S_5, S_6\} C={S1​,S2​,S5​,S6​}是一个set cover. 当给定每个集合$S_j$的权重时, 我们的目标是要找到一个总成本最低的set cover.

整数规划

令 S = { S 1 , S 2 , … , S m } \mathcal{S} = \{S_1, S_2, \ldots, S_m\} S={S1​,S2​,…,Sm​}. $\forall S\in \mathcal{S}$, 定义决策变量 x S ∈ { 0 , 1 } x_S\in \{0, 1\} xS​∈{0,1}, 代表是否选择$S$. 上述问题可以描述成如下整数规划.
min ⁡   ∑ S ∈ S c S x S s.t.  ∑ S ∋ e x S ≥ 1 , ∀ e ∈ E x S ∈ { 0 , 1 } \begin{aligned} \min\ & \sum_{S\in \mathcal{S}}c_Sx_S \\ \text{s.t. } & \sum_{S\ni e} x_S \geq 1, \quad \forall e\in E \\ & x_S \in \{0, 1\} \end{aligned} min s.t. ​S∈S∑​cS​xS​S∋e∑​xS​≥1,∀e∈ExS​∈{0,1}​

其LP relaxation的对偶问题如下:
$$\begin{array}{rl}\max & \sum _{e\in E}{y}_e\\ \text{s.t. }& \sum _{e\in S}{y}_e\le c_S,\forall S\in \mathcal{S}\\ & y_e\ge 0,\forall e\in E.\end{array}$$

互补松弛条件

• Primal complementary slackness
  $$x_S>0\Rightarrow \sum _{e\in S}{y}_e=c_S$$
• Dual complementary slackness
  $$y_e>0\Rightarrow \sum _{S\ni e}{x}_S=1$$

思路: Relax对偶互补条件: 存在$\beta >1$, 满足
$$y_e>0\Rightarrow \sum _{S\ni e}{x}_S\le \beta .$$

Primal-Dual算法

1. 初始化$x_S=0,\forall S\in \mathcal{S}$, $y_e=0,\forall e\in E$.
2. 如果存在$e\in E$未被覆盖, 则增加$y_e$直到存在$S\in \mathcal{S}$使得${\sum }_{e\in S}{y}_e=c_S$, 然后把所有满足条件的$S$添加到结果中(即, $x_S=1$), 并把相应的元素标记未已覆盖.

说明 上述$y_e$的取值可以在区间$[0,\max c_S]$通过二分查找得到.

近似比

令 f e = ∣ { S j ∣ e ∈ S j } ∣ f_e = |\{S_j | e \in S_j\}| fe​=∣{Sj​∣e∈Sj​}∣, 代表元素$e$在所有$S_j$中出现的次数. 令$f=\max f_e$, 我们有${\sum }_{S\ni e}{x}_S\le f$. 因此, 该算法的近似比为$f$, 即$\text{ALG}\le \beta \cdot \text{OPT}$.

Python实现

class SetCoverPD(object):
    """ 基于Primal-Dual的近似算法求解Set Cover问题.
    """

    def __init__(self, m, sets, costs):
        """
        注意: 输出参数的数值必须是整数.
        :param m: 元素个数. 元素的编号从0开始, e.g. 0, 1, ..., m-1
        :param sets: 元素(下标)的集合, list of sets, e.g. [{0, 2}, {1, 2, 3, 6}, {3, 4, 5}, {2, 4, 6}]
        :param costs: 每个元素集合的cost, e.g. [2, 3, 4, 1]
        """
        self._m = m
        self._sets = sets
        self._costs = costs
        self._result = []  # 计算结果: 保存结果集合在sets中的编号
        self._y = [0] * self._m  # 对偶决策变量

    def solve(self):
        """
        Primal-Dual算法.
        """
        uncovered_elements = set(range(self._m))
        ub = max(self._costs)
        while uncovered_elements:
            # 如果存在未被覆盖的元素i, 通过二分查找y[i]的值使得它满足条件:
            # y(S)=c(S), for some S 包含元素i. (称S为tight set.)
            i = uncovered_elements.pop()
            tight_sets = self._binary_search(i, 0, ub)
            # 把关于i的所有tight sets添加到结果集
            for ts in tight_sets:
                uncovered_elements -= self._sets[ts]
                self._result.append(ts)

    def _is_some_set_satisfied(self, i):
        """ 给定i(已知y[i]), 判断是否存在集合S满足y(S)>=c(S),
        其中y(S)代表S中元素对应y值之和, c(S)代表集合S的cost.
        """
        for j in range(len(self._sets)):
            set_j = self._sets[j]
            if i not in set_j:
                continue
            sum_y = sum([self._y[k] for k in set_j])
            if sum_y >= self._costs[j]:
                return True
        return False

    def _binary_search(self, i, lb, ub):
        """
        用二分法查找y[i]使得y(S) = c(S) for some S (S称为tight set)
        :param i: 元素的下标
        :param lb: y_i的下界
        :param ub: y_i的上界
        :return: list of tight sets
        """
        if ub - lb <= 1:
            self._y[i] = ub
            return self._get_tight_sets(i)
        mid = (lb + ub) // 2
        self._y[i] = mid
        if not self._is_some_set_satisfied(i):
            return self._binary_search(i, mid, ub)
        else:
            return self._binary_search(i, lb, mid)

    def _get_tight_sets(self, i):
        """
        已知y[i], 找到所有tight sets(满足y(S)=c(S)).
        """
        tight_sets = []
        for j in range(len(self._sets)):
            set_j = self._sets[j]
            if i not in set_j:
                continue
            sum_y = sum([self._y[k] for k in set_j])
            if abs(sum_y - self._costs[j]) < 1e-6:
                tight_sets.append(j)

        return tight_sets

    def print_result(self):
        print("solution:", [self._sets[i] for i in self._result])
        print("total cost:", sum([self._costs[i] for i in self._result]))


if __name__ == '__main__':
    sc = SetCoverPD(9,  # m
                    [{0, 1}, {2, 5, 8}, {3, 5}, {4, 6}, {1, 2, 3, 4}, {6, 7, 8}],  # sets
                    [4, 2, 5, 7, 8, 1]  # costs
                    )
    sc.solve()
    sc.print_result()