Implementation of a rank aggregation algorithm using the paper http://www2007.org/papers/paper286.pdf The solver cvxopt is used with the frontend picos.

rankaggregation.py

from cvxopt import matrix, sparse, spmatrix
import cvxopt.lapack
import picos


class MarkovChain():

    def __init__(self, columns, target):
        n = len(target)
        l = len(columns)

        H, m = compute_H(target, n=n)
        Ps = build_transition_matrices(columns, n=n)
        A = compute_A(Ps, n=n, l=l)
        Hs = compute_Hs(A=A, H=H, l=l, n=n, m=m)

        # we convert the constants
        Hs = map(lambda (Hi, Hi_str): picos.new_param(Hi_str, Hi), zip(Hs, ["H0", "H1", "H2", "H3"]))

        # we now start defining the problem
        prob = picos.Problem()
        gamma = prob.add_variable("gamma")
        lambdas = map(lambda (name, val): prob.add_variable(name, val),
                      zip(map("Lambda{0}".format, range(4)),
                          [1, l + m, 2, n + m]))

        map(lambda l: prob.add_constraint(l > 0), lambdas)

        V, (D, U) = build_V(Hs=Hs, lambdas=lambdas, l=l, m=m, n=n, gamma=gamma)

        prob.add_constraint(V >> 0)

        e2 = picos.new_param("e2", matrix([1, 1]))
        prob.set_objective("max", gamma - lambdas[2].T * e2 - 2 * lambdas[0])

        self.problem = prob
        self.n = n
        self.m = m
        self.l = l
        self.gamma = gamma
        self.lambdas = lambdas
        self.V = V
        self.D = D
        self.U = U
        self.Hs = Hs
        self.A = A
        self.H = H
        self.Ps = Ps

    def solve(self):
        print(self.problem)
        self.problem.solve()

        beta = compute_beta((self.Hs[0] + self.lambdas[0] * self.D).value, - self.U.T.value / 2)

        l, m, n = self.l, self.m, self.n
        alpha = beta[0:l]
        x = beta[l: l + n]
        t = beta[l + n: l + n + m]

        return alpha, x, t


def compute_beta(must_inverse, and_multiply):
    """
    This holds a tricky part: inverting H0 + lambda0*D and applying it to Ut
    Problem: it is semi-definite positive, and may not be invertible
    see http://abel.ee.ucla.edu/cvxopt/userguide/lapack.html for available operations
    `potri` neglects this lack, so maybe `getri` (or a related function) should be used
    """
    inverse = must_inverse.real()
    ipiv = matrix(0, (must_inverse.size[0], 1))
    beta = +and_multiply
    cvxopt.lapack.gesv(inverse, beta, ipiv)
    return beta


# Helpers for matrices construction
def ident(q):
    return spmatrix(1, range(q), range(q))


def null(q, r=None):
    if r is None:
        r = q
    return matrix(0, (q, r))


def build_V(Hs, lambdas, l, m, n, gamma):
    U = sum(map(lambda (l, h): l.T * h,
                zip(lambdas, Hs)[1:]))

    D = sparse([[ident(l + n), null(m, l + n)],
                [null(l + n, m), null(m, m)]])
    D = picos.new_param("D", D)

    # we now build the SDP constraint
    V11 = Hs[0] + lambdas[0] * D
    V12 = U.T / 2
    V21 = U / 2
    V22 = - gamma

    V = (V11 & V12) // (V21 & V22)

    # also return elements needed to compute  the answer
    return V, (D, U)


def build_transition_matrices(columns, n):
    return map(lambda c: build_transition_matrix(c, n), columns)


def build_transition_matrix(column, n):
    mat = matrix(0., (n, n))
    for i in range(n):
        for j in range(n):
            if column[j] > column[i] or i == j:
                mat[i, j] = 1. / (1 + n - column[i])
    return mat


def compute_H(column, n):
    h = matrix(0, (n - 1, n))
    for i, v in enumerate(column):
        if v > 1:
            h[v - 2, i] = -1
        if v < n:
            h[v - 1, i] = 1
    return h, n - 1


def compute_A(P_matrices, n, l):
    A = matrix(0., (l, n))
    for i in range(l):
        for j in range(n):
            A[i, j] = P_matrices[i][j, j]
    return A


def compute_Hs(A, H, l, n, m):

    Im = ident(m)
    Il = ident(l)
    In = ident(n)

    H0 = sparse([
        [null(l, l), - A.T, null(m, l)],
        [- A, null(n + m, n)],
        [null(l + n, m), Im]
    ])
    H1 = sparse([
        [- Il, null(m, l)],
        [null(l + m, n)],
        [null(l, m), -Im]
    ])
    H2 = sparse([
        [1] * l + [0] * (n + m),
        [0] * l + [1] * n + [0] * m
    ]).T
    H3 = sparse([
        [null(n + m, l)],
        [- In, H],
        [null(n, m), - Im]
    ])
    return H0, H1, H2, H3