diff --git a/doc/development/deprecations.rst b/doc/development/deprecations.rst index f869302f449..b4d29bd1049 100644 --- a/doc/development/deprecations.rst +++ b/doc/development/deprecations.rst @@ -11,6 +11,11 @@ Pending deprecations - Still accessible in v0.31 - Removed in v0.32 +* ``qml.math.reduced_dm`` has been deprecated. Please use ``qml.math.reduce_dm`` or ``qml.math.reduce_statevector`` instead. + + - Still accessible in v0.31 + - Removed in v0.32 + * The ``qml.specs`` dictionary will no longer support direct key access to certain keys. Instead these quantities can be accessed as fields of the new ``Resources`` object saved under ``specs_dict["resources"]``: diff --git a/doc/releases/changelog-dev.md b/doc/releases/changelog-dev.md index 6bdcbfe4519..18741780b9a 100644 --- a/doc/releases/changelog-dev.md +++ b/doc/releases/changelog-dev.md @@ -142,6 +142,9 @@ and now inherits from `qml.ops.op_math.ControlledOp`. [(#4116)](https://github.com/PennyLaneAI/pennylane/pull/4116/) +* Added `qml.math.reduce_dm` and `qml.math.reduce_statevector` to produce reduced density matrices. + Both functions have broadcasting support. + [(#4173)](https://github.com/PennyLaneAI/pennylane/pull/4173)

Breaking changes 💔

@@ -177,6 +180,9 @@ * `qml.grouping` module is removed. The functionality has been reorganized in the `qml.pauli` module. +* `qml.math.reduced_dm` has been deprecated. Please use `qml.math.reduce_dm` or `qml.math.reduce_statevector` instead. + [(#4173)](https://github.com/PennyLaneAI/pennylane/pull/4173) +

Documentation 📝

* The description of `mult` in the `qchem.Molecule` docstring now correctly states the value diff --git a/pennylane/_qubit_device.py b/pennylane/_qubit_device.py index 03fb02597be..91c08f1bd48 100644 --- a/pennylane/_qubit_device.py +++ b/pennylane/_qubit_device.py @@ -1281,7 +1281,7 @@ def density_matrix(self, wires): """ state = getattr(self, "state", None) wires = self.map_wires(wires) - return qml.math.reduced_dm(state, indices=wires, c_dtype=self.C_DTYPE) + return qml.math.reduce_statevector(state, indices=wires, c_dtype=self.C_DTYPE) def vn_entropy(self, wires, log_base): r"""Returns the Von Neumann entropy prior to measurement. diff --git a/pennylane/devices/default_mixed.py b/pennylane/devices/default_mixed.py index cd61c4e6b79..018b30076bd 100644 --- a/pennylane/devices/default_mixed.py +++ b/pennylane/devices/default_mixed.py @@ -241,6 +241,20 @@ def state(self): # User obtains state as a matrix return qnp.reshape(self._pre_rotated_state, (dim, dim)) + def density_matrix(self, wires): + """Returns the reduced density matrix over the given wires. + + Args: + wires (Wires): wires of the reduced system + + Returns: + array[complex]: complex array of shape ``(2 ** len(wires), 2 ** len(wires))`` + representing the reduced density matrix of the state prior to measurement. + """ + state = getattr(self, "state", None) + wires = self.map_wires(wires) + return qml.math.reduce_dm(state, indices=wires, c_dtype=self.C_DTYPE) + def reset(self): """Resets the device""" super().reset() diff --git a/pennylane/math/__init__.py b/pennylane/math/__init__.py index 4cfbadb7918..51d6dae65ce 100644 --- a/pennylane/math/__init__.py +++ b/pennylane/math/__init__.py @@ -71,6 +71,8 @@ mutual_info, purity, reduced_dm, + reduce_dm, + reduce_statevector, relative_entropy, sqrt_matrix, vn_entropy, diff --git a/pennylane/math/quantum.py b/pennylane/math/quantum.py index 4bd20da4611..e3e3faa7f90 100644 --- a/pennylane/math/quantum.py +++ b/pennylane/math/quantum.py @@ -15,6 +15,7 @@ # pylint: disable=import-outside-toplevel import itertools import functools +import warnings from string import ascii_letters as ABC from autoray import numpy as np @@ -172,129 +173,157 @@ def marginal_prob(prob, axis): return np.flatten(prob) -def _density_matrix_from_matrix(density_matrix, indices, check_state=False): +def reduce_dm(density_matrix, indices, check_state=False, c_dtype="complex128"): """Compute the density matrix from a state represented with a density matrix. - Args: - density_matrix (tensor_like): 2D density matrix tensor. This tensor should be of size ``(2**N, 2**N)`` for some - integer number of wires``N``. + density_matrix (tensor_like): 2D or 3D density matrix tensor. This tensor should be of size ``(2**N, 2**N)`` or + ``(batch_dim, 2**N, 2**N)``, for some integer number of wires``N``. indices (list(int)): List of indices in the considered subsystem. check_state (bool): If True, the function will check the state validity (shape and norm). + c_dtype (str): Complex floating point precision type. Returns: - tensor_like: Density matrix of size ``(2**len(wires), 2**len(wires))`` + tensor_like: Density matrix of size ``(2**len(indices), 2**len(indices))`` or ``(batch_dim, 2**len(indices), 2**len(indices))`` + + .. seealso:: :func:`pennylane.math.reduce_statevector`, :func:`pennylane.qinfo.transforms.reduced_dm`, + and :func:`pennylane.density_matrix` **Example** >>> x = np.array([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]) - >>> _density_matrix_from_matrix(x, indices=[0]) + >>> reduce_dm(x, indices=[0]) [[1.+0.j 0.+0.j] [0.+0.j 0.+0.j]] >>> y = [[0.5, 0, 0.5, 0], [0, 0, 0, 0], [0.5, 0, 0.5, 0], [0, 0, 0, 0]] - >>> _density_matrix_from_matrix(y, indices=[0]) + >>> reduce_dm(y, indices=[0]) [[0.5+0.j 0.5+0.j] [0.5+0.j 0.5+0.j]] - >>> _density_matrix_from_matrix(y, indices=[1]) + >>> reduce_dm(y, indices=[1]) [[1.+0.j 0.+0.j] [0.+0.j 0.+0.j]] >>> z = tf.Variable([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=tf.complex128) - >>> _density_matrix_from_matrix(x, indices=[1]) + >>> reduce_dm(x, indices=[1]) tf.Tensor( [[1.+0.j 0.+0.j] [0.+0.j 0.+0.j]], shape=(2, 2), dtype=complex128) + >>> x = np.array([[[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], + ... [[0, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]]) + >>> reduce_dm(x, indices=[1]) + array([[[1.+0.j, 0.+0.j], + [0.+0.j, 0.+0.j]], + [[0.+0.j, 0.+0.j], + [0.+0.j, 1.+0.j]]]) """ - dim = density_matrix.shape[0] - num_indices = int(np.log2(dim)) + density_matrix = cast(density_matrix, dtype=c_dtype) if check_state: _check_density_matrix(density_matrix) + if len(np.shape(density_matrix)) == 2: + batch_dim, dim = None, density_matrix.shape[0] + else: + batch_dim, dim = density_matrix.shape[:2] + + num_indices = int(np.log2(dim)) consecutive_indices = list(range(num_indices)) # Return the full density matrix if all the wires are given, potentially permuted if len(indices) == num_indices: - return _permute_dense_matrix(density_matrix, consecutive_indices, indices, None) + return _permute_dense_matrix(density_matrix, consecutive_indices, indices, batch_dim) + + if batch_dim is None: + density_matrix = qml.math.stack([density_matrix]) # Compute the partial trace traced_wires = [x for x in consecutive_indices if x not in indices] - density_matrix = _partial_trace(density_matrix, traced_wires) + density_matrix = _batched_partial_trace(density_matrix, traced_wires) + + if batch_dim is None: + density_matrix = density_matrix[0] + # Permute the remaining indices of the density matrix - return _permute_dense_matrix(density_matrix, sorted(indices), indices, None) + return _permute_dense_matrix(density_matrix, sorted(indices), indices, batch_dim) -def _partial_trace(density_matrix, indices): +def _batched_partial_trace(density_matrix, indices): """Compute the reduced density matrix by tracing out the provided indices. Args: - density_matrix (tensor_like): 2D density matrix tensor. This tensor should be of size ``(2**N, 2**N)`` for some - integer number of wires ``N``. + density_matrix (tensor_like): 3D density matrix tensor. This tensor should be of size + ``(batch_dim, 2**N, 2**N)``, for some integer number of wires``N``. indices (list(int)): List of indices to be traced. Returns: - tensor_like: (reduced) Density matrix of size ``(2**len(wires), 2**len(wires))`` + tensor_like: (reduced) Density matrix of size ``(batch_dim, 2**len(wires), 2**len(wires))`` **Example** - >>> x = np.array([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]) - >>> _partial_trace(x, indices=[0]) - [[1.+0.j 0.+0.j] - [0.+0.j 0.+0.j]] + >>> x = np.array([[[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], + ... [[0, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]]) + >>> _batched_partial_trace(x, indices=[0]) + array([[[1, 0], + [0, 0]], + [[0, 0], + [0, 1]]]) - >>> x = tf.Variable([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=tf.complex128) - >>> _partial_trace(x, indices=[1]) - tf.Tensor( - [[1.+0.j 0.+0.j] - [0.+0.j 0.+0.j]], shape=(2, 2), dtype=complex128) + >>> x = tf.Variable([[[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], + ... [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]]], dtype=tf.complex128) + >>> _batched_partial_trace(x, indices=[1]) + """ # Autograd does not support same indices sum in backprop if get_interface(density_matrix) == "autograd": - density_matrix = _partial_trace_autograd(density_matrix, indices) - return density_matrix + return _batched_partial_trace_autograd(density_matrix, indices) # Dimension and reshape - shape = density_matrix.shape[0] - num_indices = int(np.log2(shape)) + batch_dim, dim = density_matrix.shape[:2] + num_indices = int(np.log2(dim)) rho_dim = 2 * num_indices - density_matrix = np.reshape(density_matrix, [2] * 2 * num_indices) + density_matrix = np.reshape(density_matrix, [batch_dim] + [2] * 2 * num_indices) indices = np.sort(indices) # For loop over wires for i, target_index in enumerate(indices): target_index = target_index - i - state_indices = ABC[: rho_dim - 2 * i] + state_indices = ABC[1 : rho_dim - 2 * i + 1] state_indices = list(state_indices) target_letter = state_indices[target_index] state_indices[target_index + num_indices - i] = target_letter state_indices = "".join(state_indices) - einsum_indices = f"{state_indices}" + einsum_indices = f"a{state_indices}" density_matrix = einsum(einsum_indices, density_matrix) number_wires_sub = num_indices - len(indices) reduced_density_matrix = np.reshape( - density_matrix, (2**number_wires_sub, 2**number_wires_sub) + density_matrix, (batch_dim, 2**number_wires_sub, 2**number_wires_sub) ) return reduced_density_matrix -def _partial_trace_autograd(density_matrix, indices): +def _batched_partial_trace_autograd(density_matrix, indices): """Compute the reduced density matrix for autograd interface by tracing out the provided indices with the use of projectors as same subscripts indices are not supported in autograd backprop. """ # Dimension and reshape - shape = density_matrix.shape[0] - num_indices = int(np.log2(shape)) + batch_dim, dim = density_matrix.shape[:2] + num_indices = int(np.log2(dim)) rho_dim = 2 * num_indices - density_matrix = np.reshape(density_matrix, [2] * 2 * num_indices) + density_matrix = np.reshape(density_matrix, [batch_dim] + [2] * 2 * num_indices) kraus = cast(np.eye(2), density_matrix.dtype) @@ -306,22 +335,22 @@ def _partial_trace_autograd(density_matrix, indices): # For loop over wires for target_wire in indices: # Tensor indices of density matrix - state_indices = ABC[:rho_dim] + state_indices = ABC[1 : rho_dim + 1] # row indices of the quantum state affected by this operation - row_wires_list = [target_wire] + row_wires_list = [target_wire + 1] row_indices = "".join(ABC_ARRAY[row_wires_list].tolist()) # column indices are shifted by the number of wires col_wires_list = [w + num_indices for w in row_wires_list] col_indices = "".join(ABC_ARRAY[col_wires_list].tolist()) # indices in einsum must be replaced with new ones num_partial_trace_wires = 1 - new_row_indices = ABC[rho_dim : rho_dim + num_partial_trace_wires] + new_row_indices = ABC[rho_dim + 1 : rho_dim + num_partial_trace_wires + 1] new_col_indices = ABC[ - rho_dim + num_partial_trace_wires : rho_dim + 2 * num_partial_trace_wires + rho_dim + num_partial_trace_wires + 1 : rho_dim + 2 * num_partial_trace_wires + 1 ] # index for summation over Kraus operators kraus_index = ABC[ - rho_dim + 2 * num_partial_trace_wires : rho_dim + 2 * num_partial_trace_wires + 1 + rho_dim + 2 * num_partial_trace_wires + 1 : rho_dim + 2 * num_partial_trace_wires + 2 ] # new state indices replace row and column indices with new ones new_state_indices = functools.reduce( @@ -331,71 +360,108 @@ def _partial_trace_autograd(density_matrix, indices): ) # index mapping for einsum, e.g., 'iga,abcdef,idh->gbchef' einsum_indices = ( - f"{kraus_index}{new_row_indices}{row_indices}, {state_indices}," - f"{kraus_index}{col_indices}{new_col_indices}->{new_state_indices}" + f"{kraus_index}{new_row_indices}{row_indices}, a{state_indices}," + f"{kraus_index}{col_indices}{new_col_indices}->a{new_state_indices}" ) density_matrix = einsum(einsum_indices, kraus, density_matrix, kraus_dagger) number_wires_sub = num_indices - len(indices) reduced_density_matrix = np.reshape( - density_matrix, (2**number_wires_sub, 2**number_wires_sub) + density_matrix, (batch_dim, 2**number_wires_sub, 2**number_wires_sub) ) return reduced_density_matrix -def _density_matrix_from_state_vector(state, indices, check_state=False): +def reduce_statevector(state, indices, check_state=False, c_dtype="complex128"): """Compute the density matrix from a state vector. Args: - state (tensor_like): 1D tensor state vector. This tensor should of size ``(2**N,)`` for some integer value ``N``. + state (tensor_like): 1D or 2D tensor state vector. This tensor should of size ``(2**N,)`` + or ``(batch_dim, 2**N)``, for some integer value ``N``. indices (list(int)): List of indices in the considered subsystem. check_state (bool): If True, the function will check the state validity (shape and norm). + c_dtype (str): Complex floating point precision type. Returns: - tensor_like: Density matrix of size ``(2**len(indices), 2**len(indices))`` + tensor_like: Density matrix of size ``(2**len(indices), 2**len(indices))`` or ``(batch_dim, 2**len(indices), 2**len(indices))`` + + .. seealso:: :func:`pennylane.math.reduce_dm`, :func:`pennylane.qinfo.transforms.reduced_dm`, + and :func:`pennylane.density_matrix` **Example** >>> x = np.array([1, 0, 0, 0]) - >>> _density_matrix_from_state_vector(x, indices=[0]) + >>> reduce_statevector(x, indices=[0]) [[1.+0.j 0.+0.j] [0.+0.j 0.+0.j]] >>> y = [1, 0, 1, 0] / np.sqrt(2) - >>> _density_matrix_from_state_vector(y, indices=[0]) + >>> reduce_statevector(y, indices=[0]) [[0.5+0.j 0.5+0.j] [0.5+0.j 0.5+0.j]] - >>> _density_matrix_from_state_vector(y, indices=[1]) + >>> reduce_statevector(y, indices=[1]) [[1.+0.j 0.+0.j] [0.+0.j 0.+0.j]] >>> z = tf.Variable([1, 0, 0, 0], dtype=tf.complex128) - >>> _density_matrix_from_state_vector(z, indices=[1]) + >>> reduce_statevector(z, indices=[1]) tf.Tensor( [[1.+0.j 0.+0.j] [0.+0.j 0.+0.j]], shape=(2, 2), dtype=complex128) + >>> x = np.array([[1, 0, 0, 0], [0, 1, 0, 0]]) + >>> reduce_statevector(x, indices=[1]) + array([[[1.+0.j, 0.+0.j], + [0.+0.j, 0.+0.j]], + + [[0.+0.j, 0.+0.j], + [0.+0.j, 1.+0.j]]]) """ - dim = np.shape(state)[0] + state = cast(state, dtype=c_dtype) # Check the format and norm of the state vector if check_state: _check_state_vector(state) + if len(np.shape(state)) == 1: + batch_dim, dim = None, np.shape(state)[0] + else: + batch_dim, dim = np.shape(state)[:2] + # Get dimension of the quantum system and reshape num_wires = int(np.log2(dim)) consecutive_wires = list(range(num_wires)) - state = np.reshape(state, [2] * num_wires) + + if batch_dim is None: + state = qml.math.stack([state]) + + state = np.reshape(state, [batch_dim if batch_dim is not None else 1] + [2] * num_wires) # Get the system to be traced - traced_system = [x for x in consecutive_wires if x not in indices] + # traced_system = [x + 1 for x in consecutive_wires if x not in indices] + + # trace out the subsystem + indices1 = ABC[1 : num_wires + 1] + indices2 = "".join( + [ABC[num_wires + i + 1] if i in indices else ABC[i + 1] for i in consecutive_wires] + ) + target = "".join( + [ABC[i + 1] for i in sorted(indices)] + [ABC[num_wires + i + 1] for i in sorted(indices)] + ) + density_matrix = einsum(f"a{indices1},a{indices2}->a{target}", state, np.conj(state)) # Return the reduced density matrix by using numpy tensor product - density_matrix = np.tensordot(state, np.conj(state), axes=(traced_system, traced_system)) - density_matrix = np.reshape(density_matrix, (2 ** len(indices), 2 ** len(indices))) + # density_matrix = np.tensordot(state, np.conj(state), axes=(traced_system, traced_system)) + + if batch_dim is None: + density_matrix = np.reshape(density_matrix, (2 ** len(indices), 2 ** len(indices))) + else: + density_matrix = np.reshape( + density_matrix, (batch_dim, 2 ** len(indices), 2 ** len(indices)) + ) - return _permute_dense_matrix(density_matrix, sorted(indices), indices, None) + return _permute_dense_matrix(density_matrix, sorted(indices), indices, batch_dim) def reduced_dm(state, indices, check_state=False, c_dtype="complex128"): @@ -445,14 +511,19 @@ def reduced_dm(state, indices, check_state=False, c_dtype="complex128"): .. seealso:: :func:`pennylane.qinfo.transforms.reduced_dm` and :func:`pennylane.density_matrix` """ + warnings.warn( + "reduced_dm has been deprecated. Please use reduce_statevector or reduce_dm instead", + UserWarning, + ) + # Cast as a c_dtype array state = cast(state, dtype=c_dtype) len_state = state.shape[0] # State vector if state.shape == (len_state,): - return _density_matrix_from_state_vector(state, indices, check_state) + return reduce_statevector(state, indices, check_state) - return _density_matrix_from_matrix(state, indices, check_state) + return reduce_dm(state, indices, check_state) def purity(state, indices, check_state=False, c_dtype="complex128"): @@ -516,11 +587,11 @@ def purity(state, indices, check_state=False, c_dtype="complex128"): # If the state is a state vector but the system in question is a sub-system of the # overall state, then the purity of the sub-system still needs to be computed. if state.shape == (len_state,): - density_matrix = _density_matrix_from_state_vector(state, indices, check_state) + density_matrix = reduce_statevector(state, indices, check_state) return _compute_purity(density_matrix) # If the state is a density matrix, compute the purity. - density_matrix = _density_matrix_from_matrix(state, indices, check_state) + density_matrix = reduce_dm(state, indices, check_state) return _compute_purity(density_matrix) @@ -996,41 +1067,52 @@ def relative_entropy(state0, state1, base=None, check_state=False, c_dtype="comp def _check_density_matrix(density_matrix): """Check the shape, the trace and the positive semi-definitiveness of a matrix.""" - shape = density_matrix.shape[0] + dim = density_matrix.shape[-1] if ( - len(density_matrix.shape) != 2 - or density_matrix.shape[0] != density_matrix.shape[1] - or not np.log2(shape).is_integer() + len(density_matrix.shape) not in (2, 3) + or density_matrix.shape[-2] != dim + or not np.log2(dim).is_integer() ): - raise ValueError("Density matrix must be of shape (2**N, 2**N).") - # Check trace - trace = np.trace(density_matrix) - if not is_abstract(trace): - if not allclose(trace, 1.0, atol=1e-10): - raise ValueError("The trace of the density matrix should be one.") - # Check if the matrix is Hermitian - conj_trans = np.transpose(np.conj(density_matrix)) - if not allclose(density_matrix, conj_trans): - raise ValueError("The matrix is not Hermitian.") - # Check if positive semi-definite - evs, _ = qml.math.linalg.eigh(density_matrix) - evs = np.real(evs) - evs_non_negative = [ev for ev in evs if ev >= 0.0] - if len(evs) != len(evs_non_negative): - raise ValueError("The matrix is not positive semi-definite.") + raise ValueError("Density matrix must be of shape (2**N, 2**N) or (batch_dim, 2**N, 2**N).") + + if len(density_matrix.shape) == 2: + density_matrix = qml.math.stack([density_matrix]) + + if not is_abstract(density_matrix): + for dm in density_matrix: + # Check trace + trace = np.trace(dm) + if not allclose(trace, 1.0, atol=1e-10): + raise ValueError("The trace of the density matrix should be one.") + + # Check if the matrix is Hermitian + conj_trans = np.transpose(np.conj(dm)) + if not allclose(dm, conj_trans): + raise ValueError("The matrix is not Hermitian.") + + # Check if positive semi-definite + evs, _ = qml.math.linalg.eigh(dm) + evs = np.real(evs) + evs_non_negative = [ev for ev in evs if ev >= 0.0] + if len(evs) != len(evs_non_negative): + raise ValueError("The matrix is not positive semi-definite.") def _check_state_vector(state_vector): """Check the shape and the norm of a state vector.""" - len_state = state_vector.shape[0] - # Check format - if len(np.shape(state_vector)) != 1 or not np.log2(len_state).is_integer(): - raise ValueError("State vector must be of length 2**wires.") + dim = state_vector.shape[-1] + if len(np.shape(state_vector)) not in (1, 2) or not np.log2(dim).is_integer(): + raise ValueError("State vector must be of shape (2**wires,) or (batch_dim, 2**wires)") + + if len(state_vector.shape) == 1: + state_vector = qml.math.stack([state_vector]) + # Check norm - norm = np.linalg.norm(state_vector, ord=2) - if not is_abstract(norm): - if not allclose(norm, 1.0, atol=1e-10): - raise ValueError("Sum of amplitudes-squared does not equal one.") + if not is_abstract(state_vector): + for sv in state_vector: + norm = np.linalg.norm(sv, ord=2) + if not allclose(norm, 1.0, atol=1e-10): + raise ValueError("Sum of amplitudes-squared does not equal one.") def max_entropy(state, indices, base=None, check_state=False, c_dtype="complex128"): diff --git a/pennylane/measurements/state.py b/pennylane/measurements/state.py index 8f220216fec..22e89d2aae2 100644 --- a/pennylane/measurements/state.py +++ b/pennylane/measurements/state.py @@ -188,6 +188,6 @@ def process_state(self, state: Sequence[complex], wire_order: Wires): # qml.density_matrix wire_map = dict(zip(wire_order, range(len(wire_order)))) mapped_wires = [wire_map[w] for w in self.wires] - return qml.math.reduced_dm(state, indices=mapped_wires, c_dtype=state.dtype) + return qml.math.reduce_statevector(state, indices=mapped_wires, c_dtype=state.dtype) # qml.state return state diff --git a/pennylane/qinfo/transforms.py b/pennylane/qinfo/transforms.py index 6cee18ebe85..6d9a298c4fe 100644 --- a/pennylane/qinfo/transforms.py +++ b/pennylane/qinfo/transforms.py @@ -58,11 +58,15 @@ def wrapper(*args, **kwargs): if len(measurements) != 1 or not isinstance(measurements[0], StateMP): raise ValueError("The qfunc measurement needs to be State.") + # determine if the measurement is a state vector or a density matrix + # TODO: once we separate StateMP and DensityMatrixMP, we can replace this + # line with isinstance checks + dm_measurement = measurements[0].wires or "mixed" in qnode.device.name + dm_func = qml.math.reduce_statevector if not dm_measurement else qml.math.reduce_dm + # TODO: optimize given the wires by creating a tape with relevant operations state_built = qnode(*args, **kwargs) - density_matrix = qml.math.reduced_dm( - state_built, indices=wires, c_dtype=qnode.device.C_DTYPE - ) + density_matrix = dm_func(state_built, indices=wires, c_dtype=qnode.device.C_DTYPE) return density_matrix return wrapper diff --git a/tests/math/test_density_matrices.py b/tests/math/test_density_matrices.py index 2d43871e0e2..263bfa3410d 100644 --- a/tests/math/test_density_matrices.py +++ b/tests/math/test_density_matrices.py @@ -104,25 +104,23 @@ class TestDensityMatrixFromStateVectors: @pytest.mark.parametrize("array_func", array_funcs) @pytest.mark.parametrize("state_vector, expected_density_matrix", state_vectors) @pytest.mark.parametrize("wires", single_wires_list) - def test_density_matrix_from_state_vector_single_wires( + def test_reduce_statevector_single_wires( self, state_vector, wires, expected_density_matrix, array_func ): """Test the density matrix from state vectors for single wires.""" - # pylint: disable=protected-access state_vector = array_func(state_vector) - density_matrix = fn.quantum._density_matrix_from_state_vector(state_vector, indices=wires) + density_matrix = fn.quantum.reduce_statevector(state_vector, indices=wires) assert np.allclose(density_matrix, expected_density_matrix[wires[0]]) @pytest.mark.parametrize("array_func", array_funcs) @pytest.mark.parametrize("state_vector, expected_density_matrix", state_vectors) @pytest.mark.parametrize("wires", multiple_wires_list) - def test_density_matrix_from_state_vector_full_wires( + def test_reduce_statevector_full_wires( self, state_vector, wires, expected_density_matrix, array_func ): """Test the density matrix from state vectors for full wires.""" - # pylint: disable=protected-access state_vector = array_func(state_vector) - density_matrix = fn.quantum._density_matrix_from_state_vector(state_vector, indices=wires) + density_matrix = fn.quantum.reduce_statevector(state_vector, indices=wires) expected = expected_density_matrix[2] if wires == [1, 0]: expected = permute_two_qubit_dm(expected) @@ -136,7 +134,7 @@ def test_reduced_dm_with_state_vector_single_wires( ): """Test the reduced_dm with state vectors for single wires.""" state_vector = array_func(state_vector) - density_matrix = fn.reduced_dm(state_vector, indices=wires) + density_matrix = fn.reduce_statevector(state_vector, indices=wires) assert np.allclose(density_matrix, expected_density_matrix[wires[0]]) @pytest.mark.parametrize("array_func", array_funcs) @@ -147,7 +145,7 @@ def test_reduced_dm_with_state_vector_full_wires( ): """Test the reduced_dm with state vectors for full wires.""" state_vector = array_func(state_vector) - density_matrix = fn.reduced_dm(state_vector, indices=wires) + density_matrix = fn.reduce_statevector(state_vector, indices=wires) expected = expected_density_matrix[2] if wires == [1, 0]: expected = permute_two_qubit_dm(expected) @@ -156,12 +154,14 @@ def test_reduced_dm_with_state_vector_full_wires( @pytest.mark.parametrize("array_func", array_funcs) @pytest.mark.parametrize("state_vector, expected_density_matrix", state_vectors) @pytest.mark.parametrize("wires", multiple_wires_list) - def test_density_matrix_from_state_vector_check_state( + def test_reduce_statevector_check_state( self, state_vector, wires, expected_density_matrix, array_func ): """Test the density matrix from state vectors for single wires with state checking""" state_vector = array_func(state_vector) - density_matrix = fn.quantum.reduced_dm(state_vector, indices=wires, check_state=True) + density_matrix = fn.quantum.reduce_statevector( + state_vector, indices=wires, check_state=True + ) expected = expected_density_matrix[2] if wires == [1, 0]: expected = permute_two_qubit_dm(expected) @@ -173,25 +173,23 @@ def test_state_vector_wrong_shape(self): state_vector = [1, 0, 0] with pytest.raises(ValueError, match="State vector must be"): - fn.quantum.reduced_dm(state_vector, indices=[0], check_state=True) + fn.quantum.reduce_statevector(state_vector, indices=[0], check_state=True) def test_state_vector_wrong_norm(self): """Test that state vector with wrong norm raises an error with check_state=True""" state_vector = [0.1, 0, 0, 0] with pytest.raises(ValueError, match="Sum of amplitudes-squared does not equal one."): - fn.quantum.reduced_dm(state_vector, indices=[0], check_state=True) + fn.quantum.reduce_statevector(state_vector, indices=[0], check_state=True) - def test_density_matrix_from_state_vector_jax_jit(self): + def test_reduce_statevector_jax_jit(self): """Test jitting the density matrix from state vector function.""" # pylint: disable=protected-access from jax import jit state_vector = jnp.array([1, 0, 0, 0]) - jitted_dens_matrix_func = jit( - fn.quantum._density_matrix_from_state_vector, static_argnums=[1, 2] - ) + jitted_dens_matrix_func = jit(fn.quantum.reduce_statevector, static_argnums=[1, 2]) density_matrix = jitted_dens_matrix_func(state_vector, indices=(0, 1), check_state=True) assert np.allclose(density_matrix, [[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]) @@ -203,9 +201,7 @@ def test_wrong_shape_jax_jit(self): state_vector = jnp.array([1, 0, 0]) - jitted_dens_matrix_func = jit( - fn.quantum._density_matrix_from_state_vector, static_argnums=[1, 2] - ) + jitted_dens_matrix_func = jit(fn.quantum.reduce_statevector, static_argnums=[1, 2]) with pytest.raises(ValueError, match="State vector must be"): jitted_dens_matrix_func(state_vector, indices=(0, 1), check_state=True) @@ -216,7 +212,7 @@ def test_density_matrix_tf_jit(self): state_vector = tf.Variable([1, 0, 0, 0], dtype=tf.complex128) - density_matrix = partial(fn.reduced_dm, indices=[0]) + density_matrix = partial(fn.reduce_statevector, indices=[0]) density_matrix = tf.function( density_matrix, @@ -235,7 +231,7 @@ def test_density_matrix_c_dtype(self, array_func, state_vector, wires, c_dtype): state_vector = array_func(state_vector) if fn.get_interface(state_vector) == "jax" and c_dtype == "complex128": pytest.skip("Jax does not support complex 128") - density_matrix = fn.reduced_dm(state_vector, indices=wires, c_dtype=c_dtype) + density_matrix = fn.reduce_statevector(state_vector, indices=wires, c_dtype=c_dtype) if fn.get_interface(state_vector) == "torch": if c_dtype == "complex64": c_dtype = torch.complex64 @@ -290,14 +286,14 @@ def test_reduced_dm_with_matrix_single_wires( self, density_matrix, wires, expected_density_matrix ): """Test the reduced_dm with matrix for single wires.""" - density_matrix = fn.reduced_dm(density_matrix, indices=wires) + density_matrix = fn.reduce_dm(density_matrix, indices=wires) assert np.allclose(density_matrix, expected_density_matrix[wires[0]]) @pytest.mark.parametrize("density_matrix", list(zip(*density_matrices))[0]) @pytest.mark.parametrize("wires", multiple_wires_list) def test_reduced_dm_with_matrix_full_wires(self, density_matrix, wires): """Test the reduced_dm with matrix for full wires.""" - returned_density_matrix = fn.reduced_dm(density_matrix, indices=wires) + returned_density_matrix = fn.reduce_dm(density_matrix, indices=wires) expected = density_matrix if wires == [1, 0]: expected = permute_two_qubit_dm(expected) @@ -305,10 +301,10 @@ def test_reduced_dm_with_matrix_full_wires(self, density_matrix, wires): @pytest.mark.parametrize("density_matrix", list(zip(*density_matrices))[0]) @pytest.mark.parametrize("wires", multiple_wires_list) - def test_density_matrix_from_matrix_check(self, density_matrix, wires): + def test_reduce_dm_check(self, density_matrix, wires): """Test the density matrix from matrices for single wires with state checking""" # pylint: disable=protected-access - returned_density_matrix = fn.quantum._density_matrix_from_matrix( + returned_density_matrix = fn.quantum.reduce_dm( density_matrix, indices=wires, check_state=True ) expected = density_matrix @@ -322,36 +318,36 @@ def test_matrix_wrong_shape(self): density_matrix = [[1, 0, 0], [0, 0, 0], [0, 0, 0]] with pytest.raises(ValueError, match="Density matrix must be of shape"): - fn.quantum.reduced_dm(density_matrix, indices=[0], check_state=True) + fn.quantum.reduce_dm(density_matrix, indices=[0], check_state=True) def test_matrix_wrong_trace(self): """Test that density matrix with wrong trace raises an error with check_state=True""" density_matrix = [[0.1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]] with pytest.raises(ValueError, match="The trace of the density matrix should be one."): - fn.quantum.reduced_dm(density_matrix, indices=[0], check_state=True) + fn.quantum.reduce_dm(density_matrix, indices=[0], check_state=True) def test_matrix_not_hermitian(self): """Test that non hermitian matrix raises an error with check_state=True""" density_matrix = [[0.1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0.5, 0.9]] with pytest.raises(ValueError, match="The matrix is not Hermitian."): - fn.quantum.reduced_dm(density_matrix, indices=[0], check_state=True) + fn.quantum.reduce_dm(density_matrix, indices=[0], check_state=True) def test_matrix_not_positive_definite(self): """Test that non hermitian matrix raises an error with check_state=True""" density_matrix = [[3, 0], [0, -2]] with pytest.raises(ValueError, match="The matrix is not positive semi-definite."): - fn.quantum.reduced_dm(density_matrix, indices=[0], check_state=True) + fn.quantum.reduce_dm(density_matrix, indices=[0], check_state=True) - def test_density_matrix_from_state_vector_jax_jit(self): + def test_reduce_statevector_jax_jit(self): """Test jitting the density matrix from state vector function.""" from jax import jit state_vector = jnp.array([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]) - jitted_dens_matrix_func = jit(fn.quantum.reduced_dm, static_argnums=[1, 2]) + jitted_dens_matrix_func = jit(fn.quantum.reduce_dm, static_argnums=[1, 2]) density_matrix = jitted_dens_matrix_func(state_vector, indices=(0,), check_state=True) assert np.allclose(density_matrix, [[1, 0], [0, 0]]) @@ -363,7 +359,7 @@ def test_wrong_shape_jax_jit(self): state_vector = jnp.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]]) - jitted_dens_matrix_func = jit(fn.quantum._density_matrix_from_matrix, static_argnums=[1, 2]) + jitted_dens_matrix_func = jit(fn.quantum.reduce_dm, static_argnums=[1, 2]) with pytest.raises(ValueError, match="Density matrix must be of shape"): jitted_dens_matrix_func(state_vector, indices=(0, 1), check_state=True) @@ -381,7 +377,7 @@ def test_density_matrix_tf_jit(self): ], dtype=tf.complex128, ) - density_matrix = partial(fn.reduced_dm, indices=[0]) + density_matrix = partial(fn.reduce_dm, indices=[0]) density_matrix = tf.function( density_matrix, @@ -398,10 +394,77 @@ def test_density_matrix_c_dtype(self, density_matrix, wires, c_dtype): """Test different complex dtype.""" if fn.get_interface(density_matrix) == "jax" and c_dtype == "complex128": pytest.skip("Jax does not support complex 128") - density_matrix = fn.reduced_dm(density_matrix, indices=wires, c_dtype=c_dtype) + density_matrix = fn.reduce_dm(density_matrix, indices=wires, c_dtype=c_dtype) if fn.get_interface(density_matrix) == "torch": if c_dtype == "complex64": c_dtype = torch.complex64 elif c_dtype == "complex128": c_dtype = torch.complex128 assert density_matrix.dtype == c_dtype + + +class TestDensityMatrixBroadcasting: + """Test that broadcasting works as expected for the reduced_dm functions""" + + @pytest.mark.parametrize("array_func", array_funcs) + @pytest.mark.parametrize("wires", single_wires_list) + def test_sv_broadcast_single_wires(self, wires, array_func): + """Test that broadcasting works for state vectors and single wires""" + state_vector = array_func([sv[0] for sv in state_vectors]) + expected_density_matrices = array_func([sv[1][wires[0]] for sv in state_vectors]) + + density_matrices = fn.reduce_statevector(state_vector, indices=wires) + assert np.allclose(density_matrices, expected_density_matrices) + + @pytest.mark.parametrize("array_func", array_funcs) + @pytest.mark.parametrize("wires", multiple_wires_list) + def test_reduced_dm_with_state_vector_full_wires(self, wires, array_func): + """Test that broadcasting works for state vectors and full wires""" + state_vector = array_func([sv[0] for sv in state_vectors]) + + if wires == [0, 1]: + expected_density_matrices = array_func([sv[1][2] for sv in state_vectors]) + else: + expected_density_matrices = array_func( + [permute_two_qubit_dm(sv[1][2]) for sv in state_vectors] + ) + + density_matrices = fn.reduce_statevector(state_vector, indices=wires) + assert np.allclose(density_matrices, expected_density_matrices) + + @pytest.mark.parametrize("array_func", array_funcs) + @pytest.mark.parametrize("wires", single_wires_list) + def test_reduced_dm_with_matrix_single_wires(self, wires, array_func): + """Test that broadcasting works for density matrices and single wires""" + density_matrix = array_func([dm[0] for dm in density_matrices[:4]]) + expected_density_matrices = array_func([dm[1][wires[0]] for dm in density_matrices[:4]]) + + density_matrix = fn.reduce_dm(density_matrix, indices=wires) + assert np.allclose(density_matrix, expected_density_matrices) + + @pytest.mark.parametrize("array_func", array_funcs) + @pytest.mark.parametrize("wires", multiple_wires_list) + def test_reduced_dm_with_matrix_full_wires(self, wires, array_func): + """Test that broadcasting works for density matrices and full wires""" + density_matrix = array_func([dm[0] for dm in density_matrices[:4]]) + + if wires == [0, 1]: + expected_density_matrices = density_matrix + else: + expected_density_matrices = array_func( + [permute_two_qubit_dm(dm[0]) for dm in density_matrices[:4]] + ) + + density_matrix = fn.reduce_dm(density_matrix, indices=wires) + assert np.allclose(density_matrix, expected_density_matrices) + + +def test_deprecation_warning(): + """Test that a deprecation warning is raised when qml.math.reduced_dm is called""" + state = np.array([1, 0]) + with pytest.warns(UserWarning, match="reduced_dm has been deprecated"): + fn.reduced_dm(state, [0]) + + dm = np.array([[1, 0], [0, 0]]) + with pytest.warns(UserWarning, match="reduced_dm has been deprecated"): + fn.reduced_dm(dm, [0]) diff --git a/tests/math/test_fidelity_math.py b/tests/math/test_fidelity_math.py index 328f7212026..7bcf6b6ea84 100644 --- a/tests/math/test_fidelity_math.py +++ b/tests/math/test_fidelity_math.py @@ -91,7 +91,7 @@ def test_state_vector_wrong_amplitudes(self, state0, state1): @pytest.mark.parametrize("state0,state1", state_wrong_shape) def test_state_vector_wrong_shape(self, state0, state1): """Test that a message is raised when the state does not have the right shape.""" - with pytest.raises(ValueError, match="State vector must be of length"): + with pytest.raises(ValueError, match="State vector must be of shape"): qml.math.fidelity(state0, state1, check_state=True) d_mat_wrong_shape = [ diff --git a/tests/measurements/test_state.py b/tests/measurements/test_state.py index dc35b09d227..5198c77ecaa 100644 --- a/tests/measurements/test_state.py +++ b/tests/measurements/test_state.py @@ -19,7 +19,7 @@ from pennylane import numpy as pnp from pennylane.devices import DefaultQubit from pennylane.measurements import State, StateMP, Shots, density_matrix, expval, state -from pennylane.math.quantum import _density_matrix_from_state_vector, _density_matrix_from_matrix +from pennylane.math.quantum import reduce_statevector, reduce_dm from pennylane.math.matrix_manipulation import _permute_dense_matrix # pylint: disable=no-member, comparison-with-callable, import-outside-toplevel @@ -68,9 +68,10 @@ def test_process_state_matrix_from_vec(self, vec, wires): if len(wires) == num_wires: exp = np.outer(vec, vec.conj()) else: - exp = _density_matrix_from_state_vector(vec, wires) + exp = reduce_statevector(vec, wires) assert qml.math.allclose(processed, exp) + @pytest.mark.xfail(reason="StateMP.process_state no longer supports density matrix parameters") @pytest.mark.parametrize( "mat, wires", [ @@ -96,7 +97,7 @@ def test_process_state_matrix_from_matrix(self, mat, wires): if len(wires) == num_wires: exp = _permute_dense_matrix(mat, wires, order, None) else: - exp = _density_matrix_from_matrix(mat, wires) + exp = reduce_dm(mat, wires) assert qml.math.allclose(processed, exp) @@ -479,10 +480,11 @@ def func(): dev = func.device - assert np.allclose( - expected, - qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), - ) + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), + ) @pytest.mark.jax @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) @@ -501,10 +503,11 @@ def func(): assert np.allclose(expected, density_mat) - assert np.allclose( - expected, - qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), - ) + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), + ) @pytest.mark.tf @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) @@ -523,10 +526,11 @@ def func(): assert np.allclose(expected, density_mat) - assert np.allclose( - expected, - qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), - ) + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) def test_correct_density_matrix_product_state_first(self, dev_name): @@ -546,10 +550,11 @@ def func(): assert np.allclose(expected, density_first) - assert np.allclose( - expected, - qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), - ) + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=0).process_state(state=dev.state, wire_order=dev.wires), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) def test_correct_density_matrix_product_state_second(self, dev_name): @@ -568,10 +573,11 @@ def func(): expected = np.array([[0.5 + 0.0j, 0.5 + 0.0j], [0.5 + 0.0j, 0.5 + 0.0j]]) assert np.allclose(expected, density_second) - assert np.allclose( - expected, - qml.density_matrix(wires=1).process_state(state=dev.state, wire_order=dev.wires), - ) + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=1).process_state(state=dev.state, wire_order=dev.wires), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) @pytest.mark.parametrize("return_wire_order", ([0, 1], [1, 0])) @@ -594,12 +600,13 @@ def func(): assert np.allclose(expected, density_both) - assert np.allclose( - expected, - qml.density_matrix(wires=return_wire_order).process_state( - state=dev.state, wire_order=dev.wires - ), - ) + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=return_wire_order).process_state( + state=dev.state, wire_order=dev.wires + ), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) def test_correct_density_matrix_three_wires_first_two(self, dev_name): @@ -624,10 +631,14 @@ def func(): ] ) assert np.allclose(expected, density_full) - assert np.allclose( - expected, - qml.density_matrix(wires=[0, 1]).process_state(state=dev.state, wire_order=dev.wires), - ) + + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=[0, 1]).process_state( + state=dev.state, wire_order=dev.wires + ), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) def test_correct_density_matrix_three_wires_last_two(self, dev_name): @@ -656,10 +667,14 @@ def func(): ) assert np.allclose(expected, density) - assert np.allclose( - expected, - qml.density_matrix(wires=[1, 2]).process_state(state=dev.state, wire_order=dev.wires), - ) + + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=[1, 2]).process_state( + state=dev.state, wire_order=dev.wires + ), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) @pytest.mark.parametrize( @@ -692,12 +707,14 @@ def func(): expected = np.outer(exp_statevector.conj(), exp_statevector) assert np.allclose(expected, density_full) - assert np.allclose( - expected, - qml.density_matrix(wires=return_wire_order).process_state( - state=dev.state, wire_order=dev.wires - ), - ) + + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=return_wire_order).process_state( + state=dev.state, wire_order=dev.wires + ), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) def test_correct_density_matrix_mixed_state(self, dev_name): @@ -738,10 +755,14 @@ def func(): ) assert np.allclose(expected, density) - assert np.allclose( - expected, - qml.density_matrix(wires=[0, 1]).process_state(state=dev.state, wire_order=dev.wires), - ) + + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=[0, 1]).process_state( + state=dev.state, wire_order=dev.wires + ), + ) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"]) def test_return_with_other_types(self, dev_name): @@ -812,10 +833,14 @@ def func(): expected = np.array([[0.5 + 0.0j, 0.0 + 0.0j], [0.0 + 0.0j, 0.5 + 0.0j]]) assert np.allclose(expected, density) - assert np.allclose( - expected, - qml.density_matrix(wires=wires[1]).process_state(state=dev.state, wire_order=dev.wires), - ) + + if dev_name != "default.mixed": + assert np.allclose( + expected, + qml.density_matrix(wires=wires[1]).process_state( + state=dev.state, wire_order=dev.wires + ), + ) @pytest.mark.parametrize("wires", [[3, 1], ["b", 1000]]) @pytest.mark.parametrize("dev_name", ["default.qubit", "default.mixed"])