Fix probability Hessians with the parameter shift method.

qiskit/aqua/operators/gradients/circuit_gradients/param_shift.py

@@ -197,19 +197,24 @@ def _parameter_shift(self,
 
                 p_param = pshift_gate.params[param_index]
                 m_param = mshift_gate.params[param_index]
-
+                # For analytic gradients the circuit parameters are shifted once by +pi/2 and
+                # once by -pi/2.
                 if self.analytic:
                     shift_constant = 0.5
                     pshift_gate.params[param_index] = (p_param + (np.pi / (4 * shift_constant)))
                     mshift_gate.params[param_index] = (m_param - (np.pi / (4 * shift_constant)))
+                # For finite difference gradients the circuit parameters are shifted once by
+                # +epsilon and once by -epsilon.
                 else:
                     shift_constant = 1. / (2 * self._epsilon)
                     pshift_gate.params[param_index] = (p_param + self._epsilon)
                     mshift_gate.params[param_index] = (m_param - self._epsilon)
-
+                # The results of the shifted operators are now evaluated according the parameter
+                # shift / finite difference formula.
                 if isinstance(operator, ComposedOp):
                     shifted_op = shift_constant * (pshift_op - mshift_op)
-
+                # If the operator represents a quantum state then we apply a special combo
+                # function to evaluate probability gradients.
                 elif isinstance(operator, StateFn):
                     shifted_op = ListOp(
                         [pshift_op, mshift_op],
@@ -261,9 +266,16 @@ def get_primitives(item):
                 item = item.primitive.data
             return item
 
-        if isinstance(x, Iterable):
+        is_statefn = False
+        if isinstance(x, list):
+            # Check if all items in x are a StateFn items
+            if all(isinstance(item, StateFn) for item in x):
+                is_statefn = True
             items = [get_primitives(item) for item in x]
         else:
+            # Check if x is a StateFn item
+            if isinstance(x, StateFn):
+                is_statefn = True
             items = [get_primitives(x)]
         if isinstance(items[0], dict):
             prob_dict: Dict[str, float] = {}
@@ -273,8 +285,17 @@ def get_primitives(item):
                                      shift_constant * ((-1) ** i) * prob_counts
             return prob_dict
         elif isinstance(items[0], Iterable):
-            return shift_constant * np.subtract(np.multiply(items[0], np.conj(items[0])),
-                                                np.multiply(items[1], np.conj(items[1])))
+            # If x was given as StateFn the state amplitudes need to be multiplied in order to
+            # evaluate the sampling probabilities which are then subtracted according to the
+            # parameter shift rule.
+            if is_statefn:
+                return shift_constant * np.subtract(np.multiply(items[0], np.conj(items[0])),
+                                                    np.multiply(items[1], np.conj(items[1])))
+            # If x was not given as a StateFn the state amplitudes were already converted into
+            # sampling probabilities which are then only subtracted according to the
+            # parameter shift rule.
+            else:
+                return shift_constant * np.subtract(items[0], items[1])
         raise TypeError(
             'Probability gradients can only be evaluated from VectorStateFs or DictStateFns.')
 

qiskit/aqua/operators/gradients/circuit_qfis/lin_comb_full.py

@@ -63,34 +63,42 @@ def convert(self,
         phase_fix_observable = ~StateFn((X + 1j * Y) ^ (I ^ operator.num_qubits))
         # see https://arxiv.org/pdf/quant-ph/0108146.pdf
 
+        # Check if the given operator corresponds to a quantum state given as a circuit.
         if not isinstance(operator, CircuitStateFn):
             raise TypeError(
                 'LinCombFull is only compatible with states that are given as CircuitStateFn')
 
+        # If a single parameter is given wrap it into a list.
         if not isinstance(params, (list, np.ndarray)):
             params = [params]
         state_qc = operator.primitive
 
         # First, the operators are computed which can compensate for a potential phase-mismatch
         # between target and trained state, i.e.〈ψ|∂lψ〉
         phase_fix_states = None
+        # Add working qubit
         qr_work = QuantumRegister(1, 'work_qubit')
         work_q = qr_work[0]
         additional_qubits: Tuple[List[Qubit], List[Qubit]] = ([work_q], [])
-        # create a copy of the original state with an additional work_q register
         for param in params:
+            # Get the gates of the given quantum state which are parameterized by param
             param_gates = state_qc._parameter_table[param]
+            # Loop through the occurrences of param in the quantum state
             for m, param_occurence in enumerate(param_gates):
+                # Get the coefficients and gates for the linear combination gradient for each
+                # occurrence, see e.g. https://arxiv.org/abs/2006.06004
                 coeffs_i, gates_i = LinComb._gate_gradient_dict(param_occurence[0])[
                     param_occurence[1]]
+                # Construct the quantum states which are then evaluated for the respective QFI
+                # element.
                 for k, gate_to_insert_i in enumerate(gates_i):
                     grad_state = state_qc.copy()
                     grad_state.add_register(qr_work)
 
                     # apply Hadamard on work_q
                     LinComb.insert_gate(grad_state, param_occurence[0], HGate(),
                                         qubits=[work_q])
-                    # Fix work_q phase
+                    # Fix work_q phase such that the gradient is correct.
                     coeff_i = coeffs_i[k]
                     sign = np.sign(coeff_i)
                     is_complex = np.iscomplex(coeff_i)
@@ -141,17 +149,21 @@ def convert(self,
                                             gate_to_insert_i,
                                             additional_qubits=additional_qubits)
 
+                    # Remove unnecessary gates.
                     grad_state = self.trim_circuit(grad_state, param_occurence[0])
-
+                    # Apply the final Hadamard on the working qubit.
                     grad_state.h(work_q)
-
+                    # Add the coefficient needed for the gradient as well as the original
+                    # coefficient of the given quantum state.
                     state = np.sqrt(np.abs(coeff_i)) * operator.coeff * CircuitStateFn(grad_state)
 
-                    # Chain Rule parameter expressions
+                    # Check if the gate parameter corresponding to param is a parameter expression
                     gate_param = param_occurence[0].params[param_occurence[1]]
                     if gate_param == param:
                         state = phase_fix_observable @ state
                     else:
+                        # If the gate parameter is a parameter expressions the chain rule needs
+                        # to be taken into account.
                         if isinstance(gate_param, ParameterExpression):
                             expr_grad = DerivativeBase.parameter_expression_grad(gate_param, param)
                             state = (expr_grad * phase_fix_observable) @ state
@@ -161,38 +173,41 @@ def convert(self,
                     if m == 0 and k == 0:
                         phase_fix_state = state
                     else:
+                        # Take the product rule into account
                         phase_fix_state += state
+            # Create a list for the phase fix states
             if not phase_fix_states:
                 phase_fix_states = [phase_fix_state]
             else:
                 phase_fix_states += [phase_fix_state]
 
         # Get  4 * Re[〈∂kψ|∂lψ]
         qfi_operators = []
+        # Add a working qubit
         qr_work_qubit = QuantumRegister(1, 'work_qubit')
         work_qubit = qr_work_qubit[0]
         additional_qubits = ([work_qubit], [])
         # create a copy of the original circuit with an additional work_qubit register
         circuit = state_qc.copy()
         circuit.add_register(qr_work_qubit)
+        # Apply a Hadamard on the working qubit
         LinComb.insert_gate(circuit, state_qc._parameter_table[params[0]][0][0], HGate(),
                             qubits=[work_qubit])
 
         # Get the circuits needed to compute〈∂iψ|∂jψ〉
         for i, param_i in enumerate(params):  # loop over parameters
             qfi_ops = None
             for j, param_j in enumerate(params):
-                # Construct the circuits
+                # Get the gates of the quantum state which are parameterized by param_i
                 param_gates_i = state_qc._parameter_table[param_i]
                 for m_i, param_occurence_i in enumerate(param_gates_i):
                     coeffs_i, gates_i = LinComb._gate_gradient_dict(param_occurence_i[0])[
                         param_occurence_i[1]]
-                    # apply Hadamard on work_qubit
 
                     for k_i, gate_to_insert_i in enumerate(gates_i):
                         coeff_i = coeffs_i[k_i]
+                        # Get the gates of the quantum state which are parameterized by param_j
                         param_gates_j = state_qc._parameter_table[param_j]
-
                         for m_j, param_occurence_j in enumerate(param_gates_j):
                             coeffs_j, gates_j = \
                                 LinComb._gate_gradient_dict(param_occurence_j[0])[
@@ -204,7 +219,8 @@ def convert(self,
                                 qfi_circuit = QuantumCircuit(*circuit.qregs)
                                 qfi_circuit.data = circuit.data
 
-                                # Fix work_qubit phase
+                                # Correct the phase of the working qubit according to coefficient i
+                                # and coefficient j
                                 sign = np.sign(np.conj(coeff_i) * coeff_j)
                                 is_complex = np.iscomplex(np.conj(coeff_i) * coeff_j)
                                 if sign == -1:
@@ -230,7 +246,7 @@ def convert(self,
                                                     XGate(),
                                                     qubits=[work_qubit])
 
-                                # Insert controlled, intercepting gate - controlled by |1>
+                                # Insert controlled, intercepting gate i - controlled by |1>
                                 if isinstance(param_occurence_i[0], UGate):
                                     if param_occurence_i[1] == 0:
                                         LinComb.insert_gate(
@@ -267,7 +283,7 @@ def convert(self,
                                     qfi_circuit, gate_to_insert_i,
                                     XGate(), qubits=[work_qubit], after=True)
 
-                                # Insert controlled, intercepting gate - controlled by |0>
+                                # Insert controlled, intercepting gate j - controlled by |0>
                                 if isinstance(param_occurence_j[0], UGate):
                                     if param_occurence_j[1] == 0:
                                         LinComb.insert_gate(
@@ -310,21 +326,26 @@ def convert(self,
                                     qfi_circuit = self.trim_circuit(
                                         qfi_circuit, param_occurence_j[0]
                                     )
-
+                                # Apply final Hadamard gate
                                 qfi_circuit.h(work_qubit)
-                                # Convert the quantum circuit into a CircuitStateFn
+                                # Convert the quantum circuit into a CircuitStateFn and add the
+                                # coefficients i, j and the original operator coefficient
                                 term = np.sqrt(np.abs(coeff_i) * np.abs(coeff_j)) * operator.coeff
                                 term = term * CircuitStateFn(qfi_circuit)
 
-                                # Chain Rule Parameter Expression
+                                # Check if the gate parameters i and j are parameter expressions
                                 gate_param_i = param_occurence_i[0].params[param_occurence_i[1]]
                                 gate_param_j = param_occurence_j[0].params[param_occurence_j[1]]
 
                                 meas = deepcopy(qfi_observable)
+                                # If the gate parameter i is a parameter expression use the chain
+                                # rule.
                                 if isinstance(gate_param_i, ParameterExpression):
                                     expr_grad = DerivativeBase.parameter_expression_grad(
                                         gate_param_i, param_i)
                                     meas *= expr_grad
+                                # If the gate parameter j is a parameter expression use the chain
+                                # rule.
                                 if isinstance(gate_param_j, ParameterExpression):
                                     expr_grad = DerivativeBase.parameter_expression_grad(
                                         gate_param_j, param_j)

qiskit/aqua/operators/gradients/hessian.py

@@ -241,12 +241,7 @@ def is_coeff_c(coeff, c):
 
         elif isinstance(operator, StateFn):
             if not operator.is_measurement:
-                from .circuit_gradients import LinComb
-                if isinstance(self.hess_method, LinComb):
-                    return self.hess_method.convert(operator, params)
-                else:
-                    raise TypeError('Please set the attribute hess_method to ´lin_comb´ to compute '
-                                    'probability Hessians.')
+                return self.hess_method.convert(operator, params)
 
             else:
                 raise TypeError('The computation of Hessians is only supported for Operators which '

qiskit/aqua/operators/gradients/natural_gradient.py

@@ -91,22 +91,29 @@ def convert(self,
                             'CircuitStateFn.')
         if not isinstance(params, Iterable):
             params = [params]
+        # Instantiate the gradient
         grad = Gradient(self._grad_method, epsilon=self._epsilon).convert(operator, params)
+        # Instantiate the QFI metric which is used to re-scale the gradient
         metric = self._qfi_method.convert(operator[-1], params) * 0.25
 
+        # Define the function which compute the natural gradient from the gradient and the QFI.
         def combo_fn(x):
             c = np.real(x[0])
             a = np.real(x[1])
             if self.regularization:
+                # If a regularization method is chosen then use a regularized solver to
+                # construct the natural gradient.
                 nat_grad = NaturalGradient._regularized_sle_solver(
                     a, c, regularization=self.regularization)
             else:
                 try:
+                    # Try to solve the system of linear equations Ax = C.
                     nat_grad = np.linalg.solve(a, c)
                 except np.linalg.LinAlgError:  # singular matrix
                     nat_grad = np.linalg.lstsq(a, c)[0]
             return np.real(nat_grad)
-
+        # Define the ListOp which combines the gradient and the QFI according to the combination
+        # function defined above.
         return ListOp([grad, metric], combo_fn=combo_fn)
 
     @property

test/aqua/operators/test_gradients.py

@@ -467,7 +467,8 @@ def test_prob_grad(self, method):
                 np.testing.assert_array_almost_equal(prob_grad_result,
                                                      correct_values[i][j], decimal=1)
 
-    def test_prob_hess_lin_comb(self):
+    @data('lin_comb', 'param_shift', 'fin_diff')
+    def test_prob_hess(self, method):
         """Test the probability Hessian using linear combination of unitaries method
 
         d^2p0/da^2 = - sin(a)sin(b) / 2
@@ -488,7 +489,7 @@ def test_prob_hess_lin_comb(self):
 
         op = CircuitStateFn(primitive=qc, coeff=1.)
 
-        prob_hess = Hessian(hess_method='lin_comb').convert(operator=op, params=params)
+        prob_hess = Hessian(hess_method=method).convert(operator=op, params=params)
         values_dict = [{a: np.pi / 4, b: 0}, {a: np.pi / 4, b: np.pi / 4},
                        {a: np.pi / 2, b: np.pi}]
         correct_values = [[[0, 0], [1 / (2 * np.sqrt(2)), - 1 / (2 * np.sqrt(2))]],