生成式 AI 模型的排位算法

sigmoid 函数来自于对对数几率（log-odds）的线性建模，也就是 logistic 回归，假设样本 x 有三个特征 x1,x2,x3 ，对每条样本进行 0/1 二分类，用 logistic 回归建模的话，x 被分类为 1 的概率是：

\[ p=\frac{1}{1+e^{-(w_1x1+2_2x2+w_3x3+b)}} \]

整理后：

\[ \ln (\frac{p}{1-p})=w_1x1+w_2x2+w_3x3+b \]

左边的 p/(1-p) 就是几率，它表示分类为 1 的概率与分类为 0 概率的比值，放在更经典的概率场景里就是：赌赢的概率和输的概率的比值 – 赔率，取对数就被称为“对数几率”。

在 Bradley-Terry 模型下，右边线性方程被替换成了两个选手的实力差，分类概率 p 对应的是 i 战胜 j 的概率 \( p_{ij} \) ：

\[ \ln \frac{p_{ij}}{1 - p_{ij}} = r_i - r_j \]

注意它和 logistic 回归的公式在语义上没有任何差别，只不过这里是对单个具体样本（比赛）的关系描述。

假设只有 3 个选手，可以把它写成更一般形式：

\[ \ln \frac{p}{1 - p} = r_1x_1 + r_2x_2 + r_3x_3 + 0 \tag{1} \]

对比一般 logisitic 公式（把权重由常见的 w 符号替换成 r）：

\[ \ln \frac{p}{1-p}=r_1x1+r_2x2+r_2x3+b \tag{2} \]

可以看到 Bradley-Terry 模型只是模型权重和数据采样过程做了一些限制：

• 偏置 b 先验为 0

• 数据集里仅包括满足以下条件的数据（或者理解为采样的时候只筛选这类数据作为合法数据进行建模）：

  只有两个非 0 特征，且一个为 1 另一个为 -1, 自然也有 sum([x1,x2,…xn])=0 。

因此 Bradley-Terry 模型只是 Logistic 回归模型的一种特殊情况，完全可以放在一般 logistic 模型下进行求解。但这些限制使得公式经过变换后可以得到更具体的解释视角，也可以设计出与现实需求更贴合的优化算法，这才有后文的故事。

         ----------------------------
        /                            \
       /      Logistic Regression     \
      /      -----------------------   \
     /      /                       \   \
    /      /  Bradley-Terry Model    \   \
   /      /                           \   \
  /      ------------------------------\   \ 
 /                                          \
---------------------------------------------

比如对公式

\[ \ln \frac{p_{ij}}{1 - p_{ij}} = r_i - r_j \]

整理后：

\[ p_{ij} = \frac{1}{1 + e^{-(r_i - r_j)}} = \frac{e^{r_i - r_j}}{1 + e^{r_i - r_j}} = \frac{e^{r_i}}{e^{r_i} + e^{r_j}} \]

这可以解释成：在已知选手 i,j 的实力得分情况下，选手 i 和选手 j 对战时，获胜的概率就是对二者的实力值进行 softmax 归一化后的结果。因为有以上表达式，可以定义 \( \pi_{i} = e^{r_{i}} \) 也作为选手实力值（后文称为强度，strength），这样就有：

\[ p_{ij}=\frac{\pi_{i}}{\pi_{i}+\pi_{j}} \]

后文中，有些排位估计算法是去估计 \( \pi_{i} \) ，有些算法则是直接去估计 \( r_{i} \) ，但二者仅仅只是对同一个量在对数和指数尺度下的不同观测。

以上模型是过度参数化的，意味着如果我们对所有 \( r_i \) 加上一个相同的常数 \( c \)，模型的预测结果不会改变，因为两个队伍之间的差异 \( r_i - r_j \) 仍然保持不变。

为了避免这个问题，可以在计算出结果后将某个选手的实力值校准到特定得分上（比如最低分设置为 1000 分），或者把选手的平均水平校准到特定得分，比如取所有 \( r_{i} \) 的算术平均 \( \frac{1}{N}\sum_{i=1}^{N}r_{i} \) 为 0 ，对应到 \( \pi_{i} \) 上，就是所有 \( \pi_{i} \) 的几何平均 \( (\prod_{i=1}^{N}\pi_{i})^{\frac{1}{N}} \) 为 1.

