Programming the D-Wave One^TM System using the BlackBox compiler

Contents

Section 1

Aim, audience and required background

This material was developed to help those interested in programming the D-Wave One^TM System at a high level. This tutorial introduces the basic ideas behind the BlackBox algorithm, which is closely related to the concept of a compiler for the D-Wave processor. The term 'compiler' is used as the BlackBox code which is provided in the devloper pack is the lowest level library that allows programming of the system without an explicit knowledge of the architecture of the chip.

All of the software necessary for developers to begin learning how to program the D-Wave One^TM System can be downloaded from the DOWNLOAD DEVKITS page. As the developer portal is currently in BETA test stage, you will need to have a D-Wave developer account to download these packs. To follow this tutorial it may be helpful to have read two of the previous installments: Quantum Computing and Energy Programming Primer and Getting Started: Installing the D-Wave developer pack and an introductory 'hello multiverse' tutorial.

The high level programming language used in the included examples and source code is Python 2.6. Some experience of using the Python programming language is helpful but not required.

What you will learn

By following through the material in this white paper, you will learn:

Software requirements

In order to proceed you should install several free software packages.

SECTION 1

As with conventional processors and computing systems, it is possible to develop on a D-Wave One^TM System by programming directly in the machine language of the processor. D-Wave's developer tools, up to and including devpack 1.3, provided this level as the only option available.

In devpack 1.4, we have released a new, higher level programming interface to the D-Wave One^TM System which is called the D-Wave BlackBox compiler - BlackBox for short. Under the hood, the D-Wave One^TM System is an extremely exotic and complex machine. Programming at the machine language level is extremely difficult, even for our own internal developers. Abstracting away from this underlying complexity is critical for making the system broadly accessible.

The BlackBox compiler framework simplifies the ease of use of the D-Wave One^TM System, and dramatically extends its possible range of uses. Programming using BlackBox is so easy, users can be up and running problems on it within minutes. Here is how it works.

First, a developer codes a function which takes as input a string of bits - zeroes and ones - and outputs a real number. We'll call the function provided by the developer the generating function, and in the next section we'll provide some examples arising in high value industrial applications. Mathematically we write the generating function like this

where the \(x_k \) variables (there are a total of N of these) are binary (0 or 1) variables. The value returned by the function \(G(x_1,x_2,...,x_N ) \) has to be a real number.

Let's take a look at a simple example:

