How to solve QUBO problems in Q#?

Question

Short version:

I'm trying to solve a traveling salesman problem very similar to the traveling Santa example here: http://quantumalgorithmzoo.org/traveling_santa/, which is also included in the samples of the Microsoft Quantum samples here: https://github.com/microsoft/Quantum/tree/main/samples/simulation. In that example, they assume some parameters beta and gamma that yield favorable odds of finding the optimal route. The problem is: how do you get these parameters? But a more general question that I have is: How would you solve a QUBO problem (with a Hamiltonian of the form $H = -\sum_i h_i \sigma_i^z -\sum_{i,j} J_{i,j} \sigma_i^z \sigma_j^z$) in qsharp?

What I've tried:

1. Building upon the QAOA sample, the first thing I did was cheat: I used a classical optimizer to solve for optimal values for beta and gamma, minimizing the energy. And I calculated the energy by dumping the quantum registry to a file. With the probabilities for each state, the estimated value for the energy is simply $\sum_{states} p_{state} E_{state}$.

2. Of course, on Azure Quantum / on real quantum hardware you don't have access to the probabilities. So I tried to find ways to get precise estimates of the energy. This is where I'm struggling given the samples and documentation. I have a registry of qubits and a Hamiltonian equation that I would like to plug in, but the EstimateEnergy function in Q# either takes JordanWignerEncodingData or a statePrepUnitary and qpeUnitary. In both cases I don't really understand how I would construct them and what they do/why I need them. Efforts to estimate the energy from the phase estimation failed, but that might be due to my lack of understanding. If this is indeed a good way to solve optimization problems, are there any good resources to understand this better?

3. The final thing I tried was the principal of slowly changing the Hamiltonian from one that has an easy to prepare groundstate, to the Hamiltonian corresponding to the optimization problem you want to solve. The example and explanation are here: https://github.com/microsoft/Quantum/blob/main/samples/simulation/ising/adiabatic/AdiabaticIsing.qs#L14. Unfortunately, I seem to get stuck in different local minima depending on the rate, and none of them actually come close to the real solution. So I found this method to be not very reliable either.

I understand the question is very similar to this one, but even after reading the answer there, I'm still not sure if what I'm trying makes sense, and how to make it work in Q#. So I am hoping for more concrete answer, or literature suited for developers who followed a quantum physics course many many years ago.

Answer 1

To answer your first question, QAOA is an application of a hybrid classical-quantum algorithm, so using a classical optimizer is a perfectly valid solution here. I would suggest using e.g. scipy.optimize any other of your favorite optimization tools.

Your second question is regarding measuring the energy of the quantum state. Indeed, Q# does not allow inspection of the quantum state, because as you noted this is not possible on physical hardware either.

However, there are some ways to retrieve this info.

In an experiment on real hardware, what you would do is simply measure all qubits over and over again, let's say N times, which gives you a list of N bitstrings. If N is large this will give you a reliable probability distribution of all the qubit states. In simulation this can be done by measuring all the qubits in the qubit register in a for-loop, and counting the occurrences of each possible iteration. This should be straightforward to implement but I will suggest a solution below that uses existing library functionality which you might find helpful.

Instead of measuring the results of each qubit, we're going to implement a handy operation that lets us measure the probability of each valid state and the probability of obtaining an invalid state.

First of all, as explained in the blog post, there are the probabilities to get one of 3 valid states, and there's the probability of getting an invalid state:

// Allowed states
let state1 = [One, One, One, One, Zero, Zero];
let state2 = [Zero, One, Zero, One, One, One];
let state3 = [One, Zero, One, Zero, One, One];
let states = [state1, state2, state3];
mutable result = new Double[Length(states) + 1];

We'll keep track of the probability of getting an invalid state as the results come in by subtracting these values from 1.0.

mutable otherProb = 1.0; // Probability of any invalid state

So what we'll do is loop over all valid states, and then measure the probability of each, and return them and the probability of getting any invalid state in an array result.

for ((index, state) in Enumerated(states)) {
    let prob = MeasureProbabilityForState(state, numSegments, weights, couplings, timeX, timeZ, numMeasurements);
    set result w/= index <- prob;
    set otherProb -= prob;
}

set result w/= Length(states) <- otherProb;

The energy, or cost that we want to minimize, can then be calculated by multiplying the probabilities by the cost of each state.

(Obviously for this example this is a bit silly, as we can simply calculate the cost of each of the above allowed states and then pick the lowest value. However, in other quantum applications the space of valid states could for example be much larger, and in that case we would only need to evaluate the cost function for the states that have a nonzero probability.)

So, how to implement MeasureProbabilityForState?

We can use the operation EstimateFrequency for this. This operation measures the probability of getting a Zero result for given state prep and measurement operations. See docs here: https://docs.microsoft.com/en-us/qsharp/api/qsharp/microsoft.quantum.characterization.estimatefrequency.

The implementation can look something like this:

    operation MeasureProbabilityForState(
        state : Result[],
        numSegments: Int, 
        weights : Double[], 
        couplings : Double[], 
        timeX : Double[], 
        timeZ : Double[],
        numMeasurements: Int
    ) : Double {
        return EstimateFrequency(
            ApplyQAOA(_, numSegments, weights, couplings, timeX, timeZ), 
            MeasureRegisterIsInState(_, state),
            numSegments, 
            numMeasurements
        );
    }

where we have to check if the register is in the desired state and return Zero if it is:

    operation MeasureRegisterIsInState(register : Qubit[], state : Result[]) : Result {
        let result = MultiM(register);
        if (All<(Result, Result)>(EqualR, Zipped(result, state))) {
            return Zero;
        } else {
            return One;
        }
    }

and

    operation ApplyQAOA(
            x: Qubit[],
            numSegments: Int, 
            weights : Double[], 
            couplings : Double[], 
            timeX : Double[], 
            timeZ : Double[]
        ) : Unit {
        ApplyToEach(H, x); // prepare the uniform distribution
        for ((tz, tx) in Zipped(timeZ, timeX))
        {
            ApplyInstanceHamiltonian(numSegments, tz, weights, couplings, x); // do Exp(-i H_C tz)
            ApplyDriverHamiltonian(tx, x); // do Exp(-i H_0 tx)
        }
    }

This is just one way of using EstimateFrequency to solve this problem; I hope this example helps you figure out the best way that works for your quantum application.