线性规划 2

强对偶性的证明

接续上次的讨论,我们已经得到了弱对偶性,这导出了最大匹配不大于最小点覆盖。然而需要了解的是,在二分图中二者相等,这就对应着强对偶性。

我们首先需要考察的是原问题和对偶问题是否存在可行解。我们会逐渐从松弛形过渡到标准形。

引理 1. ARn×m,bRn,则下面两个叙述有且只有一个成立:

  1. Ax=b 有解。
  2. ATy=0, bTy=1 有解。

我们可以从几何意义上解释该引理:若一个向量 b 不存在于 A 的列张成 V 中,则其一定于一个正交于 V 的向量成钝角。

然后证明:
如果 1.2. 同时成立,则 1=bTy=xTATy=0,这矛盾。
如果 1. 不成立,则向量 b 不存在于 A 的列张成 V 中,即 rank([A  b])=rank(A)+1,从而矩阵

[ A b0T 1 ]

的秩也是 rank(A)+1。因此 [0  1]T 可以被 [A  b]T 线性表出,因此存在 y

证完。

引理 2. (Farkas) ARn×m,bRn,则下面两个叙述有且只有一个成立:

  1. Ax=b,x0 有解。
  2. ATy0, bTy0 有解。

Farkas 引理是最优化理论的基石之一。

几何意义和引理 1 相似。Ax s.t. x0 表出了 A 张成的一个凸锥。如果一个点不在凸锥内部,则这个点和凸锥一定能由一个超平面分开,且这个超平面过原点。这个超平面对应的方程是 xTy=0

然后证明:
如果 1.2. 同时成立,则 0xTATy=(Ax)Ty=bTy<0,这矛盾。
如果 1. 不成立,则不妨设 Ax=b 有解,但解不满足 x0。无解的情况由引理 1 可知。

A=(a1,a2,,an)。现在对 n 作归纳法证明。
n=1 时,条件为 x1a1=b 有解,且 x1<0。取 y=b,则有 a1Ty=bTb/x1>0, bTy=bTb<0
归纳。假设我们已经对 k<n 证明了原命题。由于 1. 不成立,所以 i=1n1xiai=b 一定没有满足 i<n, xi0 的解(否则 xn=0 即可)。由归纳,v, aiTv0(1i<n)bTv<0。假设 anTv0 则原命题成立,因此不妨假设 anTv<0
我们取

ai^=(aiTv)an(anTv)ai(i<n)b^=(bTv)an(anTv)b

则易见 i=1n1xiai^=b^ 也没有满足 i, xi0 的解。否则整理得到

i=1n1aixi+1anTv(bTvi=0n1aiTvxi)an=b

这与 1. 不成立矛盾。归纳得 w, ai^Tw0(1i<n)b^Tw<0
可以验证 y=(anTw)v(anTv)w 即为解。

证完。

推论 1 ARm×n,bRm,则下面两个叙述有且只有一个成立:

  1. Axb, x0 有解。
  2. ATy0, y0, bTy<0 有解。

证明:
如果 1.2. 同时成立,则 0xTATy=(Ax)TybTy<0,这矛盾。

如果 1. 不成立,可以通过转松弛形的方式转换等价形式。原问题等价于

[AIm][xx]=b,x,x0

无解。由 Farkas 引理,存在 y 使得

[ATIm]y0,bTy<0

证完。

我们需要的拼图都得到了,下一步就是证明强对偶性。

定理 2(强对偶性)

设线性规划的原问题和对偶如定义 1,3 形式,记为 P,D。下面的叙述有且仅有一个成立:

  1. P,D 均有最优解,且 P=D
  2. P 有可行解,D 无可行解,且 P 的目标函数在约束条件下无界。
  3. D 有可行解,P 无可行解,且 D 的目标函数在约束条件下无界。
  4. P,D 均无可行解。

假设 P,D 均有可行解。根据弱对偶性,我们只需要证明

