Add log_i() functionality to MatrixOp.

qiskit/aqua/operators/evolutions/evolved_op.py

@@ -14,7 +14,7 @@
 
 """ EvolutionOp Class """
 
-from typing import Optional, Union, Set
+from typing import Optional, Union, Set, List
 import logging
 import numpy as np
 import scipy
@@ -23,6 +23,8 @@
 
 from ..operator_base import OperatorBase
 from ..primitive_ops.primitive_op import PrimitiveOp
+from ..primitive_ops.matrix_op import MatrixOp
+from ..list_ops import ListOp
 from ..list_ops.summed_op import SummedOp
 from ..list_ops.composed_op import ComposedOp
 from ..list_ops.tensored_op import TensoredOp
@@ -111,8 +113,6 @@ def bind_parameters(self, param_dict: dict) -> OperatorBase:
         if isinstance(self.coeff, ParameterExpression):
             unrolled_dict = self._unroll_param_dict(param_dict)
             if isinstance(unrolled_dict, list):
-                # pylint: disable=import-outside-toplevel
-                from ..list_ops.list_op import ListOp
                 return ListOp([self.bind_parameters(param_dict) for param_dict in unrolled_dict])
             if self.coeff.parameters <= set(unrolled_dict.keys()):
                 binds = {param: unrolled_dict[param] for param in self.coeff.parameters}
@@ -124,11 +124,27 @@ def eval(self,
                           OperatorBase] = None) -> Union[OperatorBase, float, complex]:
         return self.to_matrix_op().eval(front=front)
 
-    def to_matrix(self, massive: bool = False) -> np.ndarray:
+    def to_matrix(self, massive: bool = False) -> Union[np.ndarray, List[np.ndarray]]:
+        if self.primitive.__class__.__name__ == ListOp.__name__:
+            return [op.exp_i().to_matrix() * self.primitive.coeff * self.coeff
+                    for op in self.primitive.oplist]
+
         prim_mat = -1.j * self.primitive.to_matrix()
         # pylint: disable=no-member
         return scipy.linalg.expm(prim_mat) * self.coeff
 
+    def to_matrix_op(self, massive: bool = False) -> OperatorBase:
+        """ Returns a ``MatrixOp`` equivalent to this Operator. """
+        if self.primitive.__class__.__name__ == ListOp.__name__:
+            return ListOp([op.exp_i().to_matrix_op() for op in self.primitive.oplist],
+                          coeff=self.primitive.coeff * self.coeff)
+
+        prim_mat = EvolvedOp(self.primitive).to_matrix(massive=massive)
+        return MatrixOp(prim_mat, coeff=self.coeff)
+
+    def log_i(self, massive: bool = False) -> OperatorBase:
+        return self.primitive * self.coeff
+
     # pylint: disable=arguments-differ
     def to_instruction(self, massive: bool = False) -> Instruction:
         return self.primitive.to_matrix_op(massive=massive).exp_i()

qiskit/aqua/operators/list_ops/list_op.py

@@ -213,7 +213,7 @@ def to_matrix(self, massive: bool = False) -> np.ndarray:
                 ' Set massive=True if you want to proceed.'.format(2**self.num_qubits))
 
         # Combination function must be able to handle classical values
