Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add density matrix and entanglement calculations for qubits #203

Merged
merged 9 commits into from
Jul 30, 2024
Merged
Show file tree
Hide file tree
Changes from 8 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 77 additions & 0 deletions unitary/alpha/quantum_world.py
Original file line number Diff line number Diff line change
Expand Up @@ -636,8 +636,85 @@ def get_binary_probabilities(
binary_probs.append(1 - one_probs[0])
return binary_probs

def density_matrix(
self, objects: Optional[Sequence[QuantumObject]] = None, count: int = 1000
) -> np.ndarray:
"""Simulates the density matrix of the given objects.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's update the description here and then add some details about how we are assuming that this is a pure state so the density matrix is |psi><psi|

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done. Thanks.


Parameters:
objects: List of QuantumObjects (currently only qubits are supported)
count: Number of measurements

Returns:
The density matrix of the specified objects.
"""
num_all_qubits = len(self.object_name_dict.values())
num_shown_qubits = len(objects) if objects is not None else num_all_qubits

specified_names = (
[obj.qubit.name for obj in objects] if objects is not None else []
)
unspecified_names = set(self.object_name_dict.keys()) - set(specified_names)

# Make sure we have all objects, starting with the specified ones in the given order.
ordered_names = specified_names + list(unspecified_names)
ordered_objects = [self.object_name_dict[name] for name in ordered_names]

# Peek the current world `count` times and get the results.
histogram = self.get_correlated_histogram(ordered_objects, count)

# Get an estimate of the state vector.
state_vector = np.array([0.0] * (2**num_all_qubits))
for key, val in histogram.items():
state_vector += self.__to_state_vector__(key) * np.sqrt(val * 1.0 / count)
density_matrix = np.outer(state_vector, state_vector)

if num_shown_qubits == num_all_qubits:
return density_matrix
else:
# We trace out the unspecified qubits.
# The reshape is required by the partial_trace method.
traced_density_matrix = cirq.partial_trace(
density_matrix.reshape((2, 2) * num_all_qubits),
range(num_shown_qubits),
)
# Reshape back to a 2-d matrix.
return traced_density_matrix.reshape(
2**num_shown_qubits, 2**num_shown_qubits
)

def measure_entanglement(self, obj1: QuantumObject, obj2: QuantumObject) -> float:
"""Measures the entanglement (i.e. quantum mutual information) of the two given objects.
See https://en.wikipedia.org/wiki/Quantum_mutual_information for the formula.

Parameters:
obj1, obj2: two quantum objects (currently only qubits are supported)

Returns:
The quantum mutual information defined as S_1 + S_2 - S_12, where S denotes (reduced)
von Neumann entropy.
"""
density_matrix_12 = self.density_matrix([obj1, obj2]).reshape(2, 2, 2, 2)
density_matrix_1 = cirq.partial_trace(density_matrix_12, [0])
density_matrix_2 = cirq.partial_trace(density_matrix_12, [1])
return (
cirq.von_neumann_entropy(density_matrix_1, validate=False)
+ cirq.von_neumann_entropy(density_matrix_2, validate=False)
- cirq.von_neumann_entropy(density_matrix_12.reshape(4, 4), validate=False)
)

def __getitem__(self, name: str) -> QuantumObject:
quantum_object = self.object_name_dict.get(name, None)
if not quantum_object:
raise KeyError(f"{name} did not exist in this world.")
return quantum_object

def __to_state_vector__(self, input: tuple) -> np.ndarray:
"""Converts the given tuple (of length N) to the corresponding state vector (of length 2**N).
e.g. (0, 1) -> [0, 1, 0, 0]
"""
num = len(input)
index = int("".join([str(i) for i in input]), 2)
state_vector = np.array([0.0] * (2**num))
state_vector[index] = 1.0
return state_vector
81 changes: 81 additions & 0 deletions unitary/alpha/quantum_world_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
import pytest

import numpy as np
import numpy.testing as testing

import cirq

Expand Down Expand Up @@ -892,3 +893,83 @@ def test_save_and_restore_snapshot(simulator, compile_to_qubits):
# Further restore would return a value error.
with pytest.raises(ValueError, match="Unable to restore any more."):
world.restore_last_snapshot()


@pytest.mark.parametrize(
("simulator", "compile_to_qubits"),
[
(cirq.Simulator, False),
(cirq.Simulator, True),
# Cannot use SparseSimulator without `compile_to_qubits` due to issue #78.
(alpha.SparseSimulator, True),
],
)
def test_density_matrix(simulator, compile_to_qubits):
rho_green = np.reshape([0, 0, 0, 1], (2, 2))
rho_red = np.reshape([1, 0, 0, 0], (2, 2))
light1 = alpha.QuantumObject("green", Light.GREEN)
light2 = alpha.QuantumObject("red1", Light.RED)
light3 = alpha.QuantumObject("red2", Light.RED)
board = alpha.QuantumWorld(
[light1, light2, light3],
sampler=simulator(),
compile_to_qubits=compile_to_qubits,
)

testing.assert_array_equal(board.density_matrix(objects=[light1]), rho_green)
testing.assert_array_equal(board.density_matrix(objects=[light2]), rho_red)
testing.assert_array_equal(board.density_matrix(objects=[light3]), rho_red)

testing.assert_array_equal(
board.density_matrix(objects=[light1, light2]), np.kron(rho_green, rho_red)
)
testing.assert_array_equal(
board.density_matrix(objects=[light2, light3]), np.kron(rho_red, rho_red)
)
testing.assert_array_equal(
board.density_matrix(objects=[light3, light1]), np.kron(rho_red, rho_green)
)

testing.assert_array_equal(
board.density_matrix(objects=[light1, light2, light3]),
np.kron(rho_green, np.kron(rho_red, rho_red)),
)


@pytest.mark.parametrize(
("simulator", "compile_to_qubits"),
[
(cirq.Simulator, False),
(cirq.Simulator, True),
# Cannot use SparseSimulator without `compile_to_qubits` due to issue #78.
(alpha.SparseSimulator, True),
],
)
def test_measure_entanglement(simulator, compile_to_qubits):
rho_green = np.reshape([0, 0, 0, 1], (2, 2))
rho_red = np.reshape([1, 0, 0, 0], (2, 2))
light1 = alpha.QuantumObject("red1", Light.RED)
light2 = alpha.QuantumObject("green", Light.GREEN)
light3 = alpha.QuantumObject("red2", Light.RED)
board = alpha.QuantumWorld(
[light1, light2, light3],
sampler=simulator(),
compile_to_qubits=compile_to_qubits,
)

# S_1 + S_2 - S_12 = 0 + 0 - 0 = 0 for all three cases.
assert round(board.measure_entanglement(light1, light2)) == 0.0
assert round(board.measure_entanglement(light1, light3)) == 0.0
assert round(board.measure_entanglement(light2, light3)) == 0.0

alpha.Superposition()(light2)
alpha.quantum_if(light2).apply(alpha.Flip())(light3)
results = board.peek([light2, light3], count=100)
assert not all(result[0] == 0 for result in results)
assert (result[0] == result[1] for result in results)
# S_1 + S_2 - S_12 = 0 + 1 - 1 = 0
assert round(board.measure_entanglement(light1, light2), 1) == 0.0
# S_1 + S_2 - S_12 = 0 + 1 - 1 = 0
assert round(board.measure_entanglement(light1, light3), 1) == 0.0
# S_1 + S_2 - S_12 = 1 + 1 - 0 = 2
assert round(board.measure_entanglement(light2, light3), 1) == 2.0
Loading