Axb,x0ATyc,y0cTxbTy0

有解;也就是系统

[A 00 ATcT bT][xy][bc0][xy]0

有解。

如果其无解,则根据推论 1,存在 z,w0,α0 满足

AwαbATzαcbTz<cTw

α0。不然,设 x^,y^ 分别是 P,D 的可行解,有

0x^TATz=(Ax^)TzbTz<cTw(ATy^)Tw=y^TAw0

这矛盾。

此时可以设 x=w/α,y=z/α。可以验证 x,y 分别是 P,D 的可行解。由弱对偶性引出

cTw=α(cTx)=α(bTy)=bTz

这与 bTz<cTw 矛盾。

因此上系统必有解,即最优解的值相同。

假设 P 无可行解,D 有可行解。由推论 1 可得 ATw0,w0,bTw<0 有解。设 yD 的可行解,能注意到对任意非负实数 λy+λw 都是可行解。D 的目标函数可以写作 bT(y+λw)=bTy+λ(bTw)。由于 bTw<0,因此取充分大的 λ 可以使目标函数取值充分小,D 无界。2. 同理。

剩余情况即为 4.

证完。

强对偶性还能导出一个很重要的性质,即互补松弛性。

推论 2 (互补松弛性) 

设线性规划的原问题和对偶如定义 1,3 形式,记为 P,D。设 x,y 分别为 P,D 的最优解,下面的叙述有且仅有一个成立:

  1. yi>0 推出 j=1nai,jxj=bi(1im)
  2. xj>0 推出 i=1mbi,jyi=cj(1jn)

具体地说,就是第 i 个分量非负的约束和对偶问题的第 i 个约束一定有一个取等。

由强对偶性,有

j=1ncjxj=i,jyiai,jxj=i=1mbiyi

这得到

i=1m(bij=1nai,jxj)yi=0

由于求和的每一项都 0,因此每一项都必须 =0。因此第 i 个分量 0 的约束和对偶问题的第 i 个约束一定有一个取等。

因此必要性得证。

1. 2. 都成立时有

cTx=j=1ncjxj=j=1n(i=1mai,jyi)xj=i=1m(j=1nai,jxj)yi=i=1mbiyi=bTy

由强对偶性得均为最优解。

因此充分性得证。

证完。


例题

Rikka With Flow

给定一个网络,每条边 (u,v) 有代价 cu,v 和费用 wu,v。对于一条边,可以将它的费用减少,不能减到负数,代价为减少量。问要让网络的最大费用循环流费用不超过 x 的最少代价。

HihoCoder 是不是早炸了?

考虑设 wi,j 为更改后的费用。然后根据例 3 的对偶形式,直接写出线性规划问题

minw,d,y(u,v)wu,vw(u,v)s.t.(u,v), du,v+yu,vwu,v(u,v), wu,vwu,v(u,v)du,v×cu,vx(u,v), du,v,wu,v0

这里需要注意的是 d 需要让 (u,v)du,v×cu,v 取到最小值。采用单纯形法解即可。

你问咋做矩阵?考虑在关联矩阵上刻画这些性质即可。由于单纯形法能求解的是最大值,因此求解 wi,j 的最大值后用 wi,j 减就行了。

在关联矩阵上做限制条件是挺常见的 trick。


CF1307G

给出一个 n 个点 m 条边的有向图,每条边有边权 wi 。有 Q 次询问,每次询问给出一个 k 。你可以把一条边修改成 wi+aiai 不一定是整数),不过需要保证 ai0aik

你要通过修改边权使得从 1n 的最短路径尽可能长,每次询问之间独立。数据保证至少存在一条从 1n 的路径,无重边自环。输出答案和标准答案的相对误差或绝对误差应不超过 106

2n50, 1mn(n1), 1q105

用关联矩阵将网络流语言转换成线性代数语言后就能显见地取对偶了,这点的直观性是配凑等方法所不具有的。

