参考 2018 年高睿泉的集训队论文《浅谈保序回归问题》。

保序回归问题

定义

保序回归问题是指，对于一个正整数 p，给定一个偏序关系（有向无环图），图中每个点 i 有权值 (w_i,y_i)，你要给每个点赋上另一个权值 f_i，使得图中任意一条有向路径上的 f_i 递增，并最小化回归代价:

我们把上述问题称作 L_p 问题。

举个简单的例子，P4331 这道题就可以转化为全序关系（链）上的 L_1 问题，可以使用单调栈 + 堆来解决。但这种方法不能拓展到偏序关系上，所以这里不深入讨论。

L_p 均值与其性质

我们定义 L_p 均值为使 \sum\limits_{i\in U}w_i|k-y_i|^p 最小的 k，也就是所有 f_i 相等时的最优取值。对目标式求导，解出其零点就是 L_p 均值，这里给出 L_1 均值为 y_i 的带权中位数，L_2 均值为 y_i 的带权平均数 （\left(\sum\limits_{i\in U}w_iy_i\right)\left/\left(\sum\limits_{i\in U}w_i\right)\right.）。可以证明，当 p>1 时，L_p 均值是唯一的；且 L_p 问题一定存在一组最优解，使得对于任意 f_i 都可以找到某个点集，使它的 L_p 均值等于 f_i。

整体二分

应用整体二分的思想，可以得到一个求问题的近似解的算法。

我们在原问题的基础上加入一些限制。定义 solve(U,l,r) 表示当前点集为 U，且 f_i 的取值在区间 [l,r] 内的 L_p 问题。设 mid=\lfloor(l+r)/2\rfloor，如果我们能够知道 U 的任意一组最优解中，有哪些点是 \leq mid 的，就可以将当前问题划分成两个规模减半的子问题了。

尝试构造另一个问题。设开区间 (a,b) 被 [l,r] 包含，且 U 的任意非空子集的 L_p 均值都不在 (a,b) 内。这个性质保证了对于任意的子问题，都存在一组最优解满足其所有元素均不在 (a,b) 中。下面给出一个引理：

若 solve(U,a,b) 的一组最优解为 \tilde z，且满足 \tilde z_{i}\notin (a,b),\forall i\in U。则存在 solve(U,l,r) 的一组最优解 z，使得 z_i\leq a 当且仅当 \tilde z_i=a。

证明：利用反证法。由于上面的性质，显然 solve(U,a,b) 的解只有 a 与 b 两种取值。下面把 \tilde z_i=a 的点称为 A 类点，\tilde z_i=b 的点称为 B 类点，假设存在 B 类点 i 满足 z_i\leq a。设点集 U 的 L_p 均值为 k，分两种情况讨论：

• 若 k\geq b，则可以发现对于一个 B 类点，它在区间 (-\infty,k] 内的回归代价递减，而在区间 [k,+\infty) 内的回归代价递增。设 t 是满足 z_i\leq a 的 B 类点中最大的 z_i，那么我们把所有 z_i=t 的 B 类点移动到 b 位置处，它们的回归代价变小，而所有点的相对大小关系不变。因此与 z 是 solve(U,l,r) 的一组最优解的假设矛盾。

• 若 k\leq a，则 z_i\leq a 的 B 类点 i 显然不会存在 A 类点 j 满足 z_i\leq z_j 这样的偏序关系。如果这样的 B 类点只有 z_i=t 这一种取值，则其所有 y_i\leq a，所以将它们改为 A 类点回归代价会减小，且所有点的相对大小仍不发生改变。因此与 \tilde z 是 solve(U,a,b) 的一组最优解的假设矛盾。对于这些 B 类点的 z_i 有多种取值的情况，则考虑设点集 V=\{i\mid z_i<t,\tilde z_i=b\}，其 L_p 均值为 k'，如果 k'\leq t 则可以将 V 中 z_i 最大的点调整为 A 类点得到 solve(U,a,b) 的一组更优解；否则类似第一种情况，将 V 中的一些 z_i 调整到 k' 处得到一组不劣于 solve(U,l,r) 的答案，如果没有比原答案更优则显然 V 中所有 z_i 的取值相同，重复上述过程即可导出矛盾。

