保序回归例子 浅谈保序回归问题

零、引言

昨天晚上，我看这篇论文，竟然看到头疼……

壹、保序回归

1.偏序关系

用 $\preceq$

如果不懂请上度娘搜，而且不懂也可以直接理解成实数的比较。

2.问题描述

给出正整数 $p$ 和代价函数 $(w_i,y_i)\;(1\leqslant i\leqslant n)$ 满足 $w_i>0$ 。给出定义在 $\{v_1,v_2,\dots,v_n\}$ 上的非严格偏序关系 $\preceq$

求实数序列 $\{x_i\}$ 满足 $v_i\preceq v_j$ 则 $x_i\leqslant x_j$，最小化回归代价
$\sum_{i=1}^{n}w_i\cdot|x_i-y_i|^p$

贰、特殊情形下的算法

1.一种贪心算法

当偏序关系为简单的链即 $v_1\preceq v_2\preceq v_3\preceq\cdots\preceq v_n$

首先，定义 $f_S(x)=\sum_{i\in S}w_i|x-y_i|^p$ 。再定义 $M(S)$ 为使得 $f_S(x)$ 最小的 $x$，即 $S$

1.1.理论推导

引理¹一：$f_S(x)$ 的导数（忽略 $x=y_i$ 处的奇点）单增，即 $f_S(x)$

引理二：$f$，即 $M(S)$ 为导数 $f$ 的零点。在 $p>1$ 时零点只有一个；在 $p=1$ 时零点有多个，不妨令 $M(S)$

定理一：若 $M(A)\leqslant M(B)$，则 $M(A)\leqslant M(A\cup B)\leqslant M(B)$

证明：考察导数即可。由 引理一 和 引理二 知 $x=M(A)$ 处 $f_{A\cup B}(x)=f_A(x)+f_B(x)$ 导数为负，故 $M(A)\leqslant M(A\cup B)$，同理可证 $M(A\cup B)\leqslant M(B)$

证明了 定理一，其实问题已经与 $p=1$

推论一：对于集族 $\mathcal F$，有 $\min_{S\in\mathcal F}M(S)\leqslant M(\bigcup_{S\in\mathcal F}S)\leqslant\max_{S\in\mathcal F}M(S)$

定理二：若最优解的 $x_i$ 的某个极长等值段为 $[l,r]$，则 $x_l=M\big([l,r]\big)$

证明：否则可以给 $x_l$ 到 $x_r$ 统一增加充分小的偏移量 $\varepsilon$ 使其逼近 $M\big([l,r]\big)$ 。由 引理一 知，这是更优的。

定理三：原问题最优解为 $x_1=x_2=x_3=\cdots=x_n$，当且仅当 $\forall i\in[1,n),\;M\big([1,i]\big)\geqslant M\big((i,n]\big)$ 。——注意 $M$

证明：必要性毋庸置疑，否则前后可以各自取到最优点 $M$ 。下证充分性。若实际上存在 $x_i\ne x_{i+1}\;(1\leqslant i<n)$，令 $i_0$ 为满足条件的最小 $i$ 。那么 $(i_0,n]$ 一定由若干等值段构成，且每段的值都是 $M$，根据 推论一 可知 $\min{M}\leqslant M\big((i_0,n]\big)\leqslant M\big([1,i_0]\big)=x_{i_0}$ 。根据限制条件，$\min M$ 就是第一个等值段，即 $x_{i_0+1}=\min M\leqslant x_{i_0}$，又由 $i_0$ 的定义 $x_{i_0}\ne x_{i_0+1}$，故其不符合题目要求，矛盾。

定理四：若 $[l,m]$ 和 $(m,r]$ 都满足 定理三 的判定条件，且 $M\big([l,m]\big)>M\big((m,r]\big)$，则 $[l,r]$ 满足 定理三 判定条件。

证明：分割点 $i\in[l,m]$ 时，右侧加上了较小的 $M\big((m,r]\big)$，仍成立；分割点 $i\in(m,r]$ 时，左侧加上了较大的 $M\big([l,m]\big)$，仍成立。

1.2.算法流程

维护单调栈，存储每一个目前得到的 $x_i$ 等值段，以及这些等值段的 $M$

