等距回归

问题

给定值$a_{1},...,a_{n}$及其对应的正权重$w_{1},...,w_{n}$，尽可能接近地用$y_{1},...y_{n}$逼近它们，这些值满足$y_{i}\leq y_{j}$类约束。

给定

值$\mathbf {a} \in {\mathbb {R}}^{n}$,

权重$\mathbf {w} \in {\mathbb {R}}^{n}$，使得对于所有$i$都有$w_{i}>0$,

约束集$E\subset \{1,...,n\}^{2}$,

最小化 $\sum _{i}w_{i}\left(y_{i}-a_{i}\right)^{2}$ 相对于 $\mathbf {y} \in {\mathbb {R}}^{n}$，受制于 $y_{i}\leq y_{j}$ 对所有 $(i,j)\in E$

如果所有权重都等于 1，则该问题称为无权单调回归，否则称为有权单调回归。

图 $G=(V=\{1,2,...n\},E=E)$ 必须是无环的。例如，约束 $x_{2}\leq x_{5}$ 和 $x_{5}\leq x_{2}$ 不能同时存在。

例子

最小化 $3(y_{1}-2)^{2}+(y_{2}-1)^{2}+(y_{3}-4)^{2}+2y_{4}^{2}$ 受制于

$y_{1}\leq y_{2}$

$y_{2}\leq y_{3}$

$y_{2}\leq y_{4}$

线性顺序

特别有趣的是线性顺序单调回归。

给定

值$\mathbf {a} \in {\mathbb {R}}^{n}$,

权重$\mathbf {w} \in {\mathbb {R}}^{n}$，使得对于所有$i$都有$w_{i}>0$,

最小化 $\sum _{i}w_{i}\left(y_{i}-a_{i}\right)^{2}$ 相对于 $\mathbf {y} \in {\mathbb {R}}^{n}$，受制于 $y_{1}\leq y_{2}\leq \;...\;\leq y_{n}$

100 个点的线性排序单调回归。

对于线性顺序单调回归，存在一种简单的线性算法，称为#相邻违反者池算法 (PAVA)。

如果所有权重都等于 1，则该问题称为无权线性排序单调回归。

线性顺序单调回归可以被认为是用非递减函数逼近给定的 1 维观测序列。它类似于平滑样条，区别在于我们使用单调性而不是平滑性来消除数据中的噪声。

例子

线性排序单调回归用于 10⁶ 个点。

• 单调回归用于 $10^{6}$ 个点。
• 左边：点 $(x_{i},a_{i})$，其中 $a_{i}=0{\text{ or }}1$。$a_{i}=1$ 的概率由逻辑函数决定。仅显示了 1000 个点。
• 右边：单调回归的结果（黑色曲线）与逻辑函数（红色曲线）对比。逻辑函数以高精度恢复。

线性排序单调回归作为模型

有时，线性排序单调回归被应用于一组观测值 $(x_{i},y_{i}),1\leq i\leq n$，其中 $x$ 是解释变量，y 是因变量。这些观测值按其 $x$ 排序，然后将单调回归应用于 $y$，同时附加约束 $x_{i}=x_{i+1}\Rightarrow y_{i}=y_{i+1}$ 对于所有 $i$。

其他变体

非欧几里得度量

有时使用其他度量代替欧几里得度量，例如 $L_{1}$ 度量

$\sum _{i}w_{i}\left|y_{i}-a_{i}\right|$

或未加权的 $L_{\infty }$ 度量

$\max _{i}\left|y_{i}-a_{i}\right|$

网格上的点

有时，值被放置在二维或更高维度的网格上。拟合值必须沿着每个维度增加，例如

最小化 $\sum _{ij}w_{ij}\left(y_{ij}-x_{ij}\right)^{2}$ 关于 y，受制于

$y_{ij}\leq y_{kl}{\text{ if }}i\leq j,\ k\leq l$

算法

相邻违规者算法

相邻违规者算法 (PAVA) 是一种用于线性阶等距回归的线性时间（和线性内存）算法。

初步考虑

该算法基于以下定理

定理：对于一个最佳解，如果 $a_{i}\geq a_{i+1}$，那么 $y_{i}=y_{i+1}$

证明：假设相反，即 $y_{i}<y_{i+1}$。然后对于足够小的 $\varepsilon$，我们可以设置

$y_{i}^{\mathrm {new} }=y_{i}+w_{i+1}\varepsilon$

$y_{i+1}^{\mathrm {new} }=y_{i+1}-w_{i}\varepsilon$

这会减少总和 $\sum _{i}w_{i}(y_{i}-a_{i})^{2}$ 而不违反约束。因此，我们最初的解不是最优解。矛盾。

由于 $y_{i}=y_{i+1}$，我们可以将两点 $(w_{i},a_{i})$ 和 $(w_{i+1},a_{i+1})$ 合并为一个新点 $\left(w_{i}+w_{i+1},{w_{i}a_{i}+w_{i+1}a_{i+1} \over w_{i}+w_{i+1}}\right)$.

