lr.py

import numpy as np
from scipy.optimize import fmin_l_bfgs_b
import time


def check_isolated_nodes(nvertices, edges):
    d = {}
    for i, j in edges:
        if i in d:
            d[i] += 1
        else:
            d[i] = 1
        if j in d:
            d[j] += 1
        else:
            d[j] = 1
    for x in xrange(nvertices):
        if x not in d:
            return True
    return False


def solve(nvertices, edges, T=10, r=None, gamma=10, eta=0.25, sigma0=10., sigmamax=1e6, verbose=False, pgtol=1e-2, factr=10., smart_init=True):


    if nvertices <= 0:
        raise ValueError("nvertices > 0 required")

    for i, j in edges:
        if i < 0 or i >= nvertices:
            raise ValueError("invalid vertex {0}".format(i))
        if j < 0 or j >= nvertices:
            raise ValueError("invalid vertex {0}".format(j))

    if check_isolated_nodes(nvertices, edges):
        raise ValueError("Isolated edge found")

    if r is None:
        m = len(edges) + 1
        for r in xrange(1, m + 1):
            if r * (r + 1) / 2. > m:
                r = r - 1
                break
    else:
        r = int(r)
        if r <= 0:
            raise ValueError("invalid value for r")
    if r > nvertices:
        r = nvertices

    if verbose:
        print "nvertices={0}, nedges={1}, r={2}".format(nvertices, len(edges), r)

    def lagrangian(R, y, sigma):
        RR_T = np.dot(R, R.T)

        # term1 = Dot(-11^T, RR^T)
        term1 = -RR_T.sum()

        # term2 = \sum_{i=1}^{m} y_i (Dot(A_i, RR^T) - b_i)
        # term3 = \sum_{i=1}^{m} (Dot(A_i, RR^T) - b_i)^2
        c0 = np.trace(RR_T) - 1.
        term2 = y[0] * c0
        term3 = c0 * c0

        for idx, (i, j) in enumerate(edges):
            cidx = 2. * RR_T[i, j]
            term2 += y[idx + 1] * cidx
            term3 += cidx * cidx

        return term1 - term2 + sigma * 0.5 * term3

    def grad(R, y, sigma):
        RR_T = np.dot(R, R.T)

        # Stilde = C - \sum_{i=1}^m ytilde_i A_i
        Stilde = -np.ones((nvertices, nvertices))
        ytilde_0 = y[0] - sigma * (np.trace(RR_T) - 1.)
        Stilde -= ytilde_0 * np.eye(nvertices)

        for idx, (i, j) in enumerate(edges):
            ytilde_idx = y[idx + 1] - 2. * sigma * RR_T[i, j]
            Stilde[i, j] -= ytilde_idx
            Stilde[j, i] -= ytilde_idx

        return 2. * np.dot(Stilde, R)

    def calcerrvec(R):
        RR_T = np.dot(R, R.T)
        evec = np.zeros(len(edges) + 1)
        evec[0] = np.trace(RR_T) - 1.
        for idx, (i, j) in enumerate(edges):
            evec[idx + 1] = 2. * RR_T[i, j]
        return evec

    if smart_init:
        nm = nvertices * (len(edges) + 1)
        Rcur = np.zeros((nvertices, r))
        for i in xrange(nvertices):
            for j in xrange(r):
                if i == j:
                    Rcur[i, j] = np.sqrt(1./r) + np.sqrt(1./nm)
                else:
                    Rcur[i, j] = np.sqrt(1./nm)
        Rflat_cur = Rcur.flatten()
        y_cur = -np.ones(len(edges) + 1)
        y_cur[0] = -nvertices
    else:
        Rflat_cur = np.random.random(nvertices * r)
        y_cur = np.random.random(len(edges) + 1)

    errvec = calcerrvec(Rflat_cur.reshape((nvertices, r)))
    vcur = np.dot(errvec, errvec)
    sigmacur = sigma0
    start = time.time()

    t = 0
    while 1:

        def slow_fn(Rflat):
            return lagrangian(Rflat.reshape((nvertices, r)), y_cur, sigmacur)

        def slow_grad(Rflat):
            G = grad(Rflat.reshape((nvertices, r)), y_cur, sigmacur)
            return G.flatten()

        Rcur = Rflat_cur.reshape((nvertices, r))
        fobj = -np.dot(Rcur, Rcur.T).sum()

        if verbose:
            print "*** starting iteration {0}, fobj={1}, sigma={2}, starting ||g||={3}".format(
                t + 1,
                fobj,
                sigmacur,
                np.linalg.norm(slow_grad(Rflat_cur)))

        Rflat_cur, _, d = fmin_l_bfgs_b(
            slow_fn, Rflat_cur, slow_grad, pgtol=pgtol, factr=factr, disp=0)
        if d['warnflag'] != 0:
            print "WARNING: could not converge inner opt"

        Rcur = Rflat_cur.reshape((nvertices, r))
        errvec = calcerrvec(Rcur)
        v = np.dot(errvec, errvec)
        fobj = -np.dot(Rcur, Rcur.T).sum()

        if v < eta * vcur:
            y_cur -= sigmacur * errvec
            vcur = v
        else:
            if verbose:
                print "Reject!"
            sigmacur = min(gamma * sigmacur, sigmamax)

        if verbose:
            print "*** Finished iteration {0}, fobj={1}, ||grad||={2}, infeas={3}, errvec[0]={4}, time={5}, ||R||={6} ***".format(
                t + 1,
                fobj,
                np.linalg.norm(d['grad']),
                v,
                errvec[0],
                time.time() - start,
                np.linalg.norm(Rcur)
            )

        t += 1
        if t >= T:
            break

def main():

    def parse(name):
        import csv
        with open(name, 'r') as fp:
            reader = csv.reader(fp, delimiter=',', skipinitialspace=True)
            for a, b in reader:
                yield int(a), int(b)

    edges = list(parse('johnson8-4-4.csv'))
    nvertices = max(map(max, edges)) + 1

    solve(nvertices, edges, T=20, r=34, gamma=np.sqrt(10.), pgtol=1e-6, factr=10., verbose=True)


if __name__ == '__main__':
    main()