PGD Max Clique Benchmarking

import numpy as np
import os
import time

def read_dimacs_clq(file_path):
    """Read a DIMACS ascii .clq graph and return the adjacency matrix A.
    A[i, j] == 1 iff (i+1, j+1) is an edge (vertices in the file are 1-indexed).
    Lines: 'c ...' comments, 'p edge <n> <m>' header, 'e <u> <v>' edges.
    """
    n = None
    edges = []
    with open(file_path, 'r') as f:
        for line in f:
            parts = line.split()
            if not parts:
                continue
            tag = parts[0]
            if tag == 'c':
                continue
            elif tag == 'p':
                # p edge n m   (some files use 'p col n m' / 'p clq n m')
                n = int(parts[2])
                m = int(parts[3])
            elif tag == 'e':
                edges.append((int(parts[1]), int(parts[2])))
    if n is None:
        raise ValueError(f"no 'p' header found in {file_path}")
    A = np.zeros((n, n), dtype=np.uint8)
    for u, v in edges:
        if u == v:
            continue  # ignore self-loops
        A[u - 1, v - 1] = 1
        A[v - 1, u - 1] = 1  # undirected -> symmetric
    return A, n, m
def build_clique_hamiltonian(A):
    A = np.asarray(A, dtype=float)
    return -A

def clique_number_from_energy(energy, sum_constraint=1.0):
    R = float(sum_constraint)
    return 1.0 / (1.0 + energy / (R * R))
def extract_clique(x, A, tol=1e-6, scale=None):
    """Read a vertex set off the support of x and check it really is a clique.
    `tol` is RELATIVE: v is in the support when x[v] > tol * scale, with scale
    defaulting to max(x). 
    
    Returns (vertices, is_clique). Vertices are 0-indexed into A; add 1 to get
    the labels used in the DIMACS file.
    """
    x = np.asarray(x, dtype=float)
    if scale is None:
        scale = x.max() if x.size else 0.0
    if scale <= 0.0:
        return np.array([], dtype=int), False
    vertices = np.flatnonzero(x > tol * scale)
    k = len(vertices)
    if k == 0:
        return vertices, False
    sub = np.asarray(A, dtype=bool)[np.ix_(vertices, vertices)]
    # A clique on k vertices has k*(k-1) ones in its (symmetric, zero-diagonal) block.
    is_clique = bool(sub.sum() == k * (k - 1))
    return vertices, is_clique
def _extend_to_maximal(chosen, cand, Ab, x):
    """Grow `chosen` while candidates remain, highest weight first.
    """
    while True:
        idx = np.flatnonzero(cand)
        if idx.size == 0:
            return chosen
        w = x[idx]
        best = w.max()
        ties = idx[w >= best - 1e-12 * max(1.0, abs(best))]
        if ties.size > 1:
            deg = Ab[np.ix_(ties, idx)].sum(1)
            v = int(ties[np.argmax(deg)])
        else:
            v = int(ties[0])
        chosen.append(v)
        cand &= Ab[v]
        cand[v] = False
def greedy_clique_from_weights(x, A):
    """Repair step: greedily grow a genuine clique, taking vertices in order of x.
    """
    Ab = np.asarray(A, dtype=bool)
    x = np.asarray(x, dtype=float)
    cand = np.ones(Ab.shape[0], dtype=bool)
    return np.array(sorted(_extend_to_maximal([], cand, Ab, x)), dtype=int)

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

instance_dir = "Instances/"
instance_name = "keller4.clq"
instance_path = os.path.join(instance_dir, instance_name)
try:
    A, n, m = read_dimacs_clq(instance_path)
    print("file loaded sucessfully")
except FileNotFoundError:
    print(f"{instance_path} does not exist")

file loaded sucessfully

name = os.path.splitext(instance_name)[0]
print(f"{name}: n = {n} vertices, m = {m} edges (declared)")
print(f"A shape: {A.shape}")
print(f"edges in A: {int(A.sum()) // 2}")
print(f"symmetric: {np.array_equal(A, A.T)}")
print(f"self-loops: {int(np.trace(A))}")
print(f"density: {A.sum() / (n * (n - 1)):.5f}")

keller4: n = 171 vertices, m = 9435 edges (declared)
A shape: (171, 171)
edges in A: 9435
symmetric: True
self-loops: 0
density: 0.64912

c = np.zeros(n)
H = build_clique_hamiltonian(A)
sum_constraint = 1
restarts = 100
learning_rate = 0.01
max_iter = 2 * 10**4

pgd_start = time.time()
solutions, energies, times, total_iters = pgd_multistart(
    Q=H,
    c=c,
    sum_constraint=sum_constraint,
    restarts=restarts,
    lr=learning_rate,
    max_iter=max_iter,
)
pgd_time = time.time() - pgd_start

best_energy = min(energies)
best_solution = solutions[int(np.argmin(energies))]
omega = clique_number_from_energy(best_energy, sum_constraint)
support, isclique = extract_clique(best_solution, A)
# Round every restart, not just the lowest-energy one.
cliques = [greedy_clique_from_weights(sol, A) for sol in solutions]
greedy = max(cliques, key=len)

print(f"restarts:{restarts}, total PGD iterations:{total_iters}")
print(f"best energy over restarts:{best_energy:.6f}")
print(f"support size:{len(support)}, support is a clique:{isclique}")
print(f"clique vertices (1-indexed):{(greedy + 1).tolist()}")
print(f"{'energy':>16}{'omega_est':>12}{'clique':>9}{'time (s)':>12}")
print(f"{best_energy:>16.6f}{omega:>12.3f}{len(greedy):>9}{pgd_time:>12.3f}")
# Never report a set that is not actually a clique.
assert extract_clique(np.isin(np.arange(n), greedy).astype(float), A)[1], \
    "greedy repair returned a non-clique"

restarts:100, total PGD iterations:77765
best energy over restarts:-0.888889
support size:15, support is a clique:False
clique vertices (1-indexed):[3, 24, 28, 37, 59, 77, 87, 105, 111]
          energy   omega_est   clique    time (s)
       -0.888889       9.000        9       0.735