Repository: Qiskit/qiskit-terra
Pull request #5617 (2021-01-12 22:31:49)

quantum_info.operator mixin interface refactor
<!--
⚠️ 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

Internal API changes to BaseOperator class to simplify making operator subclasses, and reduce copy paste of API docs.

### Details and comments

Refactors the quantum_info BaseOperator and subclasses to use mixins for defining the operator inferface. This allows operator subclasses to only include the mixin for interfaces they support. For example Operator and QuantumChannel operators support linear operations (addition, subtraction), but Clifford and Pauli operators do not.

Current mixins are:

- `GroupMixin`: this defines the basic API to be compatible with QuantumCircuits. It includes `compose`, `dot`, `tensor`, `expand`
- `AdjointMixin`: defines `conjugate`, `transpose`, `adjoint` methods.
- `MultiplyMixin`: this is not part of the BaseOperator and can be used by special operators which support scalar multiplication. It defines `__mul__` and `__truediv__` for scalars.
- `LinearMixin`: this defines linear operators which support addition, subtraction, and scalar multiplication. It includes the `MultiplyMixin`, plus defines `__add__` and `__sub__` methods.
- `ToleranceMixin`: this defines the `atol` and `rtol` properties for the operator for use in numeric comparisons

`BaseOperator` only uses `GroupMixin` for its API. Other mixins can be included by subclasses which support their interface or require their functionality. A `LinearOp` abstract class was also added that includes `LinearMixin`, `AdjointMixin` and `TolerancesMixin` with BaseOperator.

