算法原理：Iterative Luce Spectral Ranking Algorithm
(1) Luce的选择公理The Axiom of Choice
1959年，Luce在其著作Individual Choice Behavior: A Theoretical Analysis中阐述了选择公理的内容。 p ( i ∣ A ) p(i|A) p(i∣A)指的是从备选集 A A A中选择选项 i i i的概率，给定两个选项 i i i和 j j j，以及任意两个包含 i i i和 j j j的备选集 A A A和 B B B，选择公理假设
p ( i ∣ A ) p ( j ∣ A ) = p ( i ∣ B ) p ( j ∣ B ) \frac{p(i|A)}{p(j|A)}=\frac{p(i|B)}{p(j|B)} p(j∣A)p(i∣A)=p(j∣B)p(i∣B)
即，选择选项 i i i和选项 j j j的概率是独立与其他选项的。根据该假设，可以建议一个独特的参数选择模型，在两两比较（pairwise）的情况下被称为Bradley-Terry模型。假设在多集 D = { ( c l , A l ∥ l = 1 , . . . , d ) } \mathcal{D}=\{(c_l,A_l\|l=1,...,d)\} D={(cl,Al∥l=1,...,d)}中选取 d d d个独立的观测值，每个观测值由一组备选项中的一个选项 c l c_l cl组成。
当 i , j ∈ A l i,j \in A_l i,j∈Al并且 c l = i c_l=i cl=i时，记 i ≻ j i \succ j i≻j为 i i i胜于 j j j。定义有向图为 G D = ( V , E ) G_D=(V,E) GD=(V,E)，如果在 d d d中 i i i至少胜于 j j j一次，则 V = { 1 , . . . , n } V=\{1,...,n\} V={1,...,n}， ( j , i ) ∈ E (j,i) \in E (j,i)∈E。为了明确定义最大似然估计，假设有向图 G D G_D GD是强连接的。
在选择模型下的最大似然估计表示为马尔可夫链的稳态分布。用 π l \pi_l πl表示与选项 c l c_l cl相关的模型参数。对于给定的观测集 D D D，参数 π \pi π的对数似然估计是
log L ( π ∣ D ) = ∑ l = 1 d ( log π l − log ∑ j ∈ A l π j ) \log \mathcal{L}(\pi | \mathcal{D}) = \sum_{\mathcal{l}=1}^{d} \left( \log \pi_{\mathcal{l}} - \log \sum_{j \in A_{\mathcal{l}}} \pi_{\mathcal{j}} \right) logL(π∣D)=∑l=1d(logπl−log∑j∈Alπj)
对于每个选项，我们定义了两个集合。 W i ≐ { l ∣ i ∈ A l , c l = i } W_i \doteq \{l | i \in A_l, c_l = i\} Wi≐{l∣i∈Al,cl=i}为i胜于备选项的集合， L i ≐ { l ∣ i ∈ A l , c l = i } L_i \doteq \{l | i \in A_l, c_l = i\} Li≐{l∣i∈Al,cl=i}为i败于备选项的集合。对数似然函数是严格的凹函数，并且具有唯一的最大似然估计 π ^ \hat{\pi} π^。最优情况下 ▽ π ^ l o g L = 0 \triangledown_{\hat{\pi}}log\mathcal{L}=0 ▽π^logL=0，即
∂ l o g L ∂ π ^ i = ∑ l ∈ W i ( 1 π ^ i − 1 ∑ j ∈ A ) l π ^ j ) − ∑ l ∈ W i 1 ∑ j ∈ A ) l π ^ j = 0 , ∀ i \frac{\partial log\mathcal{L}}{\partial \hat{\pi}_i} = \sum_{l \in W_i}(\frac{1}{\hat{\pi}_i} - \frac{1}{\sum_{j \in A)l}\hat{\pi}_j})- \sum_{l \in W_i}\frac{1}{\sum_{j \in A)l}\hat{\pi}_j} = 0, \forall i ∂π^i∂logL=∑l∈Wi(π^i1−∑j∈A)lπ^j1)−∑l∈Wi∑j∈A)lπ^j1=0,∀i
⟺ ∑ j ≠ i ( ∑ l ∈ W i ∩ L j π ^ j ∑ t ∈ A ) l π ^ t − ∑ l ∈ W i ∩ L j π ^ i ∑ t ∈ A ) l π ^ t ) = 0 , ∀ i \iff \sum_{j \neq i}(\sum_{l \in W_i \cap L_j} \frac{\hat{\pi}_j}{\sum_{t \in A)l}\hat{\pi}_t} -\sum_{l \in W_i \cap L_j} \frac{\hat{\pi}_i}{\sum_{t \in A)l}\hat{\pi}_t} ) =0, \forall i ⟺∑j=i(∑l∈Wi∩Lj∑t∈A)lπ^tπ^j−∑l∈Wi∩Lj∑t∈A)lπ^tπ^i)=0,∀i
函数 f ( S , π ) f(S, \pi) f(S,π)能够接收观测值 S S S和模型参数 π \pi π，其中 S ⊆ D S \subseteq \mathcal{D} S⊆D，然后返回一个非负的实数。
f ( S , π ) ≐ ∑ A ∈ S 1 ∑ i ∈ A π i f(S, \pi) \doteq \sum_{A \in S}\frac{1}{\sum_{i \in A}\pi_i} f(S,π)≐∑A∈S∑i∈Aπi1
令 D i ≻ j ≐ { ( c l , A l ) ∈ D ∣ l ∈ W i ∩ L j } \mathcal{D}_{i \succ j} \doteq \{(c_l,A_l)\in \mathcal{D}|l \in W_i \cap L_j\} Di≻j≐{(cl,Al)∈D∣l∈Wi∩Lj}，即 i i i胜于 j j j的一组观测值，代入公式则为
∑ j ≠ i π ^ i ⋅ f ( D j ≻ i , π ^ ) = ∑ j ≠ i π ^ j ⋅ f ( D i ≻ j , π ^ ) , ∀ i \sum_{j \neq i} \hat{\pi}_i \cdot f(\mathcal{D}_{j \succ i}, \hat{\pi}) = \sum_{j \neq i} \hat{\pi}_j \cdot f(\mathcal{D}_{i \succ j}, \hat{\pi}), \forall i ∑j=iπ^i⋅f(Dj≻i,π^)=∑j=iπ^j⋅f(Di≻j,π^),∀i
以下为该算法的伪代码
Algorithm terative Luce Spectral Ranking
Require: observations
D
\mathcal{D}
D
1:
π
←
[
1
/
n
,
.
.
.
,
1
/
n
]
T
\pi \leftarrow [1/n,...,1/n]^T
π←[1/n,...,1/n]T
2: repeat
3:
λ
←
0
n
×
n
\lambda \leftarrow 0_{n \times n}
λ←0n×n
4: for
(
i
,
A
)
∈
D
(i, A) \in D
(i,A)∈D do
5:
\quad
for
j
∈
A
∖
{
i
}
j \in A \setminus \{i\}
j∈A∖{i} do
6:
λ
j
i
←
λ
j
i
+
1
∑
t
∈
A
π
t
\qquad \lambda_{ji} \leftarrow \lambda_{ji} + \frac{1}{\sum_{t \in A} \pi_t}
λji←λji+∑t∈Aπt1
7:
\qquad
end for
8: end for
9:
π
←
\pi \leftarrow
π← stat. dist. of Markov chain
λ
\lambda
λ
10: until convergence
具有不均匀转移率 λ j i = f ( D i ≻ j , π ) \lambda_{ji} = f(\mathcal{D}_ i\succ j,\pi) λji=f(Di≻j,π)的马尔可夫链，对于开放的概率单纯形中的任何初始分布，收敛到最大似然估计 π ^ \hat{\pi} π^。
假设所有的备选项都具有同等的强度，来近似描述公式7中的马尔可夫链。将 π \pi π固定为 [ 1 / n , . . . , 1 / n ] T [1/n,...,1/n]^T [1/n,...,1/n]T来设置转移率 λ j i ≐ f ( D i ≻ j , π ) \lambda_{ji} \doteq f(\mathcal{D}_{i\succ j},\pi) λji≐f(Di≻j,π)。对于 i ≠ j i \neq j i=j， i i i胜过 j j j对转移率 λ j i \lambda_{ji} λji的贡献是 n / ∣ A ∣ n/|A| n/∣A∣。换句话说，对于每个观察，获胜的项目会得到一个固定数量的传入速率，均匀分配给所有备选项（分配给自身的部分会被舍弃）。我们将稳态分布 π ˉ \bar{\pi} πˉ解释为对模型参数的估计。算法总结了这个过程，称为Luce谱排序（LSR）。如果我们考虑越来越多的观察，LSR会收敛到真实的模型参数 π ∗ \pi^* π∗，即使在备选项集固定的限制情况下也是如此。
令 A = { A l } \mathcal{A}=\{\mathcal{A}_l\} A={Al}为一个备选项的集合，将A划分成任意的两个非空集合 S S S和 T T T，都满足 ( ∪ A ∈ S A ) ∩ ( ∪ A ∈ T A ) ≠ ∅ (\cup_{A \in S}A)\cap(\cup_{A \in T}A) \neq \emptyset (∪A∈SA)∩(∪A∈TA)=∅。令 d l d_l dl为备选集合 A l A_l Al中选项的数量，则当 d l → ∞ d_l \to \infty dl→∞时， π ˉ → π ∗ \bar{\pi} \to \pi^* πˉ→π∗
A A A的条件确保了在渐进的情况下 G D G_{\mathcal{D}} GD是强连通的。如果选项 i i i和 j j j至少在一个备选集合中进行了比较，那么转移率的比率满足 l i m d l → ∞ λ i j / λ j i = π j ∗ / π i ∗ lim_{d_l \to \infty}\lambda_{ij}/\lambda_{ji}=\pi^*_j/\pi^*_i limdl→∞λij/λji=πj∗/πi∗
Luce选择模型广泛应用于成对比较组成的数据，由于稳态分布对时间尺度的变化是不变的，我们可以重新调整转移率，并在成对比较数据上使用LSR时设置 λ j i ≐ ∣ D i ≻ j ∣ \lambda_{ji}\doteq|\mathcal{D}_{i \succ j}| λji≐∣Di≻j∣，设 S S S为至少进行过一次比较的选项对的集合，在每对 ( i , j ) ∈ S (i,j)\in S (i,j)∈S进行了 p p p次比较的情况下，LSR严格等价于连续时间马尔可夫链形式，且Maystre et al.(2015)证明了将转移率考虑为胜利计数的比例是合理的。
2 Source Code
仅说明对pair- comparison类型数据的处理
一、导包
import functools
import numpy as np
import pdb
from choix.convergence import NormOfDifferenceTest # 用于执行两个向量或数据集之间差异的范数测试
from choix.utils import exp_transform, log_transform, statdist # 指数变换、对数变换、计算状态分布
二、定义一个函数，用于初始化最小二乘回归的马尔可夫链和权重
def _init_lsr(n_items, alpha, initial_params):
'''
:param n_items: 表示比较的项目数量
:param alpha: 表示初始化马尔可夫链的权重
:param initial_params: 表示回归模型的初始参数
'''
if initial_params is None:
# ones函数：创建一个指定形状的数组，并将其中的元素都初始化为 1
# 创建了一个长度为 n_items 的数组，并将所有元素设置为 1
weights = np.ones(n_items)
else:
weights = exp_transform(initial_params) # 进行指数变换
# 创建一个形状为 (n_items, n_items) 的二维数组 chain， 其中每个元素的值均为 alpha
# chain是马尔可夫链的转移矩阵
chain = alpha * np.ones((n_items, n_items), dtype=float)
return weights, chain
alpha实际上是一个正则化系数，当马尔可夫链不满足常返性、周期性和连通性时，就可能无法到达唯一稳态，alpha就是解决这种情况的问题的，加上正则化系数，就可以继续进行计算。
三、定义一个函数，用于迭代地细化最小二乘回归估计（Least Squares Regression estimate），直到收敛为止
def _ilsr(fun, params, max_iter, tol):
pdb.set_trace()
'''
:param fun:表示用于迭代估计的函数或方法
:param params:表示回归模型的参数，是一个可迭代对象
:param max_iter:表示最大迭代次数
:param tol:表示收敛的容差阈值-用于判断算法是否达到收敛状态的标准，当迭代过程中的参数变化小于或等于容差阈值时，可以认为算法已经收敛。
'''
# 用于判断迭代过程是否已经收敛
# order = 1 表示使用一阶范数进行差异比较。
converged = NormOfDifferenceTest(tol, order=1)
for _ in range(max_iter): # max_iter 最大迭代次数
params = fun(initial_params=params)
if converged(params):
return params
# 表示算法未能在指定次数内收敛
raise RuntimeError("Did not converge after {} iterations".format(max_iter))
三、定义一个函数，用于计算最小二乘回归（Least Squares Regression）模型参数的估计
def lsr_pairwise(n_items, data, alpha=0.0, initial_params=None):
# 调用_init_lsr进行初始化，获取初始权重和马尔可夫链
weights, chain = _init_lsr(n_items, alpha, initial_params)
# 通过遍历data中的每个比较对（winner, loser），根据比较结果更新马尔可夫链的转移率。通过增加对应链中的元素，以及使用权重进行归一化，来反映比较结果。
for winner, loser in data:
chain[loser, winner] += 1 / (weights[winner] + weights[loser])
# 从转移率矩阵中减去对角线上的元素之和，用于使转移率矩阵满足概率分布的性质
chain -= np.diag(chain.sum(axis=1))
# 通过 statdist 函数计算马尔可夫链的稳态分布，并对其进行对数变换。返回对数变换后的稳态分布作为参数估计的结果。
return log_transform(statdist(chain))
四、定义一个函数，用于使用迭代 Luce Spectral Ranking 算法（I-LSR）计算模型参数的最大似然估计（ML estimate）
def ilsr_pairwise(
n_items, data, alpha=0.0, initial_params=None, max_iter=100, tol=1e-8):
fun = functools.partial(
lsr_pairwise, n_items=n_items, data=data, alpha=alpha)
return _ilsr(fun, initial_params, max_iter, tol)
2 原理及逻辑
- 调用ilsr_pairwise函数
# 对一百条数据进行choix-ilsr_pairwise排序
list_arrange=ilsr_pairwise(100,pair_list,alpha=0.001)
这行代码调用了函数ilsr_pairwise
def ilsr_pairwise(
n_items, data, alpha=0.0, initial_params=None, max_iter=100, tol=1e-8):
fun = functools.partial(
lsr_pairwise, n_items=n_items, data=data, alpha=alpha)
return _ilsr(fun, initial_params, max_iter, tol)
ilsr_pairwise(
100, data, alpha=0.001, initial_params=None, max_iter=100, tol=1e-8)
传入参数：
进行比较的项目的数量n_item=100
这100个项目的两两比较的结果列表data
正则化参数alpha，避免马尔可夫链无法达到唯一稳态
在data列表中，每个元素都是一个元组
(i,j)，i和j是这100个item的编号，(i,j)就代表i wins over j，即i是winner，j是loser
- 调用lsr_pairwise函数
在函数ilsr_pairwise中，调用了lsr_pairwise计算最小二乘回归（Least Squares Regression）模型参数的估计
# fuctool.partial用于创建部分函数对象，相当于调用该函数
# 第一个参数是函数对象lsr_pairwise
# 其他参数传递的是lsr_pairwise的参数
fun = functools.partial(
lsr_pairwise, n_items=n_items, data=data, alpha=alpha)
部分函数是指通过固定（或部分应用）一个或多个函数的参数，从而创建一个新的函数。这个新函数可以在后续的调用中，只需要提供剩余的参数即可。
- 进入lsr_pairwise函数
def lsr_pairwise(n_items, data, alpha=0.0, initial_params=None):
# 调用_init_lsr进行初始化，获取初始权重和马尔可夫链
weights, chain = _init_lsr(n_items, alpha, initial_params)
# 通过遍历data中的每个比较对（winner, loser），根据比较结果更新马尔可夫链的转移率。通过增加对应链中的元素，以及使用权重进行归一化，来反映比较结果。
for winner, loser in data:
chain[loser, winner] += 1 / (weights[winner] + weights[loser])
# 从转移率矩阵中减去对角线上的元素之和，用于使转移率矩阵满足概率分布的性质
chain -= np.diag(chain.sum(axis=1))
# 通过 statdist 函数计算马尔可夫链的稳态分布，并对其进行对数变换。返回对数变换后的稳态分布作为参数估计的结果。
return log_transform(statdist(chain))
# 调用_init_lsr进行初始化，获取初始权重和马尔可夫链
weights, chain = _init_lsr(n_items, alpha, initial_params)
- _init_lsr函数
weights, chain = _init_lsr(100, alpha=0.001, initial_params=None)
def _init_lsr(n_items, alpha, initial_params):
if initial_params is None:
# ones函数：创建一个指定形状的数组，并将其中的元素都初始化为 1
# 创建了一个长度为 n_items 的数组，并将所有元素设置为 1
weights = np.ones(n_items)
else:
weights = exp_transform(initial_params) # 进行指数变换
# 创建一个形状为 (n_items, n_items) 的二维数组 chain， 其中每个元素的值均为 alpha
# chain是马尔可夫链的转移矩阵
chain = alpha * np.ones((n_items, n_items), dtype=float)
return weights, chain
因为initial_params=None，所以会创建一个长度为n_item=100，元素全为1的数组[1,1,...,1]作为weight
weight表示在线性回归中每个变量的权重（初始化为1）
创建一个100行100列，所有元素都是alpha=0.001的马尔可夫链的转移矩阵
[0.001,0.001,...,0.001
0.001,0.001,...,0.001
......
0.001,0.001,...,0.001]
- 创建完weight和chain后返回
lsr_pairwise继续执行
for winner, loser in data:
chain[loser, winner] += 1 / (weights[winner] + weights[loser])
loser和winner实质上都是item的编号，[loser,winner]实际上是矩阵的坐标，实际上就是更新马尔可夫链的转移概率，以考虑胜者和输者的权重值。
目的是使权重较大的项目在马尔可夫链中的转移更加概率高
# 从转移率矩阵中减去对角线上的元素之和，用于使转移率矩阵满足概率分布的性质
chain -= np.diag(chain.sum(axis=1))
# 通过 statdist 函数计算马尔可夫链的稳态分布，并对其进行对数变换。返回对数变换后的稳态分布作为参数估计的结果。
return log_transform(statdist(chain))
chain.sum(axis=1) 计算了转移矩阵每一行的元素之和，得到一个一维数组。然后，np.diag() 函数将这个一维数组转换为对角矩阵，其中对角线上的元素与原始数组的元素相同
statdist(chain) 是计算马尔可夫链的稳定分布（stationary distribution）。稳定分布是指在长时间运行后，马尔可夫链中每个状态的概率分布。它表示了马尔可夫链在长期运行后，各个状态的平稳分布情况。
log_transform() 函数对 statdist(chain) 的结果进行对数变换。这个变换可能是为了处理概率值较小的情况，或者为了将结果转换为更合适的尺度。
- 跳出
lsr_pairwise，返回ilsr_pairwise的fun
马尔可夫链达到稳态分布之后经过对数变换得到的权重值数组返回到fun函数
return _ilsr(fun, initial_params, max_iter, tol)
因为最后ilsr_pairwise返回的值是经过_ilsr计算的
- 进入_ilsr函数
def _ilsr(fun, params, max_iter, tol):
# 用于判断迭代过程是否已经收敛
# NormOfDifferenceTest 类是一个测试条件的实现，用于检查迭代过程中的差分或差分方程的误差是否达到了给定的收敛容忍度
# order = 1 表示使用一阶范数进行差异比较
converged = NormOfDifferenceTest(tol, order=1)
for _ in range(max_iter): # max_iter 最大迭代次数
# 更新params，进行迭代，使马尔可夫链不断接近稳态
params = fun(initial_params=params)
if converged(params):
return params
# 表示算法未能在指定次数内收敛
raise RuntimeError("Did not converge after {} iterations".format(max_iter))
fun返回的结果是每个元素的权重值
params是initial_params=None
最大迭代次数max_iter是默认值100
收敛的容差阈值tol是默认值1e-8，用于判断是否到达收敛状态
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
