淺談保序迴歸問題
說在前頭
本篇部落格僅僅介紹保序迴歸問題最簡單問題的解法，對於保序迴歸問題中的許多證明以及拓展延伸都略過不談，詳情請參見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\)
特殊解法
- 對於樹上的情況，可以 \(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\)；
\[f_{u,0}+= f_{v,0},f_{u,1}+=\max(f_{u,0},f_{v,1}) \]為 \(u\le v\)：
\[f_{u,1}+= f_{v,1},f_{u,0}+=\max(f_{u,0},f_{v,1}) \]直接 \(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;
}