377 simulated annealing replica exchange

[djelovina]

Issue #377

[djelovina]

Current version contains Python prototype (to me remove in the final version) and the cpp smoke test code. The algorithm need to be implemented.

[anyzelman]

This is in draft state right? (If yes, please click "convert to draft" somewhere above-right from here=)

[anyzelman]

Target set to milestone 0.9

[djelovina]

I would recommend to update unit tests in a separate PR, since tracking down test failures might be easier with clearly separated PRs @GiovaGa

[GiovaGa]

Tests are now failing because the correct behavior depends on #397

[GiovaGa]

Rebased adding fix of #398 . Smoke tests should now pass

[djelovina]

@GiovaGa
I think it might be useful to have exact solution for a randomly generated matrices, to be able to compare solution quality vs exact solutions. This python script is doing that. Is it too complicated to adopt these matrices for our tests? They need to be converted to proper format.

import numpy as np
import itertools
from collections import defaultdict

def generate_sparse_planted_qubo(
    n,
    degree,
    weight_range=(1.0, 1.0),
    seed=0
):
    """
    Generate a sparse QUBO with a planted solution.

    Returns:
        Q_diag: dict {i: Q_ii}
        Q_off: dict {(i,j): Q_ij} with i<j
        x_star: planted solution (0/1)
        E_star: known optimal energy
    """
    rng = np.random.default_rng(seed)

    # Planted solution
    x_star = rng.integers(0, 2, size=n, dtype=np.int8)

    Q_diag = defaultdict(float)
    Q_off = defaultdict(float)

    E_star = 0.0

    for i in range(n):
        # choose neighbors (avoid self-loops)
        neighbors = rng.choice(
            n - 1,
            size=degree,
            replace=False
        )
        neighbors = neighbors + (neighbors >= i)

        for j in neighbors:
            if j < i:
                continue  # enforce i < j

            w = rng.uniform(*weight_range)
            b = x_star[i] ^ x_star[j]

            if b == 0:
                # w(x_i + x_j - 2 x_i x_j)
                Q_diag[i] += w
                Q_diag[j] += w
                Q_off[(i, j)] += -2.0 * w
            else:
                # w(1 - x_i - x_j + 2 x_i x_j)
                Q_diag[i] += -w
                Q_diag[j] += -w
                Q_off[(i, j)] += 2.0 * w
                E_star += w

    return Q_diag, Q_off, x_star, E_star

def qubo_energy(Q_diag, Q_off, x):
    E = 0.0
    for i, v in Q_diag.items():
        E += v * x[i]
    for (i, j), v in Q_off.items():
        E += v * x[i] * x[j]
    return E


def brute_force_check(Q_diag, Q_off, x_star, E_star):
    n = len(x_star)

    min_energy = float("inf")
    argmins = []

    for bits in itertools.product([0, 1], repeat=n):
        x = np.array(bits, dtype=np.int8)
        E = qubo_energy(Q_diag, Q_off, x)

        if E < min_energy - 1e-9:
            min_energy = E
            argmins = [x.copy()]
        elif abs(E - min_energy) < 1e-9:
            argmins.append(x.copy())

    print(f"Planted energy   : {-E_star}")
    print(f"Minimum found    : {min_energy}")
    print(f"# ground states  : {len(argmins)}")

    planted_ok = any(np.array_equal(x, x_star) for x in argmins)

    return {
        "planted_is_optimal": planted_ok,
        "energy_matches": abs(min_energy + E_star) < 1e-9,
        "degeneracy": len(argmins),
        "ground_states": argmins,
    }

# small instance ONLY
n = 18
degree = 4

Q_diag, Q_off, x_star, E_star = generate_sparse_planted_qubo(
    n=n,
    degree=degree,
    seed=1
)

result = brute_force_check(Q_diag, Q_off, x_star, E_star)

print(result)

[GiovaGa]

That seems definitely doable, I will try to implement it and see what happens