Repository: Qiskit/qiskit-terra
Pull request #5693 (2021-01-24 18:16:39)

Exploit symmetry in QFI and Hessian construction
<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add a release note file using reno if the change needs to be
  documented in the release notes.
⚠️ If your pull request fixes an open issue, please link to the issue.

- [ ] I have added the tests to cover my changes.
- [ ] I have updated the documentation accordingly.
- [ ] I have read the CONTRIBUTING document.
-->

### Summary
resolve https://github.com/Qiskit/qiskit-terra/issues/5631


### Details and comments

1) computation of hessian was optimised
2) full_lin_comb for QFI optimised
3) overlap_block_diag for QFI optimised
4) overlap_diag **unmodified** (no obvious way for optimisation based on symmetry)

Cryoris: Thanks for the contribution @molar-volume! The goal of #5631 was to exploit the symmetry not only for construction of the Hessian/QFI but also for the evaluation (since that is the costly part). It seems your suggestion tackles the construction, is that right? 
Do you want to also improve the evaluation? Meaning that e.g. a 3x3 hessian only evaluates 6 ListOp's instead of 9 🙂 
molar-volume: Hi @Cryoris,
thanks for the review and for suggesting the `combo_fn` to optimise evaluation, I will  use it :)

Patch:
diff --git a/qiskit/opflow/gradients/circuit_qfis/lin_comb_full.py b/qiskit/opflow/gradients/circuit_qfis/lin_comb_full.py
index cc76a5254..6d7d0edc8 100644
--- a/qiskit/opflow/gradients/circuit_qfis/lin_comb_full.py
+++ b/qiskit/opflow/gradients/circuit_qfis/lin_comb_full.py
@@ -17,8 +17,8 @@ import warnings
 
 import numpy as np
 from qiskit.circuit import Gate
-from qiskit.circuit import (QuantumCircuit, QuantumRegister, ParameterVector,
-                            ParameterExpression)
+from qiskit.circuit import QuantumCircuit, QuantumRegister, ParameterVector, ParameterExpression
+from qiskit.utils.arithmetic import triu_to_dense
 
 from ...list_ops.list_op import ListOp
 from ...list_ops.summed_op import SummedOp
@@ -94,7 +94,7 @@ class LinCombFull(CircuitQFI):
         # Get the circuits needed to compute〈∂iψ|∂jψ〉
         for i, param_i in enumerate(params):  # loop over parameters
             qfi_ops = []
-            for j, param_j in enumerate(params):
+            for j, param_j in enumerate(params[i:], i):
                 # Get the gates of the quantum state which are parameterized by param_i
                 qfi_op = []
                 param_gates_i = state_qc._parameter_table[param_i]
@@ -167,7 +167,7 @@ class LinCombFull(CircuitQFI):
 
             qfi_operators.append(ListOp(qfi_ops))
         # Return the full QFI
