题解 P6621 【[2020 年联考 A 卷] 魔法商店】
灵梦
2020-06-22 20:36:30
参考 2018 年高睿泉的集训队论文《浅谈保序回归问题》。
## 保序回归问题
### 定义
保序回归问题是指，对于一个正整数 $p$，给定一个偏序关系（有向无环图），图中每个点 $i$ 有权值 $(w_i,y_i)$，你要给每个点赋上另一个权值 $f_i$，使得图中任意一条有向路径上的 $f_i$ 递增，并最小化回归代价:
$$
\sum_{i\in U} w_i|f_i-y_i|^p
$$
我们把上述问题称作 **$L_p$ 问题**。
举个简单的例子，[P4331](https://www.luogu.com.cn/problem/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$ 的代价差即可。
## 代码 ##
```cpp
#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;
}
```