This PR also refactors a lot of doc strings in the operator classes to avoid unnecessary duplication in derived classes. This is handled by the `generate_apidocs` function that must be called on a class *after* its definition (Due to however python handles abstractmethod doc strings I couldn't seem to get this function to work as a decorator when defining the class).

molar-volume: Maybe
```python
def _evolve(self, state, qargs=None):
        return SuperOp(self)._evolve(state, qargs)
```
could be moved to `QuantumChannel`, since it is same for all subclasses (except `SuperOp`).
chriseclectic: > Maybe
> 
> ```python
> def _evolve(self, state, qargs=None):
>         return SuperOp(self)._evolve(state, qargs)
> ```
> 
> could be moved to `QuantumChannel`, since it is same for all subclasses (except `SuperOp`).

That is True, however it has cyclic import issues.
ikkoham: Why do private and public methods coexist? (e.g.` _compose`, `compose`)
`_compose` is an abstract method, so it has to be defined by the user. But at the same time, if user also overrides `compose`, then `__matmul__` will become the `compose` (not `_compose`). Is this the desired behavior?
ikkoham: Why do private and public methods coexist? (e.g.` _compose`, `compose`)
`_compose` is an abstract method, so it has to be defined by the user. But at the same time, if user also overrides `compose`, then `__matmul__` will become the `compose` (not `_compose`). Is this the desired behavior?
chriseclectic: > Why do private and public methods coexist? (e.g.` _compose`, `compose`)
> `_compose` is an abstract method, so it has to be defined by the user. But at the same time, if user also overrides `compose`, then `__matmul__` will become the `compose` (not `_compose`). Is this the desired behavior?

The idea was user never overrides `compose` only `_compose`. Honestly this was mostly to prevent having to have the bit of code:
```
        if qargs is None:
            qargs = getattr(other, 'qargs', None)
```
in every subclasses compose method. But maybe thats not a good enough reason.
chriseclectic: I updated to remove the `_tensor` and `_compose` methods and go back to the original way of have `compose`, `expand`, `tensor` all abstract methods. Thanks to review feedback I think the added complexity wasn't worth the small code reduction.

Patch:
diff --git a/qiskit/quantum_info/operators/base_operator.py b/qiskit/quantum_info/operators/base_operator.py
index 54e23ef10..365ad8f55 100644
--- a/qiskit/quantum_info/operators/base_operator.py
+++ b/qiskit/quantum_info/operators/base_operator.py
@@ -15,16 +15,17 @@ Abstract BaseOperator class.
 """
 
 import copy
-from abc import abstractmethod
-
-import numpy as np
+from abc import ABC
 
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.op_shape import OpShape
 
+from .mixins import GroupMixin
+
 
-class BaseOperator:
-    """Abstract linear operator base class."""
+# pylint: disable = abstract-method
+class BaseOperator(GroupMixin, ABC):
+    """Abstract operator base class."""
 
     def __init__(self, input_dims=None, output_dims=None,
                  num_qubits=None, shape=None, op_shape=None):
@@ -132,179 +133,3 @@ class BaseOperator:
     def copy(self):
         """Make a deep copy of current operator."""
         return copy.deepcopy(self)
-
-    def adjoint(self):
-        """Return the adjoint of the operator."""
-        return self.conjugate().transpose()
-
-    @abstractmethod
-    def conjugate(self):
-        """Return the conjugate of the operator."""
-        pass
-
-    @abstractmethod
-    def transpose(self):
-        """Return the transpose of the operator."""
-        pass
-
-    @abstractmethod
-    def tensor(self, other):
-        """Return the tensor product operator self ⊗ other.
-
-        Args:
-            other (BaseOperator): a operator subclass object.
-
-        Returns:
-            BaseOperator: the tensor product operator self ⊗ other.
-
-        Raises:
-            QiskitError: if other is not an operator.
-        """
-        pass
-
-    @abstractmethod
-    def expand(self, other):
-        """Return the tensor product operator other ⊗ self.
-
-        Args:
-            other (BaseOperator): an operator object.
-
-        Returns:
-            BaseOperator: the tensor product operator other ⊗ self.
-
-        Raises:
-            QiskitError: if other is not an operator.
-        """
-        pass
-
-    @abstractmethod
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed operator.
-
-        Args:
-            other (BaseOperator): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            BaseOperator: The operator self @ other.
-
-        Raises:
-            QiskitError: if other cannot be converted to an operator, or has
-                         incompatible dimensions for specified subsystems.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            matrix operators. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
-        pass
-
-    def dot(self, other, qargs=None):
-        """Return the right multiplied operator self * other.
-
-        Args:
-            other (BaseOperator): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-
-        Returns:
-            BaseOperator: The operator self * other.
-
-        Raises:
-            QiskitError: if other cannot be converted to an operator, or has
-                         incompatible dimensions for specified subsystems.
-        """
-        return self.compose(other, qargs=qargs, front=True)
-
-    def power(self, n):
-        """Return the compose of a operator with itself n times.
-
-        Args:
-            n (int): the number of times to compose with self (n>0).
-
-        Returns:
-            BaseOperator: the n-times composed operator.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the operator
-                         are not equal, or the power is not a positive integer.
-        """
-        # NOTE: if a subclass can have negative or non-integer powers
-        # this method should be overridden in that class.
-        if not isinstance(n, (int, np.integer)) or n < 1:
-            raise QiskitError("Can only power with positive integer powers.")
-        if not self._op_shape.is_square:
-            raise QiskitError("Can only power with input_dims = output_dims.")
-        ret = self.copy()
-        for _ in range(1, n):
-            ret = ret.compose(self)
-        return ret
-
-    def _add(self, other, qargs=None):
-        """Return the linear operator self + other.
-
-        If ``qargs`` are specified the other operator will be added
-        assuming it is identity on all other subsystems.
-
-        Args:
-            other (BaseOperator): an operator object.
-            qargs (None or list): optional subsystems to add on
-                                  (Default: None)
-
-        Returns:
-            BaseOperator: the operator self + other.
-
-        Raises:
-            NotImplementedError: if subclass does not support addition.
-        """
-        raise NotImplementedError(
-            "{} does not support addition".format(type(self)))
-
-    def _multiply(self, other):
-        """Return the linear operator other * self.
-
-        Args:
-            other (complex): a complex number.
-
-        Returns:
-            BaseOperator: the linear operator other * self.
-
-        Raises:
-            NotImplementedError: if subclass does not support multiplication.
-        """
-        raise NotImplementedError(
-            "{} does not support scalar multiplication".format(type(self)))
-
-    # Overloads
-    def __matmul__(self, other):
-        return self.compose(other)
-
-    def __mul__(self, other):
-        return self.dot(other)
-
-    def __rmul__(self, other):
-        return self._multiply(other)
-
-    def __pow__(self, n):
-        return self.power(n)
-
-    def __xor__(self, other):
-        return self.tensor(other)
-
-    def __truediv__(self, other):
-        return self._multiply(1 / other)
-
-    def __add__(self, other):
-        return self._add(other)
-
-    def __sub__(self, other):
-        return self._add(-other)
-
-    def __neg__(self):
-        return self._multiply(-1)
diff --git a/qiskit/quantum_info/operators/channel/chi.py b/qiskit/quantum_info/operators/channel/chi.py
index 954e04c77..8e9fc8b75 100644
--- a/qiskit/quantum_info/operators/channel/chi.py
+++ b/qiskit/quantum_info/operators/channel/chi.py
@@ -16,6 +16,7 @@
 Chi-matrix representation of a Quantum Channel.
 """
 
+import copy
 import numpy as np
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -25,6 +26,7 @@ from qiskit.quantum_info.operators.channel.quantum_channel import QuantumChannel
 from qiskit.quantum_info.operators.channel.choi import Choi
 from qiskit.quantum_info.operators.channel.superop import SuperOp
 from qiskit.quantum_info.operators.channel.transformations import _to_chi
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class Chi(QuantumChannel):
@@ -130,43 +132,26 @@ class Chi(QuantumChannel):
         return (self._input_dim, self._output_dim, self._input_dim,
                 self._output_dim)
 
+    def _evolve(self, state, qargs=None):
+        return SuperOp(self)._evolve(state, qargs)
+
+    # ---------------------------------------------------------------------
+    # BaseOperator methods
+    # ---------------------------------------------------------------------
+
     def conjugate(self):
-        """Return the conjugate of the QuantumChannel."""
         # Since conjugation is basis dependent we transform
         # to the Choi representation to compute the
         # conjugate channel
         return Chi(Choi(self).conjugate())
 
     def transpose(self):
-        """Return the transpose of the QuantumChannel."""
-        # Since conjugation is basis dependent we transform
-        # to the Choi representation to compute the
-        # conjugate channel
         return Chi(Choi(self).transpose())
 
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            Chi: The quantum channel self @ other.
+    def adjoint(self):
+        return Chi(Choi(self).adjoint())
 
-        Raises:
-            QiskitError: if other has incompatible dimensions.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
+    def compose(self, other, qargs=None, front=False):
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
         if qargs is not None:
@@ -176,74 +161,23 @@ class Chi(QuantumChannel):
         # representation conversion to SuperOp and then convert back to Chi
         return Chi(Choi(self).compose(other, front=front))
 
-    def power(self, n):
-        """The matrix power of the channel.
-
-        Args:
-            n (int): compute the matrix power of the superoperator matrix.
-
-        Returns:
-            Chi: the matrix power of the SuperOp converted to a Chi channel.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the
-                         QuantumChannel are not equal, or the power is not an integer.
-        """
-        if n > 0:
-            return super().power(n)
-        return Chi(SuperOp(self).power(n))
-
     def tensor(self, other):
-        """Return the tensor product channel self ⊗ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            Chi: the tensor product channel self ⊗ other as a Chi object.
-
-        Raises:
-            QiskitError: if other is not a QuantumChannel subclass.
-        """
         if not isinstance(other, Chi):
             other = Chi(other)
-        input_dims = other.input_dims() + self.input_dims()
-        output_dims = other.output_dims() + self.output_dims()
-        data = np.kron(self._data, other.data)
-        return Chi(data, input_dims, output_dims)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product channel other ⊗ self.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            Chi: the tensor product channel other ⊗ self as a Chi object.
-
-        Raises:
-            QiskitError: if other is not a QuantumChannel subclass.
-        """
         if not isinstance(other, Chi):
             other = Chi(other)
-        input_dims = self.input_dims() + other.input_dims()
-        output_dims = self.output_dims() + other.output_dims()
-        data = np.kron(other.data, self._data)
-        return Chi(data, input_dims, output_dims)
-
-    def _evolve(self, state, qargs=None):
-        """Evolve a quantum state by the quantum channel.
+        return self._tensor(other, self)
 
-        Args:
-            state (DensityMatrix or Statevector): The input state.
-            qargs (list): a list of quantum state subsystem positions to apply
-                           the quantum channel on.
+    @classmethod
+    def _tensor(cls, a, b):
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
+        ret._data = np.kron(a._data, b.data)
+        return ret
 
-        Returns:
-            DensityMatrix: the output quantum state as a density matrix.
 
-        Raises:
-            QiskitError: if the quantum channel dimension does not match the
-                         specified quantum state subsystem dimensions.
-        """
-        return SuperOp(self)._evolve(state, qargs)
+# Update docstrings for API docs
+generate_apidocs(Chi)
diff --git a/qiskit/quantum_info/operators/channel/choi.py b/qiskit/quantum_info/operators/channel/choi.py
index d1177f526..c21dfcddc 100644
--- a/qiskit/quantum_info/operators/channel/choi.py
+++ b/qiskit/quantum_info/operators/channel/choi.py
@@ -16,6 +16,7 @@
 Choi-matrix representation of a Quantum Channel.
 """
 
+import copy
 import numpy as np
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -26,6 +27,7 @@ from qiskit.quantum_info.operators.op_shape import OpShape
 from qiskit.quantum_info.operators.channel.superop import SuperOp
 from qiskit.quantum_info.operators.channel.transformations import _to_choi
 from qiskit.quantum_info.operators.channel.transformations import _bipartite_tensor
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class Choi(QuantumChannel):
@@ -135,46 +137,30 @@ class Choi(QuantumChannel):
         return (self._input_dim, self._output_dim, self._input_dim,
                 self._output_dim)
 
+    def _evolve(self, state, qargs=None):
+        return SuperOp(self)._evolve(state, qargs)
+
+    # ---------------------------------------------------------------------
+    # BaseOperator methods
+    # ---------------------------------------------------------------------
+
     def conjugate(self):
-        """Return the conjugate of the QuantumChannel."""
-        return Choi(np.conj(self._data), self.input_dims(), self.output_dims())
+        ret = copy.copy(self)
+        ret._data = np.conj(self._data)
+        return ret
 
     def transpose(self):
-        """Return the transpose of the QuantumChannel."""
+        ret = copy.copy(self)
+        ret._op_shape = self._op_shape.transpose()
         # Make bipartite matrix
         d_in, d_out = self.dim
         data = np.reshape(self._data, (d_in, d_out, d_in, d_out))
         # Swap input and output indices on bipartite matrix
         data = np.transpose(data, (1, 0, 3, 2))
-        # Transpose channel has input and output dimensions swapped
-        data = np.reshape(data, (d_in * d_out, d_in * d_out))
-        return Choi(data,
-                    input_dims=self.output_dims(),
-                    output_dims=self.input_dims())
+        ret._data = np.reshape(data, (d_in * d_out, d_in * d_out))
+        return ret
 
     def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            Choi: The quantum channel self @ other.
-
-        Raises:
-            QiskitError: if other has incompatible dimensions.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
         if qargs is not None:
@@ -200,84 +186,25 @@ class Choi(QuantumChannel):
         ret._op_shape = new_shape
         return ret
 
-    def power(self, n):
-        """The matrix power of the channel.
-
-        Args:
-            n (int): compute the matrix power of the superoperator matrix.
-
-        Returns:
-            Choi: the matrix power of the SuperOp converted to a Choi channel.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the
-                         QuantumChannel are not equal, or the power is not an integer.
-        """
-        if n > 0:
-            return super().power(n)
-        return Choi(SuperOp(self).power(n))
-
     def tensor(self, other):
-        """Return the tensor product channel self ⊗ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            Choi: the tensor product channel self ⊗ other as a Choi object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        # Convert other to Choi
         if not isinstance(other, Choi):
             other = Choi(other)
-
-        input_dims = other.input_dims() + self.input_dims()
-        output_dims = other.output_dims() + self.output_dims()
-        data = _bipartite_tensor(self._data,
-                                 other.data,
-                                 shape1=self._bipartite_shape,
-                                 shape2=other._bipartite_shape)
-        return Choi(data, input_dims, output_dims)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product channel other ⊗ self.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            Choi: the tensor product channel other ⊗ self as a Choi object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        # Convert other to Choi
         if not isinstance(other, Choi):
             other = Choi(other)
+        return self._tensor(other, self)
+
+    @classmethod
+    def _tensor(cls, a, b):
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
+        ret._data = _bipartite_tensor(a._data, b.data,
+                                      shape1=a._bipartite_shape,
+                                      shape2=b._bipartite_shape)
+        return ret
 
-        input_dims = self.input_dims() + other.input_dims()
-        output_dims = self.output_dims() + other.output_dims()
-        data = _bipartite_tensor(other.data,
-                                 self._data,
-                                 shape1=other._bipartite_shape,
-                                 shape2=self._bipartite_shape)
-        return Choi(data, input_dims, output_dims)
-
-    def _evolve(self, state, qargs=None):
-        """Evolve a quantum state by the quantum channel.
-
-        Args:
-            state (DensityMatrix or Statevector): The input state.
-            qargs (list): a list of quantum state subsystem positions to apply
-                           the quantum channel on.
-
-        Returns:
-            DensityMatrix: the output quantum state as a density matrix.
 
-        Raises:
-            QiskitError: if the quantum channel dimension does not match the
-                         specified quantum state subsystem dimensions.
-        """
-        return SuperOp(self)._evolve(state, qargs)
+# Update docstrings for API docs
+generate_apidocs(Choi)
diff --git a/qiskit/quantum_info/operators/channel/kraus.py b/qiskit/quantum_info/operators/channel/kraus.py
index 37ff57dc8..e598647de 100644
--- a/qiskit/quantum_info/operators/channel/kraus.py
+++ b/qiskit/quantum_info/operators/channel/kraus.py
@@ -28,6 +28,7 @@ from qiskit.quantum_info.operators.op_shape import OpShape
 from qiskit.quantum_info.operators.channel.choi import Choi
 from qiskit.quantum_info.operators.channel.superop import SuperOp
 from qiskit.quantum_info.operators.channel.transformations import _to_kraus
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class Kraus(QuantumChannel):
@@ -190,48 +191,43 @@ class Kraus(QuantumChannel):
             accum += np.dot(np.transpose(np.conj(op)), op)
         return is_identity_matrix(accum, rtol=rtol, atol=atol)
 
+    def _evolve(self, state, qargs=None):
+        return SuperOp(self)._evolve(state, qargs)
+
+    # ---------------------------------------------------------------------
+    # BaseOperator methods
+    # ---------------------------------------------------------------------
+
     def conjugate(self):
-        """Return the conjugate of the QuantumChannel."""
+        ret = copy.copy(self)
         kraus_l, kraus_r = self._data
-        kraus_l = [k.conj() for k in kraus_l]
+        kraus_l = [np.conj(k) for k in kraus_l]
         if kraus_r is not None:
             kraus_r = [k.conj() for k in kraus_r]
-        return Kraus((kraus_l, kraus_r), self.input_dims(), self.output_dims())
+        ret._data = (kraus_l, kraus_r)
+        return ret
 
     def transpose(self):
-        """Return the transpose of the QuantumChannel."""
+        ret = copy.copy(self)
+        ret._op_shape = self._op_shape.transpose()
         kraus_l, kraus_r = self._data
-        kraus_l = [k.T for k in kraus_l]
+        kraus_l = [np.transpose(k) for k in kraus_l]
         if kraus_r is not None:
-            kraus_r = [k.T for k in kraus_r]
-        return Kraus((kraus_l, kraus_r),
-                     input_dims=self.output_dims(),
-                     output_dims=self.input_dims())
-
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            Kraus: The quantum channel self @ other.
+            kraus_r = [np.transpose(k) for k in kraus_r]
+        ret._data = (kraus_l, kraus_r)
+        return ret
 
-        Raises:
-            QiskitError: if other cannot be converted to a Kraus or has
-                         incompatible dimensions.
+    def adjoint(self):
+        ret = copy.copy(self)
+        ret._op_shape = self._op_shape.transpose()
+        kraus_l, kraus_r = self._data
+        kraus_l = [np.conj(np.transpose(k)) for k in kraus_l]
+        if kraus_r is not None:
+            kraus_r = [np.conj(np.transpose(k)) for k in kraus_r]
+        ret._data = (kraus_l, kraus_r)
+        return ret
 
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
+    def compose(self, other, qargs=None, front=False):
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
         if qargs is not None:
@@ -264,107 +260,55 @@ class Kraus(QuantumChannel):
         ret._op_shape = new_shape
         return ret
 
-    def dot(self, other, qargs=None):
-        """Return the right multiplied quantum channel self * other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-
-        Returns:
-            Kraus: The quantum channel self * other.
-
-        Raises:
-            QiskitError: if other cannot be converted to a Kraus or has
-                         incompatible dimensions.
-        """
-        return super().dot(other, qargs=qargs)
-
-    def power(self, n):
-        """The matrix power of the channel.
-
-        Args:
-            n (int): compute the matrix power of the superoperator matrix.
-
-        Returns:
-            Kraus: the matrix power of the SuperOp converted to a Kraus channel.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the
-                         QuantumChannel are not equal, or the power
-                         is not an integer.
-        """
-        if n > 0:
-            return super().power(n)
-        return Kraus(SuperOp(self).power(n))
-
     def tensor(self, other):
-        """Return the tensor product channel self ⊗ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-
-        Returns:
-            Kraus: the tensor product channel self ⊗ other as a Kraus
-            object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        return self._tensor_product(other, reverse=False)
+        if not isinstance(other, Kraus):
+            other = Kraus(other)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product channel other ⊗ self.
-
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-
-        Returns:
-            Kraus: the tensor product channel other ⊗ self as a Kraus
-            object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        return self._tensor_product(other, reverse=True)
+        if not isinstance(other, Kraus):
+            other = Kraus(other)
+        return self._tensor(other, self)
 
-    def _add(self, other, qargs=None):
-        """Return the QuantumChannel self + other.
+    @classmethod
+    def _tensor(cls, a, b):
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
 
-        If ``qargs`` are specified the other operator will be added
-        assuming it is identity on all other subsystems.
+        # Get tensor matrix
+        ka_l, ka_r = a._data
+        kb_l, kb_r = b._data
+        kab_l = [np.kron(ka, kb) for ka in ka_l for kb in kb_l]
+        if ka_r is None and kb_r is None:
+            kab_r = None
+        else:
+            if ka_r is None:
+                ka_r = ka_l
+            if kb_r is None:
+                kb_r = kb_l
+            kab_r = [np.kron(a, b) for a in ka_r for b in kb_r]
+        ret._data = (kab_l, kab_r)
+        return ret
 
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-            qargs (None or list): optional subsystems to add on
-                                  (Default: None)
+    def __add__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        if not isinstance(other, QuantumChannel):
+            other = Choi(other)
+        return self._add(other, qargs=qargs)
 
-        Returns:
-            Kraus: the linear addition channel self + other.
+    def __sub__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        if not isinstance(other, QuantumChannel):
+            other = Choi(other)
+        return self._add(-other, qargs=qargs)
 
-        Raises:
-            QiskitError: if other cannot be converted to a channel, or
-                         has incompatible dimensions.
-        """
+    def _add(self, other, qargs=None):
         # Since we cannot directly add two channels in the Kraus
         # representation we try and use the other channels method
         # or convert to the Choi representation
         return Kraus(Choi(self)._add(other, qargs=qargs))
 
     def _multiply(self, other):
-        """Return the QuantumChannel other * self.
-
-        Args:
-            other (complex): a complex number.
-
-        Returns:
-            Kraus: the scalar multiplication other * self as a Kraus object.
-
-        Raises:
-            QiskitError: if other is not a valid scalar.
-        """
         if not isinstance(other, Number):
             raise QiskitError("other is not a number")
 
@@ -385,61 +329,6 @@ class Kraus(QuantumChannel):
         ret._data = (kraus_l, kraus_r)
         return ret
 
-    def _evolve(self, state, qargs=None):
-        """Evolve a quantum state by the quantum channel.
-
-        Args:
-            state (DensityMatrix or Statevector): The input state.
-            qargs (list): a list of quantum state subsystem positions to apply
-                           the quantum channel on.
-
-        Returns:
-            DensityMatrix: the output quantum state as a density matrix.
-
-        Raises:
-            QiskitError: if the quantum channel dimension does not match the
-                         specified quantum state subsystem dimensions.
-        """
-        return SuperOp(self)._evolve(state, qargs)
 
-    def _tensor_product(self, other, reverse=False):
-        """Return the tensor product channel.
-
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-            reverse (bool): If False return self ⊗ other, if True return
-                            if True return (other ⊗ self) [Default: False
-        Returns:
-            Kraus: the tensor product channel as a Kraus object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        # Convert other to Kraus
-        if not isinstance(other, Kraus):
-            other = Kraus(other)
-
-        # Get tensor matrix
-        ka_l, ka_r = self._data
-        kb_l, kb_r = other._data
-        if reverse:
-            input_dims = self.input_dims() + other.input_dims()
-            output_dims = self.output_dims() + other.output_dims()
-            kab_l = [np.kron(b_in, a_in) for a_in in ka_l for b_in in kb_l]
-        else:
-            input_dims = other.input_dims() + self.input_dims()
-            output_dims = other.output_dims() + self.output_dims()
-            kab_l = [np.kron(a, b) for a in ka_l for b in kb_l]
-        if ka_r is None and kb_r is None:
-            kab_r = None
-        else:
-            if ka_r is None:
-                ka_r = ka_l
-            if kb_r is None:
-                kb_r = kb_l
-            if reverse:
-                kab_r = [np.kron(b_in, a_in) for a_in in ka_r for b_in in kb_r]
-            else:
-                kab_r = [np.kron(a, b) for a in ka_r for b in kb_r]
-        data = (kab_l, kab_r)
-        return Kraus(data, input_dims, output_dims)
+# Update docstrings for API docs
+generate_apidocs(Kraus)
diff --git a/qiskit/quantum_info/operators/channel/ptm.py b/qiskit/quantum_info/operators/channel/ptm.py
index c4841c439..c9241331b 100644
--- a/qiskit/quantum_info/operators/channel/ptm.py
+++ b/qiskit/quantum_info/operators/channel/ptm.py
@@ -16,6 +16,7 @@
 Pauli Transfer Matrix (PTM) representation of a Quantum Channel.
 """
 
+import copy
 import numpy as np
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -24,6 +25,7 @@ from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.channel.quantum_channel import QuantumChannel
 from qiskit.quantum_info.operators.channel.superop import SuperOp
 from qiskit.quantum_info.operators.channel.transformations import _to_ptm
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class PTM(QuantumChannel):
@@ -135,43 +137,26 @@ class PTM(QuantumChannel):
         return (self._output_dim, self._output_dim, self._input_dim,
                 self._input_dim)
 
+    def _evolve(self, state, qargs=None):
+        return SuperOp(self)._evolve(state, qargs)
+
+    # ---------------------------------------------------------------------
+    # BaseOperator methods
+    # ---------------------------------------------------------------------
+
     def conjugate(self):
-        """Return the conjugate of the QuantumChannel."""
         # Since conjugation is basis dependent we transform
         # to the SuperOp representation to compute the
         # conjugate channel
         return PTM(SuperOp(self).conjugate())
 
     def transpose(self):
-        """Return the transpose of the QuantumChannel."""
-        # Since conjugation is basis dependent we transform
-        # to the SuperOp representation to compute the
-        # conjugate channel
         return PTM(SuperOp(self).transpose())
 
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            PTM: The quantum channel self @ other.
+    def adjoint(self):
+        return PTM(SuperOp(self).adjoint())
 
-        Raises:
-            QiskitError: if other has incompatible dimensions.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
+    def compose(self, other, qargs=None, front=False):
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
         if qargs is not None:
@@ -192,75 +177,23 @@ class PTM(QuantumChannel):
         ret._op_shape = new_shape
         return ret
 
-    def power(self, n):
-        """The matrix power of the channel.
-
-        Args:
-            n (int): compute the matrix power of the superoperator matrix.
-
-        Returns:
-            PTM: the matrix power of the SuperOp converted to a PTM channel.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the
-                         QuantumChannel are not equal, or the power is not
-                         an integer.
-        """
-        if n > 0:
-            return super().power(n)
-        return PTM(SuperOp(self).power(n))
-
     def tensor(self, other):
-        """Return the tensor product channel self ⊗ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            PTM: the tensor product channel self ⊗ other as a PTM object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
         if not isinstance(other, PTM):
             other = PTM(other)
-        input_dims = other.input_dims() + self.input_dims()
-        output_dims = other.output_dims() + self.output_dims()
-        data = np.kron(self._data, other.data)
-        return PTM(data, input_dims, output_dims)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product channel other ⊗ self.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            PTM: the tensor product channel other ⊗ self as a PTM object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
         if not isinstance(other, PTM):
             other = PTM(other)
-        input_dims = self.input_dims() + other.input_dims()
-        output_dims = self.output_dims() + other.output_dims()
-        data = np.kron(other.data, self._data)
-        return PTM(data, input_dims, output_dims)
-
-    def _evolve(self, state, qargs=None):
-        """Evolve a quantum state by the quantum channel.
+        return self._tensor(other, self)
 
-        Args:
-            state (DensityMatrix or Statevector): The input state.
-            qargs (list): a list of quantum state subsystem positions to apply
-                           the quantum channel on.
+    @classmethod
+    def _tensor(cls, a, b):
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
+        ret._data = np.kron(a._data, b.data)
+        return ret
 
-        Returns:
-            DensityMatrix: the output quantum state as a density matrix.
 
-        Raises:
-            QiskitError: if the quantum channel dimension does not match the
-                         specified quantum state subsystem dimensions.
-        """
-        return SuperOp(self)._evolve(state, qargs)
+# Update docstrings for API docs
+generate_apidocs(PTM)
diff --git a/qiskit/quantum_info/operators/channel/quantum_channel.py b/qiskit/quantum_info/operators/channel/quantum_channel.py
index 882f94532..3415e528e 100644
--- a/qiskit/quantum_info/operators/channel/quantum_channel.py
+++ b/qiskit/quantum_info/operators/channel/quantum_channel.py
@@ -16,21 +16,21 @@ Abstract base class for Quantum Channels.
 
 import copy
 from abc import abstractmethod
-from numbers import Number
+from numbers import Number, Integral
 import numpy as np
 
 from qiskit.exceptions import QiskitError
-from qiskit.quantum_info.operators.base_operator import BaseOperator
-from qiskit.quantum_info.operators.tolerances import TolerancesMixin
+from qiskit.quantum_info.operators.linear_op import LinearOp
 from qiskit.quantum_info.operators.operator import Operator
 from qiskit.quantum_info.operators.predicates import is_positive_semidefinite_matrix
+from qiskit.quantum_info.operators.channel.transformations import _transform_rep
 from qiskit.quantum_info.operators.channel.transformations import _to_choi
 from qiskit.quantum_info.operators.channel.transformations import _to_kraus
 from qiskit.quantum_info.operators.channel.transformations import _to_operator
 from qiskit.quantum_info.operators.scalar_op import ScalarOp
 
 
-class QuantumChannel(BaseOperator, TolerancesMixin):
+class QuantumChannel(LinearOp):
     """Quantum channel representation base class."""
 
     def __init__(self, data, num_qubits=None, op_shape=None):
@@ -72,90 +72,103 @@ class QuantumChannel(BaseOperator, TolerancesMixin):
         """Return channel representation string"""
         return type(self).__name__
 
-    @abstractmethod
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
+    # ---------------------------------------------------------------------
+    # LinearOp methods
+    # ---------------------------------------------------------------------
 
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
+    @abstractmethod
+    def conjugate(self):
+        r"""Return the conjugate quantum channel.
+
+        .. note::
+            This is equivalent to the matrix complex conjugate in the
+            :class:`~qiskit.quantum_info.SuperOp` representation
+            ie. for a channel :math:`\mathcal{E}`, the SuperOp of
+            the conjugate channel :math:`\overline{{\mathcal{{E}}}}` is
+            :math:`S_{\overline{\mathcal{E}^\dagger}} = \overline{S_{\mathcal{E}}}`.
+        """
 
-        Returns:
-            QuantumChannel: The quantum channel self @ other.
+    @abstractmethod
+    def transpose(self):
+        r"""Return the transpose quantum channel.
+
+        .. note::
+            This is equivalent to the matrix tranpose in the
+            :class:`~qiskit.quantum_info.SuperOp` representation,
+            ie. for a channel :math:`\mathcal{E}`, the SuperOp of
+            the transpose channel :math:`\mathcal{{E}}^T` is
+            :math:`S_{mathcal{E}^T} = S_{\mathcal{E}}^T`.
+        """
 
-        Raises:
-            QiskitError: if other has incompatible dimensions.
+    def adjoint(self):
+        r"""Return the adjoint quantum channel.
 
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
+        .. note::
+            This is equivalent to the matrix Hermitian conjugate in the
+            :class:`~qiskit.quantum_info.SuperOp` representation
+            ie. for a channel :math:`\mathcal{E}`, the SuperOp of
+            the adjoint channel :math:`\mathcal{{E}}^\dagger` is
+            :math:`S_{\mathcal{E}^\dagger} = S_{\mathcal{E}}^\dagger`.
         """
-        pass
+        return self.conjugate().transpose()
 
-    def _add(self, other, qargs=None):
-        """Return the QuantumChannel self + other.
-
-        If ``qargs`` are specified the other channel will be added
-        assuming it is the identity channel on all other subsystems.
+    def power(self, n):
+        r"""Return the power of the quantum channel.
 
         Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (None or list): optional subsystems to add on
-                                  (Default: None)
+            n (float): the power exponent.
 
         Returns:
-            QuantumChannel: the linear addition self + other as a SuperOp object.
+            CLASS: the channel :math:`\mathcal{{E}} ^n`.
 
         Raises:
-            QiskitError: if other cannot be converted to a channel or
-                         has incompatible dimensions.
+            QiskitError: if the input and output dimensions of the
+                         CLASS are not equal.
+
+        .. note::
+            For non-positive or non-integer exponents the power is
+            defined as the matrix power of the
+            :class:`~qiskit.quantum_info.SuperOp` representation
+            ie. for a channel :math:`\mathcal{{E}}`, the SuperOp of
+            the powered channel :math:`\mathcal{{E}}^\n` is
+            :math:`S_{{\mathcal{{E}}^n}} = S_{{\mathcal{{E}}}}^n`.
         """
-        # NOTE: this method must be overriden for subclasses
-        # that don't have a linear matrix representation
-        # ie Kraus and Stinespring
-
-        if qargs is None:
-            qargs = getattr(other, 'qargs', None)
+        if n > 0 and isinstance(n, Integral):
+            return super().power(n)
+
+        # Conversion to superoperator
+        if self._input_dim != self._output_dim:
+            raise QiskitError("Can only take power with input_dim = output_dim.")
+        rep = self._channel_rep
+        input_dim, output_dim = self.dim
+        superop = _transform_rep(rep, 'SuperOp', self._data, input_dim, output_dim)
+        superop = np.linalg.matrix_power(superop, n)
+
+        # Convert back to original representation
+        ret = copy.copy(self)
+        ret._data = _transform_rep('SuperOp', rep, superop, input_dim, output_dim)
+        return ret
 
+    def __sub__(self, other):
+        qargs = getattr(other, 'qargs', None)
         if not isinstance(other, type(self)):
             other = type(self)(other)
+        return self._add(-other, qargs=qargs)
 
+    def _add(self, other, qargs=None):
+        # NOTE: this method must be overridden for subclasses
+        # that don't have a linear matrix representation
+        # ie Kraus and Stinespring
+        if not isinstance(other, type(self)):
+            other = type(self)(other)
         self._op_shape._validate_add(other._op_shape, qargs)
         other = ScalarOp._pad_with_identity(self, other, qargs)
-
         ret = copy.copy(self)
         ret._data = self._data + other._data
         return ret
 
-    def __sub__(self, other):
-        # Override for sub so that other is converted before being negated
-        # which can lead to problems if other is not already a quantum channel
-        if not isinstance(other, QuantumChannel):
-            qargs = getattr(other, 'qargs', None)
-            other = type(self)(other)
-            if qargs is not None:
-                other.qargs = qargs
-        return self._add(-other)
-
     def _multiply(self, other):
-        """Return the QuantumChannel other * self.
-
-        Args:
-            other (complex): a complex number.
-
-        Returns:
-            QuantumChannel: the scalar multiplication other * self.
-
-        Raises:
-            QiskitError: if other is not a valid scalar.
-        """
-        # NOTE: this method must be overriden for subclasses
+        # NOTE: this method must be overridden for subclasses
         # that don't have a linear matrix representation
         # ie Kraus and Stinespring
         if not isinstance(other, Number):
@@ -164,6 +177,10 @@ class QuantumChannel(BaseOperator, TolerancesMixin):
         ret._data = other * self._data
         return ret
 
+    # ---------------------------------------------------------------------
+    # Additional methods
+    # ---------------------------------------------------------------------
+
     def is_cptp(self, atol=None, rtol=None):
         """Return True if completely-positive trace-preserving (CPTP)."""
         choi = _to_choi(self._channel_rep, self._data, *self.dim)
diff --git a/qiskit/quantum_info/operators/channel/stinespring.py b/qiskit/quantum_info/operators/channel/stinespring.py
index 80a237a6d..33d03291f 100644
--- a/qiskit/quantum_info/operators/channel/stinespring.py
+++ b/qiskit/quantum_info/operators/channel/stinespring.py
@@ -27,6 +27,7 @@ from qiskit.quantum_info.operators.channel.kraus import Kraus
 from qiskit.quantum_info.operators.channel.choi import Choi
 from qiskit.quantum_info.operators.channel.superop import SuperOp
 from qiskit.quantum_info.operators.channel.transformations import _to_stinespring
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class Stinespring(QuantumChannel):
@@ -155,18 +156,26 @@ class Stinespring(QuantumChannel):
         check = np.dot(np.transpose(np.conj(self._data[0])), self._data[0])
         return is_identity_matrix(check, rtol=self.rtol, atol=self.atol)
 
+    def _evolve(self, state, qargs=None):
+        return SuperOp(self)._evolve(state, qargs)
+
+    # ---------------------------------------------------------------------
+    # BaseOperator methods
+    # ---------------------------------------------------------------------
+
     def conjugate(self):
-        """Return the conjugate of the QuantumChannel."""
         # pylint: disable=assignment-from-no-return
+        ret = copy.copy(self)
         stine_l = np.conjugate(self._data[0])
         stine_r = None
         if self._data[1] is not None:
             stine_r = np.conjugate(self._data[1])
-        return Stinespring((stine_l, stine_r), self.input_dims(),
-                           self.output_dims())
+        ret._data = (stine_l, stine_r)
+        return ret
 
     def transpose(self):
-        """Return the transpose of the QuantumChannel."""
+        ret = copy.copy(self)
+        ret._op_shape = self._op_shape.transpose()
         din, dout = self.dim
         dtr = self._data[0].shape[0] // dout
         stine = [None, None]
@@ -175,145 +184,85 @@ class Stinespring(QuantumChannel):
                 stine[i] = np.reshape(
                     np.transpose(np.reshape(mat, (dout, dtr, din)), (2, 1, 0)),
                     (din * dtr, dout))
-        return Stinespring(tuple(stine),
-                           input_dims=self.output_dims(),
-                           output_dims=self.input_dims())
+        ret._data = (stine[0], stine[1])
+        return ret
 
     def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            Stinespring: The quantum channel self @ other.
-
-        Raises:
-            QiskitError: if other cannot be converted to a Stinespring or has
-                         incompatible dimensions.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
         if qargs is not None:
             return Stinespring(
                 SuperOp(self).compose(other, qargs=qargs, front=front))
-
         # Otherwise we convert via Kraus representation rather than
         # superoperator to avoid unnecessary representation conversions
         return Stinespring(Kraus(self).compose(other, front=front))
 
-    def dot(self, other, qargs=None):
-        """Return the right multiplied quantum channel self * other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-
-        Returns:
-            Stinespring: The quantum channel self * other.
-
-        Raises:
-            QiskitError: if other cannot be converted to a Stinespring or has
-                         incompatible dimensions.
-        """
-        return super().dot(other, qargs=qargs)
-
-    def power(self, n):
-        """The matrix power of the channel.
-
-        Args:
-            n (int): compute the matrix power of the superoperator matrix.
-
-        Returns:
-            Stinespring: the matrix power of the SuperOp converted to a
-            Stinespring channel.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the
-                         QuantumChannel are not equal, or the power is not
-                         an integer.
-        """
-        if n > 0:
-            return super().power(n)
-        return Stinespring(SuperOp(self).power(n))
-
     def tensor(self, other):
-        """Return the tensor product channel self ⊗ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-
-        Returns:
-            Stinespring: the tensor product channel other ⊗ self as a
-            Stinespring object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        return self._tensor_product(other, reverse=False)
+        if not isinstance(other, Stinespring):
+            other = Stinespring(other)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product channel other ⊗ self.
-
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
+        if not isinstance(other, Stinespring):
+            other = Stinespring(other)
+        return self._tensor(other, self)
 
-        Returns:
-            Stinespring: the tensor product channel other ⊗ self as a
-            Stinespring object.
+    @classmethod
+    def _tensor(cls, a, b):
+        # Tensor Stinespring ops
+        sa_l, sa_r = a._data
+        sb_l, sb_r = b._data
 
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        return self._tensor_product(other, reverse=True)
+        # Reshuffle tensor dimensions
+        din_a, dout_a = a.dim
+        din_b, dout_b = b.dim
+        dtr_a = sa_l.shape[0] // dout_a
+        dtr_b = sb_l.shape[0] // dout_b
 
-    def _add(self, other, qargs=None):
-        """Return the QuantumChannel self + other.
+        shape_in = (dout_a, dtr_a, dout_b, dtr_b, din_a * din_b)
+        shape_out = (dout_a * dtr_a * dout_b * dtr_b, din_a * din_b)
+        sab_l = np.kron(sa_l, sb_l)
+        # Reravel indices
+        sab_l = np.reshape(
+            np.transpose(np.reshape(sab_l, shape_in), (0, 2, 1, 3, 4)),
+            shape_out)
 
-        If ``qargs`` are specified the other operator will be added
-        assuming it is identity on all other subsystems.
+        # Compute right Stinespring op
+        if sa_r is None and sb_r is None:
+            sab_r = None
+        else:
+            if sa_r is None:
+                sa_r = sa_l
+            elif sb_r is None:
+                sb_r = sb_l
+            sab_r = np.kron(sa_r, sb_r)
+            # Reravel indices
+            sab_r = np.reshape(
+                np.transpose(np.reshape(sab_r, shape_in), (0, 2, 1, 3, 4)),
+                shape_out)
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
+        ret._data = (sab_l, sab_r)
+        return ret
 
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-            qargs (None or list): optional subsystems to add on
-                                  (Default: None)
+    def __add__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        if not isinstance(other, QuantumChannel):
+            other = Choi(other)
+        return self._add(other, qargs=qargs)
 
-        Returns:
-            Stinespring: the linear addition channel self + other.
+    def __sub__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        if not isinstance(other, QuantumChannel):
+            other = Choi(other)
+        return self._add(-other, qargs=qargs)
 
-        Raises:
-            QiskitError: if other cannot be converted to a channel or
-                         has incompatible dimensions.
-        """
+    def _add(self, other, qargs=None):
         # Since we cannot directly add two channels in the Stinespring
         # representation we convert to the Choi representation
         return Stinespring(Choi(self)._add(other, qargs=qargs))
 
     def _multiply(self, other):
-        """Return the QuantumChannel other * self.
-
-        Args:
-            other (complex): a complex number.
-
-        Returns:
-            Stinespring: the scalar multiplication other * self.
-
-        Raises:
-            QiskitError: if other is not a valid scalar.
-        """
         if not isinstance(other, Number):
             raise QiskitError("other is not a number")
 
@@ -336,84 +285,6 @@ class Stinespring(QuantumChannel):
         ret._data = (stine_l, stine_r)
         return ret
 
-    def _evolve(self, state, qargs=None):
-        """Evolve a quantum state by the quantum channel.
-
-        Args:
-            state (DensityMatrix or Statevector): The input state.
-            qargs (list): a list of quantum state subsystem positions to apply
-                           the quantum channel on.
-
-        Returns:
-            DensityMatrix: the output quantum state as a density matrix.
 
-        Raises:
-            QiskitError: if the quantum channel dimension does not match the
-                         specified quantum state subsystem dimensions.
-        """
-        return SuperOp(self)._evolve(state, qargs)
-
-    def _tensor_product(self, other, reverse=False):
-        """Return the tensor product channel.
-
-        Args:
-            other (QuantumChannel): a quantum channel subclass.
-            reverse (bool): If False return self ⊗ other, if True return
-                            if True return (other ⊗ self) [Default: False]
-        Returns:
-            Stinespring: the tensor product channel as a Stinespring object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        # Convert other to Stinespring
-        if not isinstance(other, Stinespring):
-            other = Stinespring(other)
-
-        # Tensor Stinespring ops
-        sa_l, sa_r = self._data
-        sb_l, sb_r = other._data
-
-        # Reshuffle tensor dimensions
-        din_a, dout_a = self.dim
-        din_b, dout_b = other.dim
-        dtr_a = sa_l.shape[0] // dout_a
-        dtr_b = sb_l.shape[0] // dout_b
-        if reverse:
-            shape_in = (dout_b, dtr_b, dout_a, dtr_a, din_b * din_a)
-            shape_out = (dout_b * dtr_b * dout_a * dtr_a, din_b * din_a)
-        else:
-            shape_in = (dout_a, dtr_a, dout_b, dtr_b, din_a * din_b)
-            shape_out = (dout_a * dtr_a * dout_b * dtr_b, din_a * din_b)
-
-        # Compute left Stinespring op
-        if reverse:
-            input_dims = self.input_dims() + other.input_dims()
-            output_dims = self.output_dims() + other.output_dims()
-            sab_l = np.kron(sb_l, sa_l)
-        else:
-            input_dims = other.input_dims() + self.input_dims()
-            output_dims = other.output_dims() + self.output_dims()
-            sab_l = np.kron(sa_l, sb_l)
-        # Reravel indices
-        sab_l = np.reshape(
-            np.transpose(np.reshape(sab_l, shape_in), (0, 2, 1, 3, 4)),
-            shape_out)
-
-        # Compute right Stinespring op
-        if sa_r is None and sb_r is None:
-            sab_r = None
-        else:
-            if sa_r is None:
-                sa_r = sa_l
-            elif sb_r is None:
-                sb_r = sb_l
-            if reverse:
-                sab_r = np.kron(sb_r, sa_r)
-            else:
-                sab_r = np.kron(sa_r, sb_r)
-            # Reravel indices
-            sab_r = np.reshape(
-                np.transpose(np.reshape(sab_r, shape_in), (0, 2, 1, 3, 4)),
-                shape_out)
-        return Stinespring((sab_l, sab_r), input_dims, output_dims)
+# Update docstrings for API docs
+generate_apidocs(Stinespring)
diff --git a/qiskit/quantum_info/operators/channel/superop.py b/qiskit/quantum_info/operators/channel/superop.py
index 4211fe619..d6c991bdc 100644
--- a/qiskit/quantum_info/operators/channel/superop.py
+++ b/qiskit/quantum_info/operators/channel/superop.py
@@ -15,6 +15,7 @@
 """
 Superoperator representation of a Quantum Channel."""
 
+import copy
 import numpy as np
 
 from qiskit.circuit.quantumcircuit import QuantumCircuit
@@ -25,6 +26,7 @@ from qiskit.quantum_info.operators.operator import Operator
 from qiskit.quantum_info.operators.channel.quantum_channel import QuantumChannel
 from qiskit.quantum_info.operators.channel.transformations import _to_superop
 from qiskit.quantum_info.operators.channel.transformations import _bipartite_tensor
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class SuperOp(QuantumChannel):
@@ -129,43 +131,49 @@ class SuperOp(QuantumChannel):
         return (self._output_dim, self._output_dim, self._input_dim,
                 self._input_dim)
 
+    # ---------------------------------------------------------------------
+    # BaseOperator methods
+    # ---------------------------------------------------------------------
+
     def conjugate(self):
-        """Return the conjugate of the QuantumChannel."""
-        return SuperOp(np.conj(self._data), self.input_dims(),
-                       self.output_dims())
+        ret = copy.copy(self)
+        ret._data = np.conj(self._data)
+        return ret
 
     def transpose(self):
-        """Return the transpose of the QuantumChannel."""
-        return SuperOp(np.transpose(self._data),
-                       input_dims=self.output_dims(),
-                       output_dims=self.input_dims())
+        ret = copy.copy(self)
+        ret._data = np.transpose(self._data)
+        ret._op_shape = self._op_shape.transpose()
+        return ret
 
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed quantum channel self @ other.
+    def adjoint(self):
+        ret = copy.copy(self)
+        ret._data = np.conj(np.transpose(self._data))
+        ret._op_shape = self._op_shape.transpose()
+        return ret
 
-        Args:
-            other (QuantumChannel): a quantum channel.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
+    def tensor(self, other):
+        if not isinstance(other, SuperOp):
+            other = SuperOp(other)
+        return self._tensor(self, other)
 
-        Returns:
-            SuperOp: The quantum channel self @ other.
+    def expand(self, other):
+        if not isinstance(other, SuperOp):
+            other = SuperOp(other)
+        return self._tensor(other, self)
 
-        Raises:
-            QiskitError: if other has incompatible dimensions.
+    @classmethod
+    def _tensor(cls, a, b):
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
+        ret._data = _bipartite_tensor(a._data, b.data,
+                                      shape1=a._bipartite_shape,
+                                      shape2=b._bipartite_shape)
+        return ret
 
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            :class:`SuperOp` matrices. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
+    def compose(self, other, qargs=None, front=False):
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
-        # Convert other to SuperOp
         if not isinstance(other, SuperOp):
             other = SuperOp(other)
         # Validate dimensions are compatible and return the composed
@@ -212,78 +220,9 @@ class SuperOp(QuantumChannel):
         ret._op_shape = new_shape
         return ret
 
-    def power(self, n):
-        """Return the compose of a QuantumChannel with itself n times.
-
-        Args:
-            n (int): compute the matrix power of the superoperator matrix.
-
-        Returns:
-            SuperOp: the n-times composition channel as a SuperOp object.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the
-                         QuantumChannel are not equal, or the power is not
-                         an integer.
-        """
-        if not isinstance(n, (int, np.integer)):
-            raise QiskitError("Can only power with integer powers.")
-        if self._input_dim != self._output_dim:
-            raise QiskitError("Can only power with input_dim = output_dim.")
-        # Override base class power so we can implement more efficiently
-        # using Numpy.matrix_power
-        return SuperOp(np.linalg.matrix_power(self._data, n),
-                       self.input_dims(), self.output_dims())
-
-    def tensor(self, other):
-        """Return the tensor product channel self ⊗ other.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            SuperOp: the tensor product channel self ⊗ other as a SuperOp
-            object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        # Convert other to SuperOp
-        if not isinstance(other, SuperOp):
-            other = SuperOp(other)
-
-        input_dims = other.input_dims() + self.input_dims()
-        output_dims = other.output_dims() + self.output_dims()
-        data = _bipartite_tensor(self._data,
-                                 other.data,
-                                 shape1=self._bipartite_shape,
-                                 shape2=other._bipartite_shape)
-        return SuperOp(data, input_dims, output_dims)
-
-    def expand(self, other):
-        """Return the tensor product channel other ⊗ self.
-
-        Args:
-            other (QuantumChannel): a quantum channel.
-
-        Returns:
-            SuperOp: the tensor product channel other ⊗ self as a SuperOp
-            object.
-
-        Raises:
-            QiskitError: if other cannot be converted to a channel.
-        """
-        # Convert other to SuperOp
-        if not isinstance(other, SuperOp):
-            other = SuperOp(other)
-
-        input_dims = self.input_dims() + other.input_dims()
-        output_dims = self.output_dims() + other.output_dims()
-        data = _bipartite_tensor(other.data,
-                                 self._data,
-                                 shape1=other._bipartite_shape,
-                                 shape2=self._bipartite_shape)
-        return SuperOp(data, input_dims, output_dims)
+    # ---------------------------------------------------------------------
+    # Additional methods
+    # ---------------------------------------------------------------------
 
     def _evolve(self, state, qargs=None):
         """Evolve a quantum state by the quantum channel.
@@ -417,3 +356,7 @@ class SuperOp(QuantumChannel):
                 else:
                     new_qargs = [qargs[tup.index] for tup in qregs]
                 self._append_instruction(instr, qargs=new_qargs)
+
+
+# Update docstrings for API docs
+generate_apidocs(SuperOp)
diff --git a/qiskit/quantum_info/operators/channel/transformations.py b/qiskit/quantum_info/operators/channel/transformations.py
index a9baa5510..f129df0f4 100644
--- a/qiskit/quantum_info/operators/channel/transformations.py
+++ b/qiskit/quantum_info/operators/channel/transformations.py
@@ -25,6 +25,27 @@ from qiskit.quantum_info.operators.predicates import is_hermitian_matrix
 from qiskit.quantum_info.operators.predicates import ATOL_DEFAULT
 
 
+def _transform_rep(input_rep, output_rep, data, input_dim, output_dim):
+    """Transform a QuantumChannel between representation."""
+    if input_rep == output_rep:
+        return data
+    if output_rep == 'Choi':
+        return _to_choi(input_rep, data, input_dim, output_dim)
+    if output_rep == 'Operator':
+        return _to_operator(input_rep, data, input_dim, output_dim)
+    if output_rep == 'SuperOp':
+        return _to_superop(input_rep, data, input_dim, output_dim)
+    if output_rep == 'Kraus':
+        return _to_kraus(input_rep, data, input_dim, output_dim)
+    if output_rep == 'Chi':
+        return _to_chi(input_rep, data, input_dim, output_dim)
+    if output_rep == 'PTM':
+        return _to_ptm(input_rep, data, input_dim, output_dim)
+    if output_rep == 'Stinespring':
+        return _to_stinespring(input_rep, data, input_dim, output_dim)
+    raise QiskitError('Invalid QuantumChannel {}'.format(output_rep))
+
+
 def _to_choi(rep, data, input_dim, output_dim):
     """Transform a QuantumChannel to the Choi representation."""
     if rep == 'Choi':
diff --git a/qiskit/quantum_info/operators/linear_op.py b/qiskit/quantum_info/operators/linear_op.py
new file mode 100644
index 000000000..c77fcb606
--- /dev/null
+++ b/qiskit/quantum_info/operators/linear_op.py
@@ -0,0 +1,29 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Abstract LinearOperator class.
+"""
+
+from abc import ABC
+
+from .base_operator import BaseOperator
+from .mixins import LinearMixin, AdjointMixin, TolerancesMixin
+
+
+# pylint: disable = abstract-method
+class LinearOp(BaseOperator,
+               AdjointMixin,
+               LinearMixin,
+               TolerancesMixin,
+               ABC):
+    """Abstract linear operator base class. """
diff --git a/qiskit/quantum_info/operators/mixins/__init__.py b/qiskit/quantum_info/operators/mixins/__init__.py
new file mode 100644
index 000000000..568879538
--- /dev/null
+++ b/qiskit/quantum_info/operators/mixins/__init__.py
@@ -0,0 +1,51 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Operator Mixins
+"""
+
+from inspect import getdoc
+
+from .group import GroupMixin
+from .adjoint import AdjointMixin
+from .linear import LinearMixin
+from .multiply import MultiplyMixin
+from .tolerances import TolerancesMixin
+
+
+def generate_apidocs(cls):
+    """Decorator to format API docstrings for classes using Mixins.
+
+    This runs string replacement on the docstrings of the mixin
+    methods to replace the placeholder CLASS with the class
+    name `cls.__name__`.
+
+    Args:
+        cls (type): The class to format docstrings.
+
+    Returns:
+        cls: the original class with updated docstrings.
+    """
+    def _replace_name(mixin, methods):
+        if issubclass(cls, mixin):
+            for i in methods:
+                meth = getattr(cls, i)
+                doc = getdoc(meth)
+                if doc is not None:
+                    meth.__doc__ = doc.replace('CLASS', cls.__name__)
+
+    _replace_name(GroupMixin, ('tensor', 'expand', 'compose', 'dot', 'power'))
+    _replace_name(AdjointMixin, ('transpose', 'conjugate', 'adjoint'))
+    _replace_name(MultiplyMixin, ('_multiply',))
+    _replace_name(LinearMixin, ('_add',))
+    return cls
diff --git a/qiskit/quantum_info/operators/mixins/adjoint.py b/qiskit/quantum_info/operators/mixins/adjoint.py
new file mode 100644
index 000000000..f511c4e20
--- /dev/null
+++ b/qiskit/quantum_info/operators/mixins/adjoint.py
@@ -0,0 +1,47 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Mixin for gate operator interface.
+"""
+# pylint: disable=abstract-method
+
+from abc import ABC, abstractmethod
+
+
+class AdjointMixin(ABC):
+    """Abstract Mixin for operator adjoint and transpose operations.
+
+    This class defines the following methods
+
+        - :meth:`transpose`
+        - :meth:`conjugate`
+        - :meth:`adjoint`
+
+    The following abstract methods must be implemented by subclasses
+    using this mixin
+
+        - ``conjugate(self)``
+        - ``transpose(self)``
+    """
+
+    def adjoint(self):
+        """Return the adjoint of the CLASS."""
+        return self.conjugate().transpose()
+
+    @abstractmethod
+    def conjugate(self):
+        """Return the conjugate of the CLASS."""
+
+    @abstractmethod
+    def transpose(self):
+        """Return the transpose of the CLASS."""
diff --git a/qiskit/quantum_info/operators/mixins/group.py b/qiskit/quantum_info/operators/mixins/group.py
new file mode 100644
index 000000000..51bc618d1
--- /dev/null
+++ b/qiskit/quantum_info/operators/mixins/group.py
@@ -0,0 +1,164 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Mixin for gate operator interface.
+"""
+# pylint: disable=abstract-method
+
+from abc import ABC, abstractmethod
+from numbers import Integral
+
+from qiskit.exceptions import QiskitError
+
+
+class GroupMixin(ABC):
+    """Abstract Mixin for operator group operations.
+
+    This class defines the following methods
+
+        - :meth:`compose`
+        - :meth:`dot`
+        - :meth:`tensor`
+        - :meth:`expand`
+        - :meth:`power`
+
+    And the following operator overloads:
+
+        - ``*``, ``__mul__`` -> :meth:`dot`
+        - ``@``, ``__matmul__`` -> :meth:`compose`
+        - ``**``, ``__pow__`` -> :meth:`power`
+        - ``^``, ``__xor__`` -> `:meth:`tensor`
+
+    The following abstract methods must be implemented by subclasses
+    using this mixin
+
+        - ``compose(self, other, qargs=None, inplace=False)``
+        - ``tensor(self, other)``
+        - ``expand(self, other)``
+    """
+
+    def __mul__(self, other):
+        return self.dot(other)
+
+    def __matmul__(self, other):
+        return self.compose(other)
+
+    def __pow__(self, n):
+        return self.power(n)
+
+    def __xor__(self, other):
+        return self.tensor(other)
+
+    @abstractmethod
+    def tensor(self, other):
+        r"""Return the tensor product with another CLASS.
+
+        Args:
+            other (CLASS): a CLASS object.
+
+        Returns:
+            CLASS: the tensor product :math:`a \otimes b`, where :math:`a`
+                   is the current CLASS, and :math:`b` is the other CLASS.
+
+        .. note::
+            The tensor product can be obtained using the ``^`` binary operator.
+            Hence ``a.tensor(b)`` is equivalent to ``a ^ b``.
+
+        .. note:
+            Tensor uses reversed operator ordering to :meth:`expand`.
+            For two operators of the same type ``a.tensor(b) = b.expand(a)``.
+        """
+
+    @abstractmethod
+    def expand(self, other):
+        r"""Return the reverse-order tensor product with another CLASS.
+
+        Args:
+            other (CLASS): a CLASS object.
+
+        Returns:
+            CLASS: the tensor product :math:`b \otimes a`, where :math:`a`
+                   is the current CLASS, and :math:`b` is the other CLASS.
+
+        .. note:
+            Expand is the opposite operator ordering to :meth:`tensor`.
+            For two operators of the same type ``a.expand(b) = b.tensor(a)``.
+        """
+
+    @abstractmethod
+    def compose(self, other, qargs=None, front=False):
+        """Return the operator composition with another CLASS.
+
+        Args:
+            other (CLASS): a CLASS object.
+            qargs (list or None): Optional, a list of subsystem positions to
+                                  apply other on. If None apply on all
+                                  subsystems (default: None).
+            front (bool): If True compose using right operator multiplication,
+                          instead of left multiplication [default: False].
+
+        Returns:
+            CLASS: The composed CLASS.
+
+        Raises:
+            QiskitError: if other cannot be converted to an operator, or has
+                         incompatible dimensions for specified subsystems.
+
+        .. note::
+            Composition (``@``) is defined as `left` matrix multiplication for
+            matrix operators. That is that ``A @ B`` is equal to ``B * A``.
+            Setting ``front=True`` returns `right` matrix multiplication
+            ``A * B`` and is equivalent to the :meth:`dot` method.
+        """
+
+    def dot(self, other, qargs=None):
+        """Return the right multiplied operator self * other.
+
+        Args:
+            other (CLASS): an operator object.
+            qargs (list or None): Optional, a list of subsystem positions to
+                                  apply other on. If None apply on all
+                                  subsystems (default: None).
+
+        Returns:
+            CLASS: The operator self * other.
+
+        .. note::
+            The dot product can be obtained using the ``*`` binary operator.
+            Hence ``a.dot(b)`` is equivalent to ``a * b``. Left operator
+            multiplication can be obtained using the :meth:`compose` method.
+
+        """
+        return self.compose(other, qargs=qargs, front=True)
+
+    def power(self, n):
+        """Return the compose of a operator with itself n times.
+
+        Args:
+            n (int): the number of times to compose with self (n>0).
+
+        Returns:
+            CLASS: the n-times composed operator.
+
+        Raises:
+            QiskitError: if the input and output dimensions of the operator
+                         are not equal, or the power is not a positive integer.
+        """
+        # NOTE: if a subclass can have negative or non-integer powers
+        # this method should be overridden in that class.
+        if not isinstance(n, Integral) or n < 1:
+            raise QiskitError("Can only power with positive integer powers.")
+        ret = self
+        for _ in range(1, n):
+            ret = ret.dot(self)
+        return ret
diff --git a/qiskit/quantum_info/operators/mixins/linear.py b/qiskit/quantum_info/operators/mixins/linear.py
new file mode 100644
index 000000000..1671dd912
--- /dev/null
+++ b/qiskit/quantum_info/operators/mixins/linear.py
@@ -0,0 +1,62 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Mixin for linear operator interface.
+"""
+
+from abc import ABC, abstractmethod
+
+from .multiply import MultiplyMixin
+
+
+class LinearMixin(MultiplyMixin, ABC):
+    """Abstract Mixin for linear operator.
+
+    This class defines the following operator overloads:
+
+        - ``+`` / ``__add__``
+        - ``-`` / ``__sub__``
+        - ``*`` / ``__rmul__`
+        - ``/`` / ``__truediv__``
+        - ``__neg__``
+
+    The following abstract methods must be implemented by subclasses
+    using this mixin
+
+        - ``_add(self, other, qargs=None)``
+        - ``_multiply(self, other)``
+    """
+
+    def __add__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        return self._add(other, qargs=qargs)
+
+    def __sub__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        return self._add(-other, qargs=qargs)
+
+    @abstractmethod
+    def _add(self, other, qargs=None):
+        """Return the CLASS self + other.
+
+        If ``qargs`` are specified the other operator will be added
+        assuming it is identity on all other subsystems.
+
+        Args:
+            other (CLASS): an operator object.
+            qargs (None or list): optional subsystems to add on
+                                  (Default: None)
+
+        Returns:
+            CLASS: the CLASS self + other.
+        """
diff --git a/qiskit/quantum_info/operators/mixins/multiply.py b/qiskit/quantum_info/operators/mixins/multiply.py
new file mode 100644
index 000000000..aac7b43f5
--- /dev/null
+++ b/qiskit/quantum_info/operators/mixins/multiply.py
@@ -0,0 +1,52 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""
+Mixin for operator scalar multiplication interface.
+"""
+
+from abc import ABC, abstractmethod
+
+
+class MultiplyMixin(ABC):
+    """Abstract Mixin for scalar multiplication.
+
+    This class defines the following operator overloads:
+
+        - ``*`` / ``__rmul__`
+        - ``/`` / ``__truediv__``
+        - ``__neg__``
+
+    The following abstract methods must be implemented by subclasses
+    using this mixin
+
+        - ``_multiply(self, other)``
+    """
+    def __rmul__(self, other):
+        return self._multiply(other)
+
+    def __truediv__(self, other):
+        return self._multiply(1 / other)
+
+    def __neg__(self):
+        return self._multiply(-1)
+
+    @abstractmethod
+    def _multiply(self, other):
+        """Return the CLASS other * self.
+
+        Args:
+            other (complex): a complex number.
+
+        Returns:
+            CLASS: the CLASS other * self.
+        """
diff --git a/qiskit/quantum_info/operators/tolerances.py b/qiskit/quantum_info/operators/mixins/tolerances.py
similarity index 80%
rename from qiskit/quantum_info/operators/tolerances.py
rename to qiskit/quantum_info/operators/mixins/tolerances.py
index 4c1dd08c4..9325d6228 100644
--- a/qiskit/quantum_info/operators/tolerances.py
+++ b/qiskit/quantum_info/operators/mixins/tolerances.py
@@ -14,11 +14,12 @@
 Tolerances mixin class.
 """
 
+from abc import ABCMeta
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.predicates import ATOL_DEFAULT, RTOL_DEFAULT
 
 
-class TolerancesMeta(type):
+class TolerancesMeta(ABCMeta):
     """Metaclass to handle tolerances"""
     def __init__(cls, *args, **kwargs):
         cls._ATOL_DEFAULT = ATOL_DEFAULT
@@ -28,7 +29,7 @@ class TolerancesMeta(type):
 
     @property
     def atol(cls):
-        """The default absolute tolerance parameter for float comparisons."""
+        """Default absolute tolerance parameter for float comparisons."""
         return cls._ATOL_DEFAULT
 
     def _check_value(cls, value, value_name):
@@ -42,18 +43,18 @@ class TolerancesMeta(type):
 
     @atol.setter
     def atol(cls, value):
-        """Set the class default absolute tolerance parameter for float comparisons."""
+        """Set default absolute tolerance parameter for float comparisons."""
         cls._check_value(value, "atol")  # pylint: disable=no-value-for-parameter
         cls._ATOL_DEFAULT = value
 
     @property
     def rtol(cls):
-        """The relative tolerance parameter for float comparisons."""
+        """Default relative tolerance parameter for float comparisons."""
         return cls._RTOL_DEFAULT
 
     @rtol.setter
     def rtol(cls, value):
-        """Set the class default relative tolerance parameter for float comparisons."""
+        """Set default relative tolerance parameter for float comparisons."""
         cls._check_value(value, "rtol")  # pylint: disable=no-value-for-parameter
         cls._RTOL_DEFAULT = value
 
@@ -63,10 +64,10 @@ class TolerancesMixin(metaclass=TolerancesMeta):
 
     @property
     def atol(self):
-        """The default absolute tolerance parameter for float comparisons."""
+        """Default absolute tolerance parameter for float comparisons."""
         return self.__class__.atol
 
     @property
     def rtol(self):
-        """The relative tolerance parameter for float comparisons."""
+        """Default relative tolerance parameter for float comparisons."""
         return self.__class__.rtol
diff --git a/qiskit/quantum_info/operators/operator.py b/qiskit/quantum_info/operators/operator.py
index 9714081f3..5d7f18181 100644
--- a/qiskit/quantum_info/operators/operator.py
+++ b/qiskit/quantum_info/operators/operator.py
@@ -25,11 +25,11 @@ from qiskit.circuit.instruction import Instruction
 from qiskit.circuit.library.standard_gates import IGate, XGate, YGate, ZGate, HGate, SGate, TGate
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.predicates import is_unitary_matrix, matrix_equal
-from qiskit.quantum_info.operators.base_operator import BaseOperator
-from qiskit.quantum_info.operators.tolerances import TolerancesMixin
+from qiskit.quantum_info.operators.linear_op import LinearOp
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
-class Operator(BaseOperator, TolerancesMixin):
+class Operator(LinearOp):
     r"""Matrix operator class
 
     This represents a matrix operator :math:`M` that will
@@ -205,14 +205,12 @@ class Operator(BaseOperator, TolerancesMixin):
         return UnitaryGate(self.data)
 
     def conjugate(self):
-        """Return the conjugate of the operator."""
         # Make a shallow copy and update array
         ret = copy.copy(self)
         ret._data = np.conj(self._data)
         return ret
 
     def transpose(self):
-        """Return the transpose of the operator."""
         # Make a shallow copy and update array
         ret = copy.copy(self)
         ret._data = np.transpose(self._data)
@@ -220,29 +218,6 @@ class Operator(BaseOperator, TolerancesMixin):
         return ret
 
     def compose(self, other, qargs=None, front=False):
-        """Return the composed operator.
-
-        Args:
-            other (Operator): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            Operator: The operator self @ other.
-
-        Raise:
-            QiskitError: if operators have incompatible dimensions for
-                         composition.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            matrix operators. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
         if not isinstance(other, Operator):
@@ -292,83 +267,40 @@ class Operator(BaseOperator, TolerancesMixin):
         ret._op_shape = new_shape
         return ret
 
-    def dot(self, other, qargs=None):
-        """Return the right multiplied operator self * other.
-
-        Args:
-            other (Operator): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-
-        Returns:
-            Operator: The operator self * other.
-
-        Raises:
-            QiskitError: if other cannot be converted to an Operator or has
-                         incompatible dimensions.
-        """
-        return super().dot(other, qargs=qargs)
-
     def power(self, n):
         """Return the matrix power of the operator.
 
         Args:
-            n (int): the power to raise the matrix to.
+            n (float): the power to raise the matrix to.
 
         Returns:
-            BaseOperator: the n-times composed operator.
+            Operator: the resulting operator ``O ** n``.
 
         Raises:
             QiskitError: if the input and output dimensions of the operator
-                         are not equal, or the power is not a positive integer.
+                         are not equal.
         """
-        if not isinstance(n, int):
-            raise QiskitError("Can only take integer powers of Operator.")
         if self.input_dims() != self.output_dims():
             raise QiskitError("Can only power with input_dims = output_dims.")
-        # Override base class power so we can implement more efficiently
-        # using Numpy.matrix_power
         ret = copy.copy(self)
         ret._data = np.linalg.matrix_power(self.data, n)
         return ret
 
     def tensor(self, other):
-        """Return the tensor product operator self ⊗ other.
-
-        Args:
-            other (Operator): a operator subclass object.
-
-        Returns:
-            Operator: the tensor product operator self ⊗ other.
-
-        Raises:
-            QiskitError: if other cannot be converted to an operator.
-        """
         if not isinstance(other, Operator):
             other = Operator(other)
-        ret = copy.copy(self)
-        ret._data = np.kron(self._data, other._data)
-        ret._op_shape = self._op_shape.tensor(other._op_shape)
-        return ret
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product operator other ⊗ self.
-
-        Args:
-            other (Operator): an operator object.
-
-        Returns:
-            Operator: the tensor product operator other ⊗ self.
-
-        Raises:
-            QiskitError: if other cannot be converted to an operator.
-        """
         if not isinstance(other, Operator):
             other = Operator(other)
-        ret = copy.copy(self)
-        ret._data = np.kron(other._data, self._data)
-        ret._op_shape = self._op_shape.expand(other._op_shape)
+        return self._tensor(other, self)
+
+    @classmethod
+    def _tensor(cls, a, b):
+        ret = copy.copy(a)
+        ret._op_shape = a._op_shape.tensor(b._op_shape)
+        ret._data = np.kron(a.data, b.data)
         return ret
 
     def _add(self, other, qargs=None):
@@ -555,3 +487,7 @@ class Operator(BaseOperator, TolerancesMixin):
                 else:
                     new_qargs = [qargs[tup.index] for tup in qregs]
                 self._append_instruction(instr, qargs=new_qargs)
+
+
+# Update docstrings for API docs
+generate_apidocs(Operator)
diff --git a/qiskit/quantum_info/operators/scalar_op.py b/qiskit/quantum_info/operators/scalar_op.py
index 697318093..d8ee0d2a9 100644
--- a/qiskit/quantum_info/operators/scalar_op.py
+++ b/qiskit/quantum_info/operators/scalar_op.py
@@ -20,11 +20,12 @@ import numpy as np
 
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.base_operator import BaseOperator
-from qiskit.quantum_info.operators.tolerances import TolerancesMixin
+from qiskit.quantum_info.operators.linear_op import LinearOp
 from qiskit.quantum_info.operators.operator import Operator
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
-class ScalarOp(BaseOperator, TolerancesMixin):
+class ScalarOp(LinearOp):
     """Scalar identity operator class.
 
     This is a symbolic representation of an scalar identity operator on
@@ -64,13 +65,11 @@ class ScalarOp(BaseOperator, TolerancesMixin):
         return self._coeff
 
     def conjugate(self):
-        """Return the conjugate of the operator."""
         ret = self.copy()
         ret._coeff = np.conjugate(self.coeff)
         return ret
 
     def transpose(self):
-        """Return the transpose of the operator."""
         return self.copy()
 
     def is_unitary(self, atol=None, rtol=None):
@@ -93,71 +92,9 @@ class ScalarOp(BaseOperator, TolerancesMixin):
                         input_dims=self.input_dims(),
                         output_dims=self.output_dims())
 
