diff --git a/.pylintdict b/.pylintdict
index 07e036844..95b983f48 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -188,6 +188,7 @@ wavefunction
 wecker
 williamson
 xs
+ys
 zemlin
 zi
 zj
diff --git a/qiskit_optimization/problems/quadratic_program.py b/qiskit_optimization/problems/quadratic_program.py
index e93b1fa81..2d07e8663 100644
--- a/qiskit_optimization/problems/quadratic_program.py
+++ b/qiskit_optimization/problems/quadratic_program.py
@@ -22,12 +22,10 @@
 import numpy as np
 from docplex.mp.model import Model
 from docplex.mp.model_reader import ModelReader
-from numpy import ndarray, zeros
-from scipy.sparse import spmatrix
-
+from numpy import ndarray
 from qiskit.exceptions import MissingOptionalLibraryError
-from qiskit.opflow import I, ListOp, OperatorBase, PauliOp, PauliSumOp, SummedOp
-from qiskit.quantum_info import Pauli
+from qiskit.opflow import OperatorBase, PauliSumOp
+from scipy.sparse import spmatrix
 
 from ..exceptions import QiskitOptimizationError
 from ..infinity import INFINITY
@@ -1015,94 +1013,10 @@ def to_ising(self) -> Tuple[OperatorBase, float]:
             QiskitOptimizationError: If a variable type is not binary.
             QiskitOptimizationError: If constraints exist in the problem.
         """
-        # if problem has variables that are not binary, raise an error
-        if self.get_num_vars() > self.get_num_binary_vars():
-            raise QiskitOptimizationError(
-                "The type of variable must be a binary variable. "
-                "Use a QuadraticProgramToQubo converter to convert "
-                "integer variables to binary variables. "
-                "If the problem contains continuous variables, "
-                "currently we can not apply VQE/QAOA directly. "
-                "you might want to use an ADMM optimizer "
-                "for the problem. "
-            )
-
-        # if constraints exist, raise an error
-        if self.linear_constraints or self.quadratic_constraints:
-            raise QiskitOptimizationError(
-                "An constraint exists. "
-                "The method supports only model with no constraints. "
-                "Use a QuadraticProgramToQubo converter. "
-                "It converts inequality constraints to equality "
-                "constraints, and then, it converters equality "
-                "constraints to penalty terms of the object function."
-            )
-
-        # initialize Hamiltonian.
-        num_nodes = self.get_num_vars()
-        pauli_list = []
-        offset = 0.0
-        zero = zeros(num_nodes, dtype=bool)
-
-        # set a sign corresponding to a maximized or minimized problem.
-        # sign == 1 is for minimized problem. sign == -1 is for maximized problem.
-        sense = self.objective.sense.value
-
-        # convert a constant part of the object function into Hamiltonian.
-        offset += self.objective.constant * sense
-
-        # convert linear parts of the object function into Hamiltonian.
-        for idx, coef in self.objective.linear.to_dict().items():
-            z_p = zeros(num_nodes, dtype=bool)
-            weight = coef * sense / 2
-            z_p[idx] = True
-
-            pauli_list.append([-weight, Pauli((z_p, zero))])
-            offset += weight
-
-        # convert quadratic parts of the object function into Hamiltonian.
-        # first merge coefficients (i, j) and (j, i)
-        coeffs = {}  # type: Dict
-        for (i, j), coeff in self.objective.quadratic.to_dict().items():
-            if j < i:  # type: ignore
-                coeffs[(j, i)] = coeffs.get((j, i), 0.0) + coeff
-            else:
-                coeffs[(i, j)] = coeffs.get((i, j), 0.0) + coeff
-
-        # create Pauli terms
-        for (i, j), coeff in coeffs.items():
-
-            weight = coeff * sense / 4
-
-            if i == j:
-                offset += weight
-            else:
-                z_p = zeros(num_nodes, dtype=bool)
-                z_p[i] = True
-                z_p[j] = True
-                pauli_list.append([weight, Pauli((z_p, zero))])
-
-            z_p = zeros(num_nodes, dtype=bool)
-            z_p[i] = True
-            pauli_list.append([-weight, Pauli((z_p, zero))])
-
-            z_p = zeros(num_nodes, dtype=bool)
-            z_p[j] = True
-            pauli_list.append([-weight, Pauli((z_p, zero))])
-
-            offset += weight
-
-        # Remove paulis whose coefficients are zeros.
-        qubit_op = sum(PauliOp(pauli, coeff=coeff) for coeff, pauli in pauli_list)
-
-        # qubit_op could be the integer 0, in this case return an identity operator of
-        # appropriate size
-        if isinstance(qubit_op, OperatorBase):
-            qubit_op = qubit_op.reduce()
-        else:
-            qubit_op = I ^ num_nodes
+        # pylint: disable=cyclic-import
+        from ..translators.ising import to_ising
 
-        return qubit_op, offset
+        return to_ising(self)
 
     def from_ising(
         self,
@@ -1129,102 +1043,12 @@ def from_ising(
             QiskitOptimizationError: If there are more than 2 Pauli Zs in any Pauli term
             NotImplementedError: If the input operator is a ListOp
         """
