浅谈保序回归问题
说在前头
本篇博客仅仅介绍保序回归问题最简单问题的解法，对于保序回归问题中的许多证明以及拓展延伸都略过不谈，详情请参见2018国家集训队论文《浅谈保序回归问题》(高睿泉)。
保序回归问题
定义
给定正整数 \(k\)，以及一张点集为 \(V\)，边集为 \(E\) 的 \(DAG\) G，点 \(i\) 有权值 \(a_i,w_i\)。
请为每个点赋一个点值 \(f_i\)，使得对于任意边 \((u,v)\in E\) 有 \(f_u\le f_v\)，且 \(\sum_{i}w_i|f_i-a_i|^k\) 最小。
对于相同的 \(k\)，我们称之为 \(L_k\) 问题。
一般解法
对于树上情况有特殊的做法，但与一般解法关系不大，在此不作讨论。
考虑整体二分，同时二分出所有点的 \(f_i\)，设当前正在进行 \(solve(V,l,r)\)，表示正在求点集 \(V\) 的权值，且已经确定它们 \(f_i\) 的取值为 \([l,r]\)，令 \(mid=\dfrac{l+r}{2}\)，找一个极小值 \(\theta\)（\(\theta\) 的大小根据题目要求而定）。考虑求出哪些点的 \(f\) \(\in [l,mid]\)，哪些点的 \(f\in [mid+\theta,r]\)，然后对两边递归处理。
构造另一个问题，在这个问题中，我们只考虑 \(V\) 中的点，且 \(V\) 中的点的 \(f\) 取值只能为 \(\{mid,mid+\theta\}\) ，在此基础上满足原题限制时求出最小回归代价。
此时不加证明给出一个引理（证明请参考开头的论文）：
若 \(solve(V,l,r)\) 的最优解为\(\{b_i\}\)，那么将 \(b_i\) 中所有 \(\le mid\) 的都替换为 \(mid\)，\(>mid\) 的都替换为 \(mid+\theta\) ，得到的新解 \(\{c_i\}\) 是新问题的最优解。由此我们只需要解决新问题，就能将 \(f\le mid\) 与 \(>mid\) 的部分分离开来了。
对于新问题的求解，对于原图中的边 \(u\rightarrow v\)，如果 \(u\) 取 \(mid+\theta\)，那么 \(v\) 必须取 \(mid+\theta\)。在此基础上最小化回归代价，这就是经典的最小闭合子图问题，可以使用网络流算法来完成。
特殊解法
- 对于树上的情况，可以 \(DP\) 维护分段函数解决，具体参见论文。
- 当 \(p=1\) 时，最终答案的取值一定为整数，因此 \(\theta\) 可以取 \(1\)。
- 当图为树/基环树时有 \(\mathcal O(n\log n)\) 做法，具体请参见下方 \(loj 2470\) 的题解。
例题
洛谷P6621 [省选联考 2020 A 卷] 魔法商店
Description
有 \(n\) 件物品，有权值 \(c_i\)，价格 \(v_i\)。定义一个物品集合合法当且仅当且非空、其中所含的物品权值线性无关且集合大小尽量的大。给定两个符合条件的集合 \(A,B\)，你可以花费 \(\sum_i (x_i-v_i)^2\) 的代价，将物品价格改为 \(\{x_i\}\) 且 \(x_i\) 必须为整数。请你花费最小代价使得 \(A\) 为所有合法集合中价格总和最小的，\(B\) 为价格总和最大的。
\(n\le 1000,|A|=|B|\le 64,c_i<2^{64},v_i\le 10^6\)。
Solution
考虑找到所有的合法集合，这是非常困难的，但是注意到题目中已经给出了一个合法集合 \(A\)，考虑在 \(A\) 的基础上删掉一些物品再加入一些物品得到合法集合。注意到，这样所有在 \(A\) 的基础上只修改一个物品得到的合法集合权值都比 \(A\) 大，那么 \(A\) 就是价格最小的。
因此，我们只考虑改变哪些物品能使得集合依然合法，枚举改变的元素 \(i\)，再枚举新元素 \(j\) 插入 \(A\) 其他元素组成的线性基中，若 \(j\) 能被表示出来，那么 \(i\) 的权值必须 \(\le j\)。对 \(B\) 同样处理。最终我们就得到了一个偏序关系的 \(DAG\)，直接使用上面的保序回归问题算法即可。由于需要保证修改价格为整数，因此令 \(\theta=1\)。
Code
#include<bits/stdc++.h>
using namespace std;
const int N=1010;
typedef unsigned long long ull;
typedef long long ll;
const ll inf=1e16;
int n,m,pd[N],v[N],a[N],b[N],mx,f[N],w[N];
ull c[N],d[N];
vector<int> to[N];
inline bool insert(ull x,bool tp){
for(int i=63;i>=0;--i){
if(x&(1ull<<i)){
if(d[i]) x^=d[i];
else{tp?d[i]=x:0;return false;}
}
}
return true;
}
inline void init(){
for(int i=1;i<=m;++i){
memset(d,0,sizeof(ull)*(64));
for(int i=1;i<=n;++i) pd[i]=0;
for(int j=1;j<=m;++j){
pd[a[j]]=1;
if(i!=j) insert(c[a[j]],1);
}
for(int j=1;j<=n;++j)
if(!pd[j])
if(!insert(c[j],0)) to[a[i]].push_back(j);
}
for(int i=1;i<=m;++i){
memset(d,0,sizeof(ull)*(64));
for(int i=1;i<=n;++i) pd[i]=0;
for(int j=1;j<=m;++j){
pd[b[j]]=1;
if(i!=j) insert(c[b[j]],1);
}
for(int j=1;j<=n;++j)
if(!pd[j])
if(!insert(c[j],0)) to[j].push_back(b[i]);
}
}
namespace flow{
struct node{
int v,nxt,f;
}e[N*N];
int first[N],cnt=1,tot,s,t,cur[N],q[N],l,r,dis[N];
inline void add(int u,int v,int f){e[++cnt].v=v;e[cnt].f=f;e[cnt].nxt=first[u];first[u]=cnt;}
inline void Add(int u,int v,int f){add(u,v,f);add(v,u,0);}
inline void init(){
memset(first+1,0,sizeof(int)*(tot));
tot=0;cnt=1;
}
inline bool bfs(){
fill(dis+1,dis+tot+1,-1);
l=1;r=0;q[++r]=s;dis[s]=0;
while(l<=r){
int u=q[l++];
for(int i=first[u];i;i=e[i].nxt){
int v=e[i].v;
if(e[i].f&&dis[v]==-1){
dis[v]=dis[u]+1;
q[++r]=v;
}
}
}
return dis[t]!=-1;
}
inline int dfs(int u,int f){
if(!f||u==t) return f;
int used=0;
for(int& i=cur[u];i;i=e[i].nxt){
int v=e[i].v;
if(e[i].f&&dis[v]==dis[u]+1){
ll fl=dfs(v,min(f,e[i].f));
used+=fl;f-=fl;
e[i].f-=fl;e[i^1].f+=fl;
if(!f) break;
}
}
if(f) dis[u]=-1;
return used;
}
inline ll dinic(int S,int T){
ll f=0;s=S;t=T;
while(bfs()){
memcpy(cur,first,sizeof(int)*(tot+1));
f+=dfs(s,inf);
}
return f;
}
}
using flow::Add;
using flow::tot;
int q[N],st[N];
inline ll calc(int i,int x){return 1ll*(x-v[i])*(x-v[i]);}
inline void solve(int ql,int qr,int l,int r){
if(l==r){
for(int i=ql;i<=qr;++i) f[q[i]]=l;
return ;
}
if(ql>qr) return ;
int mid=(l+r)>>1;
flow::init();tot=qr-ql+1;
int s=++tot,t=++tot;
for(int i=ql;i<=qr;++i) pd[q[i]]=i-ql+1;
for(int i=ql;i<=qr;++i){
int u=q[i];
for(int i=0;i<to[u].size();++i){
int v=to[u][i];
if(pd[v]) Add(pd[u],pd[v],inf);
}
ll w=calc(u,mid)-calc(u,mid+1);
if(w>0) Add(s,pd[u],w);
else Add(pd[u],t,-w);
}
flow::dinic(s,t);
int L=ql,R=qr;
for(int i=ql;i<=qr;++i){
int u=q[i];
if(flow::dis[pd[u]]!=-1) st[R--]=u;
else st[L++]=u;
pd[u]=0;
}
for(int i=ql;i<=qr;++i) q[i]=st[i];
solve(ql,L-1,l,mid);solve(L,qr,mid+1,r);
}
int main(){
scanf("%d%d",&n,&m);
for(int i=1;i<=n;++i) scanf("%llu",&c[i]),w[i]=1;
for(int i=1;i<=n;++i) scanf("%d",&v[i]),mx=max(mx,v[i]);
for(int i=1;i<=m;++i) scanf("%d",&a[i]);
for(int i=1;i<=m;++i) scanf("%d",&b[i]);
init();
for(int i=1;i<=n;++i) q[i]=i;
solve(1,n,0,mx);
ll ans=0;
for(int i=1;i<=n;++i) ans+=calc(i,f[i]);
printf("%lld\n",ans);
return 0;
}
loj#2470. 「2018 集训队互测 Day 2」有向图
Description
树或基环树上的 \(L_1\) 问题。
\(n\le 30000,m\in \{n-1,n\}\)。
Solution
首先直接保序回归显然可行，但复杂度太高，考虑优化网络流求解最小闭合子图的过程。
考虑 \(DP\)，首先对于树的情况，将图看成无向图，设 \(f_{u,0/1}\) 表示只考虑 \(u\) 的子树，\(u\) 权值为 \(mid+0/1\) 时的最小回归代价，那么有转移：
转移边为 \(u\ge v\)；
为 \(u\le v\)：
直接 \(DP\) 即可做到 \(\mathcal O(n)\)。
对于基环树的情况，先随机钦定一个环上的点作为根进行 \(DP\)，但这样会漏考虑环上的最后一条边，不妨设根为 \(rt\)，最后一条边为 \((rt,u)\)。因此再开一维维护根的状态，当遍历到 \(u\)，若这条边为 \(rt<v\)，那么将 \(f_{v,1,0}\) 设为 \(inf\)，否则将 \(f_{v,0,1}\) 设为 \(inf\)。然后照常 \(DP\) 即可。最终单次操作复杂度为 \(\mathcal O(n)\),总复杂度为 \(\mathcal O(n\log n)\)。
Code
#include<bits/stdc++.h>
using namespace std;
const int N=3e5+10;
typedef long long ll;
namespace iobuff{
const int LEN=1000000;
char in[LEN+5],out[LEN+5];
char *pin=in,*pout=out,*ed=in,*eout=out+LEN;
inline char gc(void){
return pin==ed&&(ed=(pin=in)+fread(in,1,LEN,stdin),ed==in)?EOF:*pin++;
}
inline void pc(char c){
pout==eout&&(fwrite(out,1,LEN,stdout),pout=out);
(*pout++)=c;
}
inline void flush(){fwrite(out,1,pout-out,stdout),pout=out;}
template<typename T> inline void read(T &x){
static int f;
static char c;
c=gc(),f=1,x=0;
while(c<'0'||c>'9') f=(c=='-'?-1:1),c=gc();
while(c>='0'&&c<='9') x=10*x+c-'0',c=gc();
x*=f;
}
template<typename T> inline void putint(T x,char div){
static char s[15];
static int top;
top=0;
x<0?pc('-'),x=-x:0;
while(x) s[top++]=x%10,x/=10;
!top?pc('0'),0:0;
while(top--) pc(s[top]+'0');
pc(div);
}
}
using namespace iobuff;
const ll inf=4e18;
int n,m,k,pd[N],v[N],mx,f[N],w[N],first[N],cnt=1;
struct node{
int v,nxt;
}e[N<<1];
inline void add(int u,int v){e[++cnt].v=v;e[cnt].nxt=first[u];first[u]=cnt;}
char s[10];bool col[N],vis[N];
int rt,q[N],st[N],fr[N];
ll g[N][2][2];
inline ll calc(int i,int x){return 1ll*w[i]*(k==1?abs(x-v[i]):1ll*(x-v[i])*(x-v[i]));}
inline void dfs(int u,int f){
vis[u]=1;bool fl=0,tmp=0;fr[u]=f;
for(int i=first[u];i;i=e[i].nxt){
int v=e[i].v;
if((i^(f^1))&&pd[v]){
if(!vis[v]){
dfs(v,i);
int x=~i&1;
for(int t=0;t<2;++t){
g[u][t][x]=min(g[u][t][x]+g[v][t][x],inf);
g[u][t][x^1]=min(g[u][t][x^1]+min(g[v][t][0],g[v][t][1]),inf);
}
}
else if(!rt) fl=1,rt=v,tmp=(~i&1);
}
}
if(fl) g[u][tmp^1][tmp]=inf;
if(u==rt) g[u][0][1]=g[u][1][0]=inf;
}
inline void cover(int u,bool x,bool y){
col[u]=y;
for(int i=first[u];i;i=e[i].nxt){
int v=e[i].v;
if((i^fr[v])||!pd[v]) continue;
int t=~i&1;
if(y==t) cover(v,x,y);
else cover(v,x,g[v][x][0]<g[v][x][1]?0:1);
}
}
inline void solve(int ql,int qr,int l,int r){
if(l==r){
for(int i=ql;i<=qr;++i) f[q[i]]=l;
return ;
}
if(ql>qr) return ;
int mid=(l+r)>>1;
for(int i=ql;i<=qr;++i){
pd[q[i]]=1;
g[q[i]][0][0]=g[q[i]][1][0]=calc(q[i],mid);
g[q[i]][0][1]=g[q[i]][1][1]=calc(q[i],mid+1);
}
rt=0;
for(int i=ql;i<=qr;++i){
int u=q[i];
if(!vis[u]){
dfs(u,0);
int x=0,y=0;
for(int j=0;j<2;++j)for(int k=0;k<2;++k)if(g[u][j][k]<g[u][x][y]) x=j,y=k;
cover(u,x,y);
}
}
int L=ql,R=qr;
for(int i=ql;i<=qr;++i){
int u=q[i];
if(!col[u]) st[L++]=u;
else st[R--]=u;
pd[u]=0;vis[u]=0;
}
for(int i=ql;i<=qr;++i) q[i]=st[i];
solve(ql,L-1,l,mid);solve(L,qr,mid+1,r);
}
int main(){
read(n);read(m);k=1;
for(int i=1;i<=n;++i) read(v[i]),mx=max(mx,v[i]);
for(int i=1;i<=n;++i) read(w[i]);
for(int i=1,u,v;i<=m;++i){
read(u);read(v);
add(u,v);add(v,u);
}
for(int i=1;i<=n;++i) q[i]=i;
solve(1,n,0,mx);
ll ans=0;
for(int i=1;i<=n;++i) ans+=calc(i,f[i]);
printf("%lld\n",ans);
return 0;
}