-    def tensor(self, other):
-        """Return the tensor product operator self ⊗ other.
-
-        Args:
-            other (BaseOperator): an operator object.
-
-        Returns:
-            ScalarOp: if other is an ScalarOp.
-            BaseOperator: if other is not an ScalarOp.
-        """
-        if not isinstance(other, BaseOperator):
-            other = Operator(other)
-        if isinstance(other, ScalarOp):
-            ret = copy.copy(self)
-            ret._coeff = self.coeff * other.coeff
-            ret._op_shape = self._op_shape.tensor(other._op_shape)
-            return ret
-        return other.expand(self)
-
-    def expand(self, other):
-        """Return the tensor product operator other ⊗ self.
-
-        Args:
-            other (BaseOperator): an operator object.
-
-        Returns:
-            ScalarOp: if other is an ScalarOp.
-            BaseOperator: if other is not an ScalarOp.
-        """
-        if not isinstance(other, BaseOperator):
-            other = Operator(other)
-        if isinstance(other, ScalarOp):
-            ret = copy.copy(self)
-            ret._coeff = self.coeff * other.coeff
-            ret._op_shape = self._op_shape.expand(other._op_shape)
-            return ret
-        return other.tensor(self)
-
     def compose(self, other, qargs=None, front=False):
-        """Return the composed operator.
-
-        Args:
-            other (BaseOperator): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
-
-        Returns:
-            BaseOperator: The operator self @ other.
-
-        Raises:
-            QiskitError: if other has incompatible dimensions for specified
-                         subsystems.
-
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            matrix operators. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
-
         if not isinstance(other, BaseOperator):
             other = Operator(other)
 
@@ -197,19 +134,39 @@ class ScalarOp(BaseOperator, TolerancesMixin):
         """Return the power of the ScalarOp.
 
         Args:
