Contents
Bayesian Optimization of Combinatorial Structures¶
Reference: https://altema.is.tohoku.ac.jp/QA4U/
Black-Box function¶
import numpy as np
import dimod
class QUBOBlackBox:
def __init__(self, num_variables):
qubo = np.random.normal(size=(num_variables, num_variables))
self.bqm = dimod.BQM(qubo, 'BINARY')
def __call__(self, x):
if len(np.array(x).shape) == 1:
return self.bqm.energy(x)
else:
return self.bqm.energies(x)
class PolyBlackBox:
def __init__(self, num_variables, degrees=3):
poly = {key: np.random.normal()
for key in combinations(range(num_variables), degrees)}
self.bpm = dimod.BinaryPolynomial(poly, 'BINARY')
def __call__(self, x):
if len(np.array(x).shape) == 1:
return self.bpm.energy(x)
else:
return self.bpm.energies(x)
def my_black_box(x):
energy = 0
return energy
from itertools import combinations
num_variables = 16
blackbox = QUBOBlackBox(num_variables)
# blackbox = PolyBlackBox(num_variables)
# blackbox = my_black_box
Initial data¶
def quadratic_feature(x):
num_variables = len(x)
features = [1] + list(x)
for k in range(num_variables - 1):
for l in range(k + 1, num_variables):
features.append(x[k] * x[l])
return np.array(features)
def generate_init_data(blackbox, num_variables, num_init_data):
X_init, y_init = [], []
for k in range(num_init_data):
x = np.random.choice([0, 1], num_variables)
y = blackbox(x)
X_init.append(quadratic_feature(x))
y_init.append(y)
return np.array(X_init), np.array(y_init)
X, y = generate_init_data(blackbox, num_variables, num_init_data=5)
Acquisition function with QUBO¶
def acquisitions_to_bqm(acquisitions, num_variables):
qubo = np.diag(acquisitions[1: num_variables + 1]) # 0: const
n = 0
for k in range(num_variables - 1):
for l in range(k + 1, num_variables):
qubo[k, l] = acquisitions[num_variables + 1 + n]
n += 1
return dimod.BQM(qubo, 'BINARY')
def fit_acquisitions(X, y, lam=0.01):
XX_lamI_inv = np.linalg.inv(np.dot(X.T, X) + lam * np.eye(X.shape[1]))
avg = np.dot(XX_lamI_inv, np.dot(X.T, y))
var = 0.5 * XX_lamI_inv
return np.random.multivariate_normal(avg, var)
Search by annealing¶
from neal import SimulatedAnnealingSampler
from tqdm import tqdm
sampler = SimulatedAnnealingSampler()
sampling_params = dict(num_reads=10, num_sweeps=1000)
y_min_hists = []
energy_mean_hists = []
similarity_hists = [] # only QUBOBlackBox
num_iters = 200
for _ in tqdm(np.arange(num_iters)):
acquisitions = fit_acquisitions(X, y)
bqm = acquisitions_to_bqm(acquisitions, num_variables)
sampleset = sampler.sample(bqm, **sampling_params)
x_new = sampleset.lowest().record[0].sample
y_new = blackbox(x_new)
y_min_hists.append(min(y_new, np.min(y)))
X = np.vstack([X, quadratic_feature(x_new)])
y = np.append(y, y_new)
energy_mean = np.mean((bqm.energies(X[:, 1: num_variables + 1]) - blackbox(X[:, 1: num_variables + 1]))**2)
energy_mean_hists.append(energy_mean)
# only QUBOBlackBox
similarity = np.sum((bqm.to_numpy_matrix() - blackbox.bqm.to_numpy_matrix())**2) / X.shape[1]
similarity_hists.append(similarity)
100%|██████████| 200/200 [00:02<00:00, 77.19it/s]
import matplotlib.pyplot as plt
plt.plot(y_min_hists)
plt.xlabel('Iterations')
plt.ylabel('Energy (min)')
plt.show()
plt.plot(energy_mean_hists)
plt.xlabel('Iterations')
plt.ylabel('Energy (mean)')
plt.show()
plt.plot(similarity_hists)
plt.xlabel('Iterations')
plt.ylabel('MSE of acquistion')
plt.show()
plt.imshow(bqm.to_numpy_matrix().astype(float))
plt.colorbar()
plt.show()
plt.imshow(blackbox.bqm.to_numpy_matrix().astype(float))
plt.colorbar()
plt.show()
