等距迴歸

問題

給定值$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