In this case, we have N=2 variables, and the function \(G(x_1,x_2 ) \) evaluates to real numbers (well actually integers, but that's OK too!). Specifically, plugging in the 2N = 4 possible inputs, these are \(G(0,0)=0 \), \(G(1,0)=1 \), \(G(0,1)=2 \), and \(G(1,1)=3 \).

This particular example is really simple. However you as a developer can make G as gnarly and complicated as you want, as long as it takes as input binary variables and outputs real numbers. For example it could be something bizarre like:

However it doesn't even need to arise from a closed form mathematical expression. For example, it could be related to the outcome of a series of settings of control switches on a complex system, such as a flight control system for an airplane.

This function can be coded up using any programming language (here we will use Python), and will be evaluated entirely using a conventional computing system. For some problems, the conventional computing system you'd want to use to do this is a high performance computing system, such as a cluster, a supercomputer or a cloud-based solution; however it could just as easily be your laptop - it's up to you, and depends on how hard it is to evaluate the value of the generating function given an input.

Choosing a generating function matched to your application

A developer chooses \(G(x_1,x_2,...,x_N ) \) so that the lower the real number returned by the function, the 'better' the input bit string \((x_1,x_2,...,x_N ) \) was. The input bit string represents a 'potential solution', and the number returned by the function given that input gives a measure of the goodness of the potential solution - the lower \(G(x_1,x_2,...,x_N )\) is, the better the potential solution \((x_1,x_2,...,x_N ) \) is.

Let's take a look at our first example, where \(G(x_1,x_2 )=x_1+2x_2 \). In this case if we look at the returned values for all possible input combinations \(G(0,0)=0\), \(G(1,0)=1\), \(G(0,1)=2\), and \(G(1,1)=3\), we see that the choice of \(G(x_1=0,x_2=0)=0\) gives the lowest (and therefore 'best') value of G. The generating function is sometimes called an optimization objective function, although as a developer you can use the output of BlackBox for more than just optimization (more on this shortly).

So where does G come from in the first place? It comes from your needs as an applications developer. If you are building an application where the above concept is useful (or even critical) to your work, the D Wave One provides a new tool to attack that problem.

An example: the Subset Sum Problem

For example, let's say that you wanted to solve the Subset Sum Problem (SSP). In this problem, you are given a set of integers, and you want to determine whether there is a non-empty subset whose sum is zero. For example, the set might be \({-7, -3, -2, 5, 8}\) - this is the instance provided in the Wikipedia page on SSP. In this case the answer is 'yes', because the subset \( \{ 3, -2, 5\} \) sums to zero. As an applications developer, you could code up the SSP in the following way (note that there are many ways you could try to do this). Define five variables \( (x_1,x_2,x_3,x_4,x_5 ) \) where if a variable is equal to 1 it means "please include me in the sum". In this case the generating function for the above example could be

where the values of the coefficients are just the values in the original set. The first term fixes the minimum possible value of that part of G to be zero, and the second term penalizes the solution with all variables equal to zero (which otherwise would show up as an allowed solution!) but not any of the others. If we can find a configuration of variables for which \(G(x_1,x_2,x_3,x_4,x_5 )=0 \) we've got a solution to our SSP instance.

Up to this point everything has been done in a conventional framework, and all we've done is set up some machinery defining a function. Now we add the D-Wave One^TM System in the following way.

Imagine the entire computing system as consisting of two complimentary parts - think of these as the 'right brain' and 'left brain' of a hybrid computing system (see Figure 1 below). The 'left brain' is a conventional computing system whose rational, tedious job it is to compute the value of the generating function given 'guesses' at potential solutions. In general this involves a large amount of computation of the sort conventional computing systems excel at. The 'right brain' is the D-Wave One, which suggests 'creative' potential solutions, using the results obtained by the conventional computing system to quickly hone in on better and better solutions.

This split allows hybrid systems to be built that can deal with enormous amounts of data and extremely complex generating functions. The information that passes between the conventional computing system and the D-Wave One^TM System is extremely small - it's just the bit strings coming from the system representing its guesses at good answers flowing from the D-Wave One^TM System to the conventional system, and the real numbers characterizing how good those guesses were, flowing from the conventional system to the D-Wave One^TM System. The amount of data that might be required to compute the value of the generating function could be (and often is) enormous - but we can use all of the standard tactics for dealing with this using conventional systems architecture.

Figure 1. The BlackBox programming paradigm separates the evaluation of the generating function, and the process of generating potential solutions. The evaluation of the function happens in a conventional computing system, while the solution generation happens in the D-Wave One. Information flow between the two is extremely low bandwidth and is restricted to bit strings representing potential solutions flowing from the D-Wave One^TM System to the conventional system, and real numbers representing the values of the generating function evaluated for those guesses flowing from the conventional system to the D-Wave One^TM System.

Using this hybrid system is simple. The developer provides the generating function, initiates the computation, and then the system starts 'thinking' about the generating function in the following, iterative way:

When the D-Wave One^TM System generates a series of guesses at Step 4, it is doing something rather complicated. Here we won't go into any detail as to specifically what's going on in this step, but we'll outline the process. The following section is rather technical; you can skip it if you'd just like to start using BlackBox.

In order to give some intuition for what the D-Wave One^TM System is doing in Step 4, we have to take a brief detour into statistical physics.

In physics, often a physical system is described by writing down a function whose inputs can be thought of as specific physical states of the system, and whose outputs are the energies of those states.

If the physical system of interest (call it the "central system") is in a state called thermal equilibrium, the probability of being in any particular state has a known functional form called the Boltzmann distribution. This distribution has the form

which means, "if the system whose allowed energies are given by \(G(x_1,x_2,...,x_N ) \) is in thermal equilibrium at temperature T, the probability of the system being in state \((x_1,x_2,...,x_N ) \) is \( P(x_1,x_2,...,x_N ) \)". Here k is the Boltzmann constant (don't worry too much about its specific value, it won't matter for us here) and Z is something called the partition function, which is

Note that actually calculating \(P(x_1,x_2,...,x_N ) \) given \(G(x_1,x_2,...,x_N ) \) is not easy at all. But let's not worry about that right now. Let's just imagine that we had a physical system where the probabilities of measuring the bits \( (x_1,x_2,...,x_N ) \) were somehow magically given by the Boltzmann distribution. What would this look like, and what could we do with it?

Let's take the specific example of the Subset Sum Problem instance we used earlier:

Because we only have five bits here, we can explicitly write down all of the 25 terms we need to compute all of the terms in the Boltzmann distribution. It's a little tedious, but let's do it anyway!

Figure 2. Enumerating the function.

Notice that there is a single solution (highlighted in orange) that gives \( G(x_1,x_2,x_3,x_4,x_5 )=0 \): the solution with \((x_1=0,x_2=1,x_3=1,x_4=1,x_5=0) \).

Let's assume that the probability distribution over these results is a Boltzmann, and let's for the moment just arbitrarily set kT=1. We can then compute the actual probabilities for each state. These are the values in the rightmost column.

Now every time we measure the state of the system, we get one of the 32 states above with the probabilities in the righthand column (as long as the system is in thermal equilibrium at kT=1). If we wanted to find a solution to the Subset Sum Problem, given a physical system with this probability distribution, we would repeatedly sample from the probability distribution until we either found a solution or got bored of drawing samples. In this case, we see that there is a \(P_{SS}=39.3 % \) (SS stands for "single shot") chance of drawing the interesting state \( (x_1=0,x_2=1,x_3=1,x_4=1,x_5=0) \) per sample.

So if we draw a total of K samples, the probability of seeing the state we want at least once is

So if we drew say K=10 samples in this case, the probability of seeing the state we want at least once is

Back to BlackBox

The algorithms underlying the BlackBox procedure attempt to sculpt the probability distribution returned by the D-Wave One^TM System to be close to the Boltzmann distribution you'd get over the generating function. As the iteration proceeds, the probability distribution changes and, if every goes well, starts converging to something close to the Boltzmann distribution described in the previous section.

What BlackBox returns is the best solutions found during this iterative procedure. These solutions can be used in two modes.

Often a developer is looking for the best solution found. In this optimization mode, the output of this procedure is the bit string BlackBox found that returned the lowest value of the generating function.

Sometimes a larger set of samples is desirable; we call this the sampling mode. For example, if the user wants to generate not only the best solution, but a set of 'good' alternate solutions, the system natively provides this capability - the user records the guesses generated at Step 4, the number of times these guesses were generated, and the quality of the guesses. BlackBox can be configured to return a set of the best answers found, and how many times each was seen.

Let's code up our example to see how BlackBox works in practice. You should be running devpack 1.4. If you haven't downloaded and installed it yet, go to the Developer Portal and do it!

First we import some things.

from dwave_sapi import LocalConnection, BlackBoxSolver
from numpy import dot, array, prod

Next, we define a function object. This is the place where we define the generating function.

There are two parts to the function object. The first part allows a user to pass problem-specific variables into the function object. This is the section headed by the def __init__ call. This is where you pass in any variables you'll need to compute the value of the generating function.

The second part is the section headed by the def __call__ call. This is where you actually calculate the value of the generating function. The way BlackBox works is that it sends in a list, called states, and a number, called numStates. numStates is an integer that represents how many states are being passed in at that time. The states list is a concatenation of numStates total states. For example, if your problem had 3 variables, then numStates could be 4, in which case states would be a list consisting of four possible solutions concatenated together, for example states = [-1, 1, 1, 1, -1, -1, -1, -1, -1, 1, 1, 1] (this passes in the four states [-1, 1, 1], [1, -1, -1], [-1, -1, -1] and [1, 1, 1] all at once for evaluation). As a developer you won't change the syntax of the def __call__(self, states, numStates): line - just leave it like this!

OK here goes. First we create the ObjClass class.

class ObjClass(object):
#
# We send in any parameters we want to. In this case this is the values of our particular subset sum
# instance.
#
    def __init__(self, subset_sum_array):
        self.subset_sum_array = subset_sum_array
#
# Now we define and compute the value of the generating function.
#
    def __call__(self, states, numStates):
#
# First we get the length of each state; this is just the number of variables in our problem. In our
# example, this will be five.
#
        stateLen = len(states)/numStates
#
# An important point is that BlackBox natively solves problems assuming that the variables are +1/-1
# variables. In cases where you'd prefer to use 0/1 variables (such as in this example), you need to
# explicitly convert. Here we convert the states list into a new list called states_bin which stores the
# suggested answers in 0/1 variables.
#
        states_bin  = [(item+1)/2 for item in states] # converting to 0/1 from -1/+1
#
# We now create a list called ret where we'll store all the results of computing the value of the
# generating function on the states sent in.
#
        ret = []
#
# Now cycle over all the states sent in.
#
        for state_number in range(numStates):
#
# The w array stores the current state of interest.
#
            w = array(states_bin[state_number*stateLen:(state_number+1)*stateLen])
#
# Now we compute the value of the generating function, given the bit string w.
#
            result = dot(self.subset_sum_array, w)**2+prod(1-w)
#
# We then append that result to the ret list.
#
            ret.append(result)
#
# Once we've done this for all the values sent in, the result is returned as a tuple.
#
        return tuple(ret)
#
##########################################################################################################]

Now that we've coded that up, we are now ready to use BlackBox to solve the problem!

##########################################################################################################
#
# This is the main program for using BlackBox to solve a subset sum instance.
#
# First we establish a connection to a local solver.
#
conn = LocalConnection()
#
# Next we choose a solver type. Here we'll use the optimization solver defined over the full 128 variable
# graph representing a fully functional Rainier processor.
#
solver = conn.get_solver('c4-sw_optimize')
#
# The blackbox_parameter is an integer that sets a bunch of solver parameters within the BlackBox
# algorithm. As a developer you actually have a lot more flexibility than just setting a single parameter.
# But for the time being you can just use the rule of thumb that the bigger this is, the more work
# BlackBox will put into finding the best possible solution. Here we'll set it to 10. You can experiment
# with decreasing or increasing it.
#
blackbox_parameter = 10
#
# This is an array giving the values of our initial set.
#
subset_sum_array = array ([-7, -3, -2, 5, 8])
#
# num_vars is just the length of the set.
#
num_vars = len(subset_sum_array)
#
# obj is an object of the class ObjClass. Pass variables in here (in this case subset_sum_array).
#
obj = ObjClass(subset_sum_array)
#
# Set up an instance of BlackBoxSolver
#
blackbox_solver = BlackBoxSolver(solver)
#
# Now run BlackBox. Here you have to pass in obj, num_vars, and cluster_num, but everything else is
# optional. For most users you can just do what I'm doing here and tie all of these parameters to a
# single common value.
#
blackbox_answer = blackbox_solver.solve(obj, num_vars, cluster_num = 10, \
    min_iter_inner = blackbox_parameter, max_iter_outer= blackbox_parameter, \
    unchanged_threshold=blackbox_parameter, max_unchanged_objective_outer=blackbox_parameter, \
    max_unchanged_objective_inner = blackbox_parameter, \
    unchanged_best_threshold = blackbox_parameter, verbose=0)
#
# blackbox_answer returns a list of 1/-1 variables denoting the best solution it found. Since we want
# 0/1 variables we convert it to that, and cast it as an array.
#
blackbox_answer_bin = array([(item+1)/2 for item in blackbox_answer])
#
# Now we output the answer!
# 
print 'The best bit string we found was:',blackbox_answer_bin
#
# This stores the subset found in a list.
#
subset_list = []
for k in range(num_vars):
    if blackbox_answer_bin[k] ==1:
        subset_list.append(subset_sum_array[k])
#
print 'The subset this corresponds to is:', subset_list
#
# This computes the value of the generating function at the best answer found.
#
energy_of_best_solution_found = dot(subset_sum_array, blackbox_answer_bin)**2+prod(1-blackbox_answer_bin)
print 'Its energy is:', energy_of_best_solution_found
#
if energy_of_best_solution_found==0:
    print 'We found a solution... this set has a subset that sums to zero!'
#
# And that's it!
#
##########################################################################################################

Running this code should produce the following output:

>> The best bit string we found was: [0 1 1 1 0]
>> The subset this corresponds to is: [-3, -2, 5]
>> Its energy is: 0
>> We found a solution... this set has a subset that sums to zero

To change this code to suit your own particular generating function, all you need to do is pass in whatever variables you need to ObjClass, and then modify the definition of the result variable in ObjClass. It's that simple!

The above code finds a subset that sums to zero. How would you modify the code so that it looks for a subset that sums to any given number?

As a test case, you can try the set {267, 493, 869, 961, 1000, 1153, 1246, 1598, 1766, 1922} and see if you can find a subset that sums to 5842. This example is from here: https://netfiles.uiuc.edu/c-blair/www/papers/second.pdf.