因为 a<b，所以如果 U 中任意一点 i 选择了 b，则其后继节点也都只能选 b。那么问题就转化为了一个最小权闭合子图问题，节点 i 的权值为 w_i[(b-y_i)^p-(a-y_i)^p]。将权值取反可以转化为最大权闭合子图问题，这是一个经典的网络流模型，且可以根据残量网络构造方案。

回到原问题，考虑找到满足性质的区间 (mid,mid+\epsilon)，那么就可以由 solve(U,mid,mid+\epsilon) 的最优解来确定 U 中的哪些点在原问题的最优解中是 \leq mid 的。但直接找到这样的区间比较麻烦，可能还需要跑浮点网络流，会遇到精度等一系列问题。

由于 U 的任意非空子集的 L_p 的并集是有限集，所以一定存在足够小的 \epsilon 使得区间中没有这些元素。当 \epsilon 无限趋近于 0 时，节点 i 的权值为 w_i[(mid+\epsilon-y_i)^p-(mid-y_i)^p]。把所有点的权值除以 \epsilon，最优方案不变，权值变为 w_i(mid-y_i)^p 在 mid 处的导数。当 p=2 时，权值为 2w_i(mid-y_i)，于是问题得到了完美的解决。

本题题解

题目中要求了大小关系的子集都是给定的 n 个 01 向量的极大线性无关组。所以从某一组向量中拿走一个，线性空间维数会减一，也就能再加入一个向量。因此任意两组向量之间是可以由替换元素这个操作来互相转化的，并且替换的顺序没有限制。于是我们只需要从向量组 A 的每个元素向能替换它的元素连边，向量组 B 反之，即可构造出偏序关系。然后就可以跑 L_2 问题了。

然后有个问题就是这题要求修改后的价格为整数，可以这样处理：当 r-l=1 时，仍然是答案只有两种取值的情况，因此仍可以用转化成最小闭合子图问题，相应地每个点的权值设为选 r 与选 l 的代价差即可。

代码

