Gate Set Tomography handling of initial state and measurement operator (

#408)

qiskit/ignis/verification/tomography/fitters/gateset_fitter.py

@@ -114,13 +114,17 @@ def linear_inversion(self) -> Dict[str, PTM]:
         n = len(self.gateset_basis.spam_labels)
         m = len(self.gateset_basis.gate_labels)
         gram_matrix = np.zeros((n, n))
+        E = np.zeros((1, n))
+        rho = np.zeros((n, 1))
         gate_matrices = []
         for i in range(m):
             gate_matrices.append(np.zeros((n, n)))
 
         for i in range(n):  # row
+            F_i = self.gateset_basis.spam_labels[i]
+            E[0][i] = self.probs[(F_i,)]
+            rho[i][0] = self.probs[(F_i,)]
             for j in range(n):  # column
-                F_i = self.gateset_basis.spam_labels[i]
                 F_j = self.gateset_basis.spam_labels[j]
                 gram_matrix[i][j] = self.probs[(F_j, F_i)]
 
@@ -131,7 +135,10 @@ def linear_inversion(self) -> Dict[str, PTM]:
         gram_inverse = np.linalg.inv(gram_matrix)
 
         gates = [PTM(gram_inverse @ gate_matrix) for gate_matrix in gate_matrices]
-        return dict(zip(self.gateset_basis.gate_labels, gates))
+        result = dict(zip(self.gateset_basis.gate_labels, gates))
+        result['E'] = E
+        result['rho'] = gram_inverse @ rho
+        return result
 
     def _default_init_state(self, size):
         """Returns the PTM representation of the usual ground state"""
@@ -145,6 +152,13 @@ def _default_measurement_op(self, size):
             return np.array([[np.sqrt(0.5), 0, 0, np.sqrt(0.5)]])
         raise RuntimeError("No default measurement op for more than 1 qubit")
 
+    def _ideal_gateset(self, size):
+        ideal_gateset = {label: PTM(self.gateset_basis.gate_matrices[label])
+                         for label in self.gateset_basis.gate_labels}
+        ideal_gateset['E'] = self._default_measurement_op(size)
+        ideal_gateset['rho'] = self._default_init_state(size)
+        return ideal_gateset
+
     def fit(self) -> Dict:
         """
         Reconstruct a gate set from measurement data using optimization.
@@ -163,39 +177,31 @@ def fit(self) -> Dict:
         """
         linear_inversion_results = self.linear_inversion()
         n = len(self.gateset_basis.spam_labels)
-        E = self._default_measurement_op(n)
-        rho = self._default_init_state(n)
-        Gs = [self.gateset_basis.gate_matrices[label]
-              for label in self.gateset_basis.gate_labels]
-        Fs = [self.gateset_basis.spam_matrix(label)
-              for label in self.gateset_basis.spam_labels]
-        Gs_E = [linear_inversion_results[label].data
-                for label in self.gateset_basis.gate_labels]
-        gauge_opt = GaugeOptimize(Gs, Gs_E, Fs, rho)
-        Gs_E = gauge_opt.optimize()
+        gauge_opt = GaugeOptimize(self._ideal_gateset(n),
+                                  linear_inversion_results,
+                                  self.gateset_basis)
+        past_gauge_gateset = gauge_opt.optimize()
         optimizer = GST_Optimize(self.gateset_basis.gate_labels,
                                  self.gateset_basis.spam_labels,
                                  self.gateset_basis.spam_spec,
                                  self.probs)
-        optimizer.set_initial_value(E, rho, Gs_E)
+        optimizer.set_initial_value(past_gauge_gateset)
         optimization_results = optimizer.optimize()
         return optimization_results
 
 
 class GaugeOptimize():
     def __init__(self,
-                 Gs: List[np.array],
-                 initial_Gs_E: List[np.array],
-                 Fs: List[np.array],
-                 rho: np.array
+                 ideal_gateset: Dict[str, PTM],
+                 initial_gateset: Dict[str, PTM],
+                 gateset_basis: GateSetBasis,
                  ):
         """Initialize gauge optimizer fitter with the ideal and expected
             outcomes.
         Args:
-            Gs: The ideal expected gate matrices
-            initial_Gs_E: The experimentally-obtained gate approximations.
-            Fs: The SPAM circuit matrices
-            rho: The system's initial value (usually the |0> qubits)
+            ideal_gateset: The ideal expected gate matrices
+            initial_gateset: The experimentally-obtained gate approximations.
+            gateset_basis: The gateset data
 
         Additional information:
             Gauge optimization aims to find a basis in which the tomography
@@ -212,14 +218,16 @@ def __init__(self,
             difference between the gates found by experiment
             and the "expected" gates in the ideal (noiseless) case.
         """
-        self.Gs = Gs
-        self.Fs = Fs
-        self.initial_Gs_E = initial_Gs_E
-        self.d = np.shape(Gs[0])[0]
-        self.n = len(Gs)
-        self.rho = rho
-
-    def _x_to_Gs_E(self, x: np.array) -> List[np.array]:
+        self.gateset_basis = gateset_basis
+        self.ideal_gateset = ideal_gateset
+        self.initial_gateset = initial_gateset
+        self.Fs = [self.gateset_basis.spam_matrix(label)
+                   for label in self.gateset_basis.spam_labels]
+        self.d = np.shape(ideal_gateset['rho'])[0]
+        self.n = len(gateset_basis.gate_labels)
+        self.rho = ideal_gateset['rho']
+
+    def _x_to_gateset(self, x: np.array) -> Dict[str, PTM]:
         """Converts the gauge to the gateset defined by it
         Args:
             x: An array representation of the B matrix
@@ -236,9 +244,13 @@ def _x_to_Gs_E(self, x: np.array) -> List[np.array]:
         B = np.array(x).reshape((self.d, self.d))
         try:
             BB = np.linalg.inv(B)
-            return [BB @ G @ B for G in self.initial_Gs_E]
         except np.linalg.LinAlgError:
-            return np.inf
+            return None
+        gateset = {label: PTM(BB @ self.initial_gateset[label].data @ B)
+                   for label in self.gateset_basis.gate_labels}
+        gateset['E'] = self.initial_gateset['E'] @ B
+        gateset['rho'] = BB @ self.initial_gateset['rho']
+        return gateset
 
     def _obj_fn(self, x: np.array) -> float:
         """The norm-based score function for the gauge optimizer
@@ -249,10 +261,15 @@ def _obj_fn(self, x: np.array) -> float:
             The sum of norm differences between the ideal gateset
             and the one corresponding to B
         """
-        Gs_E = self._x_to_Gs_E(x)
-        return sum([np.linalg.norm(G - G_E)
-                    for (G, G_E)
-                    in zip(self.Gs, Gs_E)])
+        gateset = self._x_to_gateset(x)
+        result = sum([np.linalg.norm(gateset[label].data -
+                                     self.ideal_gateset[label].data)
+                      for label in self.gateset_basis.gate_labels])
+        result = result + np.linalg.norm(gateset['E'] -
+                                         self.ideal_gateset['E'])
+        result = result + np.linalg.norm(gateset['rho'] -
+                                         self.ideal_gateset['rho'])
+        return result
 
     def optimize(self) -> List[np.array]:
         """The main optimization method
@@ -261,7 +278,7 @@ def optimize(self) -> List[np.array]:
         """
         initial_value = np.array([(F @ self.rho).T[0] for F in self.Fs]).T
         result = opt.minimize(self._obj_fn, initial_value)
-        return self._x_to_Gs_E(result.x)
+        return self._x_to_gateset(result.x)
 
 
 def get_cholesky_like_decomposition(mat: np.array) -> np.array:
@@ -441,7 +458,7 @@ def _join_input_vector(self,
         result = E_vec + rho_vec
         for G_T in Gs_T:
             result += self._complex_matrix_to_vec(G_T)
-        return result
+        return np.array(result)
 
     def _obj_fn(self, x: np.array) -> float:
         """The MLE objective function
@@ -508,7 +525,7 @@ def _rho_trace(self, x: np.array) -> Tuple[float]:
         """
         _, rho, _ = self._split_input_vector(x)
         d = (2 ** self.qubits)  # rho is dxd and starts at variable d^2
-        rho = rho.reshape((d, d))
+        rho = self._convert_from_ptm(rho.reshape((d, d)))
         trace = sum([rho[i][i] for i in range(d)])
         return (np.real(trace), np.imag(trace))
 
@@ -644,17 +661,14 @@ def _process_result(self, x: np.array) -> Dict:
             result[self.Gs[i]] = PTM(G_matrices[i])
         return result
 
-    def set_initial_value(self,
-                          E: np.array,
-                          rho: np.array,
-                          Gs: List[np.array]
-                          ):
+    def set_initial_value(self, initial_value: Dict[str, PTM]):
         """Sets the initial value for the MLE optimization
         Args:
-            E: The POVM measurement operator.
-            rho: The inital state.
-            Gs: A list of the gate matrices.
+            initial_value: The dictionary of the initial gateset
         """
+        E = initial_value['E']
+        rho = initial_value['rho']
+        Gs = [initial_value[label] for label in self.Gs]
         self.initial_value = self._join_input_vector(E, rho, Gs)
 
     def optimize(self, initial_value: Optional[np.array] = None) -> Dict:

test/tomography/test_gateset_tomography.py

@@ -41,14 +41,29 @@ def collect_tomography_data(shots=10000,
 
     @staticmethod
     def expected_linear_inversion_gates(Gs, Fs):
-        rho = np.array([[0.5], [0], [0], [0.5]])
+        rho = Gs['rho']
+        E = Gs['E']
         B = np.array([(F @ rho).T[0] for F in Fs]).T
-        return {label: np.linalg.inv(B) @ G @ B for (label, G) in Gs.items()}
+        BB = np.linalg.inv(B)
+        gates = {label: BB @ G @ B for (label, G) in Gs.items()
+                 if label not in ['E', 'rho']}
+        gates['E'] = E @ B
+        gates['rho'] = BB @ rho
+        return gates
 
     @staticmethod
     def hs_distance(A, B):
         return sum([np.abs(x) ** 2 for x in np.nditer(A-B)])
 
+    @staticmethod
+    def convert_from_ptm(vector):
+        Id = np.sqrt(0.5) * np.array([[1, 0], [0, 1]])
+        X = np.sqrt(0.5) * np.array([[0, 1], [1, 0]])
+        Y = np.sqrt(0.5) * np.array([[0, -1j], [1j, 0]])
+        Z = np.sqrt(0.5) * np.array([[1, 0], [0, -1]])
+        v = vector.reshape(4)
+        return v[0] * Id + v[1] * X + v[2] * Y + v[3] * Z
+
     def compare_gates(self, expected_gates, result_gates, labels, delta=0.2):
         for label in labels:
             expected_gate = expected_gates[label]
@@ -69,7 +84,8 @@ def run_test_on_basis_and_noise(self,
 
         labels = gateset_basis.gate_labels
         gates = gateset_basis.gate_matrices
-
+        gates['rho'] = np.array([[np.sqrt(0.5)], [0], [0], [np.sqrt(0.5)]])
+        gates['E'] = np.array([[np.sqrt(0.5), 0, 0, np.sqrt(0.5)]])
         # apply noise if given
         for label in labels:
             if label != "Id" and noise_ptm is not None:
@@ -83,14 +99,16 @@ def run_test_on_basis_and_noise(self,
                                               gateset_basis=gateset_basis)
 
         # linear inversion test
-        expected_gates = self.expected_linear_inversion_gates(gates, Fs)
         result_gates = fitter.linear_inversion()
-        self.compare_gates(expected_gates, result_gates, labels)
+        expected_gates = self.expected_linear_inversion_gates(gates, Fs)
+        self.compare_gates(expected_gates, result_gates, labels + ['E', 'rho'])
 
         # fitter optimization test
-        expected_gates = gates
         result_gates = fitter.fit()
-        self.compare_gates(expected_gates, result_gates, labels)
+        expected_gates = gates
+        expected_gates['E'] = self.convert_from_ptm(expected_gates['E'])
+        expected_gates['rho'] = self.convert_from_ptm(expected_gates['rho'])
+        self.compare_gates(expected_gates, result_gates, labels + ['E', 'rho'])
 
     def test_noiseless_standard_basis(self):
         self.run_test_on_basis_and_noise()