浅谈保序回归问题

说在前头

本篇博客仅仅介绍保序回归问题最简单问题的解法，对于保序回归问题中的许多证明以及拓展延伸都略过不谈，详情请参见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

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#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

复制代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#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;
}