-            n (Number): the exponent for the scalar op.
+            n (float): the exponent for the scalar op.
 
         Returns:
             ScalarOp: the ``coeff ** n`` ScalarOp.
-
-        Raises:
-            QiskitError: if the input and output dimensions of the operator
-                         are not equal, or the power is not a positive integer.
         """
         ret = self.copy()
         ret._coeff = self.coeff ** n
         return ret
 
+    def tensor(self, other):
+        if not isinstance(other, BaseOperator):
+            other = Operator(other)
+
+        if isinstance(other, ScalarOp):
+            ret = copy.copy(self)
+            ret._coeff = self.coeff * other.coeff
+            ret._op_shape = self._op_shape.tensor(other._op_shape)
+            return ret
+
+        return other.expand(self)
+
+    def expand(self, other):
+        if not isinstance(other, BaseOperator):
+            other = Operator(other)
+
+        if isinstance(other, ScalarOp):
+            ret = copy.copy(self)
+            ret._coeff = self.coeff * other.coeff
+            ret._op_shape = self._op_shape.expand(other._op_shape)
+            return ret
+
+        return other.tensor(self)
+
     def _add(self, other, qargs=None):
         """Return the operator self + other.
 
@@ -290,3 +247,7 @@ class ScalarOp(BaseOperator, TolerancesMixin):
         if qargs is None:
             return other
         return ScalarOp(current.input_dims()).compose(other, qargs=qargs)