-        return self.combo_fn([op.to_matrix() for op in self.oplist]) * self.coeff
+        return self.combo_fn([op.to_matrix() * self.coeff for op in self.oplist])
 
     def to_spmatrix(self) -> Union[spmatrix, List[spmatrix]]:
         """ Returns SciPy sparse matrix representation of the Operator.
@@ -278,6 +278,19 @@ def exp_i(self) -> OperatorBase:
         from qiskit.aqua.operators import EvolvedOp
         return EvolvedOp(self)
 
+    def log_i(self, massive: bool = False) -> OperatorBase:
+        """Return a ``MatrixOp`` equivalent to log(H)/-i for this operator H. This
+        function is the effective inverse of exp_i, equivalent to finding the Hermitian
+        Operator which produces self when exponentiated. For proper ListOps, applies ``log_i``
+        to all ops in oplist.
+        """
+        if self.__class__.__name__ == ListOp.__name__:
+            return ListOp([op.log_i(massive=massive) for op in self.oplist],
+                          coeff=self.coeff,
+                          abelian=False)
+
+        return self.to_matrix_op(massive=massive).log_i(massive=massive)
+
     def __str__(self) -> str:
         main_string = "{}(\n[{}])".format(self.__class__.__name__, ',\n'.join(
             [str(op) for op in self.oplist]))

qiskit/aqua/operators/primitive_ops/matrix_op.py

@@ -114,12 +114,6 @@ def compose(self, other: OperatorBase) -> OperatorBase:
         return ComposedOp([self, other])
 
     def to_matrix(self, massive: bool = False) -> np.ndarray:
-        if self.num_qubits > 16 and not massive:
-            raise ValueError(
-                'to_matrix will return an exponentially large matrix, '
-                'in this case {0}x{0} elements.'
-                ' Set massive=True if you want to proceed.'.format(2**self.num_qubits))
-
         return self.primitive.data * self.coeff
 
     def __str__(self) -> str:

qiskit/aqua/operators/primitive_ops/pauli_op.py

@@ -209,15 +209,18 @@ def exp_i(self) -> OperatorBase:
             return PauliOp(self.primitive)
         if np.sum(sig_qubits) == 1:
             sig_qubit_index = sig_qubits.tolist().index(True)
+            coeff = np.real(self.coeff) \
+                if not isinstance(self.coeff, ParameterExpression) \
+                else self.coeff
             # Y rotation
             if corrected_x[sig_qubit_index] and corrected_z[sig_qubit_index]:
-                rot_op = PrimitiveOp(RYGate(self.coeff))
+                rot_op = PrimitiveOp(RYGate(coeff))
             # Z rotation
             elif corrected_z[sig_qubit_index]:
-                rot_op = PrimitiveOp(RZGate(self.coeff))
+                rot_op = PrimitiveOp(RZGate(coeff))
             # X rotation
             elif corrected_x[sig_qubit_index]:
-                rot_op = PrimitiveOp(RXGate(self.coeff))
+                rot_op = PrimitiveOp(RXGate(coeff))
 
             from ..operator_globals import I
             left_pad = I.tensorpower(sig_qubit_index)

qiskit/aqua/operators/primitive_ops/primitive_op.py

@@ -18,6 +18,7 @@
 import logging
 import numpy as np
 from scipy.sparse import spmatrix
+import scipy.linalg
 
 from qiskit import QuantumCircuit
 from qiskit.circuit import Instruction, ParameterExpression
@@ -189,6 +190,16 @@ def exp_i(self) -> OperatorBase:
         from qiskit.aqua.operators import EvolvedOp
         return EvolvedOp(self)
 
+    def log_i(self, massive: bool = False) -> OperatorBase:
+        """Return a ``MatrixOp`` equivalent to log(H)/-i for this operator H. This
+        function is the effective inverse of exp_i, equivalent to finding the Hermitian
+        Operator which produces self when exponentiated."""
+        # pylint: disable=cyclic-import
+        from ..operator_globals import EVAL_SIG_DIGITS
+        from .matrix_op import MatrixOp
+        return MatrixOp(np.around(scipy.linalg.logm(self.to_matrix(massive=massive)) / -1j,
+                                  decimals=EVAL_SIG_DIGITS))
+
     def __str__(self) -> str:
         raise NotImplementedError
 
@@ -252,6 +263,11 @@ def to_pauli_op(self, massive: bool = False) -> OperatorBase:
         """ Returns a sum of ``PauliOp`` s equivalent to this Operator. """
         mat_op = self.to_matrix_op(massive=massive)
         sparse_pauli = SparsePauliOp.from_operator(mat_op.primitive)
+        if not sparse_pauli.to_list():
+            # pylint: disable=import-outside-toplevel
+            from ..operator_globals import I
+            return (I ^ self.num_qubits) * 0.0
+
         return sum([PrimitiveOp(Pauli.from_label(label),
                                 coeff.real if coeff == coeff.real else coeff)
                     for (label, coeff) in sparse_pauli.to_list()]) * self.coeff

test/aqua/operators/test_evolution.py

@@ -156,6 +156,37 @@ def test_matrix_op_evolution(self):
         ref_mat = scipy.linalg.expm(-1j * op.to_matrix())
         np.testing.assert_array_almost_equal(ref_mat, exp_mat)
 
+    def test_log_i(self):
+        """ MatrixOp.log_i() test """
+        op = (-1.052373245772859 * I ^ I) + \
+             (0.39793742484318045 * I ^ Z) + \
+             (0.18093119978423156 * X ^ X) + \
+             (-0.39793742484318045 * Z ^ I) + \
+             (-0.01128010425623538 * Z ^ Z) * np.pi/2
+        # Test with CircuitOp
+        log_exp_op = op.to_matrix_op().exp_i().log_i().to_pauli_op()
+        np.testing.assert_array_almost_equal(op.to_matrix(), log_exp_op.to_matrix())
+
+        # Test with MatrixOp
+        log_exp_op = op.to_matrix_op().exp_i().to_matrix_op().log_i().to_pauli_op()
+        np.testing.assert_array_almost_equal(op.to_matrix(), log_exp_op.to_matrix())
+
+        # Test with PauliOp
+        log_exp_op = op.to_matrix_op().exp_i().to_pauli_op().log_i().to_pauli_op()
+        np.testing.assert_array_almost_equal(op.to_matrix(), log_exp_op.to_matrix())
+
+        # Test with EvolvedOp
+        log_exp_op = op.exp_i().to_pauli_op().log_i().to_pauli_op()
+        np.testing.assert_array_almost_equal(op.to_matrix(), log_exp_op.to_matrix())
+
+        # Test with proper ListOp
+        op = ListOp([(0.39793742484318045 * I ^ Z),
+                     (0.18093119978423156 * X ^ X),
+                     (-0.39793742484318045 * Z ^ I),
+                     (-0.01128010425623538 * Z ^ Z) * np.pi / 2])
+        log_exp_op = op.to_matrix_op().exp_i().to_matrix_op().log_i().to_pauli_op()
+        np.testing.assert_array_almost_equal(op.to_matrix(), log_exp_op.to_matrix())
+
     def test_matrix_op_parameterized_evolution(self):
         """ parameterized MatrixOp evolution test """
         # pylint: disable=no-member