Merge pull request #7 from qDNA-yonsei/ae-feature-map

Updates the Amplitude Encoding feature map

qdna/embedding/__init__.py

@@ -4,3 +4,4 @@
 from qdna.embedding.nqe_z_feature_map import NqeZFeatureMap
 from qdna.embedding.nqe_zz_feature_map import NqeZZFeatureMap
 from qdna.embedding.nqe_iqp_feature_map import NqeIqpFeatureMap
+from qdna.embedding.nqe_ae_feature_map import NqeAeFeatureMap

qdna/embedding/ae_feature_map.py

@@ -19,6 +19,8 @@
 
 from typing import Sequence, Mapping
 
+from math import pi
+
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.quantumregister import QuantumRegister
 from qiskit.circuit import Parameter, ParameterVector, ParameterExpression
@@ -27,7 +29,7 @@
 
 from qclib.state_preparation.util.tree_utils import children
 
-from qdna.embedding.util.state_tree_preparation import state_decomposition
+from qdna.embedding.util.state_tree_preparation import state_decomposition, _sign
 from qdna.embedding.util.angle_tree_preparation import create_angles_tree
 from qdna.embedding.util.ucr import ucr
 
@@ -44,7 +46,9 @@ def __init__(
         reps: int = 1,
         insert_barriers: bool = False,
         parameter_prefix: str = "x",
-        name: str | None = "AeFeatureMap"
+        name: str | None = "AeFeatureMap",
+        normalize: bool = False,
+        set_global_phase: bool = False
     ) -> None:
 
         super().__init__(name=name)