+
+
+# Update docstrings for API docs
+generate_apidocs(ScalarOp)
diff --git a/qiskit/quantum_info/operators/symplectic/base_pauli.py b/qiskit/quantum_info/operators/symplectic/base_pauli.py
index e7279fd06..eb81492af 100644
--- a/qiskit/quantum_info/operators/symplectic/base_pauli.py
+++ b/qiskit/quantum_info/operators/symplectic/base_pauli.py
@@ -21,9 +21,10 @@ from qiskit.exceptions import QiskitError
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.barrier import Barrier
 from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.quantum_info.operators.mixins import AdjointMixin, MultiplyMixin
 
 
-class BasePauli(BaseOperator):
+class BasePauli(BaseOperator, AdjointMixin, MultiplyMixin):
     r"""Symplectic representation of a list of N-qubit Paulis.
 
     Base class for Pauli and PauliList.
@@ -62,51 +63,35 @@ class BasePauli(BaseOperator):
     # ---------------------------------------------------------------------
 
     def tensor(self, other):
-        """Return the tensor product Pauli self ⊗ other.
-
-        Args:
-            other (BasePauli): another Pauli.
-
-        Returns:
-            BasePauli: the tensor product Pauli.
-        """
-        z = np.hstack([self._stack(other._z, self._num_paulis), self._z])
-        x = np.hstack([self._stack(other._x, self._num_paulis), self._x])
-        phase = np.mod(self._phase + other._phase, 4)
-        return BasePauli(z, x, phase)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product Pauli other ⊗ self.
+        return self._tensor(other, self)
 
-        Args:
-            other (BasePauli): another Pauli.
-
-        Returns:
-            BasePauli: the tensor product Pauli.
-        """
-        z = np.hstack([self._z, self._stack(other._z, self._num_paulis)])
-        x = np.hstack([self._x, self._stack(other._x, self._num_paulis)])
-        phase = np.mod(self._phase + other._phase, 4)
+    @classmethod
+    def _tensor(cls, a, b):
+        z = np.hstack([a._stack(b._z, a._num_paulis), a._z])
+        x = np.hstack([a._stack(b._x, a._num_paulis), a._x])
+        phase = np.mod(a._phase + b._phase, 4)
         return BasePauli(z, x, phase)
 
     # pylint: disable=arguments-differ
     def compose(self, other, qargs=None, front=False, inplace=False):
