PGD Planted 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 loadQ(path, n):
    tri = np.load(path, mmap_mode="r")
    Q = np.empty((n,n), dtype=tri.dtype)
    iu = np.triu_indices(n)
    Q[iu] = tri
    Q.T[iu] = tri
    return Q

# Instance details
optimal_energy = 12000
sum_constraint = 100
num_var = 2000
k = 44
ub =10
seed =100
name = f"STQP_n_{num_var}_k_{k}_R_{sum_constraint}_seed_{seed}_ub_{ub}"
instance_path = os.path.join(f"Instances/{name}.npy")
c = np.zeros(num_var)
try:
    Q= loadQ(instance_path,num_var)
    print("loaded file sucessfully.")
except FileNotFoundError:
    print(f"File {instance_path} does not exist.")

loaded file sucessfully.

# solver parameters
restarts = 10
learning_rate = 0.01
max_iter = 2 * 10**4

# solve using multistart PGD
pgd_start = time.time()
solutions, energies, times, total_iters = pgd_multistart(
    Q=Q,
    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))]
indices = np.where(best_solution>1e-6)[0]

print(f"restarts:{restarts}, total PGD iterations:{total_iters}")
print(f"best energy over restarts:{best_energy:.6f}")
print(f"support size:{len(indices)}")
print(f"time taken by PGD {pgd_time}")

restarts:10, total PGD iterations:16695
best energy over restarts:12815.673891
support size:14
time taken bt PGD 152.516615152359

# compute relative gap
tol = 1e-7
abs_gap = best_energy-optimal_energy
relative_gap = (round(abs_gap,6)*100)/round(optimal_energy,6)

# print result
if abs_gap>tol:
    print(f"Gradient descent solution is not optimal")
    print(f"Absolute Gap:{abs_gap}")
    print(f"Relative Gap(%):{abs(relative_gap)}")
else:
    print(f"Optimal solution found:{best_energy}.")

Gradient descent solution is not optimal
Absolute Gap:815.6738909034775
Relative Gap(%):6.797282425