Calculating expectation values in Strawberryfields

Hello! I have been trying to find the expectation value ofH_p = (n0 + n1 - 5)**2 using strawberryfields. The expansion gives H_p = n0**2 + 2*n0*n1 - 10*n0 + n1**2 - 10*n1 + 25. and to calculate its expectation value with a given state, I created this function.

def expval2():
    n02 = state.mean_photon(0)[1] + state.mean_photon(0)[0]**2
    n0n1 = 2*state.number_expectation([0, 1])[0]
    n0z = -10*state.mean_photon(0)[0]
    n12 = state.mean_photon(1)[1] + state.mean_photon(1)[0]**2
    n1z = -10*state.mean_photon(1)[0]
    return n02 + n0n1 + n0z + n12 + n1z + 25

where each of the variables shows respectively the expectation value of each term of the expansion.

Can anyone please verify if this is the correct way? And is there any better way to calculate the expectation values of any operators in Strawberryfields ( like qml.expval(H))?

Any help is much appreciated. Thanks!

Hi @Pranav_Chandarana, thank you for your question! We’re taking a look at it and should be back soon with an answer.

Hi @Pranav_chandarana!

Yes, this looks like the correct expression :+1:. A perfectly suitable way to achieve this.

Strawberry Fields does not have the same level of support for measuring arbitrary Hamiltonians as PennyLane does (much harder to compute arbitrary expvals in photonics picture).

If you had a Hamiltonian that was a linear combination of number operators, you could use poly_quad_expectation, but your Hamiltonian a quadratic polynomial of number operators, so the way you’ve taken would be the recommended way.


Thanks, @nathan! An extended question. I was playing around with QAOA to minimize this cost function using a Photonic setup. The circuit is as follows:

Where U_C is the Hamiltonian term and U_B is a mixer term. The cost function this the expectation value as defined above.

But this optimization is not working as it should. Any intuition why this is the case? I know this question is difficult to answer by looking at the code. Any help will be appreciated.

The optimization energies look like the following:

P.S. Ignore the “probability”.

1 Like

Hi @Pranav_Chandarana, is there a particular reason why you decided to use tapes? These were originally designed as a developer feature so I’m curious about your use of it.

Also, if you can copy-paste your code I can try to run it to see if I can replicate your problem.

Hello! @CatalinaAlbornoz. I use tapes just to find gradients in a TensorFlow setup (I wanted to use Adam optimizer). We can also use scipy optimizers if needed.

As the code is a little long, it is much easy to make a google collab so I have done it.
Here is the link

The energies might not be the same but as we are finding the minimum of the H_p, we should get 0 ( as QAOA should land into any values of n_0 and n_1 that sum to 5. )

I have used p^2 as a mixer and made some changes in the code so that it is easy. Hence, it might now replicate exactly but the question is:

Why the optimization does not converge to 0 if the “expval2()” is a correct function?

Hi @Pranav_Chandarana,

Thank you for sharing your code. I can replicate your results. I’m not sure where the issue is though. We’ll take a look into this.

1 Like

Hello. Any updates on this?

Hi @Pranav_Chandarana,

I’m very sorry for the delay in our response. There are several things that can be happening here:

1 - With your current cutoff you’re not being able to represent 100% of your state, probably due to the fact that you have non-gaussian gates such as Kerr gates. You can check this by using tf.norm(state.ket()) and noticing that the norm isn’t 1. You can try increasing the cutoff, although this is going to take more simulation time and is probably not going to completely solve your problem.
2 - There are no guarantees that the algorithm that you’re using will get you to the global minimum. You could be, for instance, arriving at a local minimum. You can reduce your chances of getting stuck in a local minimum by trying different initializations of your hyperparameters or introducing some randomness.
3 - It is possible that your implementation of QAOA doesn’t actually solve the problem. You can check this paper on A Quantum Approximate Optimization Algorithm for continuous problems, and this implementation of that paper by Jack Ceroni. You will notice that this CV version of QAOA looks very different from the discrete one. Again, there are no guarantees that this will solve your problem but it’s a route worth exploring.

I hope this can help you understand and resolve your problem. Please let us know if you have any further questions :slight_smile: .

1 Like

Hello @CatalinaAlbornoz. Thanks for the detailed explanation. I agree with all the points. In the meantime, I wrote another similar code (not the same) that solves the problem. In this case, instead of finding the expectation value of H_p, I find the mean photon number of modes n_0 and n_1 at each step and put the values in the function (n_0 + n_1 - 5)^2. This way, I was able to solve the problem. But, I don’t know if this is the right way or not.