-        """Return the composition of Paulis self∘other.
-
-        Note that this discards phases.
+        """Return the composition of Paulis.
 
         Args:
-            other (BasePauli): another Pauli.
-            qargs (None or list): qubits to apply dot product on (default: None).
-            front (bool): If True use `dot` composition method (default: False).
+            a ({cls}): an operator object.
+            b ({cls}): an operator object.
+            qargs (list or None): Optional, qubits to apply dot product
+                                  on (default: None).
             inplace (bool): If True update in-place (default: False).
 
         Returns:
-            BasePauli: the output Pauli.
+            {cls}: The operator a.compose(b)
 
         Raises:
             QiskitError: if number of qubits of other does not match qargs.
-        """
+        """.format(cls=type(self).__name__)
         # Validation
         if qargs is None and other.num_qubits != self.num_qubits:
             raise QiskitError(
@@ -156,37 +141,18 @@ class BasePauli(BaseOperator):
         ret._phase = np.mod(phase, 4)
         return ret
 
-    # pylint: disable=arguments-differ
-    def dot(self, other, qargs=None, inplace=False):
-        """Return the dot product of Paulis self∘other.
-
-        Note that this discards phases.
-
-        Args:
-            other (BasePauli): another Pauli.
-            qargs (None or list): qubits to apply dot product on (default: None).
-            inplace (bool): If True update in-place (default: False).
-
-        Returns:
-            BasePauli: the output Pauli.
-
-        Raises:
-            QiskitError: if number of qubits of other does not match qargs.
-        """
-        return self.compose(other, qargs=qargs, front=True, inplace=inplace)
-
     def _multiply(self, other):
-        """Multiply each Pauli in the table by a phase.
+        """Return the {cls} other * self.
 
         Args:
-            other (complex): a complex number in [1, -1j, -1, 1j]
+            other (complex): a complex number in ``[1, -1j, -1, 1j]``.
 
         Returns:
-            BasePauli: the Pauli table other * self.
+            {cls}: the {cls} other * self.
 
         Raises:
-            QiskitError: if the phase is not in the set [1, -1j, -1, 1j].
-        """
+            QiskitError: if the phase is not in the set ``[1, -1j, -1, 1j]``.
+        """.format(cls=type(self).__name__)
         if isinstance(other, (np.ndarray, list, tuple)):
             phase = np.array([self._phase_from_complex(phase) for phase in other])
         else:
@@ -267,7 +233,7 @@ class BasePauli(BaseOperator):
         # Evolve via Pauli
         if isinstance(other, BasePauli):
             ret = self.compose(other.adjoint(), qargs=qargs)
-            ret = ret.dot(other, qargs=qargs)
+            ret = ret.compose(other, front=True, qargs=qargs)
             return ret
 
         # Otherwise evolve by the inverse circuit to compute C^dg.P.C
@@ -278,7 +244,7 @@ class BasePauli(BaseOperator):
     # ---------------------------------------------------------------------
 
     def __imul__(self, other):
-        return self.dot(other, inplace=True)
+        return self.compose(other, front=True, inplace=True)
 
     def __neg__(self):
         ret = copy.copy(self)
diff --git a/qiskit/quantum_info/operators/symplectic/clifford.py b/qiskit/quantum_info/operators/symplectic/clifford.py
index a445f6ecc..6c42672ba 100644
--- a/qiskit/quantum_info/operators/symplectic/clifford.py
+++ b/qiskit/quantum_info/operators/symplectic/clifford.py
@@ -23,11 +23,12 @@ from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.quantum_info.operators.operator import Operator
 from qiskit.quantum_info.operators.scalar_op import ScalarOp
 from qiskit.quantum_info.synthesis.clifford_decompose import decompose_clifford
+from qiskit.quantum_info.operators.mixins import generate_apidocs, AdjointMixin
 from .stabilizer_table import StabilizerTable
 from .clifford_circuits import _append_circuit
 
 
-class Clifford(BaseOperator):
+class Clifford(BaseOperator, AdjointMixin):
     """An N-qubit unitary operator from the Clifford group.
 
     **Representation**
@@ -214,44 +215,38 @@ class Clifford(BaseOperator):
     # ---------------------------------------------------------------------
 
     def conjugate(self):
-        """Return the conjugate of the Clifford."""
         return Clifford._conjugate_transpose(self, 'C')
 
     def adjoint(self):
-        """Return the conjugate transpose of the Clifford"""
         return Clifford._conjugate_transpose(self, 'A')
 
     def transpose(self):
-        """Return the transpose of the Clifford."""
         return Clifford._conjugate_transpose(self, 'T')
 
-    def compose(self, other, qargs=None, front=False):
-        """Return the composed operator.
+    def tensor(self, other):
+        if not isinstance(other, Clifford):
+            other = Clifford(other)
+        return self._tensor(self, other)
 
-        Args:
-            other (Clifford): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
-            front (bool): If True compose using right operator multiplication,
-                          instead of left multiplication [default: False].
+    def expand(self, other):
+        if not isinstance(other, Clifford):
+            other = Clifford(other)
+        return self._tensor(other, self)
 
-        Returns:
-            Clifford: The operator self @ other.
+    @classmethod
+    def _tensor(cls, a, b):
+        # Pad stabilizers and destabilizers
+        destab = (b.destabilizer.expand(a.num_qubits * 'I') +
+                  a.destabilizer.tensor(b.num_qubits * 'I'))
+        stab = (b.stabilizer.expand(a.num_qubits * 'I') +
+                a.stabilizer.tensor(b.num_qubits * 'I'))
 
-        Raise:
-            QiskitError: if operators have incompatible dimensions for
-                         composition.
+        # Add the padded table
+        return Clifford(destab + stab, validate=False)
 
-        Additional Information:
-            Composition (``@``) is defined as `left` matrix multiplication for
-            matrix operators. That is that ``A @ B`` is equal to ``B * A``.
-            Setting ``front=True`` returns `right` matrix multiplication
-            ``A * B`` and is equivalent to the :meth:`dot` method.
-        """
+    def compose(self, other, qargs=None, front=False):
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
-
         # If other is a QuantumCircuit we can more efficiently compose
         # using the _append_circuit method to update each gate recursively
         # to the current Clifford, rather than converting to a Clifford first
@@ -269,47 +264,65 @@ class Clifford(BaseOperator):
 
         # Pad other with identities if composeing on subsystem
         other = self._pad_with_identity(other, qargs)
-        return self._compose_clifford(other, front=front)
 
-    def dot(self, other, qargs=None):
-        """Return the right multiplied operator self * other.
+        if front:
+            table1 = self.table
+            table2 = other.table
+        else:
+            table1 = other.table
+            table2 = self.table
 
-        Args:
-            other (Clifford): an operator object.
-            qargs (list or None): a list of subsystem positions to apply
-                                  other on. If None apply on all
-                                  subsystems [default: None].
+        num_qubits = self.num_qubits
 
-        Returns:
-            Clifford: The operator self * other.
+        array1 = table1.array.astype(int)
+        phase1 = table1.phase.astype(int)
 
-        Raises:
-            QiskitError: if operators have incompatible dimensions for
-                         composition.
-        """
-        return super().dot(other, qargs=qargs)
+        array2 = table2.array.astype(int)
+        phase2 = table2.phase.astype(int)
 
-    def tensor(self, other):
-        """Return the tensor product operator self ⊗ other.
+        # Update Pauli table
+        pauli = StabilizerTable(array2.dot(array1) % 2)
 
-        Args:
-            other (Clifford): a operator subclass object.
+        # Add phases
+        phase = np.mod(array2.dot(phase1) + phase2, 2)
 
-        Returns:
-            Clifford: the tensor product operator self ⊗ other.
-        """
-        return self._tensor_product(other, reverse=False)
+        # Correcting for phase due to Pauli multiplication
+        ifacts = np.zeros(2 * num_qubits, dtype=np.int)
 
-    def expand(self, other):
-        """Return the tensor product operator other ⊗ self.
+        for k in range(2 * num_qubits):
 
-        Args:
-            other (Clifford): an operator object.
+            row2 = array2[k]
+            x2 = table2.X[k]
+            z2 = table2.Z[k]
 
-        Returns:
-            Clifford: the tensor product operator other ⊗ self.
-        """
-        return self._tensor_product(other, reverse=True)
+            # Adding a factor of i for each Y in the image of an operator under the
+            # first operation, since Y=iXZ
+
+            ifacts[k] += np.sum(x2 & z2)
+
+            # Adding factors of i due to qubit-wise Pauli multiplication
+
+            for j in range(num_qubits):
+                x = 0
+                z = 0
+                for i in range(2 * num_qubits):
+                    if row2[i]:
+                        x1 = array1[i, j]
+                        z1 = array1[i, j + num_qubits]
+                        if (x | z) & (x1 | z1):
+                            val = np.mod(np.abs(3 * z1 - x1) - np.abs(3 * z - x) - 1, 3)
+                            if val == 0:
+                                ifacts[k] += 1
+                            elif val == 1:
+                                ifacts[k] -= 1
+                        x = np.mod(x + x1, 2)
+                        z = np.mod(z + z1, 2)
+
+        p = np.mod(ifacts, 4) // 2
+
+        phase = np.mod(phase + p, 2)
+
+        return Clifford(StabilizerTable(pauli, phase), validate=False)
 
     # ---------------------------------------------------------------------
     # Representation conversions
@@ -490,38 +503,6 @@ class Clifford(BaseOperator):
                 ret.table.X & ret.table.Z, axis=1), 2).astype(bool)
         return ret
 
-    def _tensor_product(self, other, reverse=False):
-        """Return the tensor product operator.
-
-        Args:
-            other (Clifford): another Clifford operator.
-            reverse (bool): If False return self ⊗ other, if True return
-                            if True return (other ⊗ self) [Default: False].
-        Returns:
-            Clifford: the tensor product operator.
-
-        Raises:
-            QiskitError: if other cannot be converted into an Clifford.
-        """
-        if not isinstance(other, Clifford):
-            other = Clifford(other)
-
-        if reverse:
-            cliff0 = self
-            cliff1 = other
-        else:
-            cliff0 = other
-            cliff1 = self
-
-        # Pad stabilizers and destabilizers
-        destab = (cliff0.destabilizer.expand(cliff1.num_qubits * 'I') +
-                  cliff1.destabilizer.tensor(cliff0.num_qubits * 'I'))
-        stab = (cliff0.stabilizer.expand(cliff1.num_qubits * 'I') +
-                cliff1.stabilizer.tensor(cliff0.num_qubits * 'I'))
-
-        # Add the padded table
-        return Clifford(destab + stab, validate=False)
-
     def _pad_with_identity(self, clifford, qargs):
         """Pad Clifford with identities on other subsystems."""
         if qargs is None:
@@ -543,63 +524,6 @@ class Clifford(BaseOperator):
 
         return padded
 
-    def _compose_clifford(self, other, front=False):
-        """Return the composition channel assume other is Clifford of same size as self."""
-        if front:
-            table1 = self.table
-            table2 = other.table
-        else:
-            table1 = other.table
-            table2 = self.table
-
-        num_qubits = self.num_qubits
-
-        array1 = table1.array.astype(int)
-        phase1 = table1.phase.astype(int)
-
-        array2 = table2.array.astype(int)
-        phase2 = table2.phase.astype(int)
-
-        # Update Pauli table
-        pauli = StabilizerTable(array2.dot(array1) % 2)
-
-        # Add phases
-        phase = np.mod(array2.dot(phase1) + phase2, 2)
-
-        # Correcting for phase due to Pauli multiplication
-        ifacts = np.zeros(2 * num_qubits, dtype=int)
-
-        for k in range(2 * num_qubits):
-
-            row2 = array2[k]
-            x2 = table2.X[k]
-            z2 = table2.Z[k]
 
-            # Adding a factor of i for each Y in the image of an operator under the
-            # first operation, since Y=iXZ
-
-            ifacts[k] += np.sum(x2 & z2)
-
-            # Adding factors of i due to qubit-wise Pauli multiplication
-
-            for j in range(num_qubits):
-                x = 0
-                z = 0
-                for i in range(2 * num_qubits):
-                    if row2[i]:
-                        x1 = array1[i, j]
-                        z1 = array1[i, j + num_qubits]
-                        if (x | z) & (x1 | z1):
-                            val = np.mod(np.abs(3 * z1 - x1) - np.abs(3 * z - x) - 1, 3)
-                            if val == 0:
-                                ifacts[k] += 1
-                            elif val == 1:
-                                ifacts[k] -= 1
-                        x = np.mod(x + x1, 2)
-                        z = np.mod(z + z1, 2)
-
-        p = np.mod(ifacts, 4) // 2
-
-        phase = np.mod(phase + p, 2)
-
-        return Clifford(StabilizerTable(pauli, phase), validate=False)
+# Update docstrings for API docs
+generate_apidocs(Clifford)
diff --git a/qiskit/quantum_info/operators/symplectic/pauli.py b/qiskit/quantum_info/operators/symplectic/pauli.py
index 058db394b..739a26381 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli.py
@@ -26,6 +26,7 @@ from qiskit.circuit import QuantumCircuit, Instruction
 from qiskit.circuit.library.standard_gates import IGate, XGate, YGate, ZGate
 from qiskit.circuit.library.generalized_gates import PauliGate
 from qiskit.circuit.barrier import Barrier
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
 class Pauli(BasePauli):
@@ -426,46 +427,30 @@ class Pauli(BasePauli):
     # BaseOperator methods
     # ---------------------------------------------------------------------
 
-    def tensor(self, other):
-        """Return the tensor product Pauli self ⊗ other.
-
-        Args:
-            other (Pauli): another Pauli.
-
-        Returns:
-            Pauli: the tensor product Pauli.
-        """
-        if not isinstance(other, Pauli):
-            other = Pauli(other)
-        return Pauli(super().tensor(other))
-
-    def expand(self, other):
-        """Return the tensor product Pauli other ⊗ self.
-
-        Args:
-            other (Pauli): another Pauli.
-
-        Returns:
-            Pauli: the tensor product Pauli.
-        """
-        if not isinstance(other, Pauli):
-            other = Pauli(other)
-        return Pauli(super().expand(other))
-
+    # pylint: disable=arguments-differ
     def compose(self, other, qargs=None, front=False, inplace=False):
-        """Return the composed Pauli self∘other.
+        """Return the operator composition with another Pauli.
 
         Args:
-            other (Pauli): another Pauli.
-            qargs (None or list): qubits to apply dot product on (default: None).
-            front (bool): If True use `dot` composition method (default: False).
+            other (Pauli): a Pauli object.
+            qargs (list or None): Optional, qubits to apply dot product
+                                  on (default: None).
+            front (bool): If True compose using right operator multiplication,
+                          instead of left multiplication [default: False].
             inplace (bool): If True update in-place (default: False).
 
         Returns:
-            Pauli: the output Pauli.
+            Pauli: The composed Pauli.
 
         Raises:
-            QiskitError: if other cannot be converted to a Pauli.
+            QiskitError: if other cannot be converted to an operator, or has
+                         incompatible dimensions for specified subsystems.
+
+        .. note::
+            Composition (``@``) is defined as `left` matrix multiplication for
+            matrix operators. That is that ``A @ B`` is equal to ``B * A``.
+            Setting ``front=True`` returns `right` matrix multiplication
+            ``A * B`` and is equivalent to the :meth:`dot` method.
         """
         # pylint: disable=unused-argument
         if qargs is None:
@@ -477,50 +462,51 @@ class Pauli(BasePauli):
                                      front=front,
                                      inplace=inplace))
 