-        if isinstance(qubit_op, PauliSumOp):
-            qubit_op = qubit_op.to_pauli_op()
-
-        # No support for ListOp yet, this can be added in future
-        # pylint: disable=unidiomatic-typecheck
-        if type(qubit_op) == ListOp:
-            raise NotImplementedError(
-                "Conversion of a ListOp is not supported, convert each "
-                "operator in the ListOp separately."
-            )
-
-        # add binary variables
-        for i in range(qubit_op.num_qubits):
-            self.binary_var(name="x_{0}".format(i))
-
-        # Create a QUBO matrix
-        # The Qubo matrix is an upper triangular matrix.
-        # Diagonal elements in the QUBO matrix are for linear terms of the qubit operator.
-        # The other elements in the QUBO matrix are for quadratic terms of the qubit operator.
-        qubo_matrix = zeros((qubit_op.num_qubits, qubit_op.num_qubits))
-
-        if not isinstance(qubit_op, SummedOp):
-            pauli_list = [qubit_op.to_pauli_op()]
-        else:
-            pauli_list = qubit_op.to_pauli_op()
-
-        for pauli_op in pauli_list:
-            pauli_op = pauli_op.to_pauli_op()
-            pauli = pauli_op.primitive
-            coeff = pauli_op.coeff
-            # Count the number of Pauli Zs in a Pauli term
-            lst_z = pauli.z.tolist()
-            z_index = [i for i, z in enumerate(lst_z) if z is True]
-            num_z = len(z_index)
-
-            # Add its weight of the Pauli term to the corresponding element of QUBO matrix
-            if num_z == 1:
-                qubo_matrix[z_index[0], z_index[0]] = coeff.real
-            elif num_z == 2:
-                qubo_matrix[z_index[0], z_index[1]] = coeff.real
-            else:
-                raise QiskitOptimizationError(
-                    "There are more than 2 Pauli Zs in the Pauli term {}".format(pauli.z)
-                )
+        # pylint: disable=cyclic-import
+        from ..translators.ising import from_ising
 
