Counting in Q#: number of solutions

Question

I have this program derived from Microsoft Quantum Katas for counting (see here):

namespace Quantum.Sample
{
    open Microsoft.Quantum.Primitive;
    open Microsoft.Quantum.Canon;
    open Microsoft.Quantum.Extensions.Convert;
    open Microsoft.Quantum.Extensions.Math;

    operation SprinklerAnc (queryRegister:  Qubit[],  target : Qubit,ancilla: Qubit[]) : Unit {
        body (...) {
            X(queryRegister[2]);
            X(queryRegister[3]);
            X(queryRegister[4]);
            X(queryRegister[5]);
            X(queryRegister[6]);
            X(ancilla[0]);
            X(ancilla[1]);
            X(ancilla[2]);

            CCNOT(queryRegister[0],queryRegister[1],ancilla[0]);
            CCNOT(queryRegister[1],queryRegister[2],ancilla[1]);
            CCNOT(queryRegister[0],queryRegister[2],ancilla[2]);
            (Controlled X)([ancilla[0],ancilla[1],ancilla[2],queryRegister[3],queryRegister[4],queryRegister[5],queryRegister[6]],target);
            CCNOT(queryRegister[0],queryRegister[2],ancilla[2]);
            CCNOT(queryRegister[1],queryRegister[2],ancilla[1]);
            CCNOT(queryRegister[0],queryRegister[1],ancilla[0]);

            X(ancilla[2]);
            X(ancilla[1]);
            X(ancilla[0]);
                        X(queryRegister[2]);
            X(queryRegister[6]);
            X(queryRegister[5]);
            X(queryRegister[4]);
            X(queryRegister[3]);


        }
        adjoint invert;
        controlled auto;
        controlled adjoint auto;
    }           
    operation OracleConverterImpl (markingOracle : ((Qubit[], Qubit) => Unit : Adjoint, Controlled), register : Qubit[]) : Unit {

        body (...) {
            using (target = Qubit()) {
                // Put the target into the |-⟩ state
                X(target);
                H(target);

                // Apply the marking oracle; since the target is in the |-⟩ state,
                // flipping the target if the register satisfies the oracle condition will apply a -1 factor to the state
                markingOracle(register, target);

                // Put the target back into |0⟩ so we can return it
                H(target);
                X(target);
            }
        }

        adjoint invert;
        controlled auto;
        adjoint controlled auto;
    }

    operation HadamardTransform (register : Qubit[]) : Unit {

        body (...) {
            //ApplyToEachA(H, register);

            // ApplyToEach is a library routine that is equivalent to the following code:
             let nQubits = Length(register);
             for (idxQubit in 0..nQubits - 1) {
                 H(register[idxQubit]);
             }
        }

        adjoint invert;
        controlled auto;
        controlled adjoint auto;
    }
    operation Oracle_ArbitraryPattern (queryRegister : Qubit[], target : Qubit, pattern : Bool[]) : Unit {

        body (...) {
            (ControlledOnBitString(pattern, X))(queryRegister, target);
        }

        adjoint invert;
        controlled auto;
        controlled adjoint auto;

    }   

    // Task 2.2. Conditional phase flip
    operation ConditionalPhaseFlip (register : Qubit[]) : Unit {

        body (...) {
            // Define a marking oracle which detects an all zero state
            let allZerosOracle = Oracle_ArbitraryPattern(_, _, new Bool[Length(register)]);

            // Convert it into a phase-flip oracle and apply it
            let flipOracle = OracleConverter(allZerosOracle);
            flipOracle(register);
        }

        adjoint self;
        controlled  auto;

        controlled adjoint auto;

    }



    // Task 2.3. The Grover iteration
    operation GroverIteration (register : Qubit[], oracle : (Qubit[] => Unit :   Adjoint, Controlled)) : Unit {

        body (...) {
            oracle(register);
            HadamardTransform(register);
            ConditionalPhaseFlip(register);
            HadamardTransform(register);
        }

        adjoint invert;
              controlled auto;
        controlled adjoint auto;
    }    
    function OracleConverter (markingOracle : ((Qubit[], Qubit) => Unit : Adjoint, Controlled)) : (Qubit[] => Unit : Adjoint, Controlled) {
        return OracleConverterImpl(markingOracle, _);
    }

    operation UnitaryPowerImpl (U : (Qubit[] => Unit : Adjoint, Controlled), power : Int, q : Qubit[]) : Unit {
        body (...) {
            for (i in 1..power) {
                U(q);
            }
        }
        adjoint auto;
        controlled auto;
        controlled adjoint auto;
    }