考虑令 du 为差分情况下的最短路,即满足 dud11u 的最短路长度。设 wu,vuv 的边的边权,xu,vuv 增加的边权。然后可以写出线性规划:

maxd,xdnd1s.t.dvdu+wu,v+xu,vu,vxu,vkxu,v0

自然地,我们考虑构造限制矩阵。我们令 M 为对应图的 m×n 关联矩阵,其元素取值满足 (u,v) 有一条边 iMi,u1,Mi,v+1。令 d 为差分最短路对应的 n 维向量,x 为每条边增加边权对应的 m 维向量,w 为原边权对应的 m 维向量。和上方的限制条件进行对照后不难可以构造出如下的限制:

maxd,x(1,0,,0,1,0,,0)[dx]s.t.[M Im0 1][dx][wk], [dx]0

其中 Imm×m 单位矩阵每个元素取相反数后得到的矩阵,l=(1,0,,0,1,0,,0) 是一个 m+n 维向量,其中 l1=1, ln=1,其余分量均为 0。得到标准形式后不难转对偶:

miny[wk]ys.t.[MT 0Im 1]yl,y0

注意这里 y 是个 m+1 维的向量。根据互补松弛性取等。自然展开,设 y=[α  β],其中 αu,v 表示 uv 的边对应的变量。得到如下的线性规划:

miny(u,v)wu,vαu,v+βks.t.1<u<n, vαu,vvαv,u=0vvαv,1α1,v=1vvαv,nαn,v=1βαu,v

若将 αu,v 看作流量,则流量上界均是 β。不妨令 αu,v[0,1],令单位流为 f=1/β,改写为

miny(u,v)wu,vαu,v+βkfs.t.1<u<n, vαu,vvαv,u=0vvαv,1α1,v=fvvαv,nαn,v=fαu,v[0,1]

由于答案对 f 单调,因此当 f[s,s+1](sZ) 时目标函数的最值肯定出现在边界,这时 αu,v 的值也是整数。于是我们只需要直接建图跑网络流,对答案枚举取值后取 max 即可。

总时间复杂度 O(Dinic(n,m)+qn)

Submission.


CF375E

给出一棵 n 个节点的树,树上的节点有红黑两种颜色。每次操作可以交换两个节点的颜色,问最少需要多少次操作可以使得树上任意一个点均存在与它距离 x 的黑点,在这里认为树上两个节点的距离为它们之间的最短路径长度。

2n500, 1x109

n 维向量 b 满足 bi=[ i ]n×n 大小的矩阵 A 满足 Ai,j=[ i,j x]。令解 x 满足 xi=[ i ],同时记黑色节点总数为 H。容易写出以下的整数线性规划:

minxi=1nbixis.t.i,j=1n Ai,j×xj1i=1nxi=Hxi{0,1}

我们将说明对 xii, xi0 的约束条件和 i, xi{0,1} 的约束条件得到的结果是相同的。我们称前一种约束条件对应的线性规划为普通线性规划,后一种约束条件对应的线性规划为整数线性规划。

首先不难发现普通线性规划的最优解一定大于等于整数线性规划。因此我们只需要证明普通线性规划的最优解一定小于等于整数线性规划即可。

现在需要证明一定存在一种最优解满足整数线性规划的约束条件。原因有二:

  1. 由于 j=1n Ai,j×xj1,我们不必取 xj>1。显然地,如果存在一个最优解使得 xi>1,我们都可以将初始为红色的点的 x 值增加,这样可以将 xj 降低至 1,新解一定不劣于原解。
  2. 对于一个所有变量都取到整数的最优解,将一些 xi 调整为非整数取值不会得到更优的解。

拆最后的等于为大于等于且小于等于。写成线性代数形式就是

minxbTxs.t.[A1T1T]x[1HH]x0

转对偶得到

maxy[1THH]ys.t.[AT11]yby0