@@ -56,6 +60,8 @@ def __init__(
         self._initial_state: QuantumCircuit | None = None
         self._initial_state_circuit: QuantumCircuit | None = None
         self._bounds: list[tuple[float | None, float | None]] | None = None
+        self._normalize = normalize
+        self._set_global_phase = set_global_phase
 
         if int(reps) != reps:
             raise TypeError("The value of reps should be int.")
@@ -270,7 +276,7 @@ def _build(self) -> None:
         circuit = QuantumCircuit(*self.qregs, name=self.name)
         params = ParameterVector(self._parameter_prefix, length=self.num_parameters_settable)
 
-        state_tree = state_decomposition(self.num_qubits, params)
+        state_tree = state_decomposition(self.num_qubits, params, normalize=self._normalize)
         angle_tree = create_angles_tree(state_tree)
 
         nodes = [angle_tree]
@@ -280,13 +286,21 @@ def _build(self) -> None:
         while len(nodes) > 0:
             angles = [node.angle_y for node in nodes]
             ucry = ucr(RYGate, angles)
-            circuit.compose(ucry, [target_qubit] + control_qubits[::-1], inplace=True)
+            circuit.compose(
+                ucry, [target_qubit] + control_qubits[::-1], inplace=True, wrap=False
+            )
             control_qubits.append(target_qubit)
             nodes = children(nodes)
             target_qubit -= 1
 
             if self._insert_barriers and target_qubit >= 0:
                 circuit.barrier()
 
+        # As we only have `pi` phases (negative sign of a real number),
+        # the global phase only depends on the first coordinate of the
+        # state vector.
+        if self._set_global_phase:
+            circuit.global_phase = ((_sign(params[0])-1) / -2) * pi
+
         for _ in range(self._reps):
-            self.compose(circuit, self.qubits, inplace=True)
+            self.compose(circuit, self.qubits, inplace=True, wrap=False)

qdna/embedding/nqe_base.py

@@ -75,16 +75,16 @@ def fit(self, X, Y, batch_size=25, iters=100, optimizer=None, loss_func=None, di
         if distance == 'fidelity': # only 'fidelity' for now.
             model = _NetFidelity(qnn, self.nn)
 
-        if optimizer == None:
+        if optimizer is None:
             optimizer = optim.Adam(model.parameters(), lr=0.01)
-        if loss_func == None:
+        if loss_func is None:
             loss_func = MSELoss()
 
         # Start training
         loss_list = []  # Store loss history
         model.train()   # Set model to training mode
 
-        for iter in range(iters):
+        for i in range(iters):
             X1_batch, X2_batch, Y_batch = self.new_data(batch_size, X, Y) # Random sampling of data.
 
             optimizer.zero_grad(set_to_none=True)  # Initialize gradient
@@ -97,11 +97,9 @@ def fit(self, X, Y, batch_size=25, iters=100, optimizer=None, loss_func=None, di
 
             if verbose:
                 print(
-                    "Training [{:.0f}%]\tLoss: {:.4f}\tAvg. Loss: {:.4f}".format(
-                        100.0 * (iter + 1) / iters,
-                        loss_list[-1],
-                        sum(loss_list) / len(loss_list)
-                    )
+                    f"Training [{100.0 * (i + 1) / iters:.0f}%]\t"
+                    f"Loss: {loss_list[-1]:.4f}\t"
+                    f"Avg. Loss: {sum(loss_list) / len(loss_list):.4f}"
                 )
 
         self.model = model

qdna/embedding/util/state_tree_preparation.py

@@ -17,9 +17,10 @@
 """
 
 from dataclasses import dataclass
+import numpy as np
 from qiskit.circuit import ParameterExpression
 from qiskit.utils import optionals as _optionals
-from sympy import sqrt as _sqrt
+from sympy import sqrt as sp_sqrt
 import symengine
 
 @dataclass
@@ -33,42 +34,89 @@ class Node:
     left: "Node"
     right: "Node"
     norm: ParameterExpression
+    sign: ParameterExpression
 
     def __str__(self):
         return (
             f"{self.level}_"
             f"{self.index}\n"
-            f"{self.norm}"
+            f"{self.sign*self.norm}"
         )
 
+# def _abs(self):
+#     # SymPy's `Abs` function doesn't work with hybrid optimization backwards.
+#     """Absolute value of a ParameterExpression"""
+#     if _optionals.HAS_SYMENGINE:
+#         import symengine
+#         return ParameterExpression(self._parameter_symbols, symengine.Abs(self._symbol_expr))
+#     else:
+#         from sympy import Abs as sp_abs
+#         return ParameterExpression(self._parameter_symbols, sp_abs(self._symbol_expr))
+
 def _sqrt(self):
     """Square root of a ParameterExpression"""
     if _optionals.HAS_SYMENGINE:
         return self._call(symengine.sqrt)
-    else:
-        return self._call(_sqrt)
 
-def state_decomposition(nqubits, data):
+    return self._call(sp_sqrt)
+
+def _sign(self):
+    # SymPy's `sign` function doesn't work with hybrid optimization backwards.
+    """Sign of a ParameterExpression"""
+    # if _optionals.HAS_SYMENGINE:
+    #     import symengine
+    #     return ParameterExpression(
+    #         self._parameter_symbols, symengine.sign(self._symbol_expr + 1e-4)
+    #     )
+    # else:
+    #     from sympy import sign
+    #     return ParameterExpression(
+    #         self._parameter_symbols, sign(self._symbol_expr + 1e-4)
+    #     )
+
+    # If self is exactly `-10^-4`, the algorithm is interrupted.
+    # This is because the `sign` function will return 0, zeroing the norms.
+    return self / _sqrt(self*self)
+
+def state_decomposition(nqubits, data, normalize=False):
     """
     :param nqubits: number of qubits required to generate a
                     state with the same length as the data vector (2^nqubits)
     :param data: list with exactly 2^nqubits pairs (index, amplitude)
     :return: root of the state tree
     """
+
     new_nodes = []
 
     # leafs
+    r = 10**-6
     for i, k in enumerate(data):
+        # If one of the coordinates of the state vector
+        # is exactly `-r*10^-6`, the algorithm breaks.
+        # This is because one of the norms will be zero,
+        # causing division by zero.
+        k = k + r # prevents division by zero.
+        sign = _sign(k)
         new_nodes.append(
             Node(
                 i,
                 nqubits,
                 None,
                 None,
-                k
+                k,
+                sign
             )
         )
 
+    # normalization
+    if normalize:
+        norm = 0.0
+        for i, k in enumerate(data):
+            norm += k*k
+        norm = _sqrt(norm)
+        for node in new_nodes:
+            node.norm = node.norm/norm
+
     # build state tree
     while nqubits > 0:
         nodes = new_nodes
@@ -81,16 +129,30 @@ def state_decomposition(nqubits, data):
             # The value of the first node is always 1, so there is no need to
             # calculate it.
             new_nodes.append(
-                Node(nodes[k].index // 2, nqubits, nodes[k], nodes[k + 1], 1.0)
+                Node(
+                    nodes[k].index // 2,
+                    nqubits,
+                    nodes[k],
+                    nodes[k + 1],
+                    nodes[k].sign,
+                    nodes[k].sign
+                )
             )
         else:
             while k < n_nodes:
-                norm = _sqrt(
+                norm = nodes[k].sign * _sqrt(
                     nodes[k].norm*nodes[k].norm + nodes[k + 1].norm*nodes[k + 1].norm
                 )
 
                 new_nodes.append(
-                    Node(nodes[k].index // 2, nqubits, nodes[k], nodes[k + 1], norm)
+                    Node(
+                        nodes[k].index // 2,
+                        nqubits,
+                        nodes[k],
+                        nodes[k + 1],
+                        norm,
+                        nodes[k].sign
+                    )
                 )
                 k = k + 2
 

test/embedding/test_ae_feature_map.py

@@ -28,17 +28,88 @@
 
 class TestAeFeatureMap(TestCase):
 
-    def test_state(self):
+    def test_rnd_state(self):
         n_qubits = 6
-        state_vector = np.random.rand(2**n_qubits)
+        state_vector = (np.random.rand(2**n_qubits)-0.5)*2
         state_vector = state_vector / np.linalg.norm(state_vector)
 
-        feature_map = AeFeatureMap(n_qubits)
+        feature_map = AeFeatureMap(n_qubits, normalize=False, set_global_phase=True)
 
         circuit = QuantumCircuit(n_qubits)
         circuit.compose(feature_map, list(range(n_qubits)), inplace=True)
         circuit.assign_parameters(state_vector, inplace=True)
 
         state = Statevector(circuit)
 
-        self.assertTrue(np.allclose(state_vector, state))
+        self.assertTrue(np.allclose(state_vector, state, rtol=1e-03, atol=1e-05))
+
+    def test_non_normalized_state(self):
+        # This test is important because during NQE training
+        # the parameter vector is not normalized.
+        n_qubits = 4
+        state_vector = [
+            -3.1416, 0.0000, -0.7667, 0.5835, 0.3468, 0.5320, 0.6706, 0.6756,
+            0.5270, 0.7115, 0.5165, 0.7194, -0.5392, 0.3878, 0.5628, 0.7304
+        ]
+
+        feature_map = AeFeatureMap(n_qubits, normalize=True, set_global_phase=True)
+
+        circuit = QuantumCircuit(n_qubits)
+        circuit.compose(feature_map, list(range(n_qubits)), inplace=True)
+        circuit.assign_parameters(state_vector, inplace=True)
+
+        state = Statevector(circuit)
+
+        state_vector = state_vector / np.linalg.norm(state_vector)
+        self.assertTrue(np.allclose(state_vector, state, rtol=1e-03, atol=1e-05))
+
+    def test_sparse_state(self):
+        # Ensure that embedding is preventing division by zero.
+        n_qubits = 4
+        state_vector = [
+            0.0000, 0.0000, 0.0000, 0.5835, 0.0000, 0.5320, 0.6706, 0.6756,
+            0.5270, 0.0000, 0.0000, 0.7194, 0.0000, 0.0000, 0.0000, 0.7304
+        ]
+
+        feature_map = AeFeatureMap(n_qubits, normalize=True, set_global_phase=False)
+
+        circuit = QuantumCircuit(n_qubits)
+        circuit.compose(feature_map, list(range(n_qubits)), inplace=True)
+        circuit.assign_parameters(state_vector, inplace=True)
+
+        state = Statevector(circuit)
+
+        state_vector = state_vector / np.linalg.norm(state_vector)
+        self.assertTrue(np.allclose(state_vector, state, rtol=1e-03, atol=1e-05))
+
+    def test_fixed_state(self):
+        n_qubits = 6
+        state_vector = [-0.00750879, -0.19509541,  0.13329143, -0.05016399,
+                         0.15207381,  0.07559872,  0.06584141,  0.02580133,
+                         0.07838591,  0.1797185,  -0.08671011,  0.12130839,
+                        -0.01520751, -0.16350927, -0.14340019,  0.06187,
+                        -0.15696974, -0.07004782, -0.10200592, -0.13275576,
+                         0.03975269, -0.06866376, -0.19666519, -0.19759358,
+                        -0.06338929,  0.01957956, -0.04101603, -0.17654797,
+                         0.13379771,  0.13126602, -0.06316672, -0.16629322,
+                         0.19740418, -0.01161936,  0.00660316,  0.18850766,
+                         0.05001657,  0.05713368, -0.10440029, -0.19293127,
+                        -0.21320381, -0.15065956,  0.20290586, -0.03929146,
+                         0.0254895,   0.03343766, -0.17821551, -0.1903655,
+                         0.13099242, -0.11116633,  0.15833651, -0.02708352,
+                        -0.03610143,  0.20473973,  0.03146476, -0.03039286,
+                        -0.14796486,  0.00971746, -0.15238775,  0.12996466,
+                        -0.19630254, -0.08859373, -0.12830557,  0.12446281]
+
+#        state_vector = state_vector / np.linalg.norm(state_vector)
+
+        feature_map = AeFeatureMap(n_qubits, normalize=False, set_global_phase=True)
+
+        circuit = QuantumCircuit(n_qubits)
+        circuit.compose(feature_map, list(range(n_qubits)), inplace=True)
+        circuit.assign_parameters(state_vector, inplace=True)
+
+        state = Statevector(circuit)
+
+        state_vector = state_vector / np.linalg.norm(state_vector)
+        self.assertTrue(np.allclose(state_vector, state, rtol=1e-03, atol=1e-05))