Essentially, a variational quantum algorithm needs to have an output from which we can decode the solution. This seems to be true in this case since we can find optimal circuit evolution for which the mean photon number information will give the solution.

The problem however is that since the mean photon number is essentially <n_i>, the function that QAOA evaluates at each iteration is not exactly <H_p>. Since <(n_0 + n_1 - 5)^2> \neq (<n_0> + <n_1> - 5)^2> in general. The code can be found here. If you see the last cell, the sum of the mean photon number does equal 5 so we solved the problem.

The question here is: Is this algorithm correct physics-wise? If the answer is no, then why? The algorithm is not QAOA anymore but it can be seen as a variational algorithm that can be used to solve the problem. (since all the information that we need is n_0 and n_1 should sum to 5.)

Many thanks for the help.

Hi @Pranav_Chandarana, I’m glad that you managed to find a solution!

I cannot guarantee the physical correctness of this but it does seem to follow the approach detailed in the Continuous-QAOA implementation that I shared before.
Here the author also uses this version of QAOA (which looks very different from the discrete one) for a parabolic function with a shift: f(x) = (x − a)^2

I hope this helps. Hopefully some other members of the community can validate this too! :slight_smile:

1 Like

Many thanks for the input!!

Hi @Pranav_Chandarana,

I don’t see any issues with the approach you are following. Seems to me like you want to ensure that the sum of the mean photon numbers in both modes is equal to 5, so it’s perfectly fine to use (\langle n_0\rangle) + \langle n_1\rangle - 5)^2 as the cost function.

The mistake would be in the other direction – minimizing \langle (n_0 + n_1 -5)^2 \rangle may not give you a circuit where the sum of mean photon numbers is equal to 5, precisely for the reasons you are outlining.

Hope that helps!


Hello @jmarrazola, Thanks for the input. My goal is to extract the solution from the parameterized unitary.
I wanted to clarify that I am just doing QAOA with number states. So end goal was to get state any one of the possible ground states of the H_p with high probability. Seems like that’s hard to do. So instead of extracting the solution from the end state (usual QAOA), This cost function will extract information from the mean photon number (because for instance, <n_0>,<n_1>= 0,5 will not mean that state |0,5> has high success probability in general). I guess this algorithm can be “sold” as a VQA since the goal of any VQA is to successfully encode and decode the problem with a parameterized unitary, which seems to be the case here.

This will work for “classical” Hamiltonians (def: the solutions are pure number states (analog to classical Hamiltonians in qubits where its bitstrings)). If we go for a quantum Hamiltonian, this cost function might not work (since there will be no information from the correlations).

Sorry for drilling with the details. Could you please confirm if this is true? I can discuss more but I’ve already taken much of everyone’s time. If I take more, I will include everyone as co-authors in this project. :stuck_out_tongue:

Thanks again!

So if your goal is to find the ground state of the Hamiltonian H = (n_0 + n_1 - 5)^2, then it seems to me like your original approach to define its expectation value by expanding the operator is correct.

Why did you think your code was not working? Note that the state that minimizes the expectation value of this Hamiltonian may not be one for which \langle n_0\rangle + \langle n_1\rangle = 5.

My goal is to design a VQA to find optimal solutions for H.

In QAOA, the solution information is encoded in the GS of H. But if I do <n_0> + <n_1>=5, the information of the solution will be encoded in the <n_0> and <n_1> instead of the state. So I don’t need to find the optimal state anymore. So, finding optimal mean photons numbers can be used as an alternative cost function instead of an expectation value to retrieve the solution in a VQA.

Regarding the working part, by working, I mean, it is not reaching the ground state by using full expectation value as a cost function.


There are two algorithms 1) QAOA and 2) let’s call it CAOA ( random name where we minimize the <n_0> + <n_1> values as a cost function. ).
In both algorithms, the circuit ansatz is the same. The only difference is in the cost function that we minimize.
In QAOA, the solution information will be encoded in the GS of the H. So we minimize the expected value of H.
On the other hand, in CAOA, the solution information will not be encoded in the GS but it can be retrieved by minimizing the mean photon numbers of the modes.
QAOA is not reaching the GS, but by using CAOA, we can solve the problem. Just wanted to clarify if I am doing everything right or not in the case of CAOA. In the sense that algorithm-wise there is not any problem.

Hey @Pranav_Chandarana! Just want to check in and see if there are any more concerns or questions here. Is your code working as expected?

Hello @isaacdevlugt. More or less yes. I think I can ask (if I have something) in a separate discussion topic if you want to close this discussion. Thanks!

1 Like

Awesome! Glad we could help!