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

@@ -352,6 +352,9 @@
 
 <h3>Bug fixes 🐛</h3>
 
+* `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)
 

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)

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,15 +394,17 @@ 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]),
  qml.RX(phi, wires=wires[1]),
  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)
 

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)