PAVA 算法的一步。

然而，在将两个点 $(w_{i},a_{i})$ 和 $(w_{i+1},a_{i+1})$ 合并到新的点 $\left(w_{i}^{\prime },a_{i}^{\prime }\right)$ 后，这个新点可能会违反约束 $a_{i-1}\leq a_{i}^{\prime }$。在这种情况下，它应该与 $(i-1)$ 个点合并。如果合并后的点违反了它的约束，它应该与前一个点合并，依此类推。

算法

输入

包含 n 个值的数组：initial_values[1] ... initial_values[n]

包含 n 个权重的数组：weights[1] ... weights[n]

输出

名为 results 的数组 (results[1] ... results[n])，通过 i 对 weights[i] * (initial_values[i] - results[i]) ** 2 的总和进行最小化

算法

初始化

pooled_value[1] = initial_values[1]

pooled_weight[1] = weights[1]

num_segments = 1

segment_end[0] = 0

segment_end[1] = 1

对于 current_index = 2 到 n 执行

num_segments++

pooled_value[num_segments] = initial_values[current_index]

pooled_weight[num_segments] = weights[current_index]

当 num_segments > 1 且 pooled_value[num_segments] < pooled_value[num_segments - 1] 时执行

pooled_value[num_segments - 1] =

(pooled_weight[num_segments] * pooled_value[num_segments] + pooled_weight[num_segments - 1] * pooled_value[num_segments - 1]) /

(pooled_weight[num_segments] + pooled_weight[num_segments - 1])

pooled_weight[num_segments - 1] += pooled_weight[num_segments]

num_segments--

segment_end[num_segments] = current_index

对于 segment_index = 1 到 num_segments 执行

对于 value_index = segment_end[segment_index - 1] + 1 到 segment_end[segment_index] 执行

result[value_index] = pooled_value[segment_index]

这里 $S$ 定义了每个新点对应于哪些旧点。

任意情况算法

在任意情况下，这可以作为二次问题解决。最佳算法需要 $\Theta (n^{4})$ 时间，参见

• Maxwell, WL and Muckstadt, JA (1985), "Establishing consistent and realistic reorder intervals in production-distribution systems", Operations Research 33, pp. 1316-1341.
• Spouge J, Wan H, and Wilbur WJ (2003), "Least squares isotonic regression in two dimensions", J. Optimization Theory and Apps. 117, pp. 585-605.

实现

R

isoreg

函数 isoreg 执行非加权线性排序等距回归。它不需要任何包。对于许多简单的情况，它就足够了。

使用示例

x=sort(rnorm(10000))
y=x+rnorm(10000)
y.iso=isoreg(y)$yf
plot(x,y,cex=0.2)
lines(x,y.iso,lwd=2,col=2)

该isoreg函数还将线性排序等距回归实现为模型

x=rnorm(10000)
y=x+rnorm(10000)
y.iso=isoreg(x,y)$yf
plot(x,y,cex=0.2)
lines(sort(x),y.iso,lwd=2,col=2)

Iso

Iso 包包含三个函数

• pava - 线性排序等距回归，加权或非加权。
• biviso - 2-d 等距回归
• ufit - 单峰排序（先递增后递减）

使用示例

install.packages("Iso") # should be done only once
library("Iso") # should be done once per session
x=sort(rnorm(10000))
y=x+rnorm(10000)
y.iso=pava(y)
plot(x,y,cex=0.2); lines(x,y.iso,lwd=2,col=2)

isotone

这是最先进的包。它包含两个函数

• gpava - 线性排序等距回归，加权或非加权，适用于任何指标。类似于isoreg, gpava可以将线性排序等距回归实现为模型。
• activeSet - 适用于任何指标的一般等距回归。

使用示例

install.packages("isotone") # should be done only once
library("isotone") # should be done once per session
x=sort(rnorm(10000))
y=x+rnorm(10000)
y.iso=gpava(y)$x
plot(x,y,cex=0.2); lines(x,y.iso,lwd=2,col=2)

速度比较

由于所有三个库都以某种方式实现了 PAVA 算法，因此我们可以比较它们的速度。

从下面的图形可以看出，对于非加权线性排序等距回归 (LOIR) 和 $n>200$，isoreg应该使用。对于加权 LOIR 和非加权 LOIR 以及 $n\leq 200$，pava应该使用。至于gpava, 它应该只用于非欧几里得指标。

此外，R 上加权简单排序等距回归的实现远非完美。

Java

Weka 是一款免费的机器学习算法集合，用于数据挖掘任务，由 怀卡托大学 开发，包含一个 等距回归分类器。该分类器深深植根于 Weka 的类层次结构中，无法在没有 Weka 的情况下使用。

Python