问题是，从 logistic 回归模型层面，是什么限制使得 logistic 模型变得过参数化了？ 是权重 b=0 的先验吗？看上去不是，提前解出了一个参数（或者固定某些参数进行微调），并不会导致其他参数变得冗余。那么只可能是数据采样机制限制导致的：还是以三个选手为例，数据集的每个样本 (x1,x2,x3) 满足 sum(x1,x2,x3)=0 （当然实际限制更强，三个数只能取 -1,1 和 0 的某个排列，但对于本文问题和为 0 的限制就足够解释了）, 它如何影响到模型参数呢？ 因为 r1x1+r2x2+r3x3 可以写成 r1x1+r2x2+r3(-x1-x2)= (r1-r3)x1+(r2-r3)x2, 因此独立变量只有 (r1-r3) 和 (r2-r3), 或者说是输入数据某个特征本身是线性冗余的，放到线性模型里，对应特征也冗余了。

总之，对于 N 个选手，要估计的并不是 N 个未知参数，而是 N-1 个。

接下来介绍的两个算法是估计 \( \pi_{i} \) 的。

有了对每场比赛胜负的概率表示，那么当看到具体的 m 场比赛结果之后，可以用极大似然估计来计算出决定概率的那些参数了。

假设对任意选手 i 和选手 j (j!=i）, i 战胜 j 的次数是 \( w_{ij} \) ，那么似然函数就是：

\[ L(r_1, r_2, \dots, r_n) = \prod_{(i,j)} p_{ij}^{w_{ij}} \]

注意，如果要按一般的 logistic 回归里似然函数的形式如下（并非正确）：

\[ L(r_1, r_2, \dots, r_n) = \prod_{(i,j)} p_{ij}^{w_{ij}}(1-p_{ij})^{l_{ij}} \]

这里 \( l_{ij} \) 表示 i 输给 j 的次数，但实际 \( w_{ji} \) 就表示 i 输给 j 的次数（在对所有 i,j 遍历的时候已经覆盖了），因此我们只需要考虑对“赢”建模，这其实是 Bradley-Terry 模型是一般 logistic 模型的特例在似然函数上的体现。

由于 logistic 回归的似然函数是凸的，有唯一最优解，因此以上函数自然也有该性质。

继续计算其对数似然函数为： \[ \ln L(r_1, r_2, \dots, r_n) = \sum_{(i,j)} w_{ij} \ln p_{ij} \]

展开 \( p_{ij} \) 得到：

\[ \ln L(r_1, r_2, \dots, r_n) = \sum_{(i,j) } \left[ w_{ij} \ln \left( \frac{\pi _{i}}{\pi_{i} + \pi_{j}} \right) \right] = \sum_{(i,j) } \left[ w_{ij} \ln \pi_i-w_{ij}\ln (\pi_{i}+\pi_{j}) \right] \]

为了求最大值，可以对每个 \( \pi_k \) 求偏导数，用 \( \pi_{k} \) 表示任意某个选手得分是因为对于固定的 k, 以上式子可以拆分成两个：

\[ \ln L(r_1, r_2, \dots, r_n) = \sum_{(i=k,j) } \left[ w_{kj} \ln \pi_k-w_{kj}\ln (\pi_{k}+\pi_{j}) \right] + \sum_{(i,j=k) } \left[ w_{ik} \ln \pi_i-w_{ik}\ln (\pi_{i}+\pi_{k}) \right] \]

因此可以更清晰求导：

\[ \frac{\partial \ln L}{\partial \pi_k} = \sum_{j} \left[ w_{kj} \left( \frac{1}{\pi_{k}} \right) - w_{kj} \left( \frac{1}{\pi_k + \pi_j} \right) \right] + \sum_{i} \left[ - w_{ik} \left( \frac{1}{\pi_i + \pi_k} \right) \right] \]

此时式子的 i 可以写成 j, 与前一个式子合并：

\[ \frac{\partial \ln L}{\partial \pi_k} = \sum_{j} \left[ w_{kj} \left( \frac{1}{\pi_{k}} \right) - (w_{kj} + w_{jk}) \left( \frac{1}{\pi_k + \pi_j} \right) \right] \]

可以用随机梯度下降或其他基于梯度的优化算法来估计最优解，解出之后在用上一节提到的方法校准。

Bradley–Terry model - Wikipedia 提到了 Zermelo 迭代法，思路是先将以上梯度置为 0 ，得到一个递归的关于 \( \pi_{k} \) 的式子： \[ \pi_{k} = \frac{\sum_{j}w_{kj}}{\sum_{j}(w_{kj}+w_{jk})/(\pi_{k}+\pi_{j})} \]

然后随机初始化所有 \( \pi_{k} \), 每轮迭代用以上式子计算出一组新的 \( \pi_{k} \) ，然后几何平均归一化： \[ \pi_{k} = \frac{\pi_{k}}{(\prod_{i}\pi_{i})^{\frac{1}{n}}} \]

不断迭代就可以逐渐逼近最优解。（如何证明见下节的说明）

以上迭代算法是 1929 年 Zermelo 发现的，缺点是收敛效率低。

在近 100 年后的 2023 年，Newman 在 Efficient Computation of Rankings from Pairwise Comparisons 论文中发现对求偏导后置 0 的式子做一点变换：

\[ \frac{\partial \ln L}{\partial \pi_k} = \sum_{j} \left[ w_{kj} \left( \frac{1}{\pi_{k}} \right) - (w_{kj} + w_{jk}) \left( \frac{1}{\pi_k + \pi_j} \right) \right] = 0 \]

\[ \sum_{j} w_{kj} \frac{1}{\pi_{k}} - \sum_{j} (w_{kj} + w_{jk}) \left( \frac{1}{\pi_k + \pi_j} \right) = 0 \]

等价于：

\[ \frac{1}{\pi_{k}} \sum_{j} w_{kj}\left( \frac{\pi_k + \pi_j}{\pi_k + \pi_j} \right) - \sum_{j} (w_{kj} + w_{jk}) \left( \frac{1}{\pi_k + \pi_j} \right) = 0 \]

\[ \frac{1}{\pi_{k}} \sum_{j} w_{kj}\left( \frac{\pi_k}{\pi_k + \pi_j} \right) + \frac{1}{\pi_{k}} \sum_{j} w_{kj}\left( \frac{\pi_j}{\pi_k + \pi_j} \right) - \sum_{j} (w_{kj} + w_{jk}) \left( \frac{1}{\pi_k + \pi_j} \right) = 0 \]

\[ \sum_{j} \left( \frac{w_{kj}}{\pi_k + \pi_j} \right) + \frac{1}{\pi_{k}} \sum_{j} w_{kj}\left( \frac{\pi_j}{\pi_k + \pi_j} \right) - \sum_{j} (w_{kj} + w_{jk}) \left( \frac{1}{\pi_k + \pi_j} \right) = 0 \]

\[ \frac{1}{\pi_{k}} \sum_{j} w_{kj}\left( \frac{\pi_j}{\pi_k + \pi_j} \right) - \sum_{j} w_{jk} \left( \frac{1}{\pi_k + \pi_j} \right) = 0 \]

最终得到

\[ \pi_{k} = \frac{\sum_{j}w_{kj}\pi_{j}/(\pi_{k}+\pi_{j})}{\sum_{j}(w_{jk})/(\pi_{k}+\pi_{j})} \]

然后同样对所有 \( \pi_{k} \) 的几何平均归一到 1

\[ \pi_{k}=\frac{\pi_{k}}{(\prod_{j}^{n}\pi_{k})^{\frac{1}{n}}} \]

论文证明这个迭代方式会收敛到最优解（实际证明了一系列迭代算法的收敛性，包括 Zermelo 算法，并在经验上展示该算法在一些迭代算法中收敛速度是最快的）

前文只是说选手 i,j 的得分差和胜率的对数几率是线性的，这和 logistic 回归类似，但没有解释直觉上这意味着什么， Why does the Elo rating system work? - Mathematics Stack Exchange 中有一个回答对此做了解释，个人觉得可以作为参考：

Elo 的假定是：

• 如果选手 i 击败选手 j 的概率是选手 j 击败选手 i 的概率的 10 倍，那么选手 i 的 Elo 排名应该比选手 j 高 400 分。
• 同理，假设选手 i 击败选手 j 的概率是选手 j 击败选手 k 的概率的 10 倍，那么选手 i 击败选手 k 的概率将是选手 k 击败选手 i 的 100 倍。因此，选手 i 的 Elo 排名应该比选手 k 高 800 分。

具体地，如果两个玩家的 Elo 排名差距为 400 分，那么这会导致玩家间的胜负概率比例为 10：1 ，写成公式就是：

\[ \text{Odds(i > j)} = 10^{(r_i - r_j) / 400} \]

\[ p_{ij} = \frac{1}{1 + 10^{(r_j - r_i) / 400}} - \frac{1}{1 + e^{-(r_i - r_j)\ln 10 / 400}} \]

因此这和 Bradley-Terry 模型的假定是一样的，只不过有更具体的 ln10/400 这种缩放因子（或者可以说是 temperature）。

前文说到，Elo 算法的关键特点是局部的随机梯度下降，这种算法下，我们不需要考虑多场比赛的胜负集合，而只需要关心当前比赛，也就是以下一个式子：

\[ p_{ij} = \frac{1}{1 + 10^{(r_j - r_i) / 400}} \]

这个函数就是单次比赛结果下的似然函数，为了更方便计算其对数似然的梯度，先引入一些符号：

定义 \( \sigma(x)=\frac{1}{1+e^{-x}} \) , 这是标准的 sigmoid 函数，它有一些比较好的方便计算梯度的性质：

• \( \sigma \) 对 x 的的导数是 \( \sigma (1-\sigma) \),
• \( 1-\sigma(x)=\sigma(-x) \)

再定义： \[ \sigma'(x)=\frac{1}{1+10^{-x/400}}=\sigma(\frac{\ln10}{400}x)=\sigma(Zx) \]

这里用 Z 来表示 (ln10)/400 ， \( \sigma ' \) 的性质是：

• \( \sigma' \) 函数对 x 的导数是 \( Z\sigma' (1-\sigma') \),
• \( 1-\sigma'(x)=\sigma'(-x) \)

那么单场比赛的概率可以写成： \( p_{ij}=\sigma'(r_{i}-r_{j}) \)

如果 i 赢了， \(\ln L= \log( \sigma'(r_{i}-r_{j})) \) 就是这场比赛结果的对数似然，计算梯度有：

\[\frac{\partial \ln L}{\partial r_i}= Z(1-\sigma'(r_{i}-r_{j}))\]

\[ \frac{\partial \ln L}{\partial r_j}=-Z\sigma'(r_{j}-r_{i}) \]

利用梯度提升（也就是对负对数似然 -L 做梯度下降），选手 i 的得分一次更新之后就是： \[ r_{i} = r_{i} + \eta Z(1- \sigma'(r_{i}-r_{j})) \]

选手 j 的得分一次更新之后就是：

\[ r_{j} = r_{j}+\eta Z(0- \sigma' (r_{j}-r_{i}))\]

以上区分了 i 和 j 的，但如果令 \( S_{i}=1 \) 表示 i 胜利， \( S_{i}=0 \) 表示 i 失败（j 胜利），并且 \( K=\eta Z \), 更新规则就可以只用一个行紧凑地表示：

\[ r_{i} = r_{i} + K(S_{i}- \sigma'(r_{i}-r_{j})) \]

注意如果去做多步梯度下降，那么就纯粹在单个样本上优化，它不等同于对多场比赛的全局负对数似然进行优化，但如果只朝着以上目标走一小步，那么这一步的方向大致上也是指向全局最优的方向，所以它实质就是随机梯度下降法。

和神经网络的随机梯度下降的一个不同点在于，Elo 中一个样本只能更新两个参数： \( r_i \) 和 \( r_j \) ，这是因为每个样本的特征是稀疏的，即只有 \( r_{i} \) 和 \( r_{j} \) 对应了 1 和 -1, 其他特征都是 0，但神经网络即便只有一个样本，也会更新所有参数（除模型内故意为之的随机 dropout 或者 lora 等微调部分参数的方案）

以下是其他参考资料里更常见的分数调整公式： \[ R' = R + K (S - E) \]

其中：

• \( E = \sigma'(r_{i}-r_{j}) \)
• \( K \) 是动态调整的常数，对于参加比赛少的新手，K 的值会更大，比如 32, 而参加的比赛越多，分数应该更趋于稳定，因此更新幅度会更低，比如 16，这类似于神经网络训练里的学习率衰减策略。
• \( S \) 是实际的比赛结果，1 代表胜利，0.5 代表平局，0 代表失败，这里把平局的情况也用线性插值给加进来了。

由于 K 的选择完全是经验性的，而且 ELo 本身就是对全局似然最大化的近似算法，按比赛的不同顺序去更新计算得分得到的结果也会不同， 因此这个算法的稳定性比 Bradley-Terry 批量迭代算法更差。

def count_rank(results):
    win_count = np.sum(results, axis=1)
    fail_count = np.sum(results, axis=0)
    win_rate = win_count / (win_count + fail_count)
    return win_rate
count_rank(results)

array([0.3       , 0.61538462, 0.33333333, 0.77777778])

这里的排序是： D>B>C>A

这种做法可能遇到类似田忌赛马的问题，比如假设 B 是最弱的，A 比 B 强，C,D 实力远高于 A,B，那么如果 A 同 B 对战了 100 次且每次都是 A 胜利，A 同 C,D 分别比赛了 3 次且都输了，除此之外都是 C 和 D 的比赛，胜负比例是 1:1, 那么这种情况 C 和 D 的胜率都是 50% 多一点，而 A 的胜率接近 100% 。 换句话说，这种算法使得 A 对战 C,D 失败的信息没法传递到 A 对和 B 的对战信息上，使得 A 仅靠着欺负弱小的胜率而排到了第一。（这也有点类似辛普森悖论，在局部统计上，C,D 强于 A, 但综合统计 A 强于 C,D）

当然，如果模型之间比赛的匹配是相对均匀的，那么这种传递性问题在数据层面就被消除了（类似随机控制试验消除因果混淆？），用该方法预估得分并排序是可行的，在一些足球联赛中，每个球队和其他球队都有固定的交手次数，最后统计胜率就可以了（平局后再统计净胜球）。 但在 AI 评估场景下，由于每天搜集到很多数据，并且新模型频繁发布，旧模型也可能下线， 那么模型之间的交战次数无法保证均匀（考虑到用户体验，也会尽可能少地匹配较差的模型给用户），而且这种方式要求每次构造一个大矩阵重新计算，在榜单经常更新的情况下计算量也很大。

def elo_rank(results, k=4, reversed=False):
    n = len(results)
    order = []
    for i in range(n):
        for j in range(n):
            order.extend([(i, j)] * results[i][j])
    if reversed:
        order = order[::-1]
    ratings = np.ones(n) * 1000
    for i, j in order:
        e = 1 / (1 + 10 ** ((ratings[j] - ratings[i]) / 400))
        ratings[i] += k * (1 - e)
        ratings[j] += k * (0 - (1 - e))
    print(ratings)
    delta = ratings[0] - 1000
    ratings -= delta
    return ratings

elo_rank(results, k=4, reversed=False)

[ 991.90359033 1005.53773269  992.7547053  1009.80397168]
array([1000.        , 1013.63414236, 1000.85111497, 1017.90038135])

以上算法更新中固定 K 的值，并把 A 的得分校准到了 1000, 得到的排序是 D>B>C>A, 这和上一种计算方法得到结果是一样的（说明之前对战安排相对是公平的）

Elo 的不稳定性在于 K 的选择和比赛顺序：

比如修改 k 值，得分就不会不一样（尽管有时候总体排序还是一样）

elo_rank(results, k=10, reversed=False)

array([1000.        , 1032.7481022 , 1004.99167616, 1044.40869968])

以下修改比赛顺序，排序就发生了变化： D>B>A>C

elo_rank(results, k=4, reversed=True)

array([1000.        , 1012.61767187,  999.27397405, 1016.50208449])

不过 A,C 水平本身是接近的。

参数 K 和比赛顺序对 Elo 的稳定性详细分析可以参考： NeurIPS Poster Elo Uncovered: Robustness and Best Practices in Language Model Evaluation

更新公式为：

\[ \pi_{k} = \frac{\sum_{j}w_{kj}}{\sum_{j}(w_{kj}+w_{jk})/(\pi_{k}+\pi_{j})} \]

然后每一轮更细后要对所有 \( \pi_{k} \) 的几何平均归一到 1 \[ \pi_{k}=\frac{\pi_{k}}{(\prod_{j}^{n}\pi_{j})^{\frac{1}{n}}} \]

其实现为：

def bt_rank_1929(results, epochs=20):
    n = len(results)
    ratings = np.ones(n)
    for step in range(epochs):
        for k in range(n):
            numerator = sum(results[k])
            denominator = 0
            for j in range(n):
                if j == k:
                    continue
                denominator += (results[k][j] + results[j][k]) / (
                    ratings[k] + ratings[j]
                )
            ratings[k] = numerator / denominator
        ratings = ratings / np.exp(np.mean(np.log(ratings)))
    return ratings

print(bt_rank_1929(results, epochs=3))
print(bt_rank_1929(results, epochs=20))
print(bt_rank_1929(results, epochs=30))
print(bt_rank_1929(results, epochs=40))

[0.63693151 1.21910418 0.69658672 1.84880562]
[0.63987084 1.0440904  0.65996953 2.26801384]
[0.63983646 1.04334968 0.65981745 2.27026909]
[0.63983489 1.04331601 0.65981053 2.27037173]

排序结果是 D>B>C>A, 且 30 轮以上收敛到最优结果

引入 normalize2elo 函数，根据 \( \sigma'(x)=\frac{1}{1+10^{-x/400}}=\sigma(\frac{\ln10}{400}x) \) 且 \( \pi=e^{r} \) 这层关系将得分转到 Elo 尺度，并且将 A 的得分固定为 1000:

def normalize2elo(ratings):
    elo_scores = 400 / np.log(10) * ratings
    delta = elo_scores[0] - 1000
    elo_scores -= delta
    return elo_scores

print(normalize2elo(np.log(bt_rank_1929(results, epochs=40))))

[1000.         1084.9391811  1005.34052848 1220.01162127]

\[ \pi_{k} = \frac{\sum_{j}w_{kj}\pi_{j}/(\pi_{k}+\pi_{j})}{\sum_{j}(w_{jk})/(\pi_{k}+\pi_{j})} \]

\[ \pi_{k}=\frac{\pi_{k}}{(\prod_{j}^{n}\pi_{k})^{\frac{1}{n}}} \]

def bt_rank_2023(results, epochs=20):
    n = len(results)
    ratings = np.ones(n)
    for step in range(epochs):
        for k in range(n):
            numerator = 0
            denominator = 0
            for j in range(n):
                if j == k:
                    continue
                numerator += results[k][j] * ratings[j] / (ratings[k] + ratings[j])
                denominator += results[j][k] / (ratings[k] + ratings[j])

            ratings[k] = numerator / denominator
        ratings = ratings / np.exp(np.mean(np.log(ratings)))
    return ratings

print(bt_rank_2023(results, epochs=2))
print(bt_rank_2023(results, epochs=3))
print(bt_rank_2023(results, epochs=4))
print(bt_rank_2023(results, epochs=10))

[0.67734685 1.03368552 0.62447816 2.28708969]
[0.64996562 1.02624863 0.65621883 2.28459032]
[0.63867658 1.04248289 0.661457   2.27064091]
[0.63983459 1.04331496 0.65981021 2.2703762 ]

这里第 3-4 轮就接近收敛了，转成 Elo 得分

print(normalize2elo(np.log(bt_rank_2023(results, epochs=5))))

[1000.         1085.24912639 1005.59298743 1220.12534776]

再尝试一种近似计算方法，以前文提到的例子说明，假设 DeepSeekR1 vs GPT-4o 55 胜，45 负的情况下，针对 \( Y_{ij} \) 的极大似然估计直接就就可以得到 DeepSeekR1 和 GPT-4o 的 \( p_{ij}=0.55 \) 。

那么我们有： log(0.55/0.45)= \( r_{i}-r_{j} \) （对应的是 DeepSeekR1 和 GPT-4o 的得分差）

通过统计所有模型两两对战的结果，N 个模型会搜集到 C(N,2)=N(N-1)/2 对，那么就有 N(N-1)/2 组线性方程， 要解 N-1 个变量，方程数远大于变量数，那么这线性方程几乎是无解的，不过可以用最小二乘求近似。

这种做法相当于把似然函数多个乘积拆分成不同因子，然后单独对各个子因子进行最大似然估计，再找一个近似的可能最大化各个子因子的解，它和以上提到的优化全局似然是不同的，但这种拆分成多个子因子的思路和 Elo 算法有类似之处。

def local_bt_rank(results):
    n = len(results)
    rows = n * (n - 1) // 2
    X = np.zeros((rows, n))
    y = np.zeros(rows)
    row = 0
    for i in range(n):
        for j in range(i + 1, n):
            if (results[i][j] + results[j][i]) == 0 or results[i][j] == 0:
                continue
            p_ij = results[i][j] / (results[i][j] + results[j][i])
            log_odds = np.log(p_ij / (1 - p_ij))
            X[row, i] = 1
            X[row, j] = -1
            y[row] = log_odds
            row += 1
    # solve the least squares problem
    ratings = np.linalg.lstsq(X, y, rcond=None)[0]
    return ratings

res = local_bt_rank(results)
print(res)
print(normalize2elo(res))

[-0.47073006  0.03299569 -0.37956928  0.81730365]
[1000.         1087.50612634 1015.83624921 1223.75437381]

这也得到了 D>B>C>A 的排序，并且结果和上一节的迭代算法结果 [1000., 1085.25, 1005.6, 1220.13] 也较为接近。

上一节的算法是把每对选手的胜率统计出来，联立多个 log(p/(1-p))=ri-rj 方程，这种情况下如果 p 估计出来很不精确，比如两个选手比赛很少，就会增大误差。

然而，由于每场比赛都遵循 \( p_{ij}=\sigma (r_{i}-r_{j}) \) ，并且可以观察到胜或者负的标签（对应 1 和 0），那么直接将其看作二分类问题，用封装好的 logistic 回归的计算方法（比如 lbfgs）即可：

from sklearn.linear_model import LogisticRegression

def logistic_mle_rank(results, epoches=10):
    n = results.shape[0]  # 模型数量

    X, Y, sample_weights = [], [], []

    for i in range(n):
        for j in range(n):
            row = np.zeros(n)
            row[i], row[j] = 1, -1
            X.append(row)
            Y.append(1.0)
            sample_weights.append(results[i][j])

            row = np.zeros(n)
            row[i], row[j] = -1, 1
            X.append(row)
            Y.append(0.0)
            sample_weights.append(results[i][j])

    X = np.array(X)
    Y = np.array(Y)
    sample_weights = np.array(sample_weights)

    # 训练逻辑回归模型
    lr = LogisticRegression(
        fit_intercept=False, penalty='none', tol=1e-6,
        max_iter=epoches, solver="lbfgs"
    )
    lr.fit(X, Y, sample_weight=sample_weights)

    return lr.coef_[0]

print(logistic_mle_rank(results, 8))
print(normalize2elo(logistic_mle_rank(results, 8)))

[-0.44654166  0.04240325 -0.41580551  0.81994392]
[1000.         1084.93843084 1005.33941574 1220.0110788 ]

chatbot-arena-leaderboard 所采用的是这种计算方式，参考 Chatbot Arena Leaderboard Calculation (Bradley-Terry model) - Colab

有了先验分布，我们实际可以用 MonteCarlo 方式生成样本来模拟出对战情况，然后再用以上提到的算法去估计，对比一下高斯分布和 logistic 分布的情况：

def monte_carlo_bt_rank(dist):
    dev = 400 / np.log(10)
    if dist == np.random.logistic:
        dev *= np.sqrt(3)/np.pi
    n = 4
    scores = np.array([1000, 1085, 1005., 1220])
    results = np.zeros((n, n))
    for _ in range(10000):
        i, j = np.random.choice(n, 2, replace=False)
        si, sj = dist(scores[i], dev), dist(scores[j], dev)
        if si > sj:
            results[i][j] += 1
        else:
            results[j][i] += 1
    return normalize2elo(logistic_mle_rank(results, 8))

注意对于 np.random.logistic(mu, scale), 它的第二个参数是 scale 它和标准差的关系是 \( \sigma=\frac{s\pi}{\sqrt{3}} \)

print(monte_carlo_bt_rank(np.random.normal))
print(monte_carlo_bt_rank(np.random.logistic))

[1000.         1103.80268873 1011.547248   1260.18306746]
[1000.         1105.27856832 1001.77999902 1271.48526012]

但这里似乎把分值作为高斯分布，用 Bradley-Terry 模型中方法估计出来的分数更接近原始值，为什么？