#include <cstdio>
#include <queue>
using namespace std;
typedef unsigned long long ull;
const int MAXN=1005;
const int INF=1E9;
namespace MF
{
    struct Edge
    {
        int from, to, cap, flow;
        Edge(int u, int v, int c, int f):
            from(u), to(v), cap(c), flow(f) {}
    };
    int n, s, t, flow;
    vector<Edge> edges;
    vector<int> g[MAXN];
    int h[MAXN], cur[MAXN];
    void init(int v, int a, int b)
    {
        n=v, s=a, t=b, flow=0;
        edges.clear();
        for (int i=1; i<=n; i++) g[i].clear();
    }
    void addEdge(int from, int to, int cap)
    {
        edges.push_back(Edge(from, to, cap, 0));
        edges.push_back(Edge(to, from, 0, 0));
        g[from].push_back(edges.size()-2);
        g[to].push_back(edges.size()-1);
    }
    bool bfs()
    {
        queue<int> q;
        for (int i=1; i<=n; i++) cur[i]=h[i]=0;
        h[s]=1, q.push(s);
        while (!q.empty())
        {
            int u=q.front(); q.pop();
            for (int i=0; i<g[u].size(); i++)
            {
                Edge e=edges[g[u][i]];
                if (e.cap>e.flow&&!h[e.to])
                    h[e.to]=h[u]+1, q.push(e.to);
            }
        }
        return h[t];
    }
    int dfs(int u, int f)
    {
        if (u==t) return f;
        for (int &i=cur[u]; i<g[u].size(); i++)
        {
            Edge &e=edges[g[u][i]];
            if (e.cap>e.flow&&h[e.to]==h[u]+1)
            {
                int d=dfs(e.to, min(f, e.cap-e.flow));
                if (d>0)
                {
                    e.flow+=d;
                    edges[g[u][i]^1].flow-=d;
                    return d;
                }
            }
        }
        return 0;
    }
    void Dinic()
    {
        while (bfs())
            while (int f=dfs(s, INF)) flow+=f;
    }
}
namespace XB
{
    ull d[64];
    void clear()
    {
        for (int i=0; i<64; i++) d[i]=0;
    }
    int insert(ull x)
    {
        for (int i=63; i>=0; i--)
            if (x&1ll<<i)
                if (d[i]) x^=d[i];
                else return d[i]=x, i;
        return -1;
    }
}
ull c[MAXN];
int w[MAXN], a[66], b[66];
int p[MAXN], q[MAXN], t1[MAXN], t2[MAXN];
int f[MAXN];
vector<int> G[MAXN];
void solve(int l, int r, int L, int R)
{
    if (l>r) return;
    if (R-L==1)
    {
        MF::init(r-l+3, r-l+2, r-l+3);
        for (int i=l; i<=r; i++) q[p[i]]=i-l+1;
        for (int i=l; i<=r; i++)
        {
            int u=p[i], cost=2*(w[u]-L)-1;
            if (cost>0) MF::addEdge(MF::s, q[u], cost);
            else MF::addEdge(q[u], MF::t, -cost);
            for (int v: G[u]) if (q[v])
                MF::addEdge(q[u], q[v], INF);
        }
        for (int i=l; i<=r; i++) q[p[i]]=0;
        MF::Dinic();
        for (int i=l; i<=r; i++)
            if (!MF::h[i-l+1]) f[p[i]]=L;
            else f[p[i]]=R;
        return;
    }
    int mid=(L+R)/2;
    MF::init(r-l+3, r-l+2, r-l+3);
    for (int i=l; i<=r; i++) q[p[i]]=i-l+1;
    for (int i=l; i<=r; i++)
    {
        int u=p[i];
        if (mid<w[u]) MF::addEdge(MF::s, q[u], w[u]-mid);
        else MF::addEdge(q[u], MF::t, mid-w[u]);
        for (int v: G[u]) if (q[v])
            MF::addEdge(q[u], q[v], INF);
    }
    for (int i=l; i<=r; i++) q[p[i]]=0;
    MF::Dinic();
    int c1=0, c2=0;
    for (int i=l; i<=r; i++)
        if (!MF::h[i-l+1]) t1[++c1]=p[i];
        else t2[++c2]=p[i];
    for (int i=1; i<=c1; i++) p[l+i-1]=t1[i];
    for (int i=1; i<=c2; i++) p[l+c1+i-1]=t2[i];
    solve(l, l+c1-1, L, mid);
    solve(l+c1, r, mid, R);
}
int main()
{
//  freopen("shop.in", "r", stdin);
//  freopen("shop.out", "w", stdout);
    int n, m;
    scanf("%d%d", &n, &m);
    for (int i=1; i<=n; i++) scanf("%llu", &c[i]);
    for (int i=1; i<=n; i++) scanf("%d", &w[i]);
    for (int i=1; i<=m; i++) scanf("%d", &a[i]);
    for (int i=1; i<=m; i++) scanf("%d", &b[i]);
    for (int i=1; i<=m; i++)
    {
        XB::clear();
        for (int j=1; j<=m; j++)
            if (j!=i) XB::insert(c[a[j]]);
        for (int j=1; j<=n; j++)
            if (j!=a[i])
            {
                int t=XB::insert(c[j]);
                if (~t) G[a[i]].push_back(j), XB::d[t]=0;
            }
    }
    for (int i=1; i<=m; i++)
    {
        XB::clear();
        for (int j=1; j<=m; j++)
            if (j!=i) XB::insert(c[b[j]]);
        for (int j=1; j<=n; j++)
            if (j!=b[i])
            {
                int t=XB::insert(c[j]);
                if (~t) G[j].push_back(b[i]), XB::d[t]=0;
            }
    }
    for (int i=1; i<=n; i++) p[i]=i;
    solve(1, n, 0, 1E6);
    ull ans=0;
    for (int i=1; i<=n; i++) ans+=1ll*(f[i]-w[i])*(f[i]-w[i]);
    printf("%llu\n", ans);
    return 0;
}