虽然 scikit-learn 实现等距回归，但 Andrew Tulloch 为线性排序等距回归制作了该算法的 Cython 实现，该实现速度提高了 14 到 5000 倍，具体取决于数据大小。参见 [Speeding up isotonic regression in scikit-learn by 5,000x https://tullo.ch/articles/speeding-up-isotonic-regression/]。如果您只需要代码，请点击 [这里 https://gist.github.com/ajtulloch/9447845#file-_isotonic-pyx]。

用法

校准分类概率的输出

有关更多信息，请参阅 Alexandru Niculescu-Mizil 和 Rich Caruana 的“使用监督学习预测良好概率”。

大多数监督学习算法，例如随机森林 [1]、增强树、SVM 等，擅长预测最可能的类别，但不擅长预测概率。其中一些算法倾向于高估高概率并低估低概率。（一个值得注意的例外是神经网络，它们本身会产生经过良好校准的预测。）

为了获得正确的概率，可以使用未加权线性排序等距回归

${\hat {y}}^{\mathrm {final} }=f({\hat {y}})$

其中

${\hat {y}}^{\mathrm {final} }$ 是概率的最终预测

${\hat {y}}$ 是模型给出的预测

$f$ 是一个非递减函数

为了找到 f，模型在验证集上的预测与输出变量匹配，并将等距回归应用于这些对。

另一种方法是将 ${\hat {y}}$ 通过 sigmoid 函数传递

${\hat {y}}^{\mathrm {final} }={1 \over 1+e^{a-b{\hat {y}}}}$

这被称为 Platt 校准。

对于小型数据集，Platt 校准优于等距回归。从 200 到 5000 个观测值开始，等距回归略微优于 Platt 校准。还要注意，这种等距回归比 Platt 校准所需的逻辑回归更简单、更快。

推荐模型的校准

问题

以下示例适用于推荐引擎，其通用目的是预测用户 $u$ 给项目 $i$ 的评分。

我们有一个模型 M，它可以预测在训练集上训练的另一个模型的残差 r。对于用户或项目数量较少的案例，模型会过拟合，因此需要放松预测

${\hat {r}}_{ui}^{\mathrm {final} }=f({\hat {S}}_{ui}){\hat {r}}_{ui}$

其中

${\hat {r}}_{ui}^{\mathrm {final} }$ 是对用户 $u$、物品 $i$ 的最终预测。

${\hat {r}}_{ui}$ 是模型给出的预测。

$f$ 是一个非递减函数

$S_{ui}$ 可以是 $S_{u}$（训练集中包含用户 $u$ 的观察次数），或 $S_{i}$（训练集中包含物品 $i$ 的观察次数），或 $\min(S_{u},S_{i})$，这取决于模型的性质。字母 $S$ 代表支持。

解决方案

给定来自验证数据库的三元组集合 $({\hat {r}}_{k},r_{k},S_{k})$，我们有

$\min _{f}\sum _{k}\left(r_{k}-{\hat {r}}_{k}f(S_{k})\right)^{2}=\min _{f}\sum _{k}{\hat {r}}_{k}^{2}\left(f(S_{k})-{\frac {r_{k}}{{\hat {r}}_{k}}}\right)^{2}+\operatorname {const}$

因此，对于第 $k$ 个观测值，我们应该设置权重 $w_{k}={\hat {y}}_{k}^{2}$ 和值 $y_{k}=r_{k}/{\hat {r}}_{k}$。然后，三元组 $(S_{k},w_{k},y_{k})$ 按照 $S$ 进行排序，然后对它们应用等距回归。

当一个函数直接预测评分或接受率，而不是预测其他模型的残差时，我们可以定义 $r_{ui}$ 作为该模型与一些更简单模型（例如，项目的平均评分加上用户的平均评分修正）之间的差异。这个模型被称为基本模型。

分析和校准输入变量

通常，我们知道输出变量单调地依赖于输入变量，但除此之外，我们对依赖关系几乎没有先验知识。在这种情况下，等距回归可以提供关于依赖关系的重要线索。

非度量多维标度

给定项目-项目相似性矩阵，多维标度 (MDS) 算法将项目放置到 N 维空间中，以便相似的项目彼此更靠近，而不同的项目则更远离。这对于可视化目的特别有用。

根据相似性和距离之间的关系以及距离的定义，有几种类型的 MDS。对于非度量 MDS，距离可以是相似性的任何单调函数。该函数通常使用等距回归找到，参见维基百科中关于 非度量多维标度 的内容。非度量多维标度在 R 的 MASS 库中作为 [2] 实现。

优点、缺点和后处理

优点

• 非参数
• 速度快（线性时间）
• 简单

缺点

• 区间末端的点有时会有噪声
• 只有当统计量足够大时（$n\gtrapprox 10^{3}$，理想情况下 $n>10^{4}$）时，与其他方法相比，它具有优势。

通过对等距回归的结果进行平滑，可以改善这些缺点。通过这种方式，我们可以兼顾两全其美（平滑性和单调性）。

外部链接

• 等距回归算法，作者：Quentin F. Stout - 对当前最佳可用算法的回顾
• 使用监督学习预测良好概率，作者：Alexandru Niculescu-Mizil、Rich Coruana
• 支持向量机的概率输出以及与正则化似然方法的比较，作者：John C. Platt