PGD QPLIB Benchmarking

import numpy as np
import os
import time

def project_simplex(v, sum_constraint):
    """Exact Euclidean projection onto {x >= 0, sum(x) == R}.
    """
    R = float(sum_constraint)
    v = np.asarray(v, dtype=float)
    u = np.sort(v)[::-1]
    cumulative = np.cumsum(u) - R
    ind = np.arange(1, v.size + 1)
    rho = ind[u - cumulative / ind > 0][-1]
    return np.maximum(v - cumulative[rho - 1] / rho, 0.0)
def pgd(Q, c, sum_constraint, lr=0.01, max_iter=10**4, tol=1e-9, x0=None):
    """Projected gradient descent on {x >= 0, sum(x) == R}.
    """
    Q = np.asarray(Q, dtype=float)
    c = np.asarray(c, dtype=float).ravel()
    n = Q.shape[0]
    R = float(sum_constraint)
    x = np.full(n, R / n) if x0 is None else np.asarray(x0, dtype=float)
    x = project_simplex(x, R)
    QT = Q + Q.T                                   # gradient of x'Qx is (Q+Q')x
    it = 0
    for it in range(max_iter):
        grad = QT @ x + c
        x_next = project_simplex(x - lr * grad, R)
        delta = np.linalg.norm(x_next - x)
        x = x_next                                 # keep the newer iterate
        if delta < tol:
            break
    return x, float(x @ Q @ x + c @ x), it + 1
def pgd_multistart(Q, c, sum_constraint, restarts=32, seed=0, **kwargs):
    """PGD from the simplex centre plus random starts.
    """
    rng = np.random.default_rng(seed)
    n = Q.shape[0]
    solutions, energies, times, total_iters = [], [], [], 0
    for r in range(restarts):
        start = time.time()
        x0 = None
        if r > 0:
            x0 = rng.random(n)
            x0 *= sum_constraint / x0.sum()
        x, energy, iters = pgd(Q, c, sum_constraint, x0=x0, **kwargs)
        pgd_t = time.time()-start
        solutions.append(x)
        energies.append(energy)
        times.append(pgd_t)
        total_iters += iters
    return solutions, energies, times, total_iters

def read_qplib(path):
    with open(path) as fh:
        toks = [t for t in (line.split("#")[0].strip() for line in fh) if t]
    pos = 0
    def nxt():
        nonlocal pos
        val = toks[pos]
        pos += 1
        return val
    def defaulted(length):
        """Read a `default value` / `number of non-defaults` / entries block."""
        arr = np.full(length, float(nxt()))
        for _ in range(int(nxt())):
            i, v = nxt().split()
            arr[int(i) - 1] = float(v)
        return arr
    name = nxt()
    probtype = nxt()
    sense = nxt().lower()
    n = int(nxt())
    m = int(nxt())
    # quadratic terms of the objective (each unordered pair listed once)
    Q = np.zeros((n, n), dtype=np.float32)
    for _ in range(int(nxt())):
        i, j, v = nxt().split()
        i, j, v = int(i) - 1, int(j) - 1, float(v)
        if i == j:
            Q[i, i] = v
        else:
            Q[i, j] = Q[j, i] = 0.5 * v   # split the pair coefficient over both triangles
    b = defaulted(n)            # linear terms of the objective
    obj_const = float(nxt())
    # linear terms of the constraints
    A = np.zeros((m, n))
    for _ in range(int(nxt())):
        k, i, v = nxt().split()
        A[int(k) - 1, int(i) - 1] = float(v)
    inf = float(nxt())
    lhs, rhs = defaulted(m), defaulted(m)
    lb, ub = defaulted(n), defaulted(n)
    return dict(name=name, type=probtype, sense=sense, n=n, m=m,
                Q=Q, b=b, obj_const=obj_const,
                A=A, lhs=lhs, rhs=rhs, lb=lb, ub=ub, inf=inf)

inst_dir = "Instances/"
inst_name = "QPLIB_2761.qplib"
inst = read_qplib(os.path.join(inst_dir, inst_name))
n, Q, b = inst["n"], inst["Q"], inst["b"]
M = 0.5 * Q
c = b.copy()
sum_constraint = float(inst["rhs"][0])   # sum_i x_i = 1

restarts =100
pgd_solutions, pgd_energies, times,total_iters = pgd_multistart(Q=M,
                                                  c=c,
                                                  sum_constraint=sum_constraint,
                                                  restarts=restarts,)

pgd_energies_arr= np.array(pgd_energies)
lo = min(pgd_energies_arr)

print(f"Multi-start PGD all energies: {pgd_energies}")
print(f"PGD best energy:{lo}")

Multi-start PGD all energies: [0.10097685334659422, 0.0010485291160228403, 0.07039351549786885, 0.06693547209452383, 0.0010485291160227453, 0.01768112059763443, 0.09006403320120959, 0.0010485291160227796, 0.0010485291160227885, 0.06693263741421218, 0.0010485291160227468, 0.0716207293298459, 0.06606244717371348, 0.028507570296855817, 0.038198279039776746, 0.017681120597634478, 0.0010485291160228102, 0.0010485291160228138, 0.0010485291160228386, 0.09859646900306727, 0.06606244717371344, 0.06726658906968412, 0.07284124272182768, 0.07404253782634901, 0.001048529116022831, 0.06672539316893514, 0.0010485291160227575, 0.003990561779633381, 0.03833167156964212, 0.020248759338627555, 0.05296478049682289, 0.06606244717371344, 0.0010485291160228234, 0.010117643197959296, 0.09641535524734914, 0.0039905617796334265, 0.05732327918611205, 0.12691332537961625, 0.07020143907506993, 0.001048529116022746, 0.09641535524734914, 0.07222367403229096, 0.003990561779633378, 0.06606244717371343, 0.04578102721482212, 0.0010485291160227629, 0.11999728202282746, 0.09006403320120959, 0.0010485291160227557, 0.001048529116022779, 0.120796649804321, 0.0010485291160227798, 0.0010485291160227787, 0.048849061347131705, 0.001048529116022755, 0.0010485291160228336, 0.001048529116022795, 0.06531440522549319, 0.003990561779633435, 0.05728837829848721, 0.0010485291160227683, 0.0010485291160227826, 0.0010485291160228294, 0.04624217095385548, 0.017681120597634464, 0.06672539316893508, 0.00399056177963337, 0.0736028973621164, 0.1114976268045212, 0.0010485291160227995, 0.0010485291160228138, 0.10097685334659422, 0.001048529116022808, 0.06325572427054121, 0.0601078862513565, 0.01768112059763439, 0.001048529116022766, 0.003990561779633373, 0.11363332831651837, 0.0010485291160228023, 0.017681120597634395, 0.0010485291160227757, 0.0010485291160228308, 0.04489618562362864, 0.060010792088498695, 0.10456220462874036, 0.001048529116022779, 0.08748024698693273, 0.0010485291160227603, 0.09246536949307427, 0.06325572427054117, 0.0010485291160228062, 0.09819787932686895, 0.0010485291160227848, 0.014390150851714398, 0.10097685334659424, 0.07222367403229091, 0.12079664980432094, 0.09489267716053129, 0.01768112059763447]
PGD best energy:0.0010485291160227453