tie-aware的检索指标

检索常用指标：P（precision）、R（recall）、F1、AP（average precision）、RR（reciprocal rank）、NDCG（normalized discounted cumulative gain）、ACG（average cumulative gain）、WAP（Weighted mean Average Precision）。
检索时采用的某种 metric（如 hashing 的 hamming 距离），对某个 query 产生的检索序列，可能会有同距离的点（如同一 hamming radius），同距离的点的集合称为一个 tie。tie 内各元素对于该 metric 等价，但它们不同的排序会影响最终指标的结果，例子见 [4, 7]。
对于这种情况，有两种评估策略：

• 随机选一个合法的排序情况来评估。简单，但不合理，非确定性，相当于引入了未控制的排序变量，见 [4]。
• 枚举所有合法的排序情况，分别评估，再取平均（macro-average）。合理，但暴力枚举复杂度是阶乘级。

为综合两者优势，可以设法直接导出第二种方法的平均，避免暴力枚举，见 [1]。这样导出的指标称为原指标的 tie-aware 版，加下标 ${}_T$ 表示。

Notations

• q。query 样本；
• $d(\cdot )$。所选 metric，值越小表明样本与 q 越相似；
• $V=<v_1,\dots ,v_n>$。用 d 对 q 的检索序列，有：$d(v_i)\le d(v_j),\forall i<j$。n 个检索样本根据 d 值相等割为 m 个 tie，第 i 个记为 $V_i$；
• $T=<t_0,t_1,\dots ,t_m>$。一个下标序列，$t_i$ 表示 $V_i$ 中最后一个元素的位置（$t_0$ 纯粹是为了表述方便而加的），于是 V i = { v t i − 1 + 1 , … , v t i } V_i=\{v_{t_{i-1}+1},\dots,v_{t_i}\} Vi​={vti−1​+1​,…,vti​​}，或用区间表示成 $(t_{i-1},t_i]$；
• $\text{rel}(v_i)$。$v_i$ 与 q 的（ground-truth 的）相似性。默认的情况是相似为 1，不相似为 0，如 P、R、F1、AP 中都是这种定义；
• $r_i,n_i$。分别表示 $V_i$ 中相似元素个数、总元素个数，$r_i={\sum }_{j=t_{i-1}+1}^{t_i}rel(v_j)$，$n_i=|V_i|$；
• $R_i$。$r_i$ 的前缀和，$R_i={\sum }_{j=1}^i{r}_i$。为了方便，定义 $R_0=0$；

P

考虑 P@k，可以兼容 P（即 P = P@n），后同。原始 P@k：$$P@k(V)=\frac{1}{k}\sum _{i=1}^{k}rel(v_i)$$ 忽略常系数之后就是求和。当 k 卡在第 c 个 tie 中，$k\in (t_{c-1},t_c]$，$V_c$ 之前的相似元素个数是 $R_{c-1}$，与顺序无关；之后的对 P@k 无贡献；$V_c$ 中平均每个位置有 $\frac{r_c}{n_c}$ 个相关的样本（从「是否相关」的角度共有 $C_{n_c}^{r_c}$ 种组合，其中第 j 个位置放相关元素的有 $C_1^1\cdot C_{n_c-1}^{r_c-1}$ 种），而处于前 k 的位置有 $k-t_{c-1}$ 个，所以 tie-aware 版 P@k(V) 是：$$P_T@k(V)=\frac{1}{k}[R_{c-1}+(k-t_{c-1})\frac{r_c}{n_c}]$$

R

recall 同 precision 几乎一样，只是系数不同：$$R_T@k(V)=\frac{1}{R_m}[R_{c-1}+(k-t_{c-1})\frac{r_c}{n_c}]$$ 其中 $R_m$ 就是 m 个 tie 中所有相关样本的总个数。

F1

F1 是 P 和 R 的调和平均：$$F1@k(V)=\frac{1}{\frac{1}{P@k(V)}+\frac{1}{R@k(V)}}=\frac{2}{k+R_m}\sum _{i=1}^{k}rel(v_i)$$ 分析、改法类似前面，变成：$$F{1}_T@k(V)=\frac{2}{k+R_m}[R_{c-1}+(k-t_{c-1})\frac{r_c}{n_c}]$$

