diff --git a/qdna/compression/schmidt.py b/qdna/compression/schmidt.py index 26833d3..2b31647 100644 --- a/qdna/compression/schmidt.py +++ b/qdna/compression/schmidt.py @@ -94,15 +94,11 @@ def __init__(self, params, label=None, opt_params=None): self.svd = "auto" if opt_params.get("svd") is None else \ opt_params.get("svd") - # The trash and latent qubits must take into account that the qiskit qubits are reversed. - complement = sorted( - set(range(self.num_qubits)).difference(set(self.partition)) - )[::-1] - self.latent_qubits = [ - self.num_qubits-i-1 for i in complement - ] - self.trash_qubits = sorted( - set(range(self.num_qubits)).difference(set(self.latent_qubits)) + # The trash and latent qubits must take into account that the qiskit + # qubits are reversed. See line `return circuit.reverse_bits()`. + self.trash_qubits = sorted([self.num_qubits-i-1 for i in self.partition]) + self.latent_qubits = sorted( + set(range(self.num_qubits)).difference(set(self.trash_qubits)) ) if label is None: @@ -120,17 +116,18 @@ def _define_initialize(self): circuit.initialize(self.params) return circuit.inverse() + # reg_a = trash register, reg_b = latent register. circuit, reg_a, reg_b = self._create_quantum_circuit() - # Schmidt decomposition + # Schmidt decomposition. rank, svd_u, _, svd_v = schmidt_decomposition( self.params, reg_a, rank=self.low_rank, svd=self.svd ) - # Schmidt measure of entanglement + # Schmidt measure of entanglement. e_bits = _to_qubits(rank) - # Phase 3 and 4 encode gates U and V.T + # Phase 3 and 4 encode gates U and V.T. self._encode(svd_u, circuit, reg_b) self._encode(svd_v.T, circuit, reg_a) diff --git a/qdna/graph.py b/qdna/graph.py new file mode 100644 index 0000000..6bb5c7b --- /dev/null +++ b/qdna/graph.py @@ -0,0 +1,209 @@ +# Copyright 2023 qdna-lib project. + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import itertools +import random +from dwave.samplers import SimulatedAnnealingSampler + +def min_cut_fixed_size_optimal(graph, size_a, size_b): + ''' + Optimal solution to the minimum cut problem with fixed sizes for the sets. + + O(n! / k!(n-k)!) x O(m), where m is the number of edges in the graph. + The total number of edges m in a complete graph with n nodes is given by + m = n(n-1) / 2 + ''' + nodes = list(graph.nodes()) + min_cut_weight = float('inf') + min_cut_partition = (set(), set()) + + # Iterate over all combinations for set A + for nodes_a in itertools.combinations(nodes, size_a): + set_a = set(nodes_a) + set_b = set(nodes) - set_a + + # Ensure the size of set B is as required + if len(set_b) != size_b: + continue + + # Calculate the sum of weights of edges between the two sets + cut_weight = sum( + graph[u][v]['weight'] for u in set_a for v in set_b if graph.has_edge(u, v) + ) + + # Update min cut if a lower weight is found + if cut_weight < min_cut_weight: + min_cut_weight = cut_weight + min_cut_partition = (set_a, set_b) + + return min_cut_partition, min_cut_weight + +def min_cut_fixed_size_heuristic(graph, size_a, size_b): + ''' + Heuristic approach for the Min-Cut problem with a fixed number of nodes in + each partition. + + This approach aims to find a partition of the graph where the total edge + weight crossing the cut is as low as possible, given the size constraints. + The algorithm iteratively attempts to reduce the total weight of the cut by + swapping nodes between sets A and B, provided the swap decreases the total + weight of the cut. + This is a heuristic approach and may not always find the globally optimal + solution, especially for complex or large graphs. + + O(k * n_a * n_b * m): + k: number of iterations (vary significantly based on the graph's + structure and the initial partitioning). + n_a: subsystem A number of qubits. + n_b: subsystem B number of qubits. + m: number of edges between a node and subsystem B (typically equal to n_b). + + Example worst-case scenario (n_a=n_b=n/2): + O(k * n^3) + ''' + nodes = list(graph.nodes()) + random.shuffle(nodes) # Shuffle nodes to randomize initial selection + + # Initialize sets A and B with random nodes + set_a = set(nodes[:size_a]) + set_b = set(nodes[size_a:size_a + size_b]) + + def swapping_weight(node, other_node, set_node, set_other_node): + # Entanglement (without - with) swapping. + # A positive result indicates that a swap should be made + # to reduce entanglement. That is, the entanglement with + # a swap is smaller than without. + return \ + sum(graph[node][neighbor]['weight'] + for neighbor in set_other_node if neighbor is not other_node) - \ + sum(graph[node][neighbor]['weight'] + for neighbor in set_node if neighbor is not node) + + # Iteratively try to improve the cut by swapping nodes between A and B + improved = True + while improved: + improved = False + for node_a in set_a: + for node_b in set_b: + weight_a = swapping_weight(node_a, node_b, set_a, set_b) + weight_b = swapping_weight(node_b, node_a, set_b, set_a) + total_weight = weight_a + weight_b + # Here, the entanglements of both nodes are + # considered simultaneously. + if total_weight > 0: + # Swap the nodes if the total entanglement with the swap is + # smaller. In other words, the entanglement without + # swapping is greater, which means that `total_weight` is + # positive. + set_a.remove(node_a) + set_b.remove(node_b) + set_a.add(node_b) + set_b.add(node_a) + improved = True + break + + if improved: + break + + # Sorts the sets when the subsystem sizes are equal + if size_a == size_b and sorted(set_a)[0] > sorted(set_b)[0]: + set_a, set_b = set_b, set_a + + # Calculate the sum of weights of edges between the two sets + cut_weight = sum( + graph[u][v]['weight'] for u in set_a for v in set_b if graph.has_edge(u, v) + ) + + return (set_a, set_b), cut_weight + + +# D-Wave functions + +def graph_to_qubo(graph): + ''' + Map the graph to a QUBO model. + ''' + qubo = {} + + max_weight = max(weight for _, _, weight in graph.edges(data='weight')) + + for i, j, weight in graph.edges(data='weight'): + # Set qubo_{ij} to be `max_weight - weight` for min-cut. + weight = max_weight - weight + qubo[(i, i)] = qubo.get((i, i), 0) - weight + qubo[(j, j)] = qubo.get((j, j), 0) - weight + qubo[(i, j)] = qubo.get((i, j), 0) + 2 * weight + + return qubo + +def graph_to_ising(graph): + ''' + Map the graph to an Ising model. + ''' + ising = {} + local_field = {} + + max_weight = max(weight for _, _, weight in graph.edges(data='weight')) + + for i, j, weight in graph.edges(data='weight'): + # Set qubo_{ij} to be `max_weight - weight` for min-cut. + weight = max_weight - weight + ising[(i, j)] = weight + + # Add interaction strengths and local fields. + for i in graph.nodes(): + local_field[i] = 0 # local fields + + return ising, local_field + +def min_cut_dwave(graph, size_a=None, sample_method='ising', sampler=SimulatedAnnealingSampler(), num_reads=100): + + if sample_method == 'qubo': + # Map the graph to a QUBO model. + qubo = graph_to_qubo(graph) + response = sampler.sample_qubo(qubo, num_reads=num_reads) + else: + # Map the graph to an Ising model. + ising, local_field = graph_to_ising(graph) + response = sampler.sample_ising(local_field, ising, num_reads=num_reads) + print(response) + # Find the sample with the lowest energy (best result), respecting the size of the block. + min_cut_weight = float('inf') + set_a = None + if size_a is not None: + for sample, energy in response.data(['sample', 'energy']): + if energy < min_cut_weight and sum(v for v in sample.values() if v==1) == size_a: + min_cut_weight = energy + set_a = sample + else: + set_a = response.first.sample + min_cut_weight = response.first.energy + + # Subsystem B is the complement of subsystem A. + nodes = set(graph.nodes()) + set_a = {k for k, v in set_a.items() if v == 1} + size_a = len(set_a) + set_b = nodes - set_a + size_b = len(set_b) + + # Sorts the sets when the subsystem sizes are equal. + if size_a == size_b and sorted(set_a)[0] > sorted(set_b)[0]: + set_a, set_b = set_b, set_a + + # Calculate the sum of weights of edges between the two sets. + cut_weight = sum( + graph[u][v]['weight'] for u in set_a for v in set_b if graph.has_edge(u, v) + ) + + return (set_a, set_b), cut_weight diff --git a/qdna/quantum_info.py b/qdna/quantum_info.py new file mode 100644 index 0000000..9eed232 --- /dev/null +++ b/qdna/quantum_info.py @@ -0,0 +1,149 @@ +# Copyright 2023 qdna-lib project. + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import itertools +import networkx as nx +from qiskit.quantum_info import Statevector, partial_trace +import numpy as np + +def von_neumann_entropy(rho): + ''' + Compute the Von Neumann entropy (entanglement measure). + + To calculate the entropies, it is convenient to calculate the + eigendecomposition of `rho` and use the eigenvalues `lambda_i` to determine + the entropy: + `S(rho) = -sum_i( lambda_i * ln(lambda_i) )` + ''' + evals = np.real(np.linalg.eigvals(rho.data)) + return -np.sum([e * np.log2(e) for e in evals if 0 < e < 1]) + +def concurrence(rho): + ''' + Concurrence specifically quantifies the quantum entanglement between the + two qubits. It does not consider classical correlations. + + https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.80.2245 + ''' + # Compute the spin-flipped state + sigma_y = np.array([[0, -1j], [1j, 0]]) + rho_star = np.conj(rho) + rho_tilde = np.kron(sigma_y, sigma_y) @ rho_star @ np.kron(sigma_y, sigma_y) + + # Calculate the eigenvalues of the product matrix + eigenvalues = np.linalg.eigvals(rho @ rho_tilde) + # Sort in decreasing order + eigenvalues = np.sort(np.sqrt(np.abs(eigenvalues)))[::-1] + + # Compute the concurrence + return max(0, eigenvalues[0] - sum(eigenvalues[1:])) + +def mutual_information(rho_a, rho_b, rho_ab): + ''' + Mutual information quantifies the total amount of correlation between two + qubits. It includes both classical and quantum correlations. + ''' + + # Compute the Von Neumann entropy for each density matrix + s_a = von_neumann_entropy(rho_a) + s_b = von_neumann_entropy(rho_b) + s_ab = von_neumann_entropy(rho_ab) + + # Calculate the mutual information + return s_a + s_b - s_ab + +def correlation(state_vector, set_a, set_b, correlation_measure=mutual_information): + ''' + Compute the correlation between subsystems A and B. + ''' + + if (len(set_a) > 1 or len(set_b) > 1) and correlation_measure is concurrence: + raise ValueError( + "The value of `correlation_measure` cannot be `concurrence` when " + "`len(set_a) > 1` or `len(set_b) > 1`. Choose, for example, " + "`mutual_information` instead." + ) + + psi = Statevector(state_vector) + + # Compute the reduced density matrix for the union of the two sets. + set_ab = set_a.union(set_b) + rho_ab = partial_trace(psi, list(set(range(psi.num_qubits)).difference(set_ab))) + + # Maintains the relative position between the qubits of the two subsystems. + new_set_a = [sum(i < item for i in set_ab) for item in set_a] + new_set_b = [sum(i < item for i in set_ab) for item in set_b] + + # Calculate the reduced density matrice for each set. + rho_a = partial_trace(rho_ab, new_set_b) + rho_b = partial_trace(rho_ab, new_set_a) + + if correlation_measure is mutual_information: + return correlation_measure(rho_a, rho_b, rho_ab) + + return correlation_measure(rho_ab) + +def correlation_graph(state_vector, + n_qubits, + min_set_size=1, + max_set_size=1, + correlation_measure=mutual_information +): + ''' + Initialize a graph where nodes represent qubits and the weights represent + the entanglement between pairs of qubits in a register of `n` qubits for a + pure state. + + O(n^2) x O(2^n) + ''' + if n_qubits <= max_set_size <= 0: + raise ValueError( + f"The value of `max_set_size` [{max_set_size}] must be greater " + f"than zero and less than `n_qubits` [{n_qubits}]." + ) + + if max_set_size < min_set_size <= 0: + raise ValueError( + f"The value of `min_set_size` [{min_set_size}] must be greater " + f"than zero and less or equal than `max_set_size` " + f"[{max_set_size}]." + ) + + # Create a graph. + graph = nx.Graph() + + # Add nodes for each set of qubits up to the size `max_set_size`. + for set_size in range(min_set_size, max_set_size + 1): + for qubit_set in itertools.combinations(range(n_qubits), set_size): + graph.add_node(qubit_set) + + # Add edges with weights representing entanglement. + for node_a in graph.nodes(): + set_a = set(node_a) + for node_b in graph.nodes(): + set_b = set(node_b) + # Ensure non-overlapping sets. + if not set_a.intersection(set_b): + # Compute the correlation betweem subsystems. + weight = correlation( + state_vector, + set_a, + set_b, + correlation_measure=correlation_measure + ) + + # Add an edge with the shared info as weight. + graph.add_edge(node_a, node_b, weight=weight) + + return graph diff --git a/requirements.txt b/requirements.txt index cefb2de..af847b3 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,3 +4,4 @@ qiskit-aer qiskit-algorithms torch qclib +dwave-ocean-sdk diff --git a/test/test_graph.py b/test/test_graph.py new file mode 100644 index 0000000..851cbcb --- /dev/null +++ b/test/test_graph.py @@ -0,0 +1,117 @@ +# Copyright 2023 qdna-lib project. + +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at + +# http://www.apache.org/licenses/LICENSE-2.0 + +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" +Tests for the entanglement.py module. +""" + +from unittest import TestCase +from math import isclose +import numpy as np +from qdna.quantum_info import correlation_graph, concurrence +from qdna.graph import min_cut_fixed_size_optimal, \ + min_cut_fixed_size_heuristic + + +# pylint: disable=missing-function-docstring +# pylint: disable=missing-class-docstring + + +class TestEntanglement(TestCase): + + def test_min_cut_optimal_product_state(self): + n_qubits = 3 + + state_vector1 = np.random.rand(2**n_qubits) + np.random.rand(2**n_qubits) * 1j + state_vector1 = state_vector1 / np.linalg.norm(state_vector1) + + state_vector2 = np.random.rand(2**n_qubits) + np.random.rand(2**n_qubits) * 1j + state_vector2 = state_vector2 / np.linalg.norm(state_vector2) + + state_vector = np.kron(state_vector1, state_vector2) + + graph = correlation_graph(state_vector, 6, correlation_measure=concurrence) + (set_a, set_b), cut_weight = min_cut_fixed_size_optimal(graph, 3, 3) + + self.assertTrue(isclose(cut_weight, 0.0)) + self.assertTrue(np.allclose(sorted(sum(set_a, ())), [0, 1, 2])) + self.assertTrue(np.allclose(sorted(sum(set_b, ())), [3, 4, 5])) + + def test_min_cut_heuristic_product_state(self): + n_qubits = 3 + + state_vector1 = np.random.rand(2**n_qubits) + np.random.rand(2**n_qubits) * 1j + state_vector1 = state_vector1 / np.linalg.norm(state_vector1) + + state_vector2 = np.random.rand(2**n_qubits) + np.random.rand(2**n_qubits) * 1j + state_vector2 = state_vector2 / np.linalg.norm(state_vector2) + + state_vector = np.kron(state_vector1, state_vector2) + + graph = correlation_graph(state_vector, 6, correlation_measure=concurrence) + (set_a, set_b), cut_weight = min_cut_fixed_size_heuristic(graph, 3, 3) + + self.assertTrue(isclose(cut_weight, 0.0)) + self.assertTrue(np.allclose(sorted(sum(set_a, ())), [0, 1, 2])) + self.assertTrue(np.allclose(sorted(sum(set_b, ())), [3, 4, 5])) + + def test_min_cut_optimal_fixed(self): + state_vector = [0.00000000e+00, 6.94991284e-04, 6.34061054e-02, 2.13286776e-01, + 2.05658826e-01, 8.15141431e-02, 8.40762648e-03, 0.00000000e+00, + 0.00000000e+00, 9.96282923e-03, 1.59787371e-01, 2.46020212e-01, + 2.40154124e-01, 1.88551405e-01, 1.90979862e-02, 0.00000000e+00, + 9.31053205e-04, 4.54430541e-02, 2.03018262e-01, 1.86283916e-01, + 1.52698965e-01, 1.86716850e-01, 3.98788754e-02, 0.00000000e+00, + 9.31053205e-04, 7.52123527e-02, 2.05815792e-01, 1.50703564e-01, + 1.29600833e-01, 1.44311684e-01, 7.06191583e-02, 0.00000000e+00, + 0.00000000e+00, 7.69075128e-02, 1.73599511e-01, 1.17960532e-01, + 1.25934109e-01, 1.33062442e-01, 8.39199837e-02, 0.00000000e+00, + 0.00000000e+00, 3.84563398e-02, 1.76000915e-01, 1.11207757e-01, + 1.36568022e-01, 1.59108898e-01, 6.19303017e-02, 0.00000000e+00, + 0.00000000e+00, 9.01723392e-03, 1.70193955e-01, 1.96598638e-01, + 2.21354420e-01, 1.95836976e-01, 4.14572396e-02, 6.39589243e-03, + 0.00000000e+00, 2.11959665e-04, 6.15966561e-02, 2.17282223e-01, + 2.49816212e-01, 1.27759885e-01, 2.73632167e-02, 1.20487999e-02] + + graph = correlation_graph(state_vector, 6, correlation_measure=concurrence) + (set_a, set_b), cut_weight = min_cut_fixed_size_optimal(graph, 3, 3) + + self.assertTrue(isclose(cut_weight, 0.0)) + self.assertTrue(np.allclose(sorted(sum(set_a, ())), [0, 1, 2])) + self.assertTrue(np.allclose(sorted(sum(set_b, ())), [3, 4, 5])) + + def test_min_cut_heuristic_fixed(self): + state_vector = [0.00000000e+00, 6.94991284e-04, 6.34061054e-02, 2.13286776e-01, + 2.05658826e-01, 8.15141431e-02, 8.40762648e-03, 0.00000000e+00, + 0.00000000e+00, 9.96282923e-03, 1.59787371e-01, 2.46020212e-01, + 2.40154124e-01, 1.88551405e-01, 1.90979862e-02, 0.00000000e+00, + 9.31053205e-04, 4.54430541e-02, 2.03018262e-01, 1.86283916e-01, + 1.52698965e-01, 1.86716850e-01, 3.98788754e-02, 0.00000000e+00, + 9.31053205e-04, 7.52123527e-02, 2.05815792e-01, 1.50703564e-01, + 1.29600833e-01, 1.44311684e-01, 7.06191583e-02, 0.00000000e+00, + 0.00000000e+00, 7.69075128e-02, 1.73599511e-01, 1.17960532e-01, + 1.25934109e-01, 1.33062442e-01, 8.39199837e-02, 0.00000000e+00, + 0.00000000e+00, 3.84563398e-02, 1.76000915e-01, 1.11207757e-01, + 1.36568022e-01, 1.59108898e-01, 6.19303017e-02, 0.00000000e+00, + 0.00000000e+00, 9.01723392e-03, 1.70193955e-01, 1.96598638e-01, + 2.21354420e-01, 1.95836976e-01, 4.14572396e-02, 6.39589243e-03, + 0.00000000e+00, 2.11959665e-04, 6.15966561e-02, 2.17282223e-01, + 2.49816212e-01, 1.27759885e-01, 2.73632167e-02, 1.20487999e-02] + + graph = correlation_graph(state_vector, 6, correlation_measure=concurrence) + (set_a, set_b), cut_weight = min_cut_fixed_size_heuristic(graph, 3, 3) + + self.assertTrue(isclose(cut_weight, 0.0)) + self.assertTrue(np.allclose(sorted(sum(set_a, ())), [0, 1, 2])) + self.assertTrue(np.allclose(sorted(sum(set_b, ())), [3, 4, 5]))