-        return ListOp(qfi_operators)
+        return ListOp(qfi_operators, combo_fn=triu_to_dense)
 
     @staticmethod
     def trim_circuit(circuit: QuantumCircuit,
diff --git a/qiskit/opflow/gradients/circuit_qfis/overlap_block_diag.py b/qiskit/opflow/gradients/circuit_qfis/overlap_block_diag.py
index 77b57babb..d452ccffe 100644
--- a/qiskit/opflow/gradients/circuit_qfis/overlap_block_diag.py
+++ b/qiskit/opflow/gradients/circuit_qfis/overlap_block_diag.py
@@ -10,13 +10,14 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""The module for Quantum the Fisher Information."""
+"""The module for the Quantum Fisher Information."""
 
 from typing import List, Union
 
 import numpy as np
 from scipy.linalg import block_diag
 from qiskit.circuit import Parameter, ParameterVector, ParameterExpression
+from qiskit.utils.arithmetic import triu_to_dense
 from ...list_ops.list_op import ListOp
 from ...primitive_ops.circuit_op import CircuitOp
 from ...expectations.pauli_expectation import PauliExpectation
@@ -54,7 +55,7 @@ class OverlapBlockDiag(CircuitQFI):
         Raises:
             NotImplementedError: If ``operator`` is neither ``CircuitOp`` nor ``CircuitStateFn``.
         """
-        if not isinstance(operator, (CircuitStateFn)):
+        if not isinstance(operator, (CircuitOp, CircuitStateFn)):
             raise NotImplementedError('operator must be a CircuitOp or CircuitStateFn')
         return self._block_diag_approx(operator=operator, params=params)
 
@@ -146,15 +147,14 @@ class OverlapBlockDiag(CircuitQFI):
             for i, p_i in enumerate(params):
                 generator_i = generators[p_i]
                 param_expr_i = get_parameter_expression(circuit, p_i)
-
-                for j, p_j in enumerate(params):
+                for j, p_j in enumerate(params[i:], i):
                     if i == j:
                         block[i][i] = ListOp([single_terms[i]], combo_fn=lambda x: 1 - x[0] ** 2)
                         if isinstance(param_expr_i, ParameterExpression) and not isinstance(
                                 param_expr_i, Parameter):
                             expr_grad_i = DerivativeBase.parameter_expression_grad(
                                 param_expr_i, p_i)
-                            block[i][j] *= (expr_grad_i) * (expr_grad_i)
+                            block[i][i] *= expr_grad_i * expr_grad_i
                         continue
 
                     generator_j = generators[p_j]
@@ -174,9 +174,11 @@ class OverlapBlockDiag(CircuitQFI):
                         expr_grad_j = DerivativeBase.parameter_expression_grad(param_expr_j, p_j)
                         block[i][j] *= expr_grad_j
 
-            wrapped_block = ListOp([ListOp(row) for row in block])
+            wrapped_block = ListOp([ListOp([
+                block[i][j]
+                for j in range(i, len(params))]) for i in range(len(params))],
+                                   combo_fn=triu_to_dense)
             blocks.append(wrapped_block)
 
-        block_diagonal_qfi = ListOp(oplist=blocks,
-                                    combo_fn=lambda x: np.real(block_diag(*x))[:, perm][perm, :])
-        return block_diagonal_qfi
+        return ListOp(oplist=blocks,
+                      combo_fn=lambda x: np.real(block_diag(*x))[:, perm][perm, :])
diff --git a/qiskit/opflow/gradients/circuit_qfis/overlap_diag.py b/qiskit/opflow/gradients/circuit_qfis/overlap_diag.py
index 1d0bc9f14..54fc19ce2 100644
--- a/qiskit/opflow/gradients/circuit_qfis/overlap_diag.py
+++ b/qiskit/opflow/gradients/circuit_qfis/overlap_diag.py
@@ -83,8 +83,8 @@ class OverlapDiag(CircuitQFI):
 
         """
 
-        if not isinstance(operator, CircuitStateFn):
-            raise NotImplementedError('operator must be a CircuitStateFn')
+        if not isinstance(operator, (CircuitOp, CircuitStateFn)):
+            raise NotImplementedError('operator must be a CircuitOp or CircuitStateFn')
 
         # If a single parameter is given wrap it into a list.
         if isinstance(params, ParameterExpression):
diff --git a/qiskit/opflow/gradients/hessian.py b/qiskit/opflow/gradients/hessian.py
index f7ea7ef08..0a82fca0b 100644
--- a/qiskit/opflow/gradients/hessian.py
+++ b/qiskit/opflow/gradients/hessian.py
@@ -29,6 +29,7 @@ from ..operator_base import OperatorBase
 from .gradient import Gradient
 from .hessian_base import HessianBase
 from ..exceptions import OpflowError
+from ...utils.arithmetic import triu_to_dense
 
 try:
     from jax import grad, jit
@@ -58,16 +59,6 @@ class Hessian(HessianBase):
         Returns:
             OperatorBase: An operator whose evaluation yields the Hessian
         """
-        # if input is a tuple instead of a list, wrap it into a list
-        if isinstance(params, (ParameterVector, list)):
-            # Case: a list of parameters were given, compute the Hessian for all param pairs
-            if all(isinstance(param, ParameterExpression) for param in params):
-                return ListOp(
-                    [ListOp([self.convert(operator, (p0, p1)) for p1 in params]) for p0 in params])
-            # Case: a list was given containing tuples of parameter pairs.
-            # Compute the Hessian entries corresponding to these pairs of parameters.
-            elif all(isinstance(param, tuple) for param in params):
-                return ListOp([self.convert(operator, param_pair) for param_pair in params])
 
         expec_op = PauliExpectation(group_paulis=False).convert(operator).reduce()
         cleaned_op = self._factor_coeffs_out_of_composed_op(expec_op)
@@ -98,6 +89,19 @@ class Hessian(HessianBase):
             Exception: Unintended code is reached
             MissingOptionalLibraryError: jax not installed
         """
+        # if input is a tuple instead of a list, wrap it into a list
+        if isinstance(params, (ParameterVector, list)):
+            # Case: a list of parameters were given, compute the Hessian for all param pairs
+            if all(isinstance(param, ParameterExpression) for param in params):
+                return ListOp([ListOp([
+                    self.get_hessian(operator, (p_i, p_j))
+                    for i, p_i in enumerate(params[j:], j)])
+                               for j, p_j in enumerate(params)],
+                              combo_fn=triu_to_dense)
+            # Case: a list was given containing tuples of parameter pairs.
+            elif all(isinstance(param, tuple) for param in params):
+                # Compute the Hessian entries corresponding to these pairs of parameters.
+                return ListOp([self.get_hessian(operator, param_pair) for param_pair in params])
 
         def is_coeff_c(coeff, c):
             if isinstance(coeff, ParameterExpression):
@@ -105,18 +109,6 @@ class Hessian(HessianBase):
                 return expr == c
             return coeff == c
 
-        if isinstance(params, (ParameterVector, list)):
-            # Case: a list of parameters were given, compute the Hessian for all param pairs
-            if all(isinstance(param, ParameterExpression) for param in params):
-                return ListOp(
-                    [ListOp([self.get_hessian(operator, (p0, p1)) for p1 in params])
-                     for p0 in params])
-            # Case: a list was given containing tuples of parameter pairs.
-            # Compute the Hessian entries corresponding to these pairs of parameters.
-            elif all(isinstance(param, tuple) for param in params):
-                return ListOp(
-                    [self.get_hessian(operator, param_pair) for param_pair in params])
-
         # If a gradient is requested w.r.t a single parameter, then call the
         # Gradient().get_gradient method.
         if isinstance(params, ParameterExpression):
diff --git a/qiskit/utils/arithmetic.py b/qiskit/utils/arithmetic.py
index a4eb20d26..cbb6c886e 100644
--- a/qiskit/utils/arithmetic.py
+++ b/qiskit/utils/arithmetic.py
@@ -131,3 +131,23 @@ def transpositions(permutation: List[int]) -> List[Tuple[int, int]]:
             for y in x[len_x - 1:0:-1]:
                 res.append((first, y))
     return res
+
+
+def triu_to_dense(triu: np.ndarray) -> np.ndarray:
+    """Converts upper triangular part of matrix to dense matrix.
+
+    Args:
+        triu: array in the form [[A, B, C], [D, E], [F]]
+
+    Returns:
+        Array [[A, B, C], [B, D, E], [C, E, F]]
+    """
+    dim = len(triu)
+    matrix = np.empty((dim, dim), dtype=complex)
+    for i in range(dim):
+        for j in range(dim - i):
+            matrix[i, i + j] = triu[i][j]
+            if j != 0:
+                matrix[i + j, i] = triu[i][j]
+
+    return matrix


Test patch:
diff --git a/test/python/test_util.py b/test/python/test_util.py
index 624378b6c..7787fa846 100644
--- a/test/python/test_util.py
+++ b/test/python/test_util.py
@@ -13,9 +13,11 @@
 """Tests for qiskit/_util.py"""
 
 from unittest import mock
+import numpy as np
 
 from qiskit import util
 from qiskit.test import QiskitTestCase
+from qiskit.utils.arithmetic import triu_to_dense
 
 
 class TestUtil(QiskitTestCase):
@@ -30,3 +32,14 @@ class TestUtil(QiskitTestCase):
         del cpu_count_mock, vmem_mock, platform_mock  # unused
         result = util.local_hardware_info()
         self.assertEqual(1, result['cpus'])
+
+    def test_triu_to_dense(self):
+        """ Test conversion of upper triangular matrix to dense matrix. """
+        np.random.seed(50)
+        n = np.random.randint(5, 15)
+        m = np.random.randint(-100, 100, size=(n, n))
+        symm = (m + m.T) / 2
+
+        triu = np.array([[symm[i, j] for i in range(j, n)] for j in range(n)])
+
+        self.assertTrue(np.array_equal(symm, triu_to_dense(triu)))