AP

• （2021.11.10）有个问题：当 @k 时，分母在 k 所在 tie 的计算是不是也应该按概率处理一下？这个问题目前未解决，所以目前只有用 @ALL 是保证没问题的。
• （2020.11.13）由于采用 [1] 公式在 @k 时会与普通 mAP 差别很大，所以改成与 [9] 一致，分母用「检索序列前 k 个结果中的相关样本数」，见代码注释。
• 计算 mAP@k 时关于分母应该是「整个检索序列的相关样本数」，还是「检索序列前 k 个结果中的相关样本数」可见 [9] 中笔记。本文暂时以 [1] 的公式为准，即采用前者，而与 [9] 中不同。 不过，当计算 $mA{P}_T$@ALL 的时候，这种分歧不造成影响。

原始 AP@k 公式是：$$AP@k(V)=\frac{{\sum }_{i=1}^{k}P@i(V)\cdot rel(v_i)}{R_m}$$ 对于每个位置 i，假设其落在第 $f(i)$ 个 tie：$i\in (t_{f(i)-1},t_{f(i)}]$；$P_T@i(V)$ 的分母部分是 i；$rel(v_i)$ 有两种情况：

• $rel(v_i)=0$，贡献为 0，加了等于没加；
• $rel(v_i)=1$，频率是 $\frac{r_{f(i)}}{n_{f(i)}}$。此时$P_T@i(V)$ 的分子部分是 $R_{f(i)-1}+(i-t_{f(i)-1}-1)\frac{r_{f(i)}-1}{n_{f(i)}-1}+1$，中间那项的分析类似前面，其中 -1 是除了第 i 个位置，而第 i 个位置的贡献就是第三项 +1；

所以综合两种情况：$$A{P}_T@k(V)=\frac{1}{R_m}\sum _{i=1}^k\frac{1}{i}[R_{f(i)-1}+(i-t_{f(i)-1}-1)\frac{r_{f(i)}-1}{n_{f(i)}-1}+1]\cdot \frac{r_{f(i)}}{n_{f(i)}}$$ 需要注意 $n_{f(i)}=1$ 的情况，即整个 tie 只有一个元素，此时 $\frac{r_{f(i)}-1}{n_{f(i)}-1}$ 无意义，需要将 $(i-t_{f(i)-1}-1)\frac{r_{f(i)}-1}{n_{f(i)}-1}$ 整项置零，而此元素是否有贡献由后面的 $r_{f(i)}$ 决定。

NDCG

介绍见 [5]，此时 $\text{rel}(\cdot )$ 可以是多值的，表示相似的程度。原始公式：$$DCG@k(V)=\sum _{i=1}^{k}G[rel(v_i)]\cdot D(i)$$ 其中 $G(\cdot )$ 是 gain 函数，一般取 $G(x)=2^x-1$；$D(\cdot )$ 是 discount 函数，取 $D(x)=\frac{1}{\log (1+x)}$，此时可重写成：$$DCG@k(V)=\sum _{i=1}^k\frac{2^{rel(v_i)}-1}{\log (1+i)}$$ 分子这个形式可以增大不同 $rel(v_i)$ 取值之间的差异，而当 r e l ( ⋅ ) ∈ { 0 , 1 } rel(\cdot)\in \{0,1\} rel(⋅)∈{0,1} 时与直接取 $G(x)=x$ 等价。总体就是与相似程度正相关，与位置负相关。假设对于 q，最优的检索序列是 $I$，则：$$NDCG@k(V)=\frac{DCG@k(V)}{DCG@k(I)}$$ 给定 q 时 $I$、$DCG@k(I)$ 都是常量，可以另外算。考虑 DCG@k(V) 的计算。每一个 tie $V_c$ 分别算，有两种情况：

