Computing objective function for QAOA from sampled bitstrings

PennyLane Help
1

Hello All,

I am trying to compute the objective function for the Max-Cut QAOA from the sampled bitstring, instead of computing the expectation value via qml.expval. However, when I run my code I get the error as:

p=1
Traceback (most recent call last):
  File "<path>\autograd\core.py", line 238, in vspace
    return VSpace.mappings[type(value)](value)
KeyError: <class 'numpy.int64'>

.
.
.
.
.

raise TypeError("Can't find vector space for value {} of type {}. "
TypeError: Can't find vector space for value 0 of type <class 'numpy.int64'>. Valid types are dict_keys([<class 'autograd.core.SparseObject'>, <class 'list'>, <class 'tuple'>, <class 'dict'>, <class 'numpy.ndarray'>, <class 'float'>, <class 'numpy.float64'>, <class 'numpy.float32'>, <class 'numpy.float16'>, <class 'complex'>, <class 'numpy.complex64'>, <class 'numpy.complex128'>, <class 'pennylane.numpy.tensor.tensor'>])

The cod is pasted below. Some help on this issue would be great. Thank you

@qml.qnode(dev)
def circuit(gammas, betas, n_layers=1):
    # apply Hadamards to get the n qubit |+> state
    for wire in range(n_wires):
        qml.Hadamard(wires=wire)
    # p instances of unitary operators
    for i in range(n_layers):
        U_C(gammas[i])
        U_B(betas[i])
    return qml.sample(wires=[0, 1, 2, 3])


def bitstring_to_int(bit_string_sample):
    bit_string = np.array([float(bits) for bits in bit_string_sample])
    return bit_string


def qaoa_maxcut(n_layers=1):
    print("\np={:d}".format(n_layers))
    # initialize the parameters
    init_params = 0.1 * np.random.rand(2, n_layers, requires_grad=True)

    def objective(params):
        strings = circuit(params[0], params[1], n_layers=n_layers)
        cost = 0.0
        for x_i in strings:
            for x_j in strings:
                x1 = qml.math.cast(x_i, np.float64)
                x2 = qml.math.cast(x_j, np.float64)
                cost += -x1 - x2 + 2.0 * x1 * x2
        return cost

    # initialize optimizer
    opt = qml.AdagradOptimizer(stepsize=0.01)

    # optimize parameters in objective
    params = init_params
    steps = 10
    for i in range(steps):
        params = opt.step(objective, params)
        if (i + 1) % 5 == 0:
            print("Objective after step {:5d}: {: .7f}".format(i + 1, objective(params)))

    return params

bitstrings1 = qaoa_maxcut(n_layers=1)[1]
2

Hi @Abhishek_Awasthi, welcome to the forum!

The issue here is that your parameters are being turned into an ArrayBox object. The way you can fix this is to add the following line in your objective function: params = qml.math.toarray(params)

This is how your function should look:

def objective(params):
        params = qml.math.toarray(params)
        strings = circuit(params[0], params[1], n_layers=n_layers)
        cost = 0.0
        for x_i in strings:
            for x_j in strings:
                x1 = qml.math.cast(x_i, np.float64)
                x2 = qml.math.cast(x_j, np.float64)
                cost += -x1 - x2 + 2.0 * x1 * x2
        return cost

That being said, I expect that you might have some optimization issues because if the value is being turned into ArrayBox, it might mean that you have some side effect within your cost function as explained here.

If you do get training issues my recommendation would be to go back to the original demo and change it little by little until you can pinpoint the line that is causing this behaviour.

Please let me know if this helps or if you’re still having the error!

3

Hi @CatalinaAlbornoz, Thank you for a prompt response. I have added you suggestion in the code. Good thing is it is now working without any errors. However, the trouble now is that the parameters are not getting updated at all, they stay exactly the same as the initialized values.

def qaoa_maxcut(n_layers=1):
print("\np={:d}".format(n_layers))
# initialize the parameters
init_params = 0.1 * np.random.rand(2, n_layers, requires_grad=True)
print('Initial', init_params, type(init_params))
def objective(params):
    params = qml.math.toarray(params)
    strings = circuit(params[0], params[1],  num_qubits, linear, quadratic, n_layers=1)
    cost = 0.0
    for x_i in strings:
        for x_j in strings:
            x1 = qml.math.cast(x_i, np.float64)
            x2 = qml.math.cast(x_j, np.float64)
            cost += -x1 - x2 + 2.0 * x1 * x2
    return cost

# initialize optimizer
opt = qml.AdagradOptimizer(stepsize=0.01)

# optimize parameters in objective
params = init_params
steps = 10
for i in range(steps):
    params = opt.step(objective, params)
    print(params, type(params))
    print('--------------------')        
    if (i + 1) % 5 == 0:
        print("Objective after step {:5d}: {: .7f}".format(i + 1, objective(params)))

Here is the output of the code.

p=1
Initial [[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
Objective after step     5: -120.0000000
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
[[0.04226406]
 [0.02070951]] <class 'pennylane.numpy.tensor.tensor'>
--------------------
Objective after step    10: -128.0000000

If I can get this working it would be really great. And the reason I am not computing the expected value from the circuit is because it is extremely slow.

One more thing, I tired to compare the data-structures in the original demo, and it turns out the params remains of the type <class 'autograd.numpy.numpy_boxes.ArrayBox'> . In my case though the params is of type <class 'pennylane.numpy.tensor.tensor'>. Is there a way around it? Thank you so much.

4

Hi @Abhishek_Awasthi, yes, I expected this might happen. The problem here becomes more complex. I’ve taken a look at it but haven’t found a solution so far. This is an issue with Autograd. I’ll see if I can find a way.

5

Hi @Abhishek_Awasthi :slightly_smiling_face:

My spidey sense tells me that the underlying issue here is not necessarily in the code, it’s conceptual.

You should be able to evaluate the output samples of a circuit in PennyLane, but if you are trying to optimize, things will fall apart. The reason is that PennyLane’s optimizers are based on gradient descent , and so we need to compute derivatives of your objective function with respect to params. In order to compute the derivatives of your objective function, we’ll also need to compute derivatives of circuit.

However, when you use qml.sample instead of qml.expval, the output of your circuit is not differentiable! This is because derivatives are nicely defined for real/floating-point values (like expectation values), but are not nicely defined for bitstrings. The notion of derivative requires that there are suitable points infinitesimally close to the point at which you’re evaluating, but for integers or bitstrings this is not true. So naturally PennyLane/Autograd raises an error message like TypeError: Can't find vector space for value 0 of type <class 'numpy.int64'> because it doesn’t know how to compute derivatives involving bitstrings

6

thank you so much. That would be great!

7

It is the cost function value I need to differentiate. It is dependent on the string but by the simple first principle of finite difference there must be a way to compute the gradient on PennyLane.

def update_parameters(objective, obj, params, n_layers, learning_rate, shift):
orig_params = params.copy()
new_params = params.copy()

for i in range(n_layers):
    for j in range(2):
        orig_params[j,i] += shift
        forward = objective(orig_params)
        orig_params[j,i] -= 2*shift
        backward = objective(orig_params)
        orig_params[j,i] += shift
        
        new_params[j,i] = orig_params[j,i] - learning_rate* (forward-backward)/ (2)
        
return new_params

But I was hoping that in PennyLane there must be better way.

Computation of expectation values is a computationally very expensive step at each iteration, specially if I simulate a 20 - 25 qubit QAOA problem.

8

Hi @Abhishek_Awasthi, I think the solution to your problem would be to change from sample() to probs() or state(). You can learn more about measurements here.

Is there any specific reason why you specifically need sample and not other kinds of measurements?

9

Hi,

Yeah the reason for avoiding the expval is that I have a QAOA circuit for 20 qubits and the computation of expval against the problem hamiltonian is consuming enormous amount of time.

Which is why I thought of computing the expectation value (aka cost function value) using the bit string and the corresponding QUBO. But this was causing troubles which we already have been discussing.

With probs, states I will again run into vectors of size 2^20 and again lose a lot of time.

10

Hi @Abhishek_Awasthi, I understand.

Unfortunately if you want to use gradients you cannot use sample. There may be other ways of finding speedups in your code while using expval. You can, for instance, use lightning.qubit as you device and diff_method='adjoint'. This can bring you significant speedups.

Please let me know if you need guidance onto how to use this diff_method.

I hope this helps.