采用单纯形法求解即可。

Submission.


拉格朗日对偶

其实想写写CF1696G这题 但是没啥时间了。先做多项式。

放一下 UOJ 板子的代码:

code
#include <bits/stdc++.h>
using namespace std; 
#define rep(i,s,t) for (register int i = (s), i##_ = (t) + 1; i < i##_; ++ i)
#define pre(i,s,t) for (register int i = (s), i##_ = (t) - 1; i > i##_; -- i)
const int inf = 1e9;
const double eps = 1e-9;
int n, m, t;

namespace Simp {
    // c_i              = - r[0][i]
    // a[i][j]          = - r[i][j]
    // b_i              =   r[i][0]
    // MAX_VALUE        = - r[0][0]
    const int N = 1000;

    double r[N][N];
    int row[N];

    void pivot(int out, int in) {
        swap(row[out + n], row[in]);
        double roi = - r[out][in];
        rep(i,0,n) r[out][i] /= roi;
        r[out][in] = -1 / roi;
        rep(i,0,m) if (i != out) {
            double rin = r[i][in];
            if (abs(rin) <= eps) continue;
            rep(j,0,n) { 
                if (j != in) r[i][j] += rin * r[out][j];
                r[i][in] = r[out][in] * rin;
            }
        }
    }

    bool init() {
        rep(i,1,n+m) row[i] = i;
        while (1) {
            int q = 1; // 找到最小的 b 相同最小所在行
            double b_min = r[1][0];
            rep(i,2,m) if (r[i][0] < b_min) 
                b_min = r[i][0], q = i;
            if (b_min + eps >= 0) return true; // 所有 b \ge 0
            int p = 0;
            rep(i,1,n) if (r[q][i] > eps and (p == 0 or row[i] > row[p])) 
                p = i;
            if (p == 0) break;
            pivot(q, p);
        } return false;
    }

    bool simp() {
        while (1) {
            int t = 1, k = 0; // 最小 t
            rep(j,2,n) if (r[0][j] < r[0][t]) t = j;
            if (r[0][t] >= - eps) return 1; // 最优解
            double rati_min = 1e18;
            rep(i,1,m) if (r[i][t] < - eps) {
                double rati = - r[i][0] / r[i][t];
                if (k == 0 or (rati < rati_min) or 
                   (abs(rati - rati_min) <= eps and row[i] > row[k]))
                    rati_min = rati, k = i;
            } 
            if (k == 0) break; 
            pivot(k, t);
        } 
        return 0;
    }

    void Simplex() {
        if (!init()) puts("Infeasible"), exit(0);
        if (!simp()) puts("Unbounded"), exit(0);
    }

    int col[N];
    void print_xi() {
        rep(i,n+1,n+m) col[row[i]] = i;
        rep(i,1,n) cout << (col[i] ? r[col[i] - n][0] : 0) << ' ';
        cout << endl;
    }
} // namespace Simplex Method
using namespace Simp;

signed main() {
	ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);
    cin >> n >> m >> t; cout << setprecision(12) << fixed;
    rep(i,1,n) cin >> r[0][i], r[0][i] = -r[0][i]; 
    rep(i,1,m) {
        rep(j,1,n) cin >> r[i][j], r[i][j] = -r[i][j];
        cin >> r[i][0];
    }

    Simplex();
    cout << - r[0][0] << '\n';
    if (t) print_xi();
}

我不知道为啥,jjdw 想让我把“自然展开”换成“光翼展开”。虽然我也想换。

posted @   joke3579  阅读(158)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· winform 绘制太阳,地球,月球 运作规律
· 超详细:普通电脑也行Windows部署deepseek R1训练数据并当服务器共享给他人
· 上周热点回顾(3.3-3.9)
· AI 智能体引爆开源社区「GitHub 热点速览」
· 写一个简单的SQL生成工具
点击右上角即可分享
微信分享提示