-            # If there are Pauli Xs in the Pauli term, raise an error
-            lst_x = pauli.x.tolist()
-            x_index = [i for i, x in enumerate(lst_x) if x is True]
-            if len(x_index) > 0:
-                raise QiskitOptimizationError("Pauli Xs exist in the Pauli {}".format(pauli.x))
-
-        # Initialize dicts for linear terms and quadratic terms
-        linear_terms = {}
-        quadratic_terms = {}
-
-        # For quadratic pauli terms of operator
-        # x_i * x_ j = (1 - Z_i - Z_j + Z_i * Z_j)/4
-        for i, row in enumerate(qubo_matrix):
-            for j, weight in enumerate(row):
-                # Focus on the upper triangular matrix
-                if j <= i:
-                    continue
-                # Add a quadratic term to the object function of `QuadraticProgram`
-                # The coefficient of the quadratic term in `QuadraticProgram` is
-                # 4 * weight of the pauli
-                coef = weight * 4
-                quadratic_terms[i, j] = coef
-                # Sub the weight of the quadratic pauli term from the QUBO matrix
-                qubo_matrix[i, j] -= weight
-                # Sub the weight of the linear pauli term from the QUBO matrix
-                qubo_matrix[i, i] += weight
-                qubo_matrix[j, j] += weight
-                # Sub the weight from offset
-                offset -= weight
-
-        # After processing quadratic pauli terms, only linear paulis are left
-        # x_i = (1 - Z_i)/2
-        for i in range(qubit_op.num_qubits):
-            weight = qubo_matrix[i, i]
-            # Add a linear term to the object function of `QuadraticProgram`
-            # The coefficient of the linear term in `QuadraticProgram` is
-            # 2 * weight of the pauli
-            coef = weight * 2
-            if linear:
-                # If the linear option is True, add it into linear_terms
-                linear_terms[i] = -coef
-            else:
-                # Else, add it into quadratic_terms as a diagonal element.
-                quadratic_terms[i, i] = -coef
-            # Sub the weight of the linear pauli term from the QUBO matrix
-            qubo_matrix[i, i] -= weight
-            offset += weight
-
-        # Set the objective function
-        self.minimize(constant=offset, linear=linear_terms, quadratic=quadratic_terms)
-        offset -= offset
+        other = from_ising(qubit_op, offset, linear)
+        for attr, val in vars(other).items():
+            setattr(self, attr, val)
 
     def get_feasibility_info(
         self, x: Union[List[float], np.ndarray]
diff --git a/qiskit_optimization/translators/__init__.py b/qiskit_optimization/translators/__init__.py
index 4f8655f0d..310681007 100644
--- a/qiskit_optimization/translators/__init__.py
+++ b/qiskit_optimization/translators/__init__.py
@@ -29,9 +29,19 @@
    to_docplex_mp
    from_gurobipy
    to_gurobipy
+   from_ising
+   to_ising
 """
 
 from .docplex_mp import from_docplex_mp, to_docplex_mp
 from .gurobipy import from_gurobipy, to_gurobipy
+from .ising import from_ising, to_ising
 
-_all = ["from_docplex_mp", "to_docplex_mp", "from_gurobipy", "to_gurobipy"]
+_all = [
+    "from_docplex_mp",
+    "to_docplex_mp",
+    "from_gurobipy",
+    "to_gurobipy",
+    "from_ising",
+    "to_ising",
+]
diff --git a/qiskit_optimization/translators/ising.py b/qiskit_optimization/translators/ising.py
new file mode 100644
index 000000000..0c5ecf264
--- /dev/null
+++ b/qiskit_optimization/translators/ising.py
@@ -0,0 +1,222 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2019, 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.
+
+"""Translator between an Ising Hamiltonian and a quadratic program"""
+
+import math
+from typing import Tuple, Union
+
+import numpy as np
+
+from qiskit.opflow import I, ListOp, OperatorBase, PauliOp, PauliSumOp, SummedOp
+from qiskit.quantum_info import Pauli
+from qiskit_optimization.exceptions import QiskitOptimizationError
+from qiskit_optimization.problems.quadratic_program import QuadraticProgram
+
+
+def to_ising(quad_prog: QuadraticProgram) -> Tuple[OperatorBase, float]:
+    """Return the Ising Hamiltonian of this problem.
+
+    Variables are mapped to qubits in the same order, i.e.,
+    i-th variable is mapped to i-th qubit.
+    See https://github.com/Qiskit/qiskit-terra/issues/1148 for details.
+
+    Returns:
+        qubit_op: The qubit operator for the problem
+        offset: The constant value in the Ising Hamiltonian.
+
+    Raises:
+        QiskitOptimizationError: If an integer variable or a continuous variable exists
+            in the problem.
+        QiskitOptimizationError: If constraints exist in the problem.
+    """
+    # if problem has variables that are not binary, raise an error
+    if quad_prog.get_num_vars() > quad_prog.get_num_binary_vars():
+        raise QiskitOptimizationError(
+            "The type of all variables must be binary. "
+            "You can use `QuadraticProgramToQubo` converter "
+            "to convert integer variables to binary variables. "
+            "If the problem contains continuous variables, `to_ising` cannot handle it. "
+            "You might be able to solve it with `ADMMOptimizer`."
+        )
+
+    # if constraints exist, raise an error
+    if quad_prog.linear_constraints or quad_prog.quadratic_constraints:
+        raise QiskitOptimizationError(
+            "There must be no constraint in the problem. "
+            "You can use `QuadraticProgramToQubo` converter "
+            "to convert constraints to penalty terms of the objective function."
+        )
+
+    # initialize Hamiltonian.
+    num_nodes = quad_prog.get_num_vars()
+    pauli_list = []
+    offset = 0.0
+    zero = np.zeros(num_nodes, dtype=bool)
+
+    # set a sign corresponding to a maximized or minimized problem.
+    # sign == 1 is for minimized problem. sign == -1 is for maximized problem.
+    sense = quad_prog.objective.sense.value
+
+    # convert a constant part of the object function into Hamiltonian.
+    offset += quad_prog.objective.constant * sense
+
+    # convert linear parts of the object function into Hamiltonian.
+    for idx, coef in quad_prog.objective.linear.to_dict().items():
+        z_p = zero.copy()
+        weight = coef * sense / 2
+        z_p[idx] = True
+
+        pauli_list.append(PauliOp(Pauli((z_p, zero)), -weight))
+        offset += weight
+
+    # create Pauli terms
+    for (i, j), coeff in quad_prog.objective.quadratic.to_dict().items():
+        weight = coeff * sense / 4
+
+        if i == j:
+            offset += weight
+        else:
+            z_p = zero.copy()
+            z_p[i] = True
+            z_p[j] = True
+            pauli_list.append(PauliOp(Pauli((z_p, zero)), weight))
+
+        z_p = zero.copy()
+        z_p[i] = True
+        pauli_list.append(PauliOp(Pauli((z_p, zero)), -weight))
+
+        z_p = zero.copy()
+        z_p[j] = True
+        pauli_list.append(PauliOp(Pauli((z_p, zero)), -weight))
+
+        offset += weight
+
+    # Remove paulis whose coefficients are zeros.
+    qubit_op = sum(pauli_list)
+
+    # qubit_op could be the integer 0, in this case return an identity operator of
+    # appropriate size
+    if isinstance(qubit_op, OperatorBase):
+        qubit_op = qubit_op.reduce()
+    else:
+        qubit_op = I ^ num_nodes
+
+    return qubit_op, offset
+
+
+def from_ising(
+    qubit_op: Union[OperatorBase, PauliSumOp],
+    offset: float = 0.0,
+    linear: bool = False,
+) -> QuadraticProgram:
+    r"""Create a quadratic program from a qubit operator and a shift value.
+
+    Variables are mapped to qubits in the same order, i.e.,
+    i-th variable is mapped to i-th qubit.
+    See https://github.com/Qiskit/qiskit-terra/issues/1148 for details.
+
+    Args:
+        qubit_op: The qubit operator of the problem.
+        offset: The constant term in the Ising Hamiltonian.
+        linear: If linear is True, :math:`x^2` is treated as a linear term
+            since :math:`x^2 = x` for :math:`x \in \{0,1\}`.
+            Otherwise, :math:`x^2` is treat as a quadratic term.
+            The default value is False.
+
+    Returns:
+        The quadratic program corresponding to the qubit operator.
+
+    Raises:
+        QiskitOptimizationError: if there are Pauli Xs or Ys in any Pauli term
+        QiskitOptimizationError: if there are more than 2 Pauli Zs in any Pauli term
+        QiskitOptimizationError: if any Pauli term has an imaginary coefficient
+        NotImplementedError: If the input operator is a ListOp
+    """
+    if isinstance(qubit_op, PauliSumOp):
+        qubit_op = qubit_op.to_pauli_op()
+
+    # No support for ListOp yet, this can be added in future
+    # pylint: disable=unidiomatic-typecheck
+    if type(qubit_op) == ListOp:
+        raise NotImplementedError(
+            "Conversion of a ListOp is not supported, convert each "
+            "operator in the ListOp separately."
+        )
+
+    quad_prog = QuadraticProgram()
+    quad_prog.binary_var_list(qubit_op.num_qubits)
+
+    if not isinstance(qubit_op, SummedOp):
+        pauli_list = [qubit_op.to_pauli_op()]
+    else:
+        pauli_list = qubit_op.to_pauli_op()
+
+    # prepare a matrix of coefficients of Pauli terms
+    # `pauli_coeffs_diag` is the diagonal part
+    # `pauli_coeffs_triu` is the upper triangular part
+    pauli_coeffs_diag = [0.0] * qubit_op.num_qubits
+    pauli_coeffs_triu = {}
+
+    for pauli_op in pauli_list:
+        pauli_op = pauli_op.to_pauli_op()
+        pauli = pauli_op.primitive
+        coeff = pauli_op.coeff
+
+        if not math.isclose(coeff.imag, 0.0, abs_tol=1e-10):
+            raise QiskitOptimizationError(f"Imaginary coefficient exists: {pauli_op}")
+
+        if np.any(pauli.x):
+            raise QiskitOptimizationError(f"Pauli X or Y exists in the Pauli term: {pauli}")
+
+        # indices of Pauli Zs in the Pauli term
+        z_index = np.where(pauli.z)[0]
+        num_z = len(z_index)
+
+        if num_z == 1:
+            pauli_coeffs_diag[z_index[0]] = coeff.real
+        elif num_z == 2:
+            pauli_coeffs_triu[z_index[0], z_index[1]] = coeff.real
+        else:
+            raise QiskitOptimizationError(
+                f"There are more than 2 Pauli Zs in the Pauli term: {pauli}"
+            )
+
+    linear_terms = {}
+    quadratic_terms = {}
+
+    # For quadratic pauli terms of operator
+    # x_i * x_j = (1 - Z_i - Z_j + Z_i * Z_j)/4
+    for (i, j), weight in pauli_coeffs_triu.items():
+        # Add a quadratic term to the object function of `QuadraticProgram`
+        # The coefficient of the quadratic term in `QuadraticProgram` is
+        # 4 * weight of the pauli
+        quadratic_terms[i, j] = 4 * weight
+        pauli_coeffs_diag[i] += weight
+        pauli_coeffs_diag[j] += weight
+        offset -= weight
+
+    # After processing quadratic pauli terms, only linear paulis are left
+    # x_i = (1 - Z_i)/2
+    for i, weight in enumerate(pauli_coeffs_diag):
+        # Add a linear term to the object function of `QuadraticProgram`
+        # The coefficient of the linear term in `QuadraticProgram` is
+        # 2 * weight of the pauli
+        if linear:
+            linear_terms[i] = -2 * weight
+        else:
+            quadratic_terms[i, i] = -2 * weight
+        offset += weight
+
+    quad_prog.minimize(constant=offset, linear=linear_terms, quadratic=quadratic_terms)
+
+    return quad_prog
diff --git a/releasenotes/notes/add-ising-translator-b93e22adda321e56.yaml b/releasenotes/notes/add-ising-translator-b93e22adda321e56.yaml
new file mode 100644
index 000000000..01e503151
--- /dev/null
+++ b/releasenotes/notes/add-ising-translator-b93e22adda321e56.yaml
@@ -0,0 +1,5 @@
+features:
+  - |
+    Adds translators between Ising Hamiltonian and ``QuadraticProgram``,
+    :meth:`~qiskit_optimization.translators.from_ising` and
+    :meth:`~qiskit_optimization.translators.to_ising`.
diff --git a/test/translators/test_ising.py b/test/translators/test_ising.py
new file mode 100644
index 000000000..7a37a0185
--- /dev/null
+++ b/test/translators/test_ising.py
@@ -0,0 +1,152 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 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.
+
+"""Test from_ising and to_ising"""
+
+from test.optimization_test_case import QiskitOptimizationTestCase
+
+import numpy as np
+
+from qiskit.opflow import PauliSumOp
+from qiskit_optimization.exceptions import QiskitOptimizationError
+from qiskit_optimization.problems import QuadraticProgram
+from qiskit_optimization.translators import from_ising, to_ising
+
+
+class TestIsingTranslator(QiskitOptimizationTestCase):
+    """Test from_ising and to_ising"""
+
+    def test_to_ising(self):
+        """test to_ising"""
+
+        with self.subTest("minimize"):
+            # minimize: x + x * y
+            # subject to: x, y \in {0, 1}
+            q_p = QuadraticProgram("test")
+            q_p.binary_var(name="x")
+            q_p.binary_var(name="y")
+            q_p.minimize(linear={"x": 1}, quadratic={("x", "y"): 1})
+            op, offset = to_ising(q_p)
+            op_ref = PauliSumOp.from_list([("ZI", -0.25), ("IZ", -0.75), ("ZZ", 0.25)])
+            np.testing.assert_array_almost_equal(op.to_matrix(), op_ref.to_matrix())
+            self.assertAlmostEqual(offset, 0.75)
+
+        with self.subTest("maximize"):
+            # maximize: x + x * y
+            # subject to: x, y \in {0, 1}
+            q_p = QuadraticProgram("test")
+            q_p.binary_var(name="x")
+            q_p.binary_var(name="y")
+            q_p.maximize(linear={"x": 1}, quadratic={("x", "y"): 1})
+            op, offset = to_ising(q_p)
+            op_ref = PauliSumOp.from_list([("ZI", 0.25), ("IZ", 0.75), ("ZZ", -0.25)])
+            np.testing.assert_array_almost_equal(op.to_matrix(), op_ref.to_matrix())
+            self.assertAlmostEqual(offset, -0.75)
+
+    def test_to_ising2(self):
+        """test to_ising 2"""
+
+        with self.subTest("minimize"):
+            # minimize: 1 - 2 * x1 - 2 * x2 + 4 * x1 * x2
+            # subject to: x, y \in {0, 1}
+            q_p = QuadraticProgram("test")
+            q_p.binary_var(name="x")
+            q_p.binary_var(name="y")
+            q_p.minimize(constant=1, linear={"x": -2, "y": -2}, quadratic={("x", "y"): 4})
+            op, offset = to_ising(q_p)
+            op_ref = PauliSumOp.from_list([("ZZ", 1.0)])
+            np.testing.assert_array_almost_equal(op.to_matrix(), op_ref.to_matrix())
+            self.assertAlmostEqual(offset, 0.0)
+
+        with self.subTest("maximize"):
+            # maximize: 1 - 2 * x1 - 2 * x2 + 4 * x1 * x2
+            # subject to: x, y \in {0, 1}
+            q_p = QuadraticProgram("test")
+            q_p.binary_var(name="x")
+            q_p.binary_var(name="y")
+            q_p.maximize(constant=1, linear={"x": -2, "y": -2}, quadratic={("x", "y"): 4})
+            op, offset = to_ising(q_p)
+            op_ref = PauliSumOp.from_list([("ZZ", -1.0)])
+            np.testing.assert_array_almost_equal(op.to_matrix(), op_ref.to_matrix())
+            self.assertAlmostEqual(offset, 0.0)
+
+    def test_from_ising(self):
+        """test to_from_ising"""
+        # minimize: x + x * y
+        # subject to: x, y \in {0, 1}
+        op = PauliSumOp.from_list([("ZI", -0.25), ("IZ", -0.75), ("ZZ", 0.25)])
+        with self.subTest("linear: True"):
+            q_p = from_ising(op, 0.75, linear=True)
+            self.assertEqual(q_p.get_num_vars(), op.num_qubits)
+            self.assertEqual(q_p.get_num_linear_constraints(), 0)
+            self.assertEqual(q_p.get_num_quadratic_constraints(), 0)
+            self.assertAlmostEqual(q_p.objective.constant, 0)
+            np.testing.assert_array_almost_equal(q_p.objective.linear.to_array(), [1, 0])
+            np.testing.assert_array_almost_equal(
+                q_p.objective.quadratic.to_array(), [[0, 1], [0, 0]]
+            )
+
+        with self.subTest("linear: False"):
+            q_p = from_ising(op, 0.75, linear=False)
+            self.assertEqual(q_p.get_num_vars(), op.num_qubits)
+            self.assertEqual(q_p.get_num_linear_constraints(), 0)
+            self.assertEqual(q_p.get_num_quadratic_constraints(), 0)
+            self.assertAlmostEqual(q_p.objective.constant, 0)
+            np.testing.assert_array_almost_equal(q_p.objective.linear.to_array(), [0, 0])
+            np.testing.assert_array_almost_equal(
+                q_p.objective.quadratic.to_array(), [[1, 1], [0, 0]]
+            )
+
+    def test_from_ising2(self):
+        """test to_from_ising 2"""
+        # minimize: 1 - 2 * x1 - 2 * x2 + 4 * x1 * x2
+        # subject to: x, y \in {0, 1}
+        op = PauliSumOp.from_list([("ZZ", 1)])
+        with self.subTest("linear: True"):
+            q_p = from_ising(op, 0, linear=True)
+            self.assertEqual(q_p.get_num_vars(), op.num_qubits)
+            self.assertEqual(q_p.get_num_linear_constraints(), 0)
+            self.assertEqual(q_p.get_num_quadratic_constraints(), 0)
+            self.assertAlmostEqual(q_p.objective.constant, 1)
+            np.testing.assert_array_almost_equal(q_p.objective.linear.to_array(), [-2, -2])
+            np.testing.assert_array_almost_equal(
+                q_p.objective.quadratic.to_array(), [[0, 4], [0, 0]]
+            )
+
+        with self.subTest("linear: False"):
+            q_p = from_ising(op, 0, linear=False)
+            self.assertEqual(q_p.get_num_vars(), op.num_qubits)
+            self.assertEqual(q_p.get_num_linear_constraints(), 0)
+            self.assertEqual(q_p.get_num_quadratic_constraints(), 0)
+            self.assertAlmostEqual(q_p.objective.constant, 1)
+            np.testing.assert_array_almost_equal(q_p.objective.linear.to_array(), [0, 0])
+            np.testing.assert_array_almost_equal(
+                q_p.objective.quadratic.to_array(), [[-2, 4], [0, -2]]
+            )
+
+    def test_from_ising_pauli_with_invalid_paulis(self):
+        """test to_from_ising with invalid Pauli terms"""
+        with self.assertRaises(QiskitOptimizationError):
+            op = PauliSumOp.from_list([("IX", 1)])
+            _ = from_ising(op, 0)
+
+        with self.assertRaises(QiskitOptimizationError):
+            op = PauliSumOp.from_list([("IY", 1)])
+            _ = from_ising(op, 0)
+
+        with self.assertRaises(QiskitOptimizationError):
+            op = PauliSumOp.from_list([("ZZZ", 1)])
+            _ = from_ising(op, 0)
+
+        with self.assertRaises(QiskitOptimizationError):
+            op = PauliSumOp.from_list([("IZ", 1j)])
+            _ = from_ising(op, 0)