From a04d1a29faec52cba28b824271ba2bded22c9d7e Mon Sep 17 00:00:00 2001 From: Christina Lee Date: Fri, 12 Apr 2024 18:46:38 -0400 Subject: [PATCH] fix discontuity at special case (#5448) **Context:** We were seeing an incorrect two-qubit-decomposition at one very particular value of matrix. **Description of the Change:** Adds a very small perturbation to one of the eigenvalues. **Benefits:** Decomposition matrix now matches input matrix. **Possible Drawbacks:** Still don't know why the discontinuity happens. It may now be happening in a slightly different place. **Related GitHub Issues:** Fixes #5308 [sc-58004] --------- Co-authored-by: Thomas R. Bromley <49409390+trbromley@users.noreply.github.com> --- doc/releases/changelog-dev.md | 3 ++ pennylane/math/single_dispatch.py | 7 ++++ .../decompositions/two_qubit_unitary.py | 8 +++-- tests/ops/op_math/test_decompositions.py | 34 +++++++++++++++++++ 4 files changed, 50 insertions(+), 2 deletions(-) diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index b4456e59708..4b450f41684 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -352,6 +352,9 @@

Bug fixes 🐛

+* `two_qubit_decomposition` no longer diverges at a special case of unitary matrix. + [(#5448)](https://github.com/PennyLaneAI/pennylane/pull/5448) + * The `qml.QNSPSAOptimizer` now correctly handles optimization for legacy devices that do not follow the new API design. [(#5497)](https://github.com/PennyLaneAI/pennylane/pull/5497) diff --git a/pennylane/math/single_dispatch.py b/pennylane/math/single_dispatch.py index 1999ce69edc..f1ed70589a3 100644 --- a/pennylane/math/single_dispatch.py +++ b/pennylane/math/single_dispatch.py @@ -265,6 +265,13 @@ def _take_autograd(tensor, indices, axis=None): ar.autoray._FUNC_ALIASES["tensorflow", "diag"] = "diag" +ar.register_function( + "tensorflow", + "finfo", + lambda *args, **kwargs: _i("tf").experimental.numpy.finfo(*args, **kwargs), +) + + def _tf_convert_to_tensor(x, **kwargs): if isinstance(x, _i("tf").Tensor) and "dtype" in kwargs: return _i("tf").cast(x, **kwargs) diff --git a/pennylane/ops/op_math/decompositions/two_qubit_unitary.py b/pennylane/ops/op_math/decompositions/two_qubit_unitary.py index 3f72ff49bbc..9baf4b7da0e 100644 --- a/pennylane/ops/op_math/decompositions/two_qubit_unitary.py +++ b/pennylane/ops/op_math/decompositions/two_qubit_unitary.py @@ -347,6 +347,7 @@ def _decomposition_2_cnots(U, wires): part has the same spectrum as U, and then we can recover A, B, C, D. """ # Compute the rotation angles + u = math.dot(Edag, math.dot(U, E)) gammaU = math.dot(u, math.T(u)) evs, _ = math.linalg.eig(gammaU) @@ -383,6 +384,7 @@ def _decomposition_2_cnots(U, wires): else: # For the non-special case, the eigenvalues come in conjugate pairs. # We need to find two non-conjugate eigenvalues to extract the angles. + x = math.angle(evs[0]) y = math.angle(evs[1]) @@ -392,7 +394,6 @@ def _decomposition_2_cnots(U, wires): delta = (x + y) / 2 phi = (x - y) / 2 - interior_decomp = [ qml.CNOT(wires=[wires[1], wires[0]]), qml.RZ(delta, wires=wires[0]), @@ -400,7 +401,10 @@ def _decomposition_2_cnots(U, wires): qml.CNOT(wires=[wires[1], wires[0]]), ] - RZd = qml.RZ(math.cast_like(delta, 1j), wires=0).matrix() + # need to perturb x by 5 precision to avoid a discontinuity at a special case. + # see https://github.com/PennyLaneAI/pennylane/issues/5308 + precision = qml.math.finfo(delta.dtype).eps + RZd = qml.RZ(math.cast_like(delta + 5 * precision, 1j), wires=0).matrix() RXp = qml.RX(phi, wires=0).matrix() inner_matrix = math.kron(RZd, RXp) diff --git a/tests/ops/op_math/test_decompositions.py b/tests/ops/op_math/test_decompositions.py index 83a27d508be..6f5e7e8bd0b 100644 --- a/tests/ops/op_math/test_decompositions.py +++ b/tests/ops/op_math/test_decompositions.py @@ -1230,3 +1230,37 @@ def wrapped_decomposition(U): jitted_matrix = jax.jit(wrapped_decomposition)(U) assert check_matrix_equivalence(U, jitted_matrix, atol=1e-7) + + +def test_two_qubit_decomposition_special_case_discontinuity(): + """Test that two_qubit_decomposition still provides accurate numbers at a special case.""" + + def make_unitary(theta1): + generator = ( + theta1 + / 2 + * ( + np.cos(0.2) + / 2 + * (np.array([[0, 0, 0, 0], [0, 0, 2, 0], [0, 2, 0, 0], [0, 0, 0, 0]])) + + np.sin(0.2) + / 2 + * (np.array([[0, 0, 0, 0], [0, 0, 2j, 0], [0, -2j, 0, 0], [0, 0, 0, 0]])) + ) + ) + + def expm(val): + d, U = np.linalg.eigh(-1.0j * val) + return np.dot(U, np.dot(np.diag(np.exp(1.0j * d)), np.conj(U).T)) + + mat = expm(-1j * generator) + + assert np.allclose( + np.dot(np.transpose(np.conj(mat)), mat), np.eye(len(mat)) + ), "mat is not unitary" + + return mat + + mat = make_unitary(np.pi / 2) + decomp_mat = qml.matrix(two_qubit_decomposition, wire_order=(0, 1))(mat, wires=(0, 1)) + assert qml.math.allclose(mat, decomp_mat)