    operation QPE() : Double {
        mutable phase = -1.0;
        let n=8;
                       using ((reg,phaseRegister,ancilla)=(Qubit[7 ], Qubit[n],Qubit[3]))
                       {
        // Construct a phase estimation oracle from the unitary
                         let phaseOracle = OracleConverter(SprinklerAnc(_,_,ancilla));

        let oracle = DiscreteOracle(UnitaryPowerImpl(GroverIteration(_, phaseOracle), _, _));

              //  let markingOracle = Sprinkler(_, _);

        // Allocate qubits to hold the eigenstate of U and the phase in a big endian register 

            let phaseRegisterBE = BigEndian(phaseRegister);
            // Prepare the eigenstate of U
                HadamardTransform(reg);

            // Call library
            QuantumPhaseEstimation(oracle, reg, phaseRegisterBE);
            // Read out the phase
            set phase = ToDouble(MeasureIntegerBE(phaseRegisterBE)) / ToDouble(1 <<< (n));

            ResetAll(reg);
            ResetAll(phaseRegister);
        }
        let angle = PI()*phase;
        let res = 128.0 *(1.0- PowD(Sin(angle),2.0));

        return res;
    }



}

with driver (see this):

using System;

using Microsoft.Quantum.Simulation.Core;
using Microsoft.Quantum.Simulation.Simulators;

namespace Quantum.Sample
{
    class Driver
    {
        static void Main(string[] args)
        {
            using (var qsim = new QuantumSimulator())
            {
                for (int i = 0; i < 100; i++)
                {
                    var res = QPE.Run(qsim).Result;
                    System.Console.WriteLine($"Res:{res}");
                }
            }
        }
    }
}

It should count the solutions of the 3 bits formula (not reg[0] or reg[2]) and (not reg[1] or reg[2]) and (not reg[0] or not reg[1]) This formula has 4 solutions. I consider an extended formula with 7 bits that is true when the first three bit satisfy the above formula and the other four bits are at 0. The extended formula this has 4 solutions as well but I can increase easily the number of solutions by excluding variables from the formula. I compute the phase $\phi$ with the linked program, then I compute $\theta/2$ as $2\pi\phi/2=\pi\phi$ Now the number of solutions should be $$128\sin(\theta/2)^2$$ right? I get the correct number of solutions if I compute $$128(1-\sin(\theta/2)^2)$$ It's as if the formula is negated but I can't find where the negation is. This same formula with GroverSearch from Microsoft Quantum Katas returns correct solutions.

Answer 1

I have a theory as for where the issue comes from (huge thanks to Robin for helping me figure it out!)

Grover iteration consists of four steps:

• Apply the oracle.
• Apply the Hadamard transform.
• Perform a conditional phase shift.
• Apply the Hadamard transform.

ConditionalPhaseFlip operation in the Q# code implements the third step: it gives a phase shift of -1 to the $|0\rangle$ state (if you follow Q# code closely, this routine takes a marking oracle which marks the $|0\rangle$ state and converts it to phase-flipping oracle which flips the phase of this state only).

This implementation differs from the description in Nielsen and Chuang, which says that each basis state except $|0\rangle$ gets a phase shift of -1 (the difference is global phase of -1). When you use this Grover iteration as part of Grover search algorithm, this extra global phase has no effect on measurement outcomes.

However, when you use this iteration as part of quantum counting algorithm, this global phase can actually be detected! Since phase estimation uses controlled versions of the iteration, the phase becomes relative instead of global, and starts affecting the computation.

Let's denote the Grover iteration as $G$ with eigenvalues $e^{i\theta}$ and $e^{i (2 \pi - \theta)}$ (with the correct number of solutions calculated as $M_c = N \sin^2 \frac{\theta}{2}$). The Grover iteration with an global phase of $-1$ $(-G)$ will have eigenvalues multiplied by $-1$: $-e^{i\theta}$ and $-e^{i (2 \pi - \theta)}$. Given that $e^{i\pi} = -1$, we can write them as $e^{i(\pi + \theta)}$ and $e^{i (\pi - \theta)}$. When you feed the unitary $(-G)$ into the quantum counting algorithm, you'll get $M_{-1} = N \sin^2 \frac{\theta + \pi}{2} = N \cos^2 \frac{\theta}{2}$ - which is exactly the answer you're getting!

The fix is a lot easier than figuring out what went wrong - just add a global phase of -1 to ConditionalPhaseFlip. You can use R operation for that: R(PauliI, 2.0 * PI(), register[0]);.