Liu et al. "Metric Learning from Relative Comparisons by Minimizing Squared Residual". ICDM 2012.

relative_metric.py

import numpy as np
import scipy.linalg as LA

def comparison_loss(metric, comparisons):
    loss = 0.
    for weight, xa, xb, xc, xd in comparisons:
        vab = xa - xb
        vcd = xc - xd
        dab = np.dot(vab.T, np.dot(metric, vab))
        dcd = np.dot(vcd.T, np.dot(metric, vcd))
        if dab > dcd:
            loss += weight * (np.sqrt(dab) - np.sqrt(dcd))**2 
    return loss

def regularization_loss(metric, prior_inv):
    return np.trace(np.dot(metric, prior_inv)) - np.log(LA.det(metric))

def total_loss(metric, prior_inv, comparisons):
    return comparison_loss(metric, comparisons) + regularization_loss(metric, prior_inv)

def gradient(metric, prior_inv, comparisons):
    dMetric = prior_inv - LA.inv(metric)
    for weight, xa, xb, xc, xd in comparisons:
        vab = xa - xb
        vcd = xc - xd
        dab = np.dot(vab.T, np.dot(metric, vab))
        dcd = np.dot(vcd.T, np.dot(metric, vcd))
        if dab <= dcd:
            continue # comparison already satisfied.
        dMetric += (1-np.sqrt(dcd/dab))*np.outer(vab, vab) + (1-np.sqrt(dab/dcd))*np.outer(vcd, vcd)
    return dMetric

def lsml(prior, comparisons, tol=1e-3, max_iter=1000):
    prior_inv = LA.inv(prior)
    metric = prior.copy()
    it = 0
    s_best = total_loss(metric, prior_inv, comparisons)
    print 'initial loss', s_best
    while True:
        if it > max_iter:
            break
        grad = gradient(metric, prior_inv, comparisons)
        grad_norm = LA.norm(grad)
        if grad_norm < tol:
            break
        print 'gradient norm', grad_norm
        M_best = None
        for step_size in np.logspace(-10, 0, 10):
            step_size /= grad_norm
            new_metric = metric - step_size * grad
            # spectral decomposition
            w, v = LA.eigh(new_metric)
            #if np.any(w < 0):
            #    print 'for step_size', step_size, 'some eigs neg'
            new_metric = np.dot(v, np.dot(np.diag(np.maximum(w, 0)), v.T))
            assert new_metric.ndim == 2
            cur_s = total_loss(new_metric, prior_inv, comparisons)
            if cur_s < s_best:
                l_best = step_size
                s_best = cur_s
                M_best = new_metric
        it += 1
        print 'iter', it, 'cost', s_best, 'best step', l_best * grad_norm
        if M_best is None:
            break
        metric = M_best
    return metric