Speed up GroverOperator and add GlobalPhase to its decomposition

[dwierichs]

Context:

1. The GroverOperator is applied by DefaultQubit by applying its matrix for up to 12 wires (incl).
   However, GroverOperator has a lot of structure as it is 2*P-1 for a projector P, and a custom application rule is faster than constructing and applying the matrix.
2. The compute_matrix method of GroverOperator in addition uses a rather contrived way to construct the matrix np.full((dim, dim), normalization_constant). Instead, we can simply use the broadcasting behaviour of numpy and do normalization_constant - np.eye(dim) (timings below).
3. The GroverOperator decomposition is only correct up to a minus sign (or phase of $\pi$).

Description of the Change:

1. Adds a dispatch registration for apply_operation in /qubit/apply_operation.py for GroverOperator that is not based on the matrix representation of GroverOperator.
2. Modifies the compute_matrix method of GroverOperator to skip repetitive calling of np.kron and np.outer to construct a trivial matrix.
3. Introduces a GlobalPhase($\pi$) gate in the decomposition.

Benefits:
Speedups and correct decomposition.

Possible Drawbacks:
N/A

Related GitHub Issues:

Timing code for the matrix construction:

def f(n_wires):
    s1 = np.array([1, 1]) / np.sqrt(2)
    # uniform superposition state |s>
    s = reduce(np.kron, list(it.repeat(s1, n_wires)))
    # Grover diffusion operator
    return 2 * np.outer(s, s) - np.identity(2**n_wires)

for N in [3, 6, 9, 12]:
    dim = 2**N

    print(f"N={N}")
    %timeit f(N)
    %timeit 2 / dim - np.eye(dim)
    print()

Timings for the different apply_operation options, for three scenarios: Number of wires for operation and state is equal, number of wires for state is that for operation +4, and number of wires for state is double of that for operation.

Given the results, we could do einsum for len(op.wires)<8 and the dispatch method for bigger wire counts.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (17503d3) 99.65% compared to head (cd34615) 99.64%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #4666      +/-   ##
==========================================
- Coverage   99.65%   99.64%   -0.01%     
==========================================
  Files         383      383              
  Lines       34559    34311     -248     
==========================================
- Hits        34439    34190     -249     
- Misses        120      121       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[albi3ro]

You may need to tweak "use matrix versus decomposition" threshold in default qubit as well.

It is currently located here:

pennylane/pennylane/devices/qubit/preprocess.py

Line 71 in ea7e1d4

if op.name == "GroverOperator" and len(op.wires) >= 13:

But I have an open PR that will be moving it. We can resolve that merge conflict when we get to it.

[dwierichs]

You may need to tweak "use matrix versus decomposition" threshold in default qubit as well.

Good call. I actually removed it now, because it is supposed to address the large wire counts, which is the regime where the dispatched apply_grover is particularly fast. Moreover, as the decomposition also creates a dense matrix (via MultiControlledX), this PR makes GroverOperator accessible for qubit counts (>=15 on my laptop) that were inaccessible before :)

But I have an open PR that will be moving it. We can resolve that merge conflict when we get to it.

👍

[dwierichs]

Updated timings for all interfaces (line style gives interface, color gives target wires): As there is an outlier in Tensorflow which leads to num_wires=16, num_op_wires=8 having a factor 2 slowdown, still, I will move the threshold for falling back to matrix-based method (einsum/tensordot) from num_op_wires<8 to num_op_wires<9.

[timmysilv]

pretty much good to go, just a few questions

[timmysilv]

great optimization, thanks for putting this all together! 🎉