• $V_c$ 包含在 [1,k] 之内：$0\le t_{c-1}<t_c\le k$。对于 $V_c$ 中的任一样本 v，gain 是 $G[rel(v)]$，可以遍历所有 $V_c$ 的位置，且处于每个位置的频率都是 $\frac{1}{n_c}$，所以 v 对 DCG@k 的贡献是：$$\frac{G[rel(v)]}{n_c}\sum _{i=t_{c-1}+1}^{t_c}D(i)$$ 此处 $\frac{G[rel(v)]}{n_c}$ 可以看作 v 的 gain 的均值，故 $V_c$ 中所有元素对 DCG@k 的总贡献是：$$\sum _{v\in V_c}\frac{G[rel(v)]}{n_c}\sum _{i=t_{c-1}+1}^{t_c}D(i)$$
• k 在 $V_c$ 内：$k\in (t_{c-1},t_c]$。此时 v 的 gain 均值不变（还是考虑 $V_c$ 所有可能的排列），只是位置截到前 k 就停，所以贡献是：$$\sum _{v\in V_c}\frac{G[rel(v)]}{n_c}\sum _{i=t_{c-1}+1}^{k}D(i)$$

综上，有： D C G T @ k ( V ) = ∑ c = 1 m ′ ∑ v ∈ V c G [ r e l ( v ) ] n c ∑ i = t c − 1 + 1 min ⁡ { t c , k } D ( i ) DCG_T@k(V)=\sum_{c=1}^{m'}\sum_{v\in V_c}\frac{G[rel(v)]}{n_c}\sum_{i=t_{c-1}+1}^{\min\{t_c,k\}}D(i) DCGT​@k(V)=c=1∑m′​v∈Vc​∑​nc​G[rel(v)]​i=tc−1​+1∑min{tc​,k}​D(i) 其中 m’ 使得 $k\in (t_{m^{\prime }-1},t_{m^{\prime }}]$。

Codes

• （2021.10.24 fix bug）NDCG_tie 更新，之前在处理 k 所在 tie 的 gain 均值时有错，现而改正。不过此 bug 只在 $k<ALL$ 时显现，不影响 $NDCG@ALL$ 的结果。
• 与 [10] 提供的程序对拍过，应该是正确的，样例还是 [9] 那个。注意 [10] 的代码是假设 hash 码是 0/1 的而不是 -1/1 的。因 [10] 的程序是没有 @k 的，所以测试通过是指在 $mA{P}_T@ALL$ 和 $NDC{G}_T@ALL$ 下与 [10] 结果相同。

# import numpy as np

def hamming(A, B=None):
    """A, B: [None, bit]
    elements in {-1, 1}
    """
    if B is None: B = A
    bit = A.shape[1]
    return (bit - A.dot(B.T)) // 2


def sim_mat(label, label_2=None, sparse=False):
    if label_2 is None:
        label_2 = label
    if sparse:
        S = label[:, np.newaxis] == label2[np.newaxis, :]
    else:
        S = np.dot(label, label_2.T) > 0
    return S.astype(label.dtype)


def P_R_F1_tie(qH, rH, qL, rL, k=-1, sparse=False):
    """tie-aware precision@k, recall@k, F1@k
    qH, rH: [n, bit] & [m, bit], hash code of query & database samples
    qL, rL: [n, #class] & [m, #class], label of query & database samples
        if sparse, shapes are [n] & [m]
    k: position threshold of `P@k`
    sparse: whether label is sparse class ID or one-hot vector
    """
    m = rH.shape[0]
    if (k < 0) or (k > m):
        k = m
    D = hamming(qH, rH)
    Rnk = np.argsort(D)
    S = sim_mat(qL, rL, sparse)
    # 找 k 所在 tie: t_{c-1} < k <= t_c
    D_sort = np.vstack([d[r] for d, r in zip(D, Rnk)])
    D_at_k = D_sort[:, k-1:k]  # `-1` 因为下标是 0-base
    mask_tie = np.equal(D, D_at_k).astype(np.int)
    # r_c, n_c
    nc = mask_tie.sum(1)
    rc = (mask_tie * S).sum(1)
    # find t_{c-1}
    pos = np.arange(m)[np.newaxis, :]  # [1, m]
    tie_pos = np.where(np.equal(D_sort, D_at_k), pos, np.inf)
    tc_1 = np.min(tie_pos, 1) # - 1 + 1  # `+1` to shift 0-base to 1-base
    # R_{c-1}
    mask_pre = (D < D_at_k).astype(np.int)
    Rc_1 = (mask_pre * S).sum(1)
    # P@k, R@k, F1@k
    _common = Rc_1 + (k - tc_1) * rc / nc  # [n]
    Rm = S.sum(1)
    P_at_k = _common / k
    R_at_k = _common / Rm
    F1_at_k = 2 * _common / (k + Rm)
    return P_at_k.mean(), R_at_k.mean(), F1_at_k.mean()


def mAP_tie(qH, rH, qL, rL, k=-1, sparse=False):
    """tie-aware mAP
    ref:
    - https://blog.csdn.net/HackerTom/article/details/107458334
    - https://github.com/kunhe/TALR/blob/master/%2Beval/tieAP.m
    """
    n, m = qH.shape[0], rH.shape[0]
    if (k < 0) or (k > m):
        k = m
    D = hamming(qH, rH)
    Rnk = np.argsort(D)
    S = sim_mat(qL, rL, sparse)
    AP = 0
    pos = np.arange(m)  # 0-base
    # t_fi_1[k]: t_{f(k) - 1}
    t_fi_1 = np.zeros([m])
    # r_fi[k]: #relevant samples in the tie where k lies in
    r_fi = np.zeros([m])
    # n_fi[k]: #samples in the tie where k lies in
    n_fi = np.zeros([m])
    # R_fi_1[k]: prefix sum of r_fi (exclude r_fi[k])
    R_fi_1 = np.zeros([m])
    for d, s, rnk in zip(D, S, Rnk):
        # Rm = s.sum()  # #rel in all
        s_sort = s[rnk]
        Rm = s_sort[:k].sum()  # #rel in top-k
        if 0 == Rm:
            continue
        d_unique = np.unique(d)  # ascending
        d_sort = d[rnk]
        # s_sort = s[rnk]
        _R_fi_1 = 0  # R_{f(i) - 1}
        for _d in d_unique:
            tie_idx = (d_sort == _d)
            t_fi_1[tie_idx] = pos[tie_idx].min() # - 1 + 1  # `+1` to shift 0-base to 1-base
            _r_fi = s_sort[tie_idx].sum()
            r_fi[tie_idx] = _r_fi
            n_fi[tie_idx] = tie_idx.astype(np.int).sum()
            R_fi_1[tie_idx] = _R_fi_1  # exclude `_r_fi`
            _R_fi_1 += _r_fi

        # deal with invalid terms
        n_fi_1, r_fi_1 = n_fi - 1, r_fi - 1
        idx_invalid = (n_fi_1 == 0)
        n_fi_1[idx_invalid] = 1
        r_fi_1[idx_invalid] = 0
        # in computing (i - t_{f(i)-1} - 1),
        # the lastest `-1` is megered: pos = i - 1
        kernel = (R_fi_1 + (pos - t_fi_1) * r_fi_1 / n_fi_1 + 1) * r_fi / n_fi / (pos + 1)
        AP += kernel[:k].sum() / Rm

    return AP / n


def NDCG_tie(qH, rH, qL, rL, k=-1, sparse=False):
    """tie-aware NDCG
    ref:
    - https://blog.csdn.net/HackerTom/article/details/107458334
    - https://github.com/kunhe/TALR/blob/master/%2Beval/tieNDCG.m
    """
    n, m = Dist.shape
    if (k < 0) or (k > m):
        k = m
    G = 2 ** Rel - 1
    pos = np.arange(m) + 1  # 1-base
    D_inv = 1 / np.log2(1 + pos)
    Rank = np.argsort(Dist)

    _NDCG = 0
    for g, d, rnk in zip(G, Dist, Rank):
        dcg_best = (np.sort(g)[::-1][:k] * D_inv[:k]).sum()
        if 0 == dcg_best:
            continue
        d_unique = np.unique(d)  # ascending
        d_sort = d[rnk]
        g_sort = g[rnk]
        dcg = 0
        for _d in d_unique:
            tie_idx = (d_sort == _d)
            tie_pos = pos[tie_idx]
            tc_1 = tie_pos[0] - 1  # i.e. tie_pos.min() - 1
            if tc_1 >= k:  # k <= t_{c-1} < t_c, out of range
                break  # continue
            n_c = tie_idx.astype(np.int).sum()
            tc = tie_pos[-1]
            sum_d = (1 / np.log2(1 + np.arange(tc_1 + 1, min(k, tc) + 1))).sum()
            tie_avg_gain = g_sort[tie_idx].sum() / n_c
            dcg += tie_avg_gain * sum_d

        _NDCG += dcg / dcg_best

    return _NDCG / n

multiple $R$

• 升级版，支持单个或多个 position thresholds，传 int 或 int tuple/list。
• multi_nDCG_tie 与 2021.10.24 debug 之后的单 k 版一致，在 $@ALL$ 时与 [10] 也一致。

def multi_mAP_tie(Dist, Sim, k=-1):
    """支持单 k、多 k
    多个 k 时传 int tuple/list
    """
    if isinstance(k, int):
        k = [k]
    n, m = Dist.shape
    for kid in range(len(k)):
        if (k[kid] < 0) or (k[kid] > m):
            k[kid] = m
    k = sorted(k)  # ascending
    assert k[0] != 0, "`@0` is meaningless and disallowed for efficiency"
    Rnk = np.argsort(Dist, axis=-1)
    # AP = 0
    AP = np.zeros([len(k)], dtype=np.float32)
    pos = np.arange(m)  # 0-base
    # t_fi_1[k]: t_{f(k) - 1}
    t_fi_1 = np.zeros([m])
    # r_fi[k]: #relevant samples in the tie where k lies in
    r_fi = np.zeros([m])
    # n_fi[k]: #samples in the tie where k lies in
    n_fi = np.zeros([m])
    # R_fi_1[k]: prefix sum of r_fi (exclude r_fi[k])
    R_fi_1 = np.zeros([m])
    for d, s, rnk in zip(Dist, Sim, Rnk):
        # Rm = s.sum()  # #rel in all
        s_sort = s[rnk]
        # Rm = s_sort[:k].sum()  # #rel in top-k
        # if 0 == Rm:
        #     continue
        Rm_list = s_sort.cumsum()
        if 0 == Rm_list[-1]:  # = Rm.max() = s_sort.sum()
            continue
        d_unique = np.unique(d)  # ascending
        d_sort = d[rnk]
        # s_sort = s[rnk]
        _R_fi_1 = 0  # R_{f(i) - 1}
        for _d in d_unique:
            tie_idx = (d_sort == _d)
            t_fi_1[tie_idx] = pos[tie_idx].min() # - 1 + 1  # `+1` to shift 0-base to 1-base
            _r_fi = s_sort[tie_idx].sum()
            r_fi[tie_idx] = _r_fi
            n_fi[tie_idx] = tie_idx.astype(np.int).sum()
            R_fi_1[tie_idx] = _R_fi_1  # exclude `_r_fi`
            _R_fi_1 += _r_fi

        # deal with invalid terms
        n_fi_1, r_fi_1 = n_fi - 1, r_fi - 1
        idx_invalid = (n_fi_1 == 0)
        n_fi_1[idx_invalid] = 1
        r_fi_1[idx_invalid] = 0
        # in computing (i - t_{f(i)-1} - 1),
        # the lastest `-1` is megered: pos = i - 1
        kernel = (R_fi_1 + (pos - t_fi_1) * r_fi_1 / n_fi_1 + 1) * r_fi / n_fi / (pos + 1)
        # AP += kernel[:k].sum() / Rm
        kernel_cumsum = np.cumsum(kernel)
        for kid, _k in enumerate(k):
            if Rm_list[_k - 1]:
                AP[kid] += kernel_cumsum[_k - 1] / Rm_list[_k - 1]

    _mAP = AP / n
    if 1 == _mAP.shape[0]:
        _mAP = _mAP[0]
    return _mAP


def multi_nDCG_tie(Dist, Rel, k=-1):
    if isinstance(k, int):
        k = [k]
    n, m = Dist.shape
    for kid in range(len(k)):
        if (k[kid] < 0) or (k[kid] > m):
            k[kid] = m
    k = sorted(k)  # ascending
    assert k[0] != 0, "`@0` is meaningless and disallowed for efficiency"
    G = 2 ** Rel - 1
    pos = np.arange(m) + 1  # 1-base
    D_inv = 1 / np.log2(1 + pos)
    Rank = np.argsort(Dist, axis=-1)

    _nDCG = np.zeros([len(k)], dtype=np.float32)
    for g, dist, rnk in zip(G, Dist, Rank):
        g_desc = np.sort(g)[::-1]
        if 0 == g_desc[0]:  # biggist DCG
            continue
        dcg_best_list = (g_desc * D_inv).cumsum()
        dist_unique = np.unique(dist)  # ascending
        dist_sort = dist[rnk]
        g_sort = g[rnk]

        _start_kid = 0
        dcg_list = np.zeros([len(k)], dtype=np.float32)
        for _d in dist_unique:
            tie_idx = (dist_sort == _d)
            tie_pos = pos[tie_idx]
            tc_1 = tie_pos[0] - 1  # i.e. tie_pos.min() - 1
            tc = tie_pos[-1]
            while _start_kid < len(k):
                if k[_start_kid] > tc_1:
                    break
                else:
                    _start_kid += 1
            if _start_kid >= len(k):
                break
            n_c = tie_idx.sum()
            tie_avg_gain = g_sort[tie_idx].sum() / n_c  # != g_sort[tie_idx].mean()
            tie_sum_d_list = (1 / np.log2(1 + np.arange(tc_1 + 1, tc + 1))).cumsum()
            for kid_offset, _k in enumerate(k[_start_kid:]):
                tie_sum_d = tie_sum_d_list[min(_k, tc) - tc_1 - 1]
                tie_dcg = tie_avg_gain * tie_sum_d
                dcg_list[kid_offset + _start_kid] += tie_dcg

        for kid, _k in enumerate(k):
            _ndcg = dcg_list[kid] / dcg_best_list[_k - 1]  # 0-base
            _nDCG[kid] += _ndcg

    _nDCG /= n
    if 1 == _nDCG.shape[0]:
        _nDCG = _nDCG[0]
    return _nDCG


if __name__ == "__main__":
    print("与原单 k 版对拍。结论：一致")
    N, M = 5, 20
    qB = np.random.randint(0, 2, size=(N, 3)) * 2 - 1
    rB = np.random.randint(0, 2, size=(M, 3)) * 2 - 1
    qL = np.random.randint(0, 2, size=(N, 7))
    rL = np.random.randint(0, 2, size=(M, 7))
    D = hamming(qB, rB)
    Rel = qL.dot(rL.T)
    S = (Rel > 0).astype(np.int32)
    k_list = [1] + list(range(0, M + 1, 5)[1:])
    print("k_list:", k_list)
    print("-> 对拍 mAP_tie")
    map1 = [mAP_tie(D, S, k=_k) for _k in k_list]
    map2 = multi_mAP_tie(D, S, k_list)
    print("mAP 1:", map1)
    print("mAP 2:", map2)
    print("-> 对拍 nDCG_tie")
    ndcg1 = [nDCG_tie(D, Rel, k=_k) for _k in k_list]
    ndcg2 = multi_nDCG_tie(D, Rel, k_list)
    print("nDCG 1:", ndcg1)
    print("nDCG 2:", ndcg2)

References

1. Computing Information Retrieval Performance Measures Efficiently in the Presence of Tied Scores
2. C# implementations
3. Hashing as Tie-Aware Learning to Rank
4. Tie-Breaking Bias: Effect of an Uncontrolled Parameter on Information Retrieval Evaluation
5. Discounted cumulative gain
6. NDCG及实现
7. numpy和pytorch的argsort结果不同
8. kunhe/TALR
9. mAP（@R）计算代码
10. kunhe/TALR/+eval