Make Pauli decompose coefficient real again (#4479)

* new `phased_pauli_decompose`

* make `black` happy

* make `pylint` happy

* rename to `pauli_decompose_with_phase`

* improve docstrings

* use `qml.math` in Hermitian check

* remove `pauli_decompose_with_phase`

* code cleanup

* improve docstrings from review

* more from review

doc/releases/changelog-dev.md

@@ -182,8 +182,9 @@ array([False, False])
 * When given a callable, `qml.ctrl` now does its custom pre-processing on all queued operators from the callable.
  [(#4370)](https://github.com/PennyLaneAI/pennylane/pull/4370)
 
-* `qml.pauli_decompose` is now differentiable and works with any non-Hermitian and non-square matrices.
+* `qml.pauli_decompose` is now exponentially faster and differentiable.
  [(#4395)](https://github.com/PennyLaneAI/pennylane/pull/4395)
+ [(#4479)](https://github.com/PennyLaneAI/pennylane/pull/4479)
 
 * `qml.interfaces.set_shots` accepts `Shots` object as well as `int`'s and tuples of `int`'s.
  [(#4388)](https://github.com/PennyLaneAI/pennylane/pull/4388)

pennylane/pauli/conversion.py

@@ -17,7 +17,7 @@
 from functools import reduce, singledispatch
 from itertools import product
 from operator import matmul
-from typing import Union
+from typing import Union, Tuple
 
 import pennylane as qml
 from pennylane.operation import Tensor
@@ -29,14 +29,10 @@
 
 
 # pylint: disable=too-many-branches
-def pauli_decompose(
- matrix,
- hide_identity=False,
- wire_order=None,
- pauli=False,
- padding=False,
-) -> Union[Hamiltonian, PauliSentence]:
- r"""Decomposes a matrix into a linear combination of Pauli operators.
+def _generalized_pauli_decompose(
+ matrix, hide_identity=False, wire_order=None, pauli=False, padding=False
+) -> Tuple[qml.typing.TensorLike, list]:
+ r"""Decomposes any matrix into a linear combination of Pauli operators.
 
  This method converts any matrix to a weighted sum of Pauli words acting on :math:`n` qubits
  in time :math:`O(n 4^n)`. The input matrix is first padded with zeros if its dimensions are not
@@ -59,83 +55,75 @@ def pauli_decompose(
  that are not of shape :math:`2^n\times 2^n` by padding them with zeros if ``True``.
 
  Returns:
- Union[~.Hamiltonian, ~.PauliSentence]: the matrix decomposed as a linear combination
- of Pauli operators, either as a :class:`~.Hamiltonian` or :class:`~.PauliSentence` instance.
+ Tuple[qml.math.array[complex], list]: the matrix decomposed as a linear combination of Pauli operators
+ as a tuple consisting of an array of complex coefficients and a list of corresponding Pauli terms.
 
  **Example:**
 
  We can use this function to compute the Pauli operator decomposition of an arbitrary matrix:
 
  >>> A = np.array(
- ... [[-2, -2+1j, -2, -2], [-2-1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]])
- >>> H = qml.pauli_decompose(A)
- >>> print(H)
- (-1.5) [I0 X1]
- + (-1.5) [X0 I1]
- + (-1.0) [I0 I1]
- + (-1.0) [I0 Z1]
- + (-1.0) [X0 X1]
- + (-0.5) [I0 Y1]
- + (-0.5) [X0 Z1]
- + (-0.5) [Z0 X1]
- + (-0.5) [Z0 Y1]
- + (1.0) [Y0 Y1]
-
- We can return a :class:`~.PauliSentence` instance by using the keyword argument ``pauli=True``:
-
- >>> ps = qml.pauli_decompose(A, pauli=True)
- >>> print(ps)
- -1.0 * I
- + -1.5 * X(1)
- + -0.5 * Y(1)
- + -1.0 * Z(1)
- + -1.5 * X(0)
- + -1.0 * X(0) @ X(1)
- + -0.5 * X(0) @ Z(1)
- + 1.0 * Y(0) @ Y(1)
- + -0.5 * Z(0) @ X(1)
- + -0.5 * Z(0) @ Y(1)
+ ... [[-2, -2+1j, -2, -2], [-2-1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 1j]])
+ >>> coeffs, obs = qml.pauli.conversion._generalized_pauli_decompose(A)
+ >>> coeffs
+ array([-1. +0.25j, -1.5+0.j , -0.5+0.j , -1. -0.25j, -1.5+0.j ,
+ -1. +0.j , -0.5+0.j , 1. -0.j , 0. -0.25j, -0.5+0.j ,
+ -0.5+0.j , 0. +0.25j])
+ >>> obs
+ [Identity(wires=[0]) @ Identity(wires=[1]),
+ Identity(wires=[0]) @ PauliX(wires=[1]),
+ Identity(wires=[0]) @ PauliY(wires=[1]),
+ Identity(wires=[0]) @ PauliZ(wires=[1]),
+ PauliX(wires=[0]) @ Identity(wires=[1]),
+ PauliX(wires=[0]) @ PauliX(wires=[1]),
+ PauliX(wires=[0]) @ PauliZ(wires=[1]),
+ PauliY(wires=[0]) @ PauliY(wires=[1]),
+ PauliZ(wires=[0]) @ Identity(wires=[1]),
+ PauliZ(wires=[0]) @ PauliX(wires=[1]),
+ PauliZ(wires=[0]) @ PauliY(wires=[1]),
+ PauliZ(wires=[0]) @ PauliZ(wires=[1])]
 
  We can also set custom wires using the ``wire_order`` argument:
 
- >>> ps = qml.pauli_decompose(A, pauli=True, wire_order=['a', 'b'])
- >>> print(ps)
- -1.0 * I
- + -1.5 * X(b)
- + -0.5 * Y(b)
- + -1.0 * Z(b)
- + -1.5 * X(a)
- + -1.0 * X(a) @ X(b)
- + -0.5 * X(a) @ Z(b)
- + 1.0 * Y(a) @ Y(b)
- + -0.5 * Z(a) @ X(b)
- + -0.5 * Z(a) @ Y(b)
+ >>> coeffs, obs = qml.pauli.conversion._generalized_pauli_decompose(A, wire_order=['a', 'b'])
+ >>> obs
+ [Identity(wires=['a']) @ Identity(wires=['b']),
+ Identity(wires=['a']) @ PauliX(wires=['b']),
+ Identity(wires=['a']) @ PauliY(wires=['b']),
+ Identity(wires=['a']) @ PauliZ(wires=['b']),
+ PauliX(wires=['a']) @ Identity(wires=['b']),
+ PauliX(wires=['a']) @ PauliX(wires=['b']),
+ PauliX(wires=['a']) @ PauliZ(wires=['b']),
+ PauliY(wires=['a']) @ PauliY(wires=['b']),
+ PauliZ(wires=['a']) @ Identity(wires=['b']),
+ PauliZ(wires=['a']) @ PauliX(wires=['b']),
+ PauliZ(wires=['a']) @ PauliY(wires=['b']),
+ PauliZ(wires=['a']) @ PauliZ(wires=['b'])]
 
  .. details::
- :title: Usage Details
+ :title: Advanced Usage Details
  :href: usage-decompose-operation
 
  For non-square matrices, we need to provide the ``padding=True`` keyword argument:
 
  >>> A = np.array([[-2, -2 + 1j]])
- >>> H = qml.pauli_decompose(A, padding=True)
- >>> print(H)
- ((-1+0j)) [I0]
- + ((-1+0.5j)) [X0]
- + ((-1+0j)) [Z0]
- + ((-0.5-1j)) [Y0]
+ >>> coeffs, obs = qml.pauli.conversion._generalized_pauli_decompose(A, padding=True)
+ >>> coeffs
+ ([-1. +0.j , -1. +0.5j, -0.5-1.j , -1. +0.j ])
+ >>> obs
+ [Identity(wires=[0]), PauliX(wires=[0]), PauliY(wires=[0]), PauliZ(wires=[0])]
 
  We can also use the method within a differentiable workflow and obtain gradients:
 
  >>> A = qml.numpy.array([[-2, -2 + 1j]], requires_grad=True)
  >>> dev = qml.device("default.qubit", wires=1)
  >>> @qml.qnode(dev)
  ... def circuit(A):
- ... decomp = qml.pauli_decompose(A, padding=True)
- ... qml.RX(decomp.coeffs[2], 0)
+ ... coeffs, _ = qml.pauli.conversion._generalized_pauli_decompose(A, padding=True)
+ ... qml.RX(qml.math.real(coeffs[2]), 0)
  ... return qml.expval(qml.PauliZ(0))
- >>> grad_numpy = qml.grad(circuit)(A)
- tensor([[-2.+0.j, -2.+1.j]], requires_grad=True)
+ >>> qml.grad(circuit)(A)
+ array([[0.+0.j , 0.+0.23971277j]])
 
  """
  # Ensuring original matrix is not manipulated and we support builtin types.
@@ -220,7 +208,114 @@ def pauli_decompose(
 
  coeffs = qml.math.stack(coeffs)
 
- # Convert to Hamiltonian and PauliSentence
+ if not pauli:
+ with qml.QueuingManager.stop_recording():
+ obs = [reduce(matmul, [op_map[o](w) for o, w in obs_term]) for obs_term in obs]
+
+ return (coeffs, obs)
+
+
+def pauli_decompose(
+ H, hide_identity=False, wire_order=None, pauli=False, check_hermitian=True
+) -> Union[Hamiltonian, PauliSentence]:
+ r"""Decomposes a Hermitian matrix into a linear combination of Pauli operators.
+
+ Args:
+ H (array[complex]): a Hermitian matrix of dimension :math:`2^n\times 2^n`.
+ hide_identity (bool): does not include the Identity observable within
+ the tensor products of the decomposition if ``True``.
+ wire_order (list[Union[int, str]]): the ordered list of wires with respect
+ to which the operator is represented as a matrix.
+ pauli (bool): return a PauliSentence instance if ``True``.
+ check_hermitian (bool): check if the provided matrix is Hermitian if ``True``.
+
+ Returns:
+ Union[~.Hamiltonian, ~.PauliSentence]: the matrix decomposed as a linear combination
+ of Pauli operators, returned either as a :class:`~.Hamiltonian` or :class:`~.PauliSentence`
+ instance.
+
+ **Example:**
+
+ We can use this function to compute the Pauli operator decomposition of an arbitrary Hermitian
+ matrix:
+
+ >>> A = np.array(
+ ... [[-2, -2+1j, -2, -2], [-2-1j, 0, 0, -1], [-2, 0, -2, -1], [-2, -1, -1, 0]])
+ >>> H = qml.pauli_decompose(A)
+ >>> print(H)
+ (-1.5) [I0 X1]
+ + (-1.5) [X0 I1]
+ + (-1.0) [I0 I1]
+ + (-1.0) [I0 Z1]
+ + (-1.0) [X0 X1]
+ + (-0.5) [I0 Y1]
+ + (-0.5) [X0 Z1]
+ + (-0.5) [Z0 X1]
+ + (-0.5) [Z0 Y1]
+ + (1.0) [Y0 Y1]
+
+ We can return a :class:`~.PauliSentence` instance by using the keyword argument ``pauli=True``:
+
+ >>> ps = qml.pauli_decompose(A, pauli=True)
+ >>> print(ps)
+ -1.0 * I
+ + -1.5 * X(1)
+ + -0.5 * Y(1)
+ + -1.0 * Z(1)
+ + -1.5 * X(0)
+ + -1.0 * X(0) @ X(1)
+ + -0.5 * X(0) @ Z(1)
+ + 1.0 * Y(0) @ Y(1)
+ + -0.5 * Z(0) @ X(1)
+ + -0.5 * Z(0) @ Y(1)
+
+ By default the wires are numbered [0, 1, ..., n], but we can also set custom wires using the
+ ``wire_order`` argument:
+
+ >>> ps = qml.pauli_decompose(A, pauli=True, wire_order=['a', 'b'])
+ >>> print(ps)
+ -1.0 * I
+ + -1.5 * X(b)
+ + -0.5 * Y(b)
+ + -1.0 * Z(b)
+ + -1.5 * X(a)
+ + -1.0 * X(a) @ X(b)
+ + -0.5 * X(a) @ Z(b)
+ + 1.0 * Y(a) @ Y(b)
+ + -0.5 * Z(a) @ X(b)
+ + -0.5 * Z(a) @ Y(b)
+
+ .. details::
+ :title: Theory
+ :href: theory
+
+ This method internally uses a generalized decomposition routine to convert the matrix to a
+ weighted sum of Pauli words acting on :math:`n` qubits in time :math:`O(n 4^n)`. The input
+ matrix is written as a quantum state in the computational basis following the
+ `channel-state duality <https://en.wikipedia.org/wiki/Channel-state_duality>`_.
+ A Bell basis transformation is then performed using the
+ `Walsh-Hadamard transform <https://en.wikipedia.org/wiki/Hadamard_transform>`_, after which
+ coefficients for each of the :math:`4^n` Pauli words are computed while accounting for the
+ phase from each ``PauliY`` term occurring in the word.
+
+ """
+ n = int(qml.math.log2(qml.math.shape(H)[0]))
+ N = 2**n
+
+ if check_hermitian:
+ if H.shape != (N, N):
+ raise ValueError("The matrix should have shape (2**n, 2**n), for any qubit number n>=1")
+
+ if not qml.math.allclose(H, qml.math.conjugate(qml.math.transpose(H))):
+ raise ValueError("The matrix is not Hermitian")
+
+ coeffs, obs = _generalized_pauli_decompose(
+ H, hide_identity=hide_identity, wire_order=wire_order, pauli=pauli, padding=True
+ )
+
+ if check_hermitian:
+ coeffs = qml.math.real(coeffs)
+
  if pauli:
  return PauliSentence(
  {
@@ -229,7 +324,6 @@ def pauli_decompose(
  }
  )
 
- obs = [reduce(matmul, [op_map[o](w) for o, w in obs_term]) for obs_term in obs]
  return Hamiltonian(coeffs, obs)