新加入一个数时，首先将其自己视为单独的 $x_i$ 等值段，放入栈中。重复检查栈顶两个元素：若序列中靠前的那个的 $M$ 大于靠后者的 $M$，根据 定理四 必然会合并成大等值段，合并即可。

时间复杂度是 $\mathcal O(n)$ 次计算 $M$ 。当 $p=1$ 时，用 $\tt splay$ 做到均摊 $\mathcal O(\log n)$，当 $p$ 为偶数时使用二分等方式求导数零点，是 $\mathcal O(p\log V)$ 或 $\mathcal O(p)$

2.维护折线算法

在 $\rm CF$ 上被称为 $\textrm{slope trick}$，比较常见，即维护类似绝对值函数的图像。

只适用于树形（递归子问题后合并）的 $p=1$

叁、一般问题的算法

1.整体二分

这又是个让人摸不着头脑的巧妙的 $\rm idea$ 。为了方便叙述，定义 $p\notin(a,b)$ 向 $[a,b]$ 取整的结果是：若 $p\leqslant a$ 则结果为 $a$，否则结果为 $b$

我们构造一个新问题，叫做 $\{a,b\}$ 型问题，表示额外添加限制 $a\leqslant x_i\leqslant b$

定理五：若 $y_i\notin(a,b)$，且存在最优解 $\{x\}$ 使得 $x_i\notin(a,b)$，则 $\{a,b\}$ 型问题的最优解是将 $x_i$ 向 $[a,b]$ 取整后的序列 $\{x_i$

证明：类似 定理二，极大等值连通块的 $x_i$ 应当取到 $M$ 。那么 $\{a,b\}$ 型问题当然要保持靠近 $M$ 的性质，即向 $[a,b]$ 取整呗。显然方案仍然是合法的。²

所以，我们求解 $\{a,b\}$ 型问题后，看看 $x_i=a$ 还是 $x_i=b$ 就知道了原问题中 $x_i\leqslant a$ 还是 $b\leqslant x_i$，实现了二分的目标。

当 $p=1$ 时，由 定理二，$x_i=M\in\{y_1,y_2,\dots,y_n\}$，即存在最优解 $\{x_i\}\subseteqq\{y_i\}$ 。设 $\{y_i\}$ 已排序，求解 $\{y_i,y_{i+1}\}$

当 $p>1$ 时，解可能在实数范围内，故只能在实数域上二分，并且采用 $\{\gamma,\;\gamma{\rm+}\varepsilon\}\;(\varepsilon\to 0^+)$

2.网络流模型

由于 $x_i$ 向 $[a,b]$ 取整后，$x_i\in\{a,b\}$，所以决策只有两种，可以用最大权闭合子图模型解决，偏序关系都可以变成 $+\infty$ 容量的边。当 $p>1$ 时，两种决策的价值差太小，精度误差过大。需要另辟蹊径：价值只在物品上，将一个物品选 $\gamma$ 和选 $(\gamma{+}\varepsilon)$ 的价值都除以 $\varepsilon$，选择方案不变。而建边时，权值是 $f(\gamma{+}\varepsilon)-f(\gamma)$，除以 $\varepsilon$

在 $\gamma=y_i$ 的奇点上，由于 $\varepsilon$ 从正方向逼近 $0$，导数也是有定义的。

3.整数解

类似 定理五，若规定 $x_i\in\Z$，则考虑 $\{a,a{\rm+}1\}$ 型问题，即可知原最优解 $x_i\leqslant a$

证明过程与 定理五 无异：靠近 $M$

肆、代码实现

一般来说，可以直接将 $v_i\preceq v_j$ 直接建成边 $j\to i$

当我们需要对点集 $S$ 进行整体二分时，任取 $i,j\in S$，则 $i,j$ 在有向图上的任意路径上的所有点都在 $S$

有向图中的边，直接变为容量 $+\infty$ 的边，表示 $x_j=\gamma$ 则 $x_i=\gamma$ 。最小割中的 $S$ 为选择 $\gamma$ 的，所以 $i\to T$ 的容量是 $w_i(|\gamma-y_i|^p-|\gamma+\varepsilon-y_i|^p)$ 。——若该值为负，则改为 $S\to i$

例题：限制 $x_i\in\Z$，核心代码如下。

# define _go(i,x) for(int i=head[x]; ~i; i=e[i].nxt)
const int MAXN = 1003, MAXM = 200005;
const int INF = 0x3fffffff;
namespace Network{
	struct Edge{ int to, nxt, val; };
	Edge e[MAXM]; int head[MAXN], cntEdge;
	void _addEdge(int a, int b, int c){
		e[cntEdge].to = b, e[cntEdge].nxt = head[a];
		e[cntEdge].val = c, head[a] = cntEdge ++;
	}
	void addEdge(int a, int b, int c = 0){
		_addEdge(a,b,c), _addEdge(b,a,0);
	}
	void clear(const int &n){
		cntEdge = 0, memset(head+1,-1,n<<2);
	}
	int dis[MAXN], cur[MAXN];
	# define _PARAM_ const int &source,\
		const int &sink, const int &n
	bool bfs(_PARAM_){
		static int que[MAXN];
		memset(dis+1,-1,n<<2); // clear
		dis[source] = 0, que[1] = source;
		for(int x,*fro=que,*bac=que+1; fro!=bac; )
			_go(i,x=*(++fro)) if(e[i].val)
				if(!(~dis[e[i].to])){
					dis[e[i].to] = dis[x]+1;
					++ bac, *bac = e[i].to;
				}
		return dis[sink] != -1;
	}
	int dfs(int x, int inFlow, const int &sink){
		int sum = 0; if(x == sink) return inFlow;
		for(int &i=cur[x]; ~i; i=e[i].nxt)
			if(e[i].val && dis[e[i].to] == dis[x]+1){
				int d = dfs(e[i].to,std::min(
					inFlow-sum,e[i].val),sink);
				e[i].val -= d, e[i^1].val += d;
				if((sum += d) == inFlow) return sum;
			}
		dis[x] = -1; return sum;
	}
	int dinic(_PARAM_){
		int res = 0;
		while(bfs(source,sink,n)){
			memcpy(cur+1,head+1,n<<2);
			res += dfs(source,INF,sink);
		}
		return res;
	}
	bool is_inS(const int &x){ return ~dis[x]; }
}

struct Edge{ int to, nxt; };
Edge e[MAXM]; int head[MAXN], cntEdge;
void addEdge(int a, int b){
	e[cntEdge].to = b, e[cntEdge].nxt = head[a];
	head[a] = cntEdge ++;
}

int v[MAXN], ans[MAXN], haxi[MAXN];
void solve(int *l, int *r, int ql, int qr){
	if(l == r || ql == qr){
		for(int *i=l; i!=r; ++i) ans[*i] = ql;
		return ; // answer found
	}
	const int mid = (ql+qr)>>1;
	for(int *i=l,p=1; i!=r; ++i,++p) haxi[*i] = p;
	const int source = int(r-l)+1, sink = source+1;
	Network::clear(sink); // recalc
	for(int *i=l,p=1; i!=r; ++i,++p){
		int w = (mid-v[*i])<<1|1; // if choose (mid+1)
		if(w > 0) Network::addEdge(source,p,w);
		else if(w) Network::addEdge(p,sink,-w);
		_go(j,*i) if(haxi[e[j].to]) // in the range
			Network::addEdge(p,haxi[e[j].to],INF);
	}
	Network::dinic(source,sink,sink);
	for(int *i=l; i!=r; ++i) haxi[*i] = 0;
	static int _tmp[MAXN]; // auxiliary
	std::copy(l,r,_tmp); // backup
	int *pl = l, *pr = r; ///< refill
	for(int *i=_tmp,p=1; pl!=pr; ++i,++p)
		if(Network::is_inS(p)) *pl = *i, ++ pl;
		else -- pr, *pr = *i; // choose mid+1
	solve(l,pl,ql,mid), solve(pr,r,mid+1,qr);
}

---

1. “引理” 的意思就是，过于显然而略去证明，只引述一下就差不多得了。 ↩︎
2. 论文中，这里用了不少的篇幅，我便猜我的证明多有不严谨之处，希望读者能够简略指出。虽说修锅的可能性不大，顶多让它变为 “参见原论文”。 ↩︎