说在前头
本篇博客仅仅介绍保序回归问题最简单问题的解法，对于保序回归问题中的许多证明以及拓展延伸都略过不谈，详情请参见2018国家集训队论文《浅谈保序回归问题》(高睿泉)。
保序回归问题
定义
给定正整数 (k)，以及一张点集为 (V)，边集为 (E) 的 (DAG) G，点 (i) 有权值 (a_i,w_i)。
请为每个点赋一个点值 (f_i)，使得对于任意边 ((u,v)in E) 有 (f_ule 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})，找一个极小值 ( heta)（( heta) 的大小根据题目要求而定）。考虑求出哪些点的 (f) (in [l,mid])，哪些点的 (fin [mid+ heta,r])，然后对两边递归处理。
构造另一个问题，在这个问题中，我们只考虑 (V) 中的点，且 (V) 中的点的 (f) 取值只能为 ({mid,mid+ heta}) ，在此基础上满足原题限制时求出最小回归代价。
此时不加证明给出一个引理（证明请参考开头的论文）：
若 (solve(V,l,r)) 的最优解为({b_i})，那么将 (b_i) 中所有 (le mid) 的都替换为 (mid)，(>mid) 的都替换为 (mid+ heta) ，得到的新解 ({c_i}) 是新问题的最优解。由此我们只需要解决新问题，就能将 (fle mid) 与 (>mid) 的部分分离开来了。
对于新问题的求解，对于原图中的边 (u ightarrow v)，如果 (u) 取 (mid+ heta)，那么 (v) 必须取 (mid+ heta)。在此基础上最小化回归代价，这就是经典的最小闭合子图问题，可以使用网络流算法来完成。
特殊解法
- 对于树上的情况，可以 (DP) 维护分段函数解决，具体参见论文。
- 当 (p=1) 时，最终答案的取值一定为整数，因此 ( heta) 可以取 (1)。
- 当图为树/基环树时有 (mathcal O(nlog 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) 为价格总和最大的。
(nle 1000,|A|=|B|le 64,c_i<2^{64},v_ile 10^6)。
Solution
考虑找到所有的合法集合，这是非常困难的，但是注意到题目中已经给出了一个合法集合 (A)，考虑在 (A) 的基础上删掉一些物品再加入一些物品得到合法集合。注意到，这样所有在 (A) 的基础上只修改一个物品得到的合法集合权值都比 (A) 大，那么 (A) 就是价格最小的。
因此，我们只考虑改变哪些物品能使得集合依然合法，枚举改变的元素 (i)，再枚举新元素 (j) 插入 (A) 其他元素组成的线性基中，若 (j) 能被表示出来，那么 (i) 的权值必须 (le j)。对 (B) 同样处理。最终我们就得到了一个偏序关系的 (DAG)，直接使用上面的保序回归问题算法即可。由于需要保证修改价格为整数，因此令 ( heta=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
",ans);
return 0;
}
loj#2470. 「2018 集训队互测 Day 2」有向图
Description
树或基环树上的 (L_1) 问题。
(nle 30000,min {n-1,n})。
Solution
首先直接保序回归显然可行，但复杂度太高，考虑优化网络流求解最小闭合子图的过程。
考虑 (DP)，首先对于树的情况，将图看成无向图，设 (f_{u,0/1}) 表示只考虑 (u) 的子树，(u) 权值为 (mid+0/1) 时的最小回归代价，那么有转移：
转移边为 (uge v)；
为 (ule 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(nlog 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
",ans);
return 0;
}