+    # pylint: disable=arguments-differ
     def dot(self, other, qargs=None, inplace=False):
-        """Return the dot product Pauli self∘other.
+        """Return the right multiplied operator self * other.
 
         Args:
-            other (Pauli): another Pauli.
-            qargs (None or list): qubits to apply dot product on (default: None).
+            other (Pauli): an operator object.
+            qargs (list or None): Optional, qubits to apply dot product
+                                  on (default: None).
             inplace (bool): If True update in-place (default: False).
 
         Returns:
-            Pauli: the dot outer product table.
+            Pauli: The operator self * other.
 
-        Raises:
-            QiskitError: if other cannot be converted to a Pauli.
-        """
-        return Pauli(super().dot(other, qargs=qargs, inplace=inplace))
+        .. note::
+            The dot product can be obtained using the ``*`` binary operator.
+            Hence ``a.dot(b)`` is equivalent to ``a * b``. Left operator
+            multiplication can be obtained using the :meth:`compose` method.
 
-    def _multiply(self, other):
-        """Multiply Pauli by a phase.
+        """
+        return self.compose(other, qargs=qargs, front=True, inplace=inplace)
 
-        Args:
-            other (complex): a complex number in [1, -1j, -1, 1j]
+    def tensor(self, other):
+        if not isinstance(other, Pauli):
+            other = Pauli(other)
+        return Pauli(super().tensor(other))
 
-        Returns:
-            Pauli: the Pauli other * self.
+    def expand(self, other):
+        if not isinstance(other, Pauli):
+            other = Pauli(other)
+        return Pauli(super().expand(other))
 
-        Raises:
-            QiskitError: if the phase is not in the set [1, -1j, -1, 1j].
-        """
+    def _multiply(self, other):
         return Pauli(super()._multiply(other))
 
     def conjugate(self):
-        """Return the conjugated Pauli."""
         return Pauli(super().conjugate())
 
     def transpose(self):
-        """Return the transposed Pauli."""
         return Pauli(super().transpose())
 
     def adjoint(self):
-        """Return the adjoint Pauli."""
         return Pauli(super().adjoint())
 
     def inverse(self):
-        """Return the inverse Pauli."""
+        """Return the inverse of the Pauli."""
         return Pauli(super().adjoint())
 
     # ---------------------------------------------------------------------
@@ -1060,3 +1046,7 @@ def _phase_from_label(label):
     if label not in phases:
         raise QiskitError("Invalid Pauli phase label '{}'".format(label))
     return phases.get(label)
+
+
+# Update docstrings for API docs
+generate_apidocs(Pauli)
diff --git a/qiskit/quantum_info/operators/symplectic/pauli_table.py b/qiskit/quantum_info/operators/symplectic/pauli_table.py
index 73b620dc4..3c37588bd 100644
--- a/qiskit/quantum_info/operators/symplectic/pauli_table.py
+++ b/qiskit/quantum_info/operators/symplectic/pauli_table.py
@@ -12,7 +12,7 @@
 """
 Symplectic Pauli Table Class
 """
-# pylint: disable=invalid-name
+# pylint: disable=invalid-name, useless-super-delegation
 
 import numpy as np
 
@@ -21,9 +21,10 @@ from qiskit.quantum_info.operators.base_operator import BaseOperator
 from qiskit.quantum_info.operators.scalar_op import ScalarOp
 from qiskit.quantum_info.operators.symplectic.pauli import Pauli
 from qiskit.quantum_info.operators.custom_iterator import CustomIterator
+from qiskit.quantum_info.operators.mixins import generate_apidocs, AdjointMixin
 
 
-class PauliTable(BaseOperator):
+class PauliTable(BaseOperator, AdjointMixin):
     r"""Symplectic representation of a list Pauli matrices.
 
     **Symplectic Representation**
@@ -509,9 +510,7 @@ class PauliTable(BaseOperator):
         """
         if not isinstance(other, PauliTable):
             other = PauliTable(other)
-        x1, x2 = self._block_stack(self.X, other.X)
-        z1, z2 = self._block_stack(self.Z, other.Z)
-        return PauliTable(np.hstack([x2, x1, z2, z1]))
+        return self._tensor(self, other)
 
     def expand(self, other):
         """Return the expand output product of two tables.
@@ -542,9 +541,7 @@ class PauliTable(BaseOperator):
         """
         if not isinstance(other, PauliTable):
             other = PauliTable(other)
-        x1, x2 = self._block_stack(self.X, other.X)
-        z1, z2 = self._block_stack(self.Z, other.Z)
-        return PauliTable(np.hstack([x1, x2, z1, z2]))
+        return self._tensor(other, self)
 
     def compose(self, other, qargs=None, front=True):
         """Return the compose output product of two tables.
@@ -631,6 +628,12 @@ class PauliTable(BaseOperator):
         """
         return self.compose(other, qargs=qargs, front=True)
 
+    @classmethod
+    def _tensor(cls, a, b):
+        x1, x2 = a._block_stack(a.X, b.X)
+        z1, z2 = a._block_stack(a.Z, b.Z)
+        return PauliTable(np.hstack([x2, x1, z2, z1]))
+
     def _add(self, other, qargs=None):
         """Append with another PauliTable.
 
@@ -663,6 +666,10 @@ class PauliTable(BaseOperator):
         padded = padded.compose(other, qargs=qargs)
         return PauliTable(np.vstack((self._array, padded._array)))
 
+    def __add__(self, other):
+        qargs = getattr(other, 'qargs', None)
+        return self._add(other, qargs=qargs)
+
     def conjugate(self):
         """Not implemented."""
         raise NotImplementedError(
@@ -1058,3 +1065,7 @@ class PauliTable(BaseOperator):
             def __getitem__(self, key):
                 return self.obj._to_matrix(self.obj.array[key], sparse=sparse)
         return MatrixIterator(self)
+
+
+# Update docstrings for API docs
+generate_apidocs(PauliTable)
diff --git a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
index 2eb2574dc..110d70115 100644
--- a/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
+++ b/qiskit/quantum_info/operators/symplectic/sparse_pauli_op.py
@@ -17,15 +17,15 @@ from numbers import Number
 import numpy as np
 
 from qiskit.exceptions import QiskitError
-from qiskit.quantum_info.operators.base_operator import BaseOperator
-from qiskit.quantum_info.operators.tolerances import TolerancesMixin
+from qiskit.quantum_info.operators.linear_op import LinearOp
 from qiskit.quantum_info.operators.operator import Operator
 from qiskit.quantum_info.operators.symplectic.pauli_table import PauliTable
 from qiskit.quantum_info.operators.symplectic.pauli_utils import pauli_basis
 from qiskit.quantum_info.operators.custom_iterator import CustomIterator
+from qiskit.quantum_info.operators.mixins import generate_apidocs
 
 
-class SparsePauliOp(BaseOperator, TolerancesMixin):
+class SparsePauliOp(LinearOp):
     """Sparse N-qubit operator in a Pauli basis representation.
 
     This is a sparse representation of an N-qubit matrix
@@ -68,7 +68,7 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
             raise QiskitError("coeff vector is incorrect shape for number"
                               " of Paulis {} != {}".format(self._coeffs.shape,
                                                            self._table.size))
-        # Initialize BaseOperator
+        # Initialize LinearOp
         super().__init__(num_qubits=self._table.num_qubits)
 
     def __array__(self, dtype=None):
@@ -141,35 +141,10 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         self.coeffs[key] = value.coeffs
 
     # ---------------------------------------------------------------------
-    # BaseOperator Methods
+    # LinearOp Methods
     # ---------------------------------------------------------------------
 
-    def is_unitary(self, atol=None, rtol=None):
-        """Return True if operator is a unitary matrix.
-
-        Args:
-            atol (float): Optional. Absolute tolerance for checking if
-                          coefficients are zero (Default: 1e-8).
-            rtol (float): Optinoal. relative tolerance for checking if
-                          coefficients are zero (Default: 1e-5).
-
-        Returns:
-            bool: True if the operator is unitary, False otherwise.
-        """
-        # Get default atol and rtol
-        if atol is None:
-            atol = self.atol
-        if rtol is None:
-            rtol = self.rtol
-
-        # Compose with adjoint
-        val = self.compose(self.adjoint()).simplify()
-        # See if the result is an identity
-        return (val.size == 1 and np.isclose(val.coeffs[0], 1.0, atol=atol, rtol=rtol)
-                and not np.any(val.table.X) and not np.any(val.table.Z))
-
     def conjugate(self):
-        """Return the conjugate of the operator."""
         # Conjugation conjugates phases and also Y.conj() = -Y
         # Hence we need to multiply conjugated coeffs by -1
         # for rows with an odd number of Y terms.
@@ -179,7 +154,6 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         return ret
 
     def transpose(self):
-        """Return the transpose of the operator."""
         # The only effect transposition has is Y.T = -Y
         # Hence we need to multiply coeffs by -1 for rows with an
         # odd number of Y terms.
@@ -189,7 +163,6 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         return ret
 
     def adjoint(self):
-        """Return the adjoint of the operator."""
         # Pauli's are self adjoint, so we only need to
         # conjugate the phases
         ret = self.copy()
@@ -197,22 +170,6 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         return ret
 
     def compose(self, other, qargs=None, front=False):
-        """Return the composition channel self∘other.
-
-        Args:
-            other (SparsePauliOp): an operator object.
-            qargs (list or None): a list of subsystem positions to compose other on.
-            front (bool or None): If False compose in standard order other(self(input))
-                          otherwise compose in reverse order self(other(input))
-                          [default: False]
-
-        Returns:
-            SparsePauliOp: The composed operator.
-
-        Raises:
-            QiskitError: if other cannot be converted to an Operator or has
-            incompatible dimensions.
-        """
         # pylint: disable=invalid-name
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
@@ -253,78 +210,23 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         coeffs *= (-1j) ** np.array(np.sum(minus_i, axis=1), dtype=int)
         return SparsePauliOp(table, coeffs)
 
-    def dot(self, other, qargs=None):
-        """Return the composition channel self∘other.
-
-        Args:
-            other (SparsePauliOp): an operator object.
-            qargs (list or None): a list of subsystem positions to compose other on.
-
-        Returns:
-            SparsePauliOp: The composed operator.
-
-        Raises:
-            QiskitError: if other cannot be converted to an Operator or has
-            incompatible dimensions.
-        """
-        return self.compose(other, qargs=qargs, front=True)
-
     def tensor(self, other):
-        """Return the tensor product operator self ⊗ other.
-
-        Args:
-            other (SparsePauliOp): a operator subclass object.
-
-        Returns:
-            SparsePauliOp: the tensor product operator self ⊗ other.
-
-        Raises:
-            QiskitError: if other cannot be converted to a SparsePauliOp
-                         operator.
-        """
         if not isinstance(other, SparsePauliOp):
             other = SparsePauliOp(other)
-        table = self.table.tensor(other.table)
-        coeffs = np.kron(self.coeffs, other.coeffs)
-        return SparsePauliOp(table, coeffs)
+        return self._tensor(self, other)
 
     def expand(self, other):
-        """Return the tensor product operator other ⊗ self.
-
-        Args:
-            other (SparsePauliOp): an operator object.
-
-        Returns:
-            SparsePauliOp: the tensor product operator other ⊗ self.
-
-        Raises:
-            QiskitError: if other cannot be converted to a SparsePauliOp
-                         operator.
-        """
         if not isinstance(other, SparsePauliOp):
             other = SparsePauliOp(other)
-        table = self.table.expand(other.table)
-        coeffs = np.kron(self.coeffs, other.coeffs)
+        return self._tensor(other, self)
+
+    @classmethod
+    def _tensor(cls, a, b):
+        table = a.table.tensor(b.table)
+        coeffs = np.kron(a.coeffs, b.coeffs)
         return SparsePauliOp(table, coeffs)
 
     def _add(self, other, qargs=None):
-        """Return the operator self + other.
-
-        If ``qargs`` are specified the other operator will be added
-        assuming it is identity on all other subsystems.
-
-        Args:
-            other (SparsePauliOp): an operator object.
-            qargs (None or list): optional subsystems to add on
-                                  (Default: None)
-
-        Returns:
-            SparsePauliOp: the operator self + other.
-
-        Raises:
-            QiskitError: if other cannot be converted to a SparsePauliOp
-                         or has incompatible dimensions.
-        """
         if qargs is None:
             qargs = getattr(other, 'qargs', None)
 
@@ -339,17 +241,6 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         return ret
 
     def _multiply(self, other):
-        """Return the operator other * self.
-
-        Args:
-            other (complex): a complex number.
-
-        Returns:
-            SparsePauliOp: the operator other * self.
-
-        Raises:
-            QiskitError: if other is not a valid complex number.
-        """
         if not isinstance(other, Number):
             raise QiskitError("other is not a number")
         if other == 0:
@@ -365,13 +256,37 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
     # Utility Methods
     # ---------------------------------------------------------------------
 
+    def is_unitary(self, atol=None, rtol=None):
+        """Return True if operator is a unitary matrix.
+
+        Args:
+            atol (float): Optional. Absolute tolerance for checking if
+                          coefficients are zero (Default: 1e-8).
+            rtol (float): Optional. relative tolerance for checking if
+                          coefficients are zero (Default: 1e-5).
+
+        Returns:
+            bool: True if the operator is unitary, False otherwise.
+        """
+        # Get default atol and rtol
+        if atol is None:
+            atol = self.atol
+        if rtol is None:
+            rtol = self.rtol
+
+        # Compose with adjoint
+        val = self.compose(self.adjoint()).simplify()
+        # See if the result is an identity
+        return (val.size == 1 and np.isclose(val.coeffs[0], 1.0, atol=atol, rtol=rtol)
+                and not np.any(val.table.X) and not np.any(val.table.Z))
+
     def simplify(self, atol=None, rtol=None):
         """Simplify PauliTable by combining duplicaties and removing zeros.
 
         Args:
             atol (float): Optional. Absolute tolerance for checking if
                           coefficients are zero (Default: 1e-8).
-            rtol (float): Optinoal. relative tolerance for checking if
+            rtol (float): Optional. relative tolerance for checking if
                           coefficients are zero (Default: 1e-5).
 
         Returns:
@@ -416,7 +331,7 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
             obj (Operator): an N-qubit operator.
             atol (float): Optional. Absolute tolerance for checking if
                           coefficients are zero (Default: 1e-8).
-            rtol (float): Optinoal. relative tolerance for checking if
+            rtol (float): Optional. relative tolerance for checking if
                           coefficients are zero (Default: 1e-5).
 
         Returns:
@@ -427,9 +342,9 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
         """
         # Get default atol and rtol
         if atol is None:
-            atol = SparsePauliOp._ATOL_DEFAULT
+            atol = SparsePauliOp.atol
         if rtol is None:
-            rtol = SparsePauliOp._RTOL_DEFAULT
+            rtol = SparsePauliOp.rtol
 
         if not isinstance(obj, Operator):
             obj = Operator(obj)
@@ -573,3 +488,7 @@ class SparsePauliOp(BaseOperator, TolerancesMixin):
                 return coeff * mat
 
         return MatrixIterator(self)
+
+
+# Update docstrings for API docs
+generate_apidocs(SparsePauliOp)
diff --git a/qiskit/quantum_info/operators/symplectic/stabilizer_table.py b/qiskit/quantum_info/operators/symplectic/stabilizer_table.py
index 73655fd6c..dd1ec3461 100644
--- a/qiskit/quantum_info/operators/symplectic/stabilizer_table.py
+++ b/qiskit/quantum_info/operators/symplectic/stabilizer_table.py
@@ -12,16 +12,17 @@
 """
 Symplectic Stabilizer Table Class
 """
-# pylint: disable=invalid-name, abstract-method
+# pylint: disable=invalid-name, abstract-method, useless-super-delegation
 
 import numpy as np
 
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.custom_iterator import CustomIterator
 from qiskit.quantum_info.operators.symplectic.pauli_table import PauliTable
+from qiskit.quantum_info.operators.mixins import generate_apidocs, AdjointMixin
 
 
-class StabilizerTable(PauliTable):
+class StabilizerTable(PauliTable, AdjointMixin):
     r"""Symplectic representation of a list Stabilizer matrices.
 
     **Symplectic Representation**
@@ -488,11 +489,7 @@ class StabilizerTable(PauliTable):
         """
         if not isinstance(other, StabilizerTable):
             other = StabilizerTable(other)
-        pauli = super().tensor(other)
-        phase1, phase2 = self._block_stack(self.phase, other.phase)
-
-        phase = np.logical_xor(phase1, phase2)
-        return StabilizerTable(pauli, phase)
+        return self._tensor(self, other)
 
     def expand(self, other):
         """Return the expand output product of two tables.
@@ -524,10 +521,7 @@ class StabilizerTable(PauliTable):
         """
         if not isinstance(other, StabilizerTable):
             other = StabilizerTable(other)
-        pauli = super().expand(other)
-        phase1, phase2 = self._block_stack(self.phase, other.phase)
-        phase = np.logical_xor(phase1, phase2)
-        return StabilizerTable(pauli, phase)
+        return self._tensor(other, self)
 
     def compose(self, other, qargs=None, front=False):
         """Return the compose output product of two tables.
@@ -576,6 +570,8 @@ class StabilizerTable(PauliTable):
         Raises:
             QiskitError: if other cannot be converted to a StabilizerTable.
         """
+        if qargs is None:
+            qargs = getattr(other, 'qargs', None)
         if not isinstance(other, StabilizerTable):
             other = StabilizerTable(other)
         if qargs is None and other.num_qubits != self.num_qubits:
@@ -650,7 +646,14 @@ class StabilizerTable(PauliTable):
         Raises:
             QiskitError: if other cannot be converted to a StabilizerTable.
         """
-        return super().dot(other, qargs=qargs)
+        return self.compose(other, qargs=qargs, front=True)
+
+    @classmethod
+    def _tensor(cls, a, b):
+        pauli = super()._tensor(a, b)
+        phase1, phase2 = a._block_stack(a.phase, b.phase)
+        phase = np.logical_xor(phase1, phase2)
+        return StabilizerTable(pauli, phase)
 
     def _add(self, other, qargs=None):
         """Append with another StabilizerTable.
@@ -1043,3 +1046,7 @@ class StabilizerTable(PauliTable):
                                            self.obj.phase[key],
                                            sparse=sparse)
         return MatrixIterator(self)
+
+
+# Update docstrings for API docs
+generate_apidocs(StabilizerTable)
diff --git a/qiskit/quantum_info/states/densitymatrix.py b/qiskit/quantum_info/states/densitymatrix.py
index f769492e1..acc65fb3c 100644
--- a/qiskit/quantum_info/states/densitymatrix.py
+++ b/qiskit/quantum_info/states/densitymatrix.py
@@ -22,7 +22,7 @@ from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.states.quantum_state import QuantumState
-from qiskit.quantum_info.operators.tolerances import TolerancesMixin
+from qiskit.quantum_info.operators.mixins.tolerances import TolerancesMixin
 from qiskit.quantum_info.operators.op_shape import OpShape
 from qiskit.quantum_info.operators.operator import Operator
 from qiskit.quantum_info.operators.scalar_op import ScalarOp
diff --git a/qiskit/quantum_info/states/statevector.py b/qiskit/quantum_info/states/statevector.py
index c2bfa48b6..0c4a143e9 100644
--- a/qiskit/quantum_info/states/statevector.py
+++ b/qiskit/quantum_info/states/statevector.py
@@ -24,7 +24,7 @@ from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.instruction import Instruction
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.states.quantum_state import QuantumState
-from qiskit.quantum_info.operators.tolerances import TolerancesMixin
+from qiskit.quantum_info.operators.mixins.tolerances import TolerancesMixin
 from qiskit.quantum_info.operators.operator import Operator
 from qiskit.quantum_info.operators.op_shape import OpShape
 from qiskit.quantum_info.operators.predicates import matrix_equal
diff --git a/releasenotes/notes/base-operator-api-8465ba5ad8856d9f.yaml b/releasenotes/notes/base-operator-api-8465ba5ad8856d9f.yaml
new file mode 100644
index 000000000..524475194
--- /dev/null
+++ b/releasenotes/notes/base-operator-api-8465ba5ad8856d9f.yaml
@@ -0,0 +1,7 @@
+---
+upgrade:
+  - |
+    The internal API of the
+    ``qiskit.quantum_info.operators.base_operator.BaseOperator`` abstract class
+    has been changed to simply making operator subclasses and to improve the
+    handling of subsystem dimensions.


Test patch:
diff --git a/test/python/quantum_info/operators/channel/test_chi.py b/test/python/quantum_info/operators/channel/test_chi.py
index 6ac2a01cc..d441e5d9c 100644
--- a/test/python/quantum_info/operators/channel/test_chi.py
+++ b/test/python/quantum_info/operators/channel/test_chi.py
@@ -243,12 +243,6 @@ class TestChi(ChannelTestCase):
         targ3 = Chi(self.depol_chi(1 - p_id3))
         self.assertEqual(chan3, targ3)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        chan = Chi(self.depol_chi(1))
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, chan.power, 0.5)
-
     def test_add(self):
         """Test add method."""
         mat1 = 0.5 * self.chiI
diff --git a/test/python/quantum_info/operators/channel/test_choi.py b/test/python/quantum_info/operators/channel/test_choi.py
index 6fbb207bc..c405f0cd6 100644
--- a/test/python/quantum_info/operators/channel/test_choi.py
+++ b/test/python/quantum_info/operators/channel/test_choi.py
@@ -335,12 +335,6 @@ class TestChoi(ChannelTestCase):
         targ3 = Choi(self.depol_choi(1 - p_id3))
         self.assertEqual(chan3, targ3)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        chan = Choi(self.depol_choi(1))
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, chan.power, 0.5)
-
     def test_add(self):
         """Test add method."""
         mat1 = 0.5 * self.choiI
diff --git a/test/python/quantum_info/operators/channel/test_kraus.py b/test/python/quantum_info/operators/channel/test_kraus.py
index 331142ac8..8cc984754 100644
--- a/test/python/quantum_info/operators/channel/test_kraus.py
+++ b/test/python/quantum_info/operators/channel/test_kraus.py
@@ -318,12 +318,6 @@ class TestKraus(ChannelTestCase):
         targ3b = rho @ Kraus(self.depol_kraus(1 - p_id3))
         self.assertEqual(rho @ chan3, targ3b)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        chan = Kraus(self.depol_kraus(0.9))
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, chan.power, 0.5)
-
     def test_add(self):
         """Test add method."""
         # Random input test state
diff --git a/test/python/quantum_info/operators/channel/test_ptm.py b/test/python/quantum_info/operators/channel/test_ptm.py
index 9f85783ad..d66fd6feb 100644
--- a/test/python/quantum_info/operators/channel/test_ptm.py
+++ b/test/python/quantum_info/operators/channel/test_ptm.py
@@ -236,12 +236,6 @@ class TestPTM(ChannelTestCase):
         targ3 = PTM(self.depol_ptm(1 - p_id3))
         self.assertEqual(chan3, targ3)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        chan = PTM(self.depol_ptm(1))
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, chan.power, 0.5)
-
     def test_add(self):
         """Test add method."""
         mat1 = 0.5 * self.ptmI
diff --git a/test/python/quantum_info/operators/channel/test_stinespring.py b/test/python/quantum_info/operators/channel/test_stinespring.py
index 55297db3e..b6bf3ff78 100644
--- a/test/python/quantum_info/operators/channel/test_stinespring.py
+++ b/test/python/quantum_info/operators/channel/test_stinespring.py
@@ -310,12 +310,6 @@ class TestStinespring(ChannelTestCase):
         rho_targ = rho_init @ Stinespring(self.depol_stine(1 - p_id3))
         self.assertEqual(rho_init @ chan, rho_targ)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        chan = Stinespring(self.depol_stine(0.9))
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, chan.power, 0.5)
-
     def test_add(self):
         """Test add method."""
         # Random input test state
diff --git a/test/python/quantum_info/operators/channel/test_superop.py b/test/python/quantum_info/operators/channel/test_superop.py
index 7e0cd957d..7d10f774f 100644
--- a/test/python/quantum_info/operators/channel/test_superop.py
+++ b/test/python/quantum_info/operators/channel/test_superop.py
@@ -541,12 +541,6 @@ class TestSuperOp(ChannelTestCase):
         targ3 = SuperOp(self.depol_sop(1 - p_id3))
         self.assertEqual(chan3, targ3)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        chan = SuperOp(self.depol_sop(1))
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, chan.power, 0.5)
-
     def test_add(self):
         """Test add method."""
         mat1 = 0.5 * self.sopI
diff --git a/test/python/quantum_info/operators/symplectic/test_pauli.py b/test/python/quantum_info/operators/symplectic/test_pauli.py
index 64f9071a3..9d48fb990 100644
--- a/test/python/quantum_info/operators/symplectic/test_pauli.py
+++ b/test/python/quantum_info/operators/symplectic/test_pauli.py
@@ -308,18 +308,6 @@ class TestPauli(QiskitTestCase):
         op = Pauli(label)
         self.assertTrue(op ** 2, iden)
 
-    def test_power_except(self):
-        """Test power method raises exceptions."""
-        op = Pauli('YXZ')
-        # Non-integer power raises error
-        self.assertRaises(QiskitError, op.power, 0.5)
-
-    def test_add_except(self):
-        """Test add method raises exceptions."""
-        P1 = Pauli('X')
-        P2 = Pauli('Z')
-        self.assertRaises(NotImplementedError, P1._add, P2)
-
     @data(1, 1.0, -1, -1.0, 1j, -1j)
     def test_multiply(self, val):
         """Test multiply method."""
diff --git a/test/python/quantum_info/operators/symplectic/test_pauli_table.py b/test/python/quantum_info/operators/symplectic/test_pauli_table.py
index 3f15e1b85..a39025eca 100644
--- a/test/python/quantum_info/operators/symplectic/test_pauli_table.py
+++ b/test/python/quantum_info/operators/symplectic/test_pauli_table.py
@@ -572,7 +572,7 @@ class TestPauliTableOperator(QiskitTestCase):
 
         value = pauli1.expand(pauli2)
         target = PauliTable.from_labels(
-            [j + i for i in labels1 for j in labels2])
+            [j + i for j in labels2 for i in labels1])
         self.assertEqual(value, target)
 
     def test_compose_1q(self):
diff --git a/test/python/quantum_info/operators/symplectic/test_stabilizer_table.py b/test/python/quantum_info/operators/symplectic/test_stabilizer_table.py
index ca76925fe..b0895db6f 100644
--- a/test/python/quantum_info/operators/symplectic/test_stabilizer_table.py
+++ b/test/python/quantum_info/operators/symplectic/test_stabilizer_table.py
@@ -848,23 +848,28 @@ class TestStabilizerTableMethods(QiskitTestCase):
                 self.assertTrue(np.all(i.toarray() == stab_mat(labels[idx])))
 
     def test_tensor(self):
-        """Test tensor and expand methods."""
+        """Test tensor method."""
         labels1 = ['-XX', 'YY']
         labels2 = ['III', '-ZZZ']
         stab1 = StabilizerTable.from_labels(labels1)
         stab2 = StabilizerTable.from_labels(labels2)
 
-        with self.subTest(msg='tensor'):
-            target = StabilizerTable.from_labels(
-                ['-XXIII', 'XXZZZ', 'YYIII', '-YYZZZ'])
-            value = stab1.tensor(stab2)
-            self.assertEqual(value, target)
+        target = StabilizerTable.from_labels(
+            ['-XXIII', 'XXZZZ', 'YYIII', '-YYZZZ'])
+        value = stab1.tensor(stab2)
+        self.assertEqual(value, target)
 
-        with self.subTest(msg='expand'):
-            target = StabilizerTable.from_labels(
-                ['-IIIXX', 'ZZZXX', 'IIIYY', '-ZZZYY'])
-            value = stab1.expand(stab2)
-            self.assertEqual(value, target)
+    def test_expand(self):
+        """Test expand method."""
+        labels1 = ['-XX', 'YY']
+        labels2 = ['III', '-ZZZ']
+        stab1 = StabilizerTable.from_labels(labels1)
+        stab2 = StabilizerTable.from_labels(labels2)
+
+        target = StabilizerTable.from_labels(
+            ['-IIIXX', 'IIIYY', 'ZZZXX', '-ZZZYY'])
+        value = stab1.expand(stab2)
+        self.assertEqual(value, target)
 
     def test_compose(self):
         """Test compose and dot methods."""
diff --git a/test/python/quantum_info/operators/test_operator.py b/test/python/quantum_info/operators/test_operator.py
index 6dad48245..ab6864ef5 100644
--- a/test/python/quantum_info/operators/test_operator.py
+++ b/test/python/quantum_info/operators/test_operator.py
@@ -528,8 +528,8 @@ class TestOperator(OperatorTestCase):
         assert_allclose(op12.data, Operator(mat12).data)
 
     def test_power_except(self):
-        """Test power method raises exceptions."""
-        op = Operator(self.rand_matrix(3, 3))
+        """Test power method raises exceptions if not square."""
+        op = Operator(self.rand_matrix(2, 3))
         # Non-integer power raises error
         self.assertRaises(QiskitError, op.power, 0.5)