Use a faster Bradley-Terry implementation

Author: dustalov

python/evalica/__init__.py

@@ -57,7 +57,7 @@ def _make_matrix(
         win_matrix: npt.NDArray[np.float64],
         tie_matrix: npt.NDArray[np.float64],
         win_weight: float = 1.,
-        tie_weight: float = 0.,
+        tie_weight: float = .5,
         nan: float = 0.0,
 ) -> npt.NDArray[np.float64]:
     with np.errstate(all="ignore"):
@@ -391,6 +391,11 @@ def bradley_terry(
         The Method of Paired Comparisons. Biometrika. 39, 324–345 (1952).
         <https://doi.org/10.2307/2334029>.
 
+    Quote:
+        Newman, M.E.J.: Efficient Computation of Rankings from Pairwise Comparisons.
+        Journal of Machine Learning Research. 24, 1&ndash;25 (2023).
+        <https://www.jmlr.org/papers/v24/22-1086.html>.
+
     Args:
         xs: The left-hand side elements.
         ys: The right-hand side elements.

python/evalica/naive.py

@@ -13,7 +13,7 @@
 
 
 def pairwise_scores(scores: npt.NDArray[np.number[Any]]) -> npt.NDArray[np.float64]:
-    if not scores.shape[0]:
+    if not scores.size:
         return np.zeros((0, 0))
 
     return np.nan_to_num(scores[:, np.newaxis] / (scores + scores[:, np.newaxis]))
@@ -63,31 +63,30 @@ def bradley_terry(
         tolerance: float = 1e-6,
         limit: int = 100,
 ) -> tuple[npt.NDArray[np.float64], int]:
-    with np.errstate(all="ignore"):
-        totals = matrix.T + matrix
-
-    active = totals > 0
+    scores = np.ones(matrix.shape[0])
 
-    wins = matrix.sum(axis=1)
+    converged, iterations = False, 0
 
-    normalized = np.zeros_like(matrix, dtype=float)
+    if not matrix.size:
+        return scores, iterations
 
-    scores = np.ones(matrix.shape[0])
     scores_new = scores.copy()
 
-    converged, iterations = False, 0
-
     while not converged and iterations < limit:
         iterations += 1
 
         with np.errstate(all="ignore"):
-            sums = np.add.outer(scores, scores)
+            for i in range(matrix.shape[0]):
+                sums = scores_new[i] + scores_new
+
+                numerator = np.sum(matrix[i] * scores_new / sums)
+                denominator = np.sum(matrix[:, i] / sums)
 
-            normalized[active] = totals[active] / sums[active]
+                scores_new[i] = numerator / denominator
 
-            scores_new[:] = wins
-            scores_new /= normalized.sum(axis=0)
-            scores_new /= scores_new.sum()
+            geometric_mean = np.exp(np.mean(np.log(scores_new)))
+
+            scores_new /= geometric_mean
 
         scores_new[:] = np.nan_to_num(scores_new, nan=tolerance)
 
@@ -105,14 +104,18 @@ def newman(
         tolerance: float = 1e-6,
         limit: int = 100,
 ) -> tuple[npt.NDArray[np.float64], float, int]:
-    win_tie_half = win_matrix + tie_matrix / 2
+    win_tie_half = np.nan_to_num(win_matrix + tie_matrix / 2, nan=tolerance)
 
     scores = np.ones(win_matrix.shape[0])
-    scores_new = scores.copy()
-    v_new = v
 
     converged, iterations = False, 0
 
+    if not win_matrix.size and not tie_matrix.size:
+        return scores, v, iterations
+
+    scores_new = scores.copy()
+    v_new = v
+
     while not converged and iterations < limit:
         iterations += 1
 
@@ -199,7 +202,7 @@ def eigen(
         tolerance: float = 1e-6,
         limit: int = 100,
 ) -> tuple[npt.NDArray[np.float64], int]:
-    if not matrix.shape[0]:
+    if not matrix.size:
         return np.zeros(0, dtype=np.float64), 0
 
     n = matrix.shape[0]
@@ -227,7 +230,7 @@ def pagerank_matrix(
         matrix: npt.NDArray[np.float64],
         damping: float,
 ) -> npt.NDArray[np.float64]:
-    if not matrix.shape[0]:
+    if not matrix.size:
         return np.zeros(0, dtype=np.float64)
 
     p = 1. / int(matrix.shape[0])

python/evalica/test_evalica.py

@@ -188,7 +188,7 @@ def test_bradley_terry(comparison: Comparison, win_weight: float, tie_weight: fl
         assert len(result.scores) == len(set(xs) | set(ys))
         assert np.isfinite(result.scores).all()
         assert result.scores.is_monotonic_decreasing
-        assert result.iterations > 0
+        assert not xs or result.iterations > 0
         assert result.limit > 0
 
     assert_series_equal(result_pyo3.scores, result_naive.scores, rtol=1e-3, check_like=True)
@@ -216,8 +216,8 @@ def test_newman(comparison: Comparison, v_init: float) -> None:
         assert len(result.scores) == len(set(xs) | set(ys))
         assert np.isfinite(result.scores).all()
         assert result.scores.is_monotonic_decreasing
-        assert np.isfinite(result.v)
-        assert result.iterations > 0
+        assert not xs or np.isfinite(result.v)
+        assert not xs or result.iterations > 0
         assert result.limit > 0
 
         if np.isfinite(v_init):
@@ -226,7 +226,8 @@ def test_newman(comparison: Comparison, v_init: float) -> None:
             assert result.v_init is v_init
 
     assert_series_equal(result_pyo3.scores, result_naive.scores, check_like=True)
-    assert result_pyo3.v == pytest.approx(result_naive.v)
+
+    assert not np.isfinite(v_init) or result_pyo3.v == pytest.approx(result_naive.v)
 
 
 @given(
@@ -445,9 +446,12 @@ def test_bradley_terry_dataset(comparison: Comparison, comparison_golden: pd.Ser
     result_pyo3 = evalica.bradley_terry(xs, ys, winners, weights=weights, solver="pyo3")
     result_naive = evalica.bradley_terry(xs, ys, winners, weights=weights, solver="naive")
 
-    assert_series_equal(result_naive.scores, comparison_golden, rtol=1e-4, check_like=True)
-    assert_series_equal(result_pyo3.scores, comparison_golden, rtol=1e-4, check_like=True)
-    assert_series_equal(result_pyo3.scores, result_naive.scores, check_like=True)
+    scores_pyo3 = result_pyo3.scores / result_pyo3.scores.sum()
+    scores_naive = result_naive.scores / result_naive.scores.sum()
+
+    assert_series_equal(scores_naive, comparison_golden, rtol=1e-4, check_like=True)
+    assert_series_equal(scores_pyo3, comparison_golden, rtol=1e-4, check_like=True)
+    assert_series_equal(scores_pyo3, scores_naive, check_like=True)
 
 
 @pytest.mark.parametrize(("algorithm", "dataset"), [

src/bradley_terry.rs

@@ -1,11 +1,11 @@
 use std::ops::{AddAssign, DivAssign};
 
-use ndarray::{Array1, Array2, ArrayView2, Axis, ErrorKind, ScalarOperand, ShapeError};
-use num_traits::Float;
+use ndarray::{Array1, ArrayView2, Axis, ErrorKind, ScalarOperand, ShapeError};
+use num_traits::{Float, FromPrimitive};
 
-use crate::utils::{nan_to_num, one_nan_to_num};
+use crate::utils::{nan_to_num, one_nan_to_num, win_plus_tie_matrix};
 
-pub fn bradley_terry<A: Float + ScalarOperand + AddAssign + DivAssign>(
+pub fn bradley_terry<A: Float + FromPrimitive + ScalarOperand + AddAssign + DivAssign>(
     matrix: &ArrayView2<A>,
     tolerance: A,
     limit: usize,
@@ -14,41 +14,35 @@ pub fn bradley_terry<A: Float + ScalarOperand + AddAssign + DivAssign>(
         return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
     }
 
-    let totals = &matrix.t().clone() + matrix;
-
-    let active = totals
-        .indexed_iter()
-        .filter(|((_, _), &total)| total > A::zero())
-        .collect::<Vec<((usize, usize), &A)>>();
-
-    let wins = matrix.sum_axis(Axis(1));
-
-    let mut normalized = Array2::zeros(matrix.raw_dim());
-
     let mut scores = Array1::ones(matrix.shape()[0]);
 
     let mut converged = false;
     let mut iterations = 0;
 
+    if matrix.is_empty() {
+        return Ok((scores, iterations));
+    }
+
     while !converged && iterations < limit {
         iterations += 1;
 
-        for ((i, j), &v) in active.iter() {
-            let i = *i;
-            let j = *j;
+        let mut scores_new = scores.clone();
 
-            normalized[[i, j]] = v / (scores[i] + scores[j]);
-        }
+        for i in 0..matrix.nrows() {
+            let mut numerator = A::zero();
+            let mut denominator = A::zero();
 
-        let mut scores_new = &wins / &normalized.sum_axis(Axis(0));
+            for j in 0..matrix.ncols() {
+                let sum_scores = scores_new[i] + scores_new[j];
+                numerator += matrix[[i, j]] * scores_new[j] / sum_scores;
+                denominator += matrix[[j, i]] / sum_scores;
+            }
 
-        // Otherwise the result is different from what is computed by NumPy
-        let mut scores_new_sum = A::zero();
-        for &score in scores_new.iter() {
-            scores_new_sum += score;
+            scores_new[i] = numerator / denominator;
         }
 
-        scores_new /= scores_new_sum;
+        let geometric_mean = scores_new.mapv(|x| x.ln()).mean().unwrap().exp();
+        scores_new /= geometric_mean;
 
         nan_to_num(&mut scores_new, tolerance);
 
@@ -72,15 +66,20 @@ pub fn newman(
         return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
     }
 
-    let win_tie_half = win_matrix + &(tie_matrix / 2.0);
-
     let mut scores = Array1::<f64>::ones(win_matrix.shape()[0]);
     let mut v = v_init;
-    let mut v_new = v;
 
     let mut converged = false;
     let mut iterations = 0;
 
+    if win_matrix.is_empty() && tie_matrix.is_empty() {
+        return Ok((scores, v, iterations));
+    }
+
+    let win_tie_half = win_plus_tie_matrix(&win_matrix, &tie_matrix, 1.0, 0.5, tolerance);
+
+    let mut v_new = v;
+
     while !converged && iterations < limit {
         iterations += 1;
 
@@ -156,7 +155,9 @@ mod tests {
         assert_eq!(actual.len(), matrix.shape()[0]);
         assert_ne!(iterations, 0);
 
-        for (left, right) in actual.iter().zip(expected.iter()) {
+        let actual_normalized = actual.clone() / actual.sum();
+
+        for (left, right) in actual_normalized.iter().zip(expected.iter()) {
             assert_abs_diff_eq!(left, right, epsilon = tolerance * 1e1);
         }
     }
@@ -191,7 +192,7 @@ mod tests {
             tolerance,
             100,
         )
-            .unwrap();
+        .unwrap();
 
         assert_eq!(actual.len(), win_matrix.shape()[0]);
         assert_eq!(actual.len(), tie_matrix.shape()[0]);

src/elo.rs

@@ -85,7 +85,7 @@ mod tests {
             1.0,
             0.5,
         )
-            .unwrap();
+        .unwrap();
 
         for (a, b) in actual.iter().zip(expected.iter()) {
             assert!((a - b).abs() < 1e-0, "a = {}, b = {}", a, b);