diff --git a/.pylintdict b/.pylintdict
index 5b1d7af55..5378ad385 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -126,6 +126,7 @@ nosignatures
 np
 num
 numpy
+numpyminimumeigensolver
 october
 optimality
 optimizationresult
diff --git a/README.md b/README.md
index 5e55d38ce..cedd16911 100644
--- a/README.md
+++ b/README.md
@@ -68,9 +68,9 @@ from docplex.mp.model import Model
 from qiskit_optimization.algorithms import MinimumEigenOptimizer
 from qiskit_optimization.translators import from_docplex_mp
 
-from qiskit.utils import algorithm_globals, QuantumInstance
-from qiskit import BasicAer
-from qiskit.algorithms import QAOA
+from qiskit.utils import algorithm_globals
+from qiskit.primitives import Sampler
+from qiskit.algorithms.minimum_eigensolvers import QAOA
 from qiskit.algorithms.optimizers import SPSA
 
 # Generate a graph of 4 nodes
@@ -97,9 +97,8 @@ seed = 1234
 algorithm_globals.random_seed = seed
 
 spsa = SPSA(maxiter=250)
-backend = BasicAer.get_backend('qasm_simulator')
-q_i = QuantumInstance(backend=backend, seed_simulator=seed, seed_transpiler=seed)
-qaoa = QAOA(optimizer=spsa, reps=5, quantum_instance=q_i)
+sampler = Sampler()
+qaoa = QAOA(sampler=sampler, optimizer=spsa, reps=5)
 algorithm = MinimumEigenOptimizer(qaoa)
 result = algorithm.solve(problem)
 print(result.prettyprint())  # prints solution, x=[1, 0, 1, 0], the cost, fval=4
diff --git a/qiskit_optimization/__init__.py b/qiskit_optimization/__init__.py
index 5753d7c90..23c9443cd 100644
--- a/qiskit_optimization/__init__.py
+++ b/qiskit_optimization/__init__.py
@@ -26,10 +26,11 @@
 A uniform interface as well as automatic conversion between different problem representations
 allows users to solve problems using a large set of algorithms, from variational quantum algorithms,
 such as the Quantum Approximate Optimization Algorithm
-(:class:`~qiskit.algorithms.QAOA`), to
+(:class:`~qiskit.algorithms.minimum_eigensolver.QAOA`), to
 `Grover Adaptive Search <https://arxiv.org/abs/quant-ph/9607014>`_
 (:class:`~algorithms.GroverOptimizer`), leveraging
-fundamental :mod:`~qiskit.algorithms` provided by Qiskit Terra. Furthermore, the modular design
+fundamental :mod:`~qiskit.algorithms.minimum_eigensolver` provided by Qiskit Terra.
+Furthermore, the modular design
 of the optimization module allows it to be easily extended and facilitates rapid development and
 testing of new algorithms. Compatible classical optimizers are also provided for testing,
 validation, and benchmarking.
diff --git a/qiskit_optimization/algorithms/admm_optimizer.py b/qiskit_optimization/algorithms/admm_optimizer.py
index 40c3911e7..881e05be7 100644
--- a/qiskit_optimization/algorithms/admm_optimizer.py
+++ b/qiskit_optimization/algorithms/admm_optimizer.py
@@ -17,21 +17,21 @@
 from typing import List, Optional, Tuple, cast
 
 import numpy as np
-from qiskit.algorithms import NumPyMinimumEigensolver
+from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
 
+from ..converters import MaximizeToMinimize
+from ..problems.constraint import Constraint
+from ..problems.linear_constraint import LinearConstraint
+from ..problems.linear_expression import LinearExpression
+from ..problems.quadratic_program import QuadraticProgram
+from ..problems.variable import Variable, VarType
 from .minimum_eigen_optimizer import MinimumEigenOptimizer
 from .optimization_algorithm import (
-    OptimizationResultStatus,
     OptimizationAlgorithm,
     OptimizationResult,
+    OptimizationResultStatus,
 )
 from .slsqp_optimizer import SlsqpOptimizer
-from ..problems.constraint import Constraint
-from ..problems.linear_constraint import LinearConstraint
-from ..problems.linear_expression import LinearExpression
-from ..problems.quadratic_program import QuadraticProgram
-from ..problems.variable import VarType, Variable
-from ..converters import MaximizeToMinimize
 
 UPDATE_RHO_BY_TEN_PERCENT = 0
 UPDATE_RHO_BY_RESIDUALS = 1
diff --git a/qiskit_optimization/algorithms/grover_optimizer.py b/qiskit_optimization/algorithms/grover_optimizer.py
index d6a75cde0..b8359c0c6 100644
--- a/qiskit_optimization/algorithms/grover_optimizer.py
+++ b/qiskit_optimization/algorithms/grover_optimizer.py
@@ -16,6 +16,7 @@
 import math
 from copy import deepcopy
 from typing import Optional, Dict, Union, List, cast
+import warnings
 
 import numpy as np
 
@@ -24,6 +25,7 @@
 from qiskit.utils import QuantumInstance, algorithm_globals
 from qiskit.algorithms.amplitude_amplifiers.grover import Grover
 from qiskit.circuit.library import QuadraticForm
+from qiskit.primitives import BaseSampler
 from qiskit.providers import Backend
 from qiskit.quantum_info import partial_trace
 from .optimization_algorithm import (
@@ -36,6 +38,7 @@
     QuadraticProgramToQubo,
     QuadraticProgramConverter,
 )
+from ..exceptions import QiskitOptimizationError
 from ..problems import Variable
 from ..problems.quadratic_program import QuadraticProgram
 
@@ -54,6 +57,7 @@ def __init__(
             Union[QuadraticProgramConverter, List[QuadraticProgramConverter]]
         ] = None,
         penalty: Optional[float] = None,
+        sampler: Optional[BaseSampler] = None,
     ) -> None:
         """
         Args:
@@ -66,20 +70,35 @@ def __init__(
                 :class:`~qiskit_optimization.converters.QuadraticProgramToQubo` will be used.
             penalty: The penalty factor used in the default
                 :class:`~qiskit_optimization.converters.QuadraticProgramToQubo` converter
+            sampler: A Sampler to use for sampling the results of the circuits.
 
         Raises:
+            ValueError: If both a quantum instance and sampler are set.
             TypeError: When there one of converters is an invalid type.
         """
         self._num_value_qubits = num_value_qubits
         self._num_key_qubits = 0
         self._n_iterations = num_iterations
-        self._quantum_instance = None  # type: Optional[QuantumInstance]
         self._circuit_results = {}  # type: dict
+        self._converters = self._prepare_converters(converters, penalty)
 
-        if quantum_instance is not None:
-            self.quantum_instance = quantum_instance
+        if quantum_instance is not None and sampler is not None:
+            raise ValueError("Only one of quantum_instance or sampler can be passed, not both!")
 
-        self._converters = self._prepare_converters(converters, penalty)
+        self._quantum_instance = None  # type: Optional[QuantumInstance]
+        if quantum_instance is not None:
+            warnings.warn(
+                "The quantum_instance argument has been superseded by the sampler argument. "
+                "This argument will be deprecated in a future release and subsequently "
+                "removed after that.",
+                category=PendingDeprecationWarning,
+                stacklevel=2,
+            )
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", category=PendingDeprecationWarning)
+                self.quantum_instance = quantum_instance
+
+        self._sampler = sampler
 
     @property
     def quantum_instance(self) -> QuantumInstance:
@@ -88,6 +107,13 @@ def quantum_instance(self) -> QuantumInstance:
         Returns:
             The quantum instance used in the algorithm.
         """
+        warnings.warn(
+            "The quantum_instance argument has been superseded by the sampler argument. "
+            "This argument will be deprecated in a future release and subsequently "
+            "removed after that.",
+            category=PendingDeprecationWarning,
+            stacklevel=2,
+        )
         return self._quantum_instance
 
     @quantum_instance.setter
@@ -97,6 +123,13 @@ def quantum_instance(self, quantum_instance: Union[Backend, QuantumInstance]) ->
         Args:
             quantum_instance: The quantum instance to be used in the algorithm.
         """
+        warnings.warn(
+            "The GroverOptimizer.quantum_instance setter is pending deprecation. "
+            "This property will be deprecated in a future release and subsequently "
+            "removed after that.",
+            category=PendingDeprecationWarning,
+            stacklevel=2,
+        )
         if isinstance(quantum_instance, Backend):
             self._quantum_instance = QuantumInstance(quantum_instance)
         else:
@@ -162,11 +195,16 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
             The result of the optimizer applied to the problem.
 
         Raises:
+            ValueError: If a quantum instance or a sampler has not been provided.
+            ValueError: If both a quantum instance and sampler are set.
             AttributeError: If the quantum instance has not been set.
             QiskitOptimizationError: If the problem is incompatible with the optimizer.
         """
-        if self.quantum_instance is None:
-            raise AttributeError("The quantum instance or backend has not been set.")
+        if self._sampler is None and self._quantum_instance is None:
+            raise ValueError("A quantum instance or sampler must be provided.")
+
+        if self._quantum_instance is not None and self._sampler is not None:
+            raise ValueError("Only one of quantum_instance or sampler can be passed, not both!")
 
         self._verify_compatibility(problem)
 
@@ -199,7 +237,7 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
         # Initialize oracle helper object.
         qr_key_value = QuantumRegister(self._num_key_qubits + self._num_value_qubits)
         orig_constant = problem_.objective.constant
-        measurement = not self.quantum_instance.is_statevector
+        measurement = self._quantum_instance is None or not self._quantum_instance.is_statevector
         oracle, is_good_state = self._get_oracle(qr_key_value)
 
         while not optimum_found:
@@ -246,15 +284,19 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
                     threshold = optimum_value
 
                     # trace out work qubits and store samples
-                    if self._quantum_instance.is_statevector:
-                        indices = list(range(n_key, len(outcome)))
-                        rho = partial_trace(self._circuit_results, indices)
-                        self._circuit_results = cast(Dict, np.diag(rho.data) ** 0.5)
-                    else:
+                    if self._sampler is not None:
                         self._circuit_results = {
                             i[-1 * n_key :]: v for i, v in self._circuit_results.items()
                         }
-
+                    else:
+                        if self._quantum_instance.is_statevector:
+                            indices = list(range(n_key, len(outcome)))
+                            rho = partial_trace(self._circuit_results, indices)
+                            self._circuit_results = cast(Dict, np.diag(rho.data) ** 0.5)
+                        else:
+                            self._circuit_results = {
+                                i[-1 * n_key :]: v for i, v in self._circuit_results.items()
+                            }
                     raw_samples = self._eigenvector_to_solutions(
                         self._circuit_results, problem_init
                     )
@@ -312,33 +354,52 @@ def solve(self, problem: QuadraticProgram) -> OptimizationResult:
 
     def _measure(self, circuit: QuantumCircuit) -> str:
         """Get probabilities from the given backend, and picks a random outcome."""
-        probs = self._get_probs(circuit)
+        probs = self._get_prob_dist(circuit)
         logger.info("Frequencies: %s", probs)
         # Pick a random outcome.
         return algorithm_globals.random.choice(list(probs.keys()), 1, p=list(probs.values()))[0]
 
-    def _get_probs(self, qc: QuantumCircuit) -> Dict[str, float]:
+    def _get_prob_dist(self, qc: QuantumCircuit) -> Dict[str, float]:
         """Gets probabilities from a given backend."""
         # Execute job and filter results.
-        result = self.quantum_instance.execute(qc)
-        if self.quantum_instance.is_statevector:
-            state = result.get_statevector(qc)
-            if not isinstance(state, np.ndarray):
-                state = state.data
-            keys = [
-                bin(i)[2::].rjust(int(np.log2(len(state))), "0")[::-1] for i in range(0, len(state))
-            ]
-            probs = [abs(a) ** 2 for a in state]
-            total = math.fsum(probs)
-            probs = [p / total for p in probs]
-            hist = {key: prob for key, prob in zip(keys, probs) if prob > 0}
-            self._circuit_results = state
+        if self._sampler is not None:
+            job = self._sampler.run([qc])
+
+            try:
+                result = job.result()
+            except Exception as exc:
+                raise QiskitOptimizationError("Sampler job failed.") from exc
+            quasi_dist = result.quasi_dists[0]
+            bit_length = (len(quasi_dist) - 1).bit_length()
+            prob_dist = {f"{i:0{bit_length}b}"[::-1]: v for i, v in quasi_dist.items()}
+            self._circuit_results = {
+                f"{i:0{bit_length}b}": v**0.5
+                for i, v in quasi_dist.items()
+                if not np.isclose(v, 0)
+            }
         else:
-            state = result.get_counts(qc)
-            shots = self.quantum_instance.run_config.shots
-            hist = {key[::-1]: val / shots for key, val in sorted(state.items()) if val > 0}
-            self._circuit_results = {b: (v / shots) ** 0.5 for (b, v) in state.items()}
-        return hist
+            result = self._quantum_instance.execute(qc)
+            if self._quantum_instance.is_statevector:
+                state = result.get_statevector(qc)
+                if not isinstance(state, np.ndarray):
+                    state = state.data
+                keys = [
+                    bin(i)[2::].rjust(int(np.log2(len(state))), "0")[::-1]
+                    for i in range(0, len(state))
+                ]
+                probs = [abs(a) ** 2 for a in state]
+                total = math.fsum(probs)
+                probs = [p / total for p in probs]
+                prob_dist = {key: prob for key, prob in zip(keys, probs) if prob > 0}
+                self._circuit_results = state
+            else:
+                state = result.get_counts(qc)
+                shots = self._quantum_instance.run_config.shots
+                prob_dist = {
+                    key[::-1]: val / shots for key, val in sorted(state.items()) if val > 0
+                }
+                self._circuit_results = {b: (v / shots) ** 0.5 for (b, v) in state.items()}
+        return prob_dist
 
     @staticmethod
     def _bin_to_int(v: str, num_value_bits: int) -> int:
diff --git a/qiskit_optimization/algorithms/minimum_eigen_optimizer.py b/qiskit_optimization/algorithms/minimum_eigen_optimizer.py
index 0a86bba28..193790bd5 100644
--- a/qiskit_optimization/algorithms/minimum_eigen_optimizer.py
+++ b/qiskit_optimization/algorithms/minimum_eigen_optimizer.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2020, 2021.
+# (C) Copyright IBM 2020, 2022.
 #
 # 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
@@ -11,24 +11,37 @@
 # that they have been altered from the originals.
 
 """A wrapper for minimum eigen solvers to be used within the optimization module."""
-from typing import Optional, Union, List, cast
+from typing import List, Optional, Union, cast
 
 import numpy as np
+from qiskit.algorithms.minimum_eigen_solvers import MinimumEigensolver as LegacyMinimumEigensolver
+from qiskit.algorithms.minimum_eigen_solvers import (
+    MinimumEigensolverResult as LegacyMinimumEigensolverResult,
+)
+from qiskit.algorithms.minimum_eigensolvers import (
+    NumPyMinimumEigensolver,
+    NumPyMinimumEigensolverResult,
+    SamplingMinimumEigensolver,
+    SamplingMinimumEigensolverResult,
+)
+from qiskit.opflow import OperatorBase, PauliOp, PauliSumOp
 
-from qiskit.algorithms import MinimumEigensolver, MinimumEigensolverResult
-from qiskit.opflow import OperatorBase
+from ..converters.quadratic_program_to_qubo import QuadraticProgramConverter, QuadraticProgramToQubo
+from ..exceptions import QiskitOptimizationError
+from ..problems.quadratic_program import QuadraticProgram, Variable
 from .optimization_algorithm import (
-    OptimizationResultStatus,
     OptimizationAlgorithm,
     OptimizationResult,
+    OptimizationResultStatus,
     SolutionSample,
 )
-from ..exceptions import QiskitOptimizationError
-from ..converters.quadratic_program_to_qubo import (
-    QuadraticProgramToQubo,
-    QuadraticProgramConverter,
-)
-from ..problems.quadratic_program import QuadraticProgram, Variable
+
+MinimumEigensolver = Union[
+    SamplingMinimumEigensolver, NumPyMinimumEigensolver, LegacyMinimumEigensolver
+]
+MinimumEigensolverResult = Union[
+    SamplingMinimumEigensolverResult, NumPyMinimumEigensolverResult, LegacyMinimumEigensolverResult
+]
 
 
 class MinimumEigenOptimizationResult(OptimizationResult):
@@ -101,7 +114,7 @@ class MinimumEigenOptimizer(OptimizationAlgorithm):
 
     .. code-block::
 
-        from qiskit.algorithms import QAOA
+        from qiskit.algorithms.minimum_eigensolver import QAOA
         from qiskit_optimization.problems import QuadraticProgram
         from qiskit_optimization.algorithms import MinimumEigenOptimizer
         problem = QuadraticProgram()
@@ -141,7 +154,7 @@ def __init__(
         if not min_eigen_solver.supports_aux_operators():
             raise QiskitOptimizationError(
                 "Given MinimumEigensolver does not return the eigenstate "
-                + "and is not supported by the MinimumEigenOptimizer."
+                "and is not supported by the MinimumEigenOptimizer."
             )
         self._min_eigen_solver = min_eigen_solver
         self._penalty = penalty
@@ -206,6 +219,9 @@ def _solve_internal(
         # only try to solve non-empty Ising Hamiltonians
         eigen_result: Optional[MinimumEigensolverResult] = None
         if operator.num_qubits > 0:
+            # NumPyEigensolver does not accept PauliOp but PauliSumOp
+            if isinstance(operator, PauliOp):
+                operator = PauliSumOp.from_list([(operator.primitive.to_label(), operator.coeff)])
             # approximate ground state of operator using min eigen solver
             eigen_result = self._min_eigen_solver.compute_minimum_eigenvalue(operator)
             # analyze results
diff --git a/qiskit_optimization/algorithms/optimization_algorithm.py b/qiskit_optimization/algorithms/optimization_algorithm.py
index 050d78ca5..4b8fb3080 100644
--- a/qiskit_optimization/algorithms/optimization_algorithm.py
+++ b/qiskit_optimization/algorithms/optimization_algorithm.py
@@ -15,14 +15,16 @@
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
 from enum import Enum
-from typing import List, Union, Any, Optional, Dict, Type, Tuple, cast
+from typing import Any, Dict, List, Optional, Tuple, Type, Union, cast
 from warnings import warn
 
 import numpy as np
+from qiskit.opflow import DictStateFn, StateFn
+from qiskit.quantum_info import Statevector
+from qiskit.result import QuasiDistribution
 
-from qiskit.opflow import StateFn, DictStateFn
+from ..converters.quadratic_program_to_qubo import QuadraticProgramConverter, QuadraticProgramToQubo
 from ..exceptions import QiskitOptimizationError
-from ..converters.quadratic_program_to_qubo import QuadraticProgramToQubo, QuadraticProgramConverter
 from ..problems.quadratic_program import QuadraticProgram, Variable
 
 
@@ -518,7 +520,7 @@ def _interpret_samples(
 
     @staticmethod
     def _eigenvector_to_solutions(
-        eigenvector: Union[dict, np.ndarray, StateFn],
+        eigenvector: Union[QuasiDistribution, Statevector, dict, np.ndarray, StateFn],
         qubo: QuadraticProgram,
         min_probability: float = 1e-6,
     ) -> List[SolutionSample]:
@@ -566,7 +568,25 @@ def generate_solution(bitstr, qubo, probability):
             )
 
         solutions = []
-        if isinstance(eigenvector, dict):
+        if isinstance(eigenvector, QuasiDistribution):
+            probabilities = eigenvector.binary_probabilities()
+            # iterate over all samples
+            for bitstr, sampling_probability in probabilities.items():
+                # add the bitstring, if the sampling probability exceeds the threshold
+                if sampling_probability >= min_probability:
+                    solutions.append(generate_solution(bitstr, qubo, sampling_probability))
+
+        elif isinstance(eigenvector, Statevector):
+            probabilities = eigenvector.probabilities()
+            num_qubits = eigenvector.num_qubits
+            # iterate over all states and their sampling probabilities
+            for i, sampling_probability in enumerate(probabilities):
+                # add the i-th state if the sampling probability exceeds the threshold
+                if sampling_probability >= min_probability:
+                    bitstr = f"{i:b}".rjust(num_qubits, "0")
+                    solutions.append(generate_solution(bitstr, qubo, sampling_probability))
+
+        elif isinstance(eigenvector, dict):
             # When eigenvector is a dict, square the values since the values are normalized.
             # See https://github.com/Qiskit/qiskit-terra/pull/5496 for more details.
             probabilities = {bitstr: val**2 for (bitstr, val) in eigenvector.items()}
@@ -579,7 +599,6 @@ def generate_solution(bitstr, qubo, probability):
         elif isinstance(eigenvector, np.ndarray):
             num_qubits = int(np.log2(eigenvector.size))
             probabilities = np.abs(eigenvector * eigenvector.conj())
-
             # iterate over all states and their sampling probabilities
             for i, sampling_probability in enumerate(probabilities):
                 # add the i-th state if the sampling probability exceeds the threshold
@@ -588,6 +607,8 @@ def generate_solution(bitstr, qubo, probability):
                     solutions.append(generate_solution(bitstr, qubo, sampling_probability))
 
         else:
-            raise TypeError("Unsupported format of eigenvector. Provide a dict or numpy.ndarray.")
-
+            raise TypeError(
+                f"Eigenvector should be QuasiDistribution, Statevector, dict or numpy.ndarray. "
+                f"But, it was {type(eigenvector)}."
+            )
         return solutions
diff --git a/qiskit_optimization/algorithms/recursive_minimum_eigen_optimizer.py b/qiskit_optimization/algorithms/recursive_minimum_eigen_optimizer.py
index 26c1184a6..c3149f3fc 100644
--- a/qiskit_optimization/algorithms/recursive_minimum_eigen_optimizer.py
+++ b/qiskit_optimization/algorithms/recursive_minimum_eigen_optimizer.py
@@ -14,28 +14,22 @@
 
 from copy import deepcopy
 from enum import Enum
-from typing import Optional, Union, List, Tuple, Dict, cast
+from typing import Dict, List, Optional, Tuple, Union, cast
 
 import numpy as np
-from qiskit.algorithms import NumPyMinimumEigensolver
+from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
 from qiskit.utils.validation import validate_min
 
-from .minimum_eigen_optimizer import (
-    MinimumEigenOptimizer,
-    MinimumEigenOptimizationResult,
-)
+from ..converters.quadratic_program_to_qubo import QuadraticProgramConverter, QuadraticProgramToQubo
+from ..exceptions import QiskitOptimizationError
+from ..problems import Variable
+from ..problems.quadratic_program import QuadraticProgram
+from .minimum_eigen_optimizer import MinimumEigenOptimizationResult, MinimumEigenOptimizer
 from .optimization_algorithm import (
-    OptimizationResultStatus,
     OptimizationAlgorithm,
     OptimizationResult,
+    OptimizationResultStatus,
 )
-from ..converters.quadratic_program_to_qubo import (
-    QuadraticProgramToQubo,
-    QuadraticProgramConverter,
-)
-from ..exceptions import QiskitOptimizationError
-from ..problems import Variable
-from ..problems.quadratic_program import QuadraticProgram
 
 
 class IntermediateResult(Enum):
@@ -123,7 +117,7 @@ class RecursiveMinimumEigenOptimizer(OptimizationAlgorithm):
 
     .. code-block:: python
 
-        from qiskit.algorithms import QAOA
+        from qiskit.algorithms.minimum_eigensolver import QAOA
         from qiskit_optimization.problems import QuadraticProgram
         from qiskit_optimization.algorithms import (
             MinimumEigenOptimizer, RecursiveMinimumEigenOptimizer
@@ -167,7 +161,7 @@ def __init__(
             min_num_vars_optimizer: This optimizer is used after the recursive scheme for the
                 problem with the remaining variables. Default value is
                 :class:`~qiskit_optimization.algorithms.MinimumEigenOptimizer` created on top of
-                :class:`~qiskit.algorithms.minimum_eigen_solver.NumPyMinimumEigensolver`.
+                :class:`~qiskit.algorithms.minimum_eigensolver.NumPyMinimumEigensolver`.
             penalty: The factor that is used to scale the penalty terms corresponding to linear
                 equality constraints.
             history: Whether the intermediate results are stored.
diff --git a/qiskit_optimization/algorithms/warm_start_qaoa_optimizer.py b/qiskit_optimization/algorithms/warm_start_qaoa_optimizer.py
index 98016cbdd..c244df1f5 100644
--- a/qiskit_optimization/algorithms/warm_start_qaoa_optimizer.py
+++ b/qiskit_optimization/algorithms/warm_start_qaoa_optimizer.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2021.
+# (C) Copyright IBM 2021, 2022.
 #
 # 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
@@ -14,26 +14,20 @@
 
 import copy
 from abc import ABC, abstractmethod
-from typing import Optional, List, Union, Dict, Tuple, cast
+from typing import Dict, List, Optional, Tuple, Union, cast
 
 import numpy as np
 from qiskit import QuantumCircuit
-from qiskit.algorithms import QAOA
+from qiskit.algorithms import QAOA as LegacyQAOA
+from qiskit.algorithms.minimum_eigensolvers import QAOA
 from qiskit.circuit import Parameter
 
-from .minimum_eigen_optimizer import (
-    MinimumEigenOptimizer,
-    MinimumEigenOptimizationResult,
-)
-from .optimization_algorithm import (
-    OptimizationAlgorithm,
-    OptimizationResultStatus,
-    SolutionSample,
-)
 from ..converters.quadratic_program_converter import QuadraticProgramConverter
 from ..exceptions import QiskitOptimizationError
 from ..problems.quadratic_program import QuadraticProgram
 from ..problems.variable import VarType
+from .minimum_eigen_optimizer import MinimumEigenOptimizationResult, MinimumEigenOptimizer
+from .optimization_algorithm import OptimizationAlgorithm, OptimizationResultStatus, SolutionSample
 
 
 class BaseAggregator(ABC):
@@ -209,7 +203,7 @@ def __init__(
         self,
         pre_solver: OptimizationAlgorithm,
         relax_for_pre_solver: bool,
-        qaoa: QAOA,
+        qaoa: Union[QAOA, LegacyQAOA],
         epsilon: float = 0.25,
         num_initial_solutions: int = 1,
         warm_start_factory: Optional[WarmStartQAOAFactory] = None,
diff --git a/qiskit_optimization/applications/optimization_application.py b/qiskit_optimization/applications/optimization_application.py
index 9680695ce..32aae008e 100644
--- a/qiskit_optimization/applications/optimization_application.py
+++ b/qiskit_optimization/applications/optimization_application.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2018, 2021.
+# (C) Copyright IBM 2018, 2022.
 #
 # 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
@@ -11,13 +11,15 @@
 # that they have been altered from the originals.
 
 """An abstract class for optimization application classes."""
-from typing import Union, Dict
-from collections import OrderedDict
 from abc import ABC, abstractmethod
+from collections import OrderedDict
+from typing import Dict, Union
 
 import numpy as np
-
 from qiskit.opflow import StateFn
+from qiskit.quantum_info import Statevector
+from qiskit.result import QuasiDistribution
+
 from qiskit_optimization.algorithms import OptimizationResult
 from qiskit_optimization.problems.quadratic_program import QuadraticProgram
 
@@ -59,25 +61,45 @@ def _result_to_x(self, result: Union[OptimizationResult, np.ndarray]) -> np.ndar
         return x
 
     @staticmethod
-    def sample_most_likely(state_vector: Union[np.ndarray, Dict]) -> np.ndarray:
+    def sample_most_likely(
+        state_vector: Union[QuasiDistribution, Statevector, np.ndarray, Dict]
+    ) -> np.ndarray:
         """Compute the most likely binary string from state vector.
 
         Args:
-            state_vector: state vector or counts.
+            state_vector: state vector or counts or quasi-probabilities.
 
         Returns:
             binary string as numpy.ndarray of ints.
+
+        Raises:
+            ValueError: if state_vector is not QuasiDistribution, Statevector,
+                np.ndarray, or dict.
         """
-        if isinstance(state_vector, (OrderedDict, dict)):
+        if isinstance(state_vector, QuasiDistribution):
+            probabilities = state_vector.binary_probabilities()
+            binary_string = max(probabilities.items(), key=lambda kv: kv[1])[0]
+            x = np.asarray([int(y) for y in reversed(list(binary_string))])
+            return x
+        elif isinstance(state_vector, Statevector):
+            probabilities = state_vector.probabilities()
+            n = state_vector.num_qubits
+            k = np.argmax(np.abs(probabilities))
+            x = np.zeros(n)
+            for i in range(n):
+                x[i] = k % 2
+                k >>= 1
+            return x
+        elif isinstance(state_vector, (OrderedDict, dict)):
             # get the binary string with the largest count
-            binary_string = sorted(state_vector.items(), key=lambda kv: kv[1])[-1][0]
+            binary_string = max(state_vector.items(), key=lambda kv: kv[1])[0]
             x = np.asarray([int(y) for y in reversed(list(binary_string))])
             return x
         elif isinstance(state_vector, StateFn):
             binary_string = list(state_vector.sample().keys())[0]
             x = np.asarray([int(y) for y in reversed(list(binary_string))])
             return x
-        else:
+        elif isinstance(state_vector, np.ndarray):
             n = int(np.log2(state_vector.shape[0]))
             k = np.argmax(np.abs(state_vector))
             x = np.zeros(n)
@@ -85,3 +107,8 @@ def sample_most_likely(state_vector: Union[np.ndarray, Dict]) -> np.ndarray:
                 x[i] = k % 2
                 k >>= 1
             return x
+        else:
+            raise ValueError(
+                "state vector should be QuasiDistribution, Statevector, ndarray, or dict. "
+                f"But it is {type(state_vector)}."
+            )
diff --git a/qiskit_optimization/deprecation.py b/qiskit_optimization/deprecation.py
index c87ab4f7c..772d9791d 100644
--- a/qiskit_optimization/deprecation.py
+++ b/qiskit_optimization/deprecation.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2021.
+# (C) Copyright IBM 2021, 2022.
 #
 # 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
diff --git a/releasenotes/notes/add-primitives-support-31af39549b5e66e3.yaml b/releasenotes/notes/add-primitives-support-31af39549b5e66e3.yaml
new file mode 100644
index 000000000..94d48700d
--- /dev/null
+++ b/releasenotes/notes/add-primitives-support-31af39549b5e66e3.yaml
@@ -0,0 +1,24 @@
+---
+prelude: >
+    Qiskit Optimization 0.5 supports the new algorithms introduced in Qiskit Terra 0.22
+    which in turn rely on the `Qiskit Primitives <https://qiskit.org/documentation/apidoc/primitives.html>`_.
+    Qiskit Optimization 0.5 still supports the former algorithms based on :class:`qiskit.utils.QuantumInstance`,
+    but they will be deprecated and then removed, along with the support here, in future releases.
+
+features:
+  - |
+    The :class:`~.MinimumEigenOptimizer` class takes the primitives-based algorithms
+    (:class:`qiskit.algorithms.minimum_eigensolvers.SamplingMinimumEigensolver` and
+    :class:`qiskit.algorithms.minimum_eigensolvers.NumPyMinimumEigensolver`)
+    as ``min_eigen_solver`` argument.
+    The former algorithm :class:`qiskit.algorithms.MinimumEigensolver`
+    is pending deprecation and will be deprecated and subsequently removed in future releases.
+    Note that :class:`qiskit.algorithms.minimum_eigensolvers.SamplingVQE` supersedes
+    :class:`qiskit.algorithms.VQE` for :class:`~.MinimumEigenOptimizer`.
+    :class:`qiskit.algorithms.minimum_eigensolvers.NumPyMinimumEigensolver` also supersedes
+    :class:`qiskit.algorithms.NumPyMinimumEigensolver`.
+  - |
+    The :class:`~.WarmStartQAOAOptimizer` class takes the primitives-based QAOA
+    (:class:`qiskit.algorithms.minimum_eigensolvers.QAOA`) as ``qaoa`` argument.
+    The former algorithm :class:`qiskit.algorithms.QAOA`
+    is pending deprecation and will be deprecated and subsequently removed in future releases.
diff --git a/releasenotes/notes/grover-opt-primitive-de82d051d6cee2e4.yaml b/releasenotes/notes/grover-opt-primitive-de82d051d6cee2e4.yaml
new file mode 100644
index 000000000..7aea9b34a
--- /dev/null
+++ b/releasenotes/notes/grover-opt-primitive-de82d051d6cee2e4.yaml
@@ -0,0 +1,9 @@
+---
+features:
+  - |
+    The :class:`~.GroverOptimizer` class has a new keyword argument, ``sampler`` which is
+    used to run the algorithm using an instance of the :class:`qiskit.primitives.BaseSampler`
+    interface to calculate the results. This new argument supersedes the
+    the ``quantum_instance`` argument and accordingly, ``quantum_instance``
+    is pending deprecation and will be deprecated and subsequently removed in
+    future releases.
diff --git a/requirements.txt b/requirements.txt
index 1cddf9a79..9b49a7012 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,4 @@
-qiskit-terra>=0.20.0
+qiskit-terra>=0.22.*
 scipy>=1.4
 numpy>=1.17
 docplex>=2.21.207,!=2.24.231
diff --git a/test/algorithms/legacy/test_min_eigen_optimizer.py b/test/algorithms/legacy/test_min_eigen_optimizer.py
new file mode 100644
index 000000000..e6e99109f
--- /dev/null
+++ b/test/algorithms/legacy/test_min_eigen_optimizer.py
@@ -0,0 +1,406 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018, 2022.
+#
+# 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 Min Eigen Optimizer with the legacy MinimumEigensolver"""
+
+import unittest
+from test.optimization_test_case import QiskitOptimizationTestCase
+
+from test.runtime.fake_vqeruntime import FakeVQERuntimeProvider, FakeQAOARuntimeProvider
+
+import numpy as np
+from ddt import data, ddt
+from qiskit import BasicAer
+from qiskit.algorithms import QAOA, VQE, NumPyMinimumEigensolver
+from qiskit.algorithms.optimizers import COBYLA, SPSA
+from qiskit.circuit.library import TwoLocal
+from qiskit.providers.basicaer import QasmSimulatorPy
+from qiskit.utils import QuantumInstance, algorithm_globals
+import qiskit_optimization.optionals as _optionals
+from qiskit_optimization.algorithms import (
+    CplexOptimizer,
+    MinimumEigenOptimizer,
+    MinimumEigenOptimizationResult,
+)
+from qiskit_optimization.algorithms.optimization_algorithm import (
+    OptimizationResultStatus,
+)
+from qiskit_optimization.converters import (
+    InequalityToEquality,
+    IntegerToBinary,
+    LinearEqualityToPenalty,
+    MaximizeToMinimize,
+    QuadraticProgramToQubo,
+)
+from qiskit_optimization.problems import QuadraticProgram
+from qiskit_optimization.runtime import VQEClient, QAOAClient
+
+
+@ddt
+class TestMinEigenOptimizer(QiskitOptimizationTestCase):
+    """Min Eigen Optimizer Tests."""
+
+    def setUp(self):
+        super().setUp()
+
+        # setup minimum eigen solvers
+        self.min_eigen_solvers = {}
+
+        # exact eigen solver
+        with self.assertWarns(PendingDeprecationWarning):
+            self.min_eigen_solvers["exact"] = NumPyMinimumEigensolver()
+
+            # QAOA
+            optimizer = COBYLA()
+            self.min_eigen_solvers["qaoa"] = QAOA(optimizer=optimizer)
+            # simulators
+            self.sv_simulator = QuantumInstance(
+                BasicAer.get_backend("statevector_simulator"),
+                seed_simulator=123,
+                seed_transpiler=123,
+            )
+            self.qasm_simulator = QuantumInstance(
+                BasicAer.get_backend("qasm_simulator"),
+                seed_simulator=123,
+                seed_transpiler=123,
+            )
+        # test minimize
+        self.op_minimize = QuadraticProgram()
+        self.op_minimize.integer_var(0, 3, "x")
+        self.op_minimize.binary_var("y")
+        self.op_minimize.minimize(linear={"x": 1, "y": 2})
+        self.op_minimize.linear_constraint(linear={"x": 1, "y": 1}, sense=">=", rhs=1, name="xy")
+
+        # test maximize
+        self.op_maximize = QuadraticProgram()
+        self.op_maximize.integer_var(0, 3, "x")
+        self.op_maximize.binary_var("y")
+        self.op_maximize.maximize(linear={"x": 1, "y": 2})
+        self.op_maximize.linear_constraint(linear={"x": 1, "y": 1}, sense="<=", rhs=1, name="xy")
+
+        # test bit ordering
+        self.op_ordering = QuadraticProgram("bit ordering")
+        self.op_ordering.binary_var("x")
+        self.op_ordering.binary_var("y")
+        self.op_ordering.minimize(linear={"x": 1, "y": -2})
+
+    @data(
+        ("exact", None, "op_ip1.lp"),
+        ("qaoa", "statevector_simulator", "op_ip1.lp"),
+        ("qaoa", "qasm_simulator", "op_ip1.lp"),
+    )
+    @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
+    def test_min_eigen_optimizer(self, config):
+        """Min Eigen Optimizer Test"""
+        try:
+            # unpack configuration
+            min_eigen_solver_name, backend, filename = config
+
+            # get minimum eigen solver
+            min_eigen_solver = self.min_eigen_solvers[min_eigen_solver_name]
+            if backend:
+                min_eigen_solver.quantum_instance = BasicAer.get_backend(backend)
+
+            # construct minimum eigen optimizer
+            min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+
+            # load optimization problem
+            problem = QuadraticProgram()
+            lp_file = self.get_resource_path(filename, "algorithms/resources")
+            problem.read_from_lp_file(lp_file)
+
+            # solve problem with cplex
+            cplex = CplexOptimizer(cplex_parameters={"threads": 1, "randomseed": 1})
+            cplex_result = cplex.solve(problem)
+
+            # solve problem
+            with self.assertWarns(PendingDeprecationWarning):
+                result = min_eigen_optimizer.solve(problem)
+            self.assertIsNotNone(result)
+
+            # analyze results
+            self.assertAlmostEqual(cplex_result.fval, result.fval)
+
+            # check that eigensolver result is present
+            self.assertIsNotNone(result.min_eigen_solver_result)
+        except RuntimeError as ex:
+            self.fail(str(ex))
+
+    @data(
+        ("op_ip1.lp", -470, 12, OptimizationResultStatus.SUCCESS),
+        ("op_ip1.lp", np.inf, None, OptimizationResultStatus.FAILURE),
+    )
+    @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
+    def test_min_eigen_optimizer_with_filter(self, config):
+        """Min Eigen Optimizer Test"""
+        try:
+            # unpack configuration
+            filename, lowerbound, fval, status = config
+
+            # get minimum eigen solver
+            with self.assertWarns(PendingDeprecationWarning):
+                min_eigen_solver = NumPyMinimumEigensolver()
+
+            # set filter
+            # pylint: disable=unused-argument
+            def filter_criterion(x, v, aux):
+                return v > lowerbound
+
+            min_eigen_solver.filter_criterion = filter_criterion
+
+            # construct minimum eigen optimizer
+            min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+
+            # load optimization problem
+            problem = QuadraticProgram()
+            lp_file = self.get_resource_path(filename, "algorithms/resources")
+            problem.read_from_lp_file(lp_file)
+
+            # solve problem
+            with self.assertWarns(PendingDeprecationWarning):
+                result = min_eigen_optimizer.solve(problem)
+            self.assertIsNotNone(result)
+
+            # analyze results
+            self.assertAlmostEqual(fval, result.fval)
+            self.assertEqual(status, result.status)
+
+            # check that eigensolver result is present
+            self.assertIsNotNone(result.min_eigen_solver_result)
+        except RuntimeError as ex:
+            self.fail(str(ex))
+
+    def test_converter_list(self):
+        """Test converter list"""
+        op = QuadraticProgram()
+        op.integer_var(0, 3, "x")
+        op.binary_var("y")
+
+        op.maximize(linear={"x": 1, "y": 2})
+        op.linear_constraint(linear={"x": 1, "y": 1}, sense="LE", rhs=3, name="xy_leq")
+        with self.assertWarns(PendingDeprecationWarning):
+            min_eigen_solver = NumPyMinimumEigensolver()
+        # a single converter
+        qp2qubo = QuadraticProgramToQubo()
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver, converters=qp2qubo)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(op)
+        self.assertEqual(result.fval, 4)
+        # a list of converters
+        ineq2eq = InequalityToEquality()
+        int2bin = IntegerToBinary()
+        penalize = LinearEqualityToPenalty()
+        max2min = MaximizeToMinimize()
+        converters = [ineq2eq, int2bin, penalize, max2min]
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver, converters=converters)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(op)
+        self.assertEqual(result.fval, 4)
+        with self.assertRaises(TypeError):
+            invalid = [qp2qubo, "invalid converter"]
+            MinimumEigenOptimizer(min_eigen_solver, converters=invalid)
+
+    def test_samples_numpy_eigen_solver(self):
+        """Test samples for NumPyMinimumEigensolver"""
+        # test minimize
+        with self.assertWarns(PendingDeprecationWarning):
+            min_eigen_solver = NumPyMinimumEigensolver()
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(self.op_minimize)
+        opt_sol = 1
+        success = OptimizationResultStatus.SUCCESS
+        self.assertEqual(result.fval, opt_sol)
+        self.assertEqual(len(result.samples), 1)
+        np.testing.assert_array_almost_equal(result.samples[0].x, [1, 0])
+        self.assertAlmostEqual(result.samples[0].fval, opt_sol)
+        self.assertAlmostEqual(result.samples[0].probability, 1.0)
+        self.assertEqual(result.samples[0].status, success)
+        self.assertEqual(len(result.raw_samples), 1)
+        np.testing.assert_array_almost_equal(result.raw_samples[0].x, [1, 0, 0, 0, 0])
+        self.assertAlmostEqual(result.raw_samples[0].fval, opt_sol)
+        self.assertAlmostEqual(result.raw_samples[0].probability, 1.0)
+        self.assertEqual(result.raw_samples[0].status, success)
+        # test maximize
+        with self.assertWarns(PendingDeprecationWarning):
+            min_eigen_solver = NumPyMinimumEigensolver()
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(self.op_maximize)
+        opt_sol = 2
+        self.assertEqual(result.fval, opt_sol)
+        self.assertEqual(len(result.samples), 1)
+        np.testing.assert_array_almost_equal(result.samples[0].x, [0, 1])
+        self.assertAlmostEqual(result.samples[0].fval, opt_sol)
+        self.assertAlmostEqual(result.samples[0].probability, 1.0)
+        self.assertEqual(result.samples[0].status, success)
+        self.assertEqual(len(result.raw_samples), 1)
+        np.testing.assert_array_almost_equal(result.raw_samples[0].x, [0, 0, 1, 0])
+        # optimizer internally deals with minimization problem
+        self.assertAlmostEqual(
+            self.op_maximize.objective.sense.value * result.raw_samples[0].fval, opt_sol
+        )
+        self.assertAlmostEqual(result.raw_samples[0].probability, 1.0)
+        self.assertEqual(result.raw_samples[0].status, success)
+
+    @data("sv", "qasm")
+    def test_samples_qaoa(self, simulator):
+        """Test samples for QAOA"""
+        # test minimize
+        algorithm_globals.random_seed = 4
+        quantum_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(optimizer=COBYLA(), quantum_instance=quantum_instance, reps=2)
+        min_eigen_optimizer = MinimumEigenOptimizer(qaoa)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(self.op_minimize)
+        success = OptimizationResultStatus.SUCCESS
+        opt_sol = 1
+        self.assertAlmostEqual(sum(s.probability for s in result.samples), 1)
+        self.assertAlmostEqual(sum(s.probability for s in result.raw_samples), 1)
+        self.assertAlmostEqual(min(s.fval for s in result.samples), 0)
+        self.assertAlmostEqual(min(s.fval for s in result.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in result.raw_samples), opt_sol)
+        for sample in result.raw_samples:
+            self.assertEqual(sample.status, success)
+        np.testing.assert_array_almost_equal(result.x, [1, 0])
+        self.assertAlmostEqual(result.fval, result.samples[0].fval)
+        self.assertEqual(result.status, result.samples[0].status)
+        self.assertAlmostEqual(result.samples[0].fval, opt_sol)
+        self.assertEqual(result.samples[0].status, success)
+        np.testing.assert_array_almost_equal(result.raw_samples[0].x, [1, 0, 0, 0, 0])
+        self.assertAlmostEqual(result.raw_samples[0].fval, opt_sol)
+        self.assertEqual(result.raw_samples[0].status, success)
+        # test maximize
+        opt_sol = 2
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(optimizer=COBYLA(), quantum_instance=quantum_instance, reps=2)
+        min_eigen_optimizer = MinimumEigenOptimizer(qaoa)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(self.op_maximize)
+        self.assertAlmostEqual(sum(s.probability for s in result.samples), 1)
+        self.assertAlmostEqual(sum(s.probability for s in result.raw_samples), 1)
+        self.assertAlmostEqual(max(s.fval for s in result.samples), 5)
+        self.assertAlmostEqual(max(s.fval for s in result.samples if s.status == success), opt_sol)
+        # optimizer internally deals with minimization problem
+        self.assertAlmostEqual(
+            max(self.op_maximize.objective.sense.value * s.fval for s in result.raw_samples),
+            opt_sol,
+        )
+        for sample in result.raw_samples:
+            self.assertEqual(sample.status, success)
+        np.testing.assert_array_almost_equal(result.x, [0, 1])
+        self.assertEqual(result.fval, opt_sol)
+        self.assertEqual(result.status, success)
+        np.testing.assert_array_almost_equal(result.samples[0].x, [0, 1])
+        self.assertAlmostEqual(result.samples[0].fval, opt_sol)
+        self.assertEqual(result.samples[0].status, success)
+        np.testing.assert_array_almost_equal(result.raw_samples[0].x, [0, 0, 1, 0])
+        # optimizer internally deals with minimization problem
+        self.assertAlmostEqual(
+            self.op_maximize.objective.sense.value * result.raw_samples[0].fval, opt_sol
+        )
+        self.assertEqual(result.raw_samples[0].status, success)
+        # test bit ordering
+        opt_sol = -2
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(optimizer=COBYLA(), quantum_instance=quantum_instance, reps=2)
+        min_eigen_optimizer = MinimumEigenOptimizer(qaoa)
+        with self.assertWarns(PendingDeprecationWarning):
+            result = min_eigen_optimizer.solve(self.op_ordering)
+        self.assertEqual(result.fval, opt_sol)
+        np.testing.assert_array_almost_equal(result.x, [0, 1])
+        self.assertEqual(result.status, success)
+        result.raw_samples.sort(key=lambda x: x.probability, reverse=True)
+        np.testing.assert_array_almost_equal(result.x, result.raw_samples[0].x)
+        self.assertAlmostEqual(sum(s.probability for s in result.samples), 1, delta=1e-5)
+        self.assertAlmostEqual(sum(s.probability for s in result.raw_samples), 1, delta=1e-5)
+        self.assertAlmostEqual(min(s.fval for s in result.samples), -2)
+        self.assertAlmostEqual(min(s.fval for s in result.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in result.raw_samples), opt_sol)
+        for sample in result.raw_samples:
+            self.assertEqual(sample.status, success)
+        np.testing.assert_array_almost_equal(result.samples[0].x, [0, 1])
+        self.assertAlmostEqual(result.samples[0].fval, opt_sol)
+        self.assertEqual(result.samples[0].status, success)
+        np.testing.assert_array_almost_equal(result.raw_samples[0].x, [0, 1])
+        self.assertAlmostEqual(result.raw_samples[0].fval, opt_sol)
+        self.assertEqual(result.raw_samples[0].status, success)
+
+    @data("sv", "qasm")
+    def test_samples_vqe(self, simulator):
+        """Test samples for VQE"""
+        # test minimize
+        algorithm_globals.random_seed = 1
+        quantum_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
+        opt_sol = -2
+        success = OptimizationResultStatus.SUCCESS
+        optimizer = SPSA(maxiter=100)
+        ry_ansatz = TwoLocal(5, "ry", "cz", reps=3, entanglement="full")
+        with self.assertWarns(PendingDeprecationWarning):
+            vqe_mes = VQE(ry_ansatz, optimizer=optimizer, quantum_instance=quantum_instance)
+        vqe = MinimumEigenOptimizer(vqe_mes)
+        with self.assertWarns(PendingDeprecationWarning):
+            results = vqe.solve(self.op_ordering)
+        self.assertEqual(results.fval, opt_sol)
+        np.testing.assert_array_almost_equal(results.x, [0, 1])
+        self.assertEqual(results.status, success)
+        results.raw_samples.sort(key=lambda x: x.probability, reverse=True)
+        np.testing.assert_array_almost_equal(results.x, results.raw_samples[0].x)
+        self.assertAlmostEqual(sum(s.probability for s in results.samples), 1, delta=1e-5)
+        self.assertAlmostEqual(sum(s.probability for s in results.raw_samples), 1, delta=1e-5)
+        self.assertAlmostEqual(min(s.fval for s in results.samples), -2)
+        self.assertAlmostEqual(min(s.fval for s in results.samples if s.status == success), opt_sol)
+        self.assertAlmostEqual(min(s.fval for s in results.raw_samples), opt_sol)
+        for sample in results.raw_samples:
+            self.assertEqual(sample.status, success)
+        np.testing.assert_array_almost_equal(results.samples[0].x, [0, 1])
+        self.assertAlmostEqual(results.samples[0].fval, opt_sol)
+        self.assertEqual(results.samples[0].status, success)
+        np.testing.assert_array_almost_equal(results.raw_samples[0].x, [0, 1])
+        self.assertAlmostEqual(results.raw_samples[0].fval, opt_sol)
+        self.assertEqual(results.raw_samples[0].status, success)
+
+    @data("vqe", "qaoa")
+    def test_runtime(self, subroutine):
+        """Test vqe and qaoa runtime"""
+        optimizer = {"name": "SPSA", "maxiter": 100}
+        backend = QasmSimulatorPy()
+
+        if subroutine == "vqe":
+            ry_ansatz = TwoLocal(5, "ry", "cz", reps=3, entanglement="full")
+            initial_point = np.random.default_rng(42).random(ry_ansatz.num_parameters)
+            solver = VQEClient(
+                ansatz=ry_ansatz,
+                optimizer=optimizer,
+                initial_point=initial_point,
+                backend=backend,
+                provider=FakeVQERuntimeProvider(),
+            )
+        else:
+            reps = 2
+            initial_point = np.random.default_rng(42).random(2 * reps)
+            solver = QAOAClient(
+                optimizer=optimizer,
+                reps=reps,
+                initial_point=initial_point,
+                backend=backend,
+                provider=FakeQAOARuntimeProvider(),
+            )
+
+        opt = MinimumEigenOptimizer(solver)
+        result = opt.solve(self.op_ordering)
+        self.assertIsInstance(result, MinimumEigenOptimizationResult)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/algorithms/legacy/test_recursive_optimization.py b/test/algorithms/legacy/test_recursive_optimization.py
new file mode 100644
index 000000000..675f8419a
--- /dev/null
+++ b/test/algorithms/legacy/test_recursive_optimization.py
@@ -0,0 +1,212 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2018, 2022.
+#
+# 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 Recursive Min Eigen Optimizer with legacy MinimumEigensolver."""
+
+import unittest
+from test import QiskitOptimizationTestCase
+
+import numpy as np
+
+from qiskit import BasicAer
+from qiskit.utils import algorithm_globals, QuantumInstance
+
+from qiskit.algorithms import NumPyMinimumEigensolver, QAOA
+
+import qiskit_optimization.optionals as _optionals
+from qiskit_optimization.algorithms import (
+    MinimumEigenOptimizer,
+    CplexOptimizer,
+    RecursiveMinimumEigenOptimizer,
+    WarmStartQAOAOptimizer,
+    SlsqpOptimizer,
+)
+from qiskit_optimization.algorithms.recursive_minimum_eigen_optimizer import (
+    IntermediateResult,
+)
+from qiskit_optimization.problems import QuadraticProgram
+from qiskit_optimization.converters import (
+    IntegerToBinary,
+    InequalityToEquality,
+    LinearEqualityToPenalty,
+    QuadraticProgramToQubo,
+)
+
+
+class TestRecursiveMinEigenOptimizer(QiskitOptimizationTestCase):
+    """Recursive Min Eigen Optimizer Tests."""
+
+    @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
+    def test_recursive_min_eigen_optimizer(self):
+        """Test the recursive minimum eigen optimizer."""
+        filename = "op_ip1.lp"
+        # get minimum eigen solver
+        with self.assertWarns(PendingDeprecationWarning):
+            min_eigen_solver = NumPyMinimumEigensolver()
+
+        # construct minimum eigen optimizer
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(
+            min_eigen_optimizer, min_num_vars=4
+        )
+
+        # load optimization problem
+        problem = QuadraticProgram()
+        lp_file = self.get_resource_path(filename, "algorithms/resources")
+        problem.read_from_lp_file(lp_file)
+
+        # solve problem with cplex
+        cplex = CplexOptimizer()
+        cplex_result = cplex.solve(problem)
+
+        # solve problem
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(problem)
+
+        # analyze results
+        np.testing.assert_array_almost_equal(cplex_result.x, result.x, 4)
+        self.assertAlmostEqual(cplex_result.fval, result.fval)
+
+    @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
+    def test_recursive_history(self):
+        """Tests different options for history."""
+        filename = "op_ip1.lp"
+        # load optimization problem
+        problem = QuadraticProgram()
+        lp_file = self.get_resource_path(filename, "algorithms/resources")
+        problem.read_from_lp_file(lp_file)
+
+        # get minimum eigen solver
+        with self.assertWarns(PendingDeprecationWarning):
+            min_eigen_solver = NumPyMinimumEigensolver()
+
+        # construct minimum eigen optimizer
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+
+        # no history
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(
+            min_eigen_optimizer,
+            min_num_vars=4,
+            history=IntermediateResult.NO_ITERATIONS,
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(problem)
+        self.assertIsNotNone(result.replacements)
+        self.assertIsNotNone(result.history)
+        self.assertIsNotNone(result.history[0])
+        self.assertEqual(len(result.history[0]), 0)
+        self.assertIsNone(result.history[1])
+
+        # only last iteration in the history
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(
+            min_eigen_optimizer,
+            min_num_vars=4,
+            history=IntermediateResult.LAST_ITERATION,
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(problem)
+        self.assertIsNotNone(result.replacements)
+        self.assertIsNotNone(result.history)
+        self.assertIsNotNone(result.history[0])
+        self.assertEqual(len(result.history[0]), 0)
+        self.assertIsNotNone(result.history[1])
+
+        # full history
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(
+            min_eigen_optimizer,
+            min_num_vars=4,
+            history=IntermediateResult.ALL_ITERATIONS,
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(problem)
+        self.assertIsNotNone(result.replacements)
+        self.assertIsNotNone(result.history)
+        self.assertIsNotNone(result.history[0])
+        self.assertGreater(len(result.history[0]), 1)
+        self.assertIsNotNone(result.history[1])
+
+    @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
+    def test_recursive_warm_qaoa(self):
+        """Test the recursive optimizer with warm start qaoa."""
+        seed = 1234
+        algorithm_globals.random_seed = seed
+        backend = BasicAer.get_backend("statevector_simulator")
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(
+                quantum_instance=QuantumInstance(
+                    backend=backend, seed_simulator=seed, seed_transpiler=seed
+                ),
+                reps=1,
+            )
+        warm_qaoa = WarmStartQAOAOptimizer(
+            pre_solver=SlsqpOptimizer(), relax_for_pre_solver=True, qaoa=qaoa
+        )
+
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(warm_qaoa, min_num_vars=4)
+
+        # load optimization problem
+        problem = QuadraticProgram()
+        lp_file = self.get_resource_path("op_ip1.lp", "algorithms/resources")
+        problem.read_from_lp_file(lp_file)
+
+        # solve problem with cplex
+        cplex = CplexOptimizer(cplex_parameters={"threads": 1, "randomseed": 1})
+        cplex_result = cplex.solve(problem)
+
+        # solve problem
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(problem)
+
+        # analyze results
+        np.testing.assert_array_almost_equal(cplex_result.x, result.x, 4)
+        self.assertAlmostEqual(cplex_result.fval, result.fval)
+
+    def test_converter_list(self):
+        """Test converter list"""
+        op = QuadraticProgram()
+        op.integer_var(0, 3, "x")
+        op.binary_var("y")
+
+        op.maximize(linear={"x": 1, "y": 2})
+        op.linear_constraint(linear={"y": 1, "x": 1}, sense="LE", rhs=3, name="xy_leq")
+
+        # construct minimum eigen optimizer
+        with self.assertWarns(PendingDeprecationWarning):
+            min_eigen_solver = NumPyMinimumEigensolver()
+        min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
+        # a single converter
+        qp2qubo = QuadraticProgramToQubo()
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(
+            min_eigen_optimizer, min_num_vars=2, converters=qp2qubo
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(op)
+        self.assertEqual(result.fval, 4)
+        # a list of converters
+        ineq2eq = InequalityToEquality()
+        int2bin = IntegerToBinary()
+        penalize = LinearEqualityToPenalty()
+        converters = [ineq2eq, int2bin, penalize]
+        recursive_min_eigen_optimizer = RecursiveMinimumEigenOptimizer(
+            min_eigen_optimizer, min_num_vars=2, converters=converters
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result = recursive_min_eigen_optimizer.solve(op)
+        self.assertEqual(result.fval, 4)
+        # invalid converters
+        with self.assertRaises(TypeError):
+            invalid = [qp2qubo, "invalid converter"]
+            RecursiveMinimumEigenOptimizer(min_eigen_optimizer, min_num_vars=2, converters=invalid)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/algorithms/legacy/test_warm_start_qaoa.py b/test/algorithms/legacy/test_warm_start_qaoa.py
new file mode 100644
index 000000000..3893f0e23
--- /dev/null
+++ b/test/algorithms/legacy/test_warm_start_qaoa.py
@@ -0,0 +1,140 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021, 2022.
+#
+# 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 warm start QAOA optimizer with legacy QAOA. """
+
+import unittest
+from test import QiskitOptimizationTestCase
+
+import numpy as np
+
+from docplex.mp.model import Model
+from qiskit import BasicAer
+from qiskit.algorithms import QAOA
+
+import qiskit_optimization.optionals as _optionals
+from qiskit_optimization.algorithms import SlsqpOptimizer
+from qiskit_optimization.algorithms.goemans_williamson_optimizer import (
+    GoemansWilliamsonOptimizer,
+)
+from qiskit_optimization.algorithms.warm_start_qaoa_optimizer import (
+    MeanAggregator,
+    WarmStartQAOAOptimizer,
+)
+from qiskit_optimization.applications.max_cut import Maxcut
+from qiskit_optimization.translators import from_docplex_mp
+
+
+class TestWarmStartQAOAOptimizer(QiskitOptimizationTestCase):
+    """Tests for the warm start QAOA optimizer."""
+
+    @unittest.skipIf(not _optionals.HAS_CVXPY, "CVXPY not available.")
+    def test_max_cut(self):
+        """Basic test on the max cut problem."""
+        graph = np.array(
+            [
+                [0.0, 1.0, 2.0, 0.0],
+                [1.0, 0.0, 1.0, 0.0],
+                [2.0, 1.0, 0.0, 1.0],
+                [0.0, 0.0, 1.0, 0.0],
+            ]
+        )
+
+        presolver = GoemansWilliamsonOptimizer(num_cuts=10)
+        problem = Maxcut(graph).to_quadratic_program()
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(quantum_instance=backend, reps=1)
+        aggregator = MeanAggregator()
+        optimizer = WarmStartQAOAOptimizer(
+            pre_solver=presolver,
+            relax_for_pre_solver=False,
+            qaoa=qaoa,
+            epsilon=0.25,
+            num_initial_solutions=10,
+            aggregator=aggregator,
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result_warm = optimizer.solve(problem)
+
+        self.assertIsNotNone(result_warm)
+        self.assertIsNotNone(result_warm.x)
+        np.testing.assert_almost_equal([0, 0, 1, 0], result_warm.x, 3)
+        self.assertIsNotNone(result_warm.fval)
+        np.testing.assert_almost_equal(4, result_warm.fval, 3)
+
+    def test_constrained_binary(self):
+        """Constrained binary optimization problem."""
+        model = Model()
+        v = model.binary_var(name="v")
+        w = model.binary_var(name="w")
+        # pylint:disable=invalid-name
+        t = model.binary_var(name="t")
+
+        model.minimize(v + w + t)
+        model.add_constraint(2 * v + 10 * w + t <= 3, "cons1")
+        model.add_constraint(v + w + t >= 2, "cons2")
+
+        problem = from_docplex_mp(model)
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(quantum_instance=backend, reps=1)
+        aggregator = MeanAggregator()
+        optimizer = WarmStartQAOAOptimizer(
+            pre_solver=SlsqpOptimizer(),
+            relax_for_pre_solver=True,
+            qaoa=qaoa,
+            epsilon=0.25,
+            aggregator=aggregator,
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result_warm = optimizer.solve(problem)
+
+        self.assertIsNotNone(result_warm)
+        self.assertIsNotNone(result_warm.x)
+        np.testing.assert_almost_equal([1, 0, 1], result_warm.x, 3)
+        self.assertIsNotNone(result_warm.fval)
+        np.testing.assert_almost_equal(2, result_warm.fval, 3)
+
+    def test_simple_qubo(self):
+        """Test on a simple QUBO problem."""
+        model = Model()
+        # pylint:disable=invalid-name
+        u = model.binary_var(name="u")
+        v = model.binary_var(name="v")
+
+        model.minimize((u - v + 2) ** 2)
+        problem = from_docplex_mp(model)
+
+        backend = BasicAer.get_backend("statevector_simulator")
+        with self.assertWarns(PendingDeprecationWarning):
+            qaoa = QAOA(quantum_instance=backend, reps=1)
+        optimizer = WarmStartQAOAOptimizer(
+            pre_solver=SlsqpOptimizer(),
+            relax_for_pre_solver=True,
+            qaoa=qaoa,
+            epsilon=0.25,
+        )
+        with self.assertWarns(PendingDeprecationWarning):
+            result_warm = optimizer.solve(problem)
+
+        self.assertIsNotNone(result_warm)
+        self.assertIsNotNone(result_warm.x)
+        np.testing.assert_almost_equal([0, 1], result_warm.x, 3)
+        self.assertIsNotNone(result_warm.fval)
+        np.testing.assert_almost_equal(1, result_warm.fval, 3)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/algorithms/test_grover_optimizer.py b/test/algorithms/test_grover_optimizer.py
index a029590fa..a39c4f202 100644
--- a/test/algorithms/test_grover_optimizer.py
+++ b/test/algorithms/test_grover_optimizer.py
@@ -19,7 +19,8 @@
 from ddt import data, ddt
 from docplex.mp.model import Model
 from qiskit.utils import QuantumInstance, algorithm_globals, optionals
-from qiskit.algorithms import NumPyMinimumEigensolver
+from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
+from qiskit.primitives import Sampler
 from qiskit_optimization.algorithms import (
     GroverOptimizer,
     MinimumEigenOptimizer,
@@ -58,6 +59,35 @@ def setUp(self):
         )
         self.n_iter = 8
 
+    def _prepare_grover_optimizer(
+        self, num_value_qubits, num_iterations, simulator, converters=None
+    ):
+        """Prepare GroverOptimizer."""
+        if simulator == "statevector":
+            with self.assertWarns(PendingDeprecationWarning):
+                grover_optimizer = GroverOptimizer(
+                    num_value_qubits=num_value_qubits,
+                    num_iterations=num_iterations,
+                    converters=converters,
+                    quantum_instance=self.sv_simulator,
+                )
+        elif simulator == "qasm":
+            with self.assertWarns(PendingDeprecationWarning):
+                grover_optimizer = GroverOptimizer(
+                    num_value_qubits=num_value_qubits,
+                    num_iterations=num_iterations,
+                    converters=converters,
+                    quantum_instance=self.qasm_simulator,
+                )
+        else:
+            grover_optimizer = GroverOptimizer(
+                num_value_qubits=num_value_qubits,
+                num_iterations=num_iterations,
+                converters=converters,
+                sampler=Sampler(),
+            )
+        return grover_optimizer
+
     def validate_results(self, problem, results):
         """Validate the results object returned by GroverOptimizer."""
         # Get expected value.
@@ -71,7 +101,8 @@ def validate_results(self, problem, results):
             results.fval, problem.objective.sense.value * results.intermediate_fval
         )
 
-    def test_qubo_gas_int_zero(self):
+    @data("statevector", "qasm", "sampler")
+    def test_qubo_gas_int_zero(self, simulator):
         """Test for when the answer is zero."""
 
         # Input.
@@ -82,13 +113,16 @@ def test_qubo_gas_int_zero(self):
         op = from_docplex_mp(model)
 
         # Will not find a negative, should return 0.
-        gmf = GroverOptimizer(1, num_iterations=1, quantum_instance=self.sv_simulator)
-        results = gmf.solve(op)
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=1, num_iterations=1, simulator=simulator
+        )
+        results = grover_optimizer.solve(op)
         np.testing.assert_array_almost_equal(results.x, [0, 0])
         self.assertEqual(results.fval, 0.0)
         self.assertAlmostEqual(results.fval, results.intermediate_fval)
 
-    def test_qubo_gas_int_simple(self):
+    @data("statevector", "qasm", "sampler")
+    def test_qubo_gas_int_simple(self, simulator):
         """Test for simple case, with 2 linear coeffs and no quadratic coeffs or constants."""
 
         # Input.
@@ -99,15 +133,18 @@ def test_qubo_gas_int_simple(self):
         op = from_docplex_mp(model)
 
         # Get the optimum key and value.
-        gmf = GroverOptimizer(4, num_iterations=self.n_iter, quantum_instance=self.sv_simulator)
-        results = gmf.solve(op)
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=4, num_iterations=self.n_iter, simulator=simulator
+        )
+        results = grover_optimizer.solve(op)
         self.validate_results(op, results)
 
         self.assertIsNotNone(results.operation_counts)
         self.assertEqual(results.n_input_qubits, 2)
         self.assertEqual(results.n_output_qubits, 4)
 
-    def test_qubo_gas_int_simple_maximize(self):
+    @data("statevector", "qasm", "sampler")
+    def test_qubo_gas_int_simple_maximize(self, simulator):
         """Test for simple case, but with maximization."""
 
         # Input.
@@ -118,11 +155,13 @@ def test_qubo_gas_int_simple_maximize(self):
         op = from_docplex_mp(model)
 
         # Get the optimum key and value.
-        gmf = GroverOptimizer(4, num_iterations=self.n_iter, quantum_instance=self.sv_simulator)
-        results = gmf.solve(op)
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=4, num_iterations=self.n_iter, simulator=simulator
+        )
+        results = grover_optimizer.solve(op)
         self.validate_results(op, results)
 
-    @data("sv", "qasm")
+    @data("statevector", "qasm", "sampler")
     def test_qubo_gas_int_paper_example(self, simulator):
         """
         Test the example from https://arxiv.org/abs/1912.04088 using the state vector simulator
@@ -138,12 +177,14 @@ def test_qubo_gas_int_paper_example(self, simulator):
         op = from_docplex_mp(model)
 
         # Get the optimum key and value.
-        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
-        gmf = GroverOptimizer(6, num_iterations=self.n_iter, quantum_instance=q_instance)
-        results = gmf.solve(op)
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=6, num_iterations=self.n_iter, simulator=simulator
+        )
+        results = grover_optimizer.solve(op)
         self.validate_results(op, results)
 
-    def test_converter_list(self):
+    @data("statevector", "qasm", "sampler")
+    def test_converter_list(self, simulator):
         """Test converters list"""
         # Input.
 
@@ -156,13 +197,10 @@ def test_converter_list(self):
         # Get the optimum key and value.
         # a single converter.
         qp2qubo = QuadraticProgramToQubo()
-        gmf = GroverOptimizer(
-            4,
-            num_iterations=self.n_iter,
-            quantum_instance=self.sv_simulator,
-            converters=qp2qubo,
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=4, num_iterations=self.n_iter, simulator=simulator
         )
-        results = gmf.solve(op)
+        results = grover_optimizer.solve(op)
         self.validate_results(op, results)
 
         # a list of converters
@@ -171,25 +209,22 @@ def test_converter_list(self):
         penalize = LinearEqualityToPenalty()
         max2min = MaximizeToMinimize()
         converters = [ineq2eq, int2bin, penalize, max2min]
-        gmf = GroverOptimizer(
-            4,
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=4,
             num_iterations=self.n_iter,
-            quantum_instance=self.sv_simulator,
+            simulator=simulator,
             converters=converters,
         )
-        results = gmf.solve(op)
+        results = grover_optimizer.solve(op)
         self.validate_results(op, results)
         # invalid converters
         with self.assertRaises(TypeError):
             invalid = [qp2qubo, "invalid converter"]
-            GroverOptimizer(
-                4,
-                num_iterations=self.n_iter,
-                quantum_instance=self.sv_simulator,
-                converters=invalid,
+            grover_optimizer = self._prepare_grover_optimizer(
+                4, num_iterations=self.n_iter, simulator=simulator, converters=invalid
             )
 
-    @data("sv", "qasm")
+    @data("statevector", "qasm", "sampler")
     def test_samples_and_raw_samples(self, simulator):
         """Test samples and raw_samples"""
         algorithm_globals.random_seed = 2
@@ -198,9 +233,8 @@ def test_samples_and_raw_samples(self, simulator):
         op.binary_var("y")
         op.minimize(linear={"x": 1, "y": 2})
         op.linear_constraint(linear={"x": 1, "y": 1}, sense=">=", rhs=1, name="xy")
-        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
-        grover_optimizer = GroverOptimizer(
-            8, num_iterations=self.n_iter, quantum_instance=q_instance
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=8, num_iterations=self.n_iter, simulator=simulator
         )
         opt_sol = 1
         success = OptimizationResultStatus.SUCCESS
@@ -221,12 +255,11 @@ def test_samples_and_raw_samples(self, simulator):
         self.assertEqual(results.status, results.raw_samples[0].status)
         np.testing.assert_array_almost_equal([1, 0, 0, 0, 0], results.raw_samples[0].x)
 
-    @data("sv", "qasm")
+    @data("statevector", "qasm", "sampler")
     def test_bit_ordering(self, simulator):
         """Test bit ordering"""
         # test minimize
         algorithm_globals.random_seed = 2
-        q_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
         mdl = Model("docplex model")
         x = mdl.binary_var("x")
         y = mdl.binary_var("y")
@@ -234,8 +267,8 @@ def test_bit_ordering(self, simulator):
         op = from_docplex_mp(mdl)
         opt_sol = -2
         success = OptimizationResultStatus.SUCCESS
-        grover_optimizer = GroverOptimizer(
-            3, num_iterations=self.n_iter, quantum_instance=q_instance
+        grover_optimizer = self._prepare_grover_optimizer(
+            num_value_qubits=3, num_iterations=self.n_iter, simulator=simulator
         )
         results = grover_optimizer.solve(op)
         self.assertEqual(results.fval, opt_sol)
diff --git a/test/algorithms/test_min_eigen_optimizer.py b/test/algorithms/test_min_eigen_optimizer.py
index 257fe7fe4..a5b518e78 100644
--- a/test/algorithms/test_min_eigen_optimizer.py
+++ b/test/algorithms/test_min_eigen_optimizer.py
@@ -10,30 +10,28 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-""" Test Min Eigen Optimizer """
+""" Test Min Eigen Optimizer with the primitive-based minimum eigensolver """
 
 import unittest
 from test.optimization_test_case import QiskitOptimizationTestCase
-
-from test.runtime.fake_vqeruntime import FakeVQERuntimeProvider, FakeQAOARuntimeProvider
+from test.runtime.fake_vqeruntime import FakeQAOARuntimeProvider, FakeVQERuntimeProvider
 
 import numpy as np
-from ddt import data, ddt
-from qiskit import BasicAer
-from qiskit.algorithms import QAOA, VQE, NumPyMinimumEigensolver
+from ddt import data, ddt, unpack
+from qiskit.algorithms.minimum_eigensolvers import QAOA, NumPyMinimumEigensolver, SamplingVQE
 from qiskit.algorithms.optimizers import COBYLA, SPSA
 from qiskit.circuit.library import TwoLocal
+from qiskit.primitives import Sampler
 from qiskit.providers.basicaer import QasmSimulatorPy
-from qiskit.utils import QuantumInstance, algorithm_globals
+from qiskit.utils import algorithm_globals
+
 import qiskit_optimization.optionals as _optionals
 from qiskit_optimization.algorithms import (
     CplexOptimizer,
-    MinimumEigenOptimizer,
     MinimumEigenOptimizationResult,
+    MinimumEigenOptimizer,
 )
-from qiskit_optimization.algorithms.optimization_algorithm import (
-    OptimizationResultStatus,
-)
+from qiskit_optimization.algorithms.optimization_algorithm import OptimizationResultStatus
 from qiskit_optimization.converters import (
     InequalityToEquality,
     IntegerToBinary,
@@ -42,7 +40,7 @@
     QuadraticProgramToQubo,
 )
 from qiskit_optimization.problems import QuadraticProgram
-from qiskit_optimization.runtime import VQEProgram, QAOAProgram
+from qiskit_optimization.runtime import QAOAClient, VQEClient
 
 
 @ddt
@@ -52,26 +50,14 @@ class TestMinEigenOptimizer(QiskitOptimizationTestCase):
     def setUp(self):
         super().setUp()
 
+        self._seed = 123
+
         # setup minimum eigen solvers
-        self.min_eigen_solvers = {}
-
-        # exact eigen solver
-        self.min_eigen_solvers["exact"] = NumPyMinimumEigensolver()
-
-        # QAOA
-        optimizer = COBYLA()
-        self.min_eigen_solvers["qaoa"] = QAOA(optimizer=optimizer)
-        # simulators
-        self.sv_simulator = QuantumInstance(
-            BasicAer.get_backend("statevector_simulator"),
-            seed_simulator=123,
-            seed_transpiler=123,
-        )
-        self.qasm_simulator = QuantumInstance(
-            BasicAer.get_backend("qasm_simulator"),
-            seed_simulator=123,
-            seed_transpiler=123,
-        )
+        self.min_eigen_solvers = {
+            "exact": NumPyMinimumEigensolver(),
+            "qaoa": QAOA(sampler=Sampler(), optimizer=COBYLA()),
+        }
+
         # test minimize
         self.op_minimize = QuadraticProgram()
         self.op_minimize.integer_var(0, 3, "x")
@@ -94,20 +80,19 @@ def setUp(self):
 
     @data(
         ("exact", None, "op_ip1.lp"),
-        ("qaoa", "statevector_simulator", "op_ip1.lp"),
-        ("qaoa", "qasm_simulator", "op_ip1.lp"),
+        ("qaoa", None, "op_ip1.lp"),
+        ("qaoa", 10000, "op_ip1.lp"),
     )
+    @unpack
     @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
-    def test_min_eigen_optimizer(self, config):
+    def test_min_eigen_optimizer(self, min_eigen_solver_name, shots, filename):
         """Min Eigen Optimizer Test"""
         try:
-            # unpack configuration
-            min_eigen_solver_name, backend, filename = config
-
             # get minimum eigen solver
             min_eigen_solver = self.min_eigen_solvers[min_eigen_solver_name]
-            if backend:
-                min_eigen_solver.quantum_instance = BasicAer.get_backend(backend)
+            if min_eigen_solver_name == "qaoa":
+                min_eigen_solver.sampler.options.shots = shots
+                min_eigen_solver.sampler.options.seed = self._seed
 
             # construct minimum eigen optimizer
             min_eigen_optimizer = MinimumEigenOptimizer(min_eigen_solver)
@@ -137,13 +122,11 @@ def test_min_eigen_optimizer(self, config):
         ("op_ip1.lp", -470, 12, OptimizationResultStatus.SUCCESS),
         ("op_ip1.lp", np.inf, None, OptimizationResultStatus.FAILURE),
     )
+    @unpack
     @unittest.skipIf(not _optionals.HAS_CPLEX, "CPLEX not available.")
-    def test_min_eigen_optimizer_with_filter(self, config):
+    def test_min_eigen_optimizer_with_filter(self, filename, lowerbound, fval, status):
         """Min Eigen Optimizer Test"""
         try:
-            # unpack configuration
-            filename, lowerbound, fval, status = config
-
             # get minimum eigen solver
             min_eigen_solver = NumPyMinimumEigensolver()
 
@@ -241,13 +224,12 @@ def test_samples_numpy_eigen_solver(self):
         self.assertAlmostEqual(result.raw_samples[0].probability, 1.0)
         self.assertEqual(result.raw_samples[0].status, success)
 
-    @data("sv", "qasm")
-    def test_samples_qaoa(self, simulator):
+    def test_samples_qaoa(self):
         """Test samples for QAOA"""
         # test minimize
         algorithm_globals.random_seed = 4
-        quantum_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
-        qaoa = QAOA(optimizer=COBYLA(), quantum_instance=quantum_instance, reps=2)
+        sampler = Sampler()
+        qaoa = QAOA(sampler=sampler, optimizer=COBYLA(), reps=2)
         min_eigen_optimizer = MinimumEigenOptimizer(qaoa)
         result = min_eigen_optimizer.solve(self.op_minimize)
         success = OptimizationResultStatus.SUCCESS
@@ -269,7 +251,7 @@ def test_samples_qaoa(self, simulator):
         self.assertEqual(result.raw_samples[0].status, success)
         # test maximize
         opt_sol = 2
-        qaoa = QAOA(optimizer=COBYLA(), quantum_instance=quantum_instance, reps=2)
+        qaoa = QAOA(sampler=sampler, optimizer=COBYLA(), reps=2)
         min_eigen_optimizer = MinimumEigenOptimizer(qaoa)
         result = min_eigen_optimizer.solve(self.op_maximize)
         self.assertAlmostEqual(sum(s.probability for s in result.samples), 1)
@@ -297,7 +279,7 @@ def test_samples_qaoa(self, simulator):
         self.assertEqual(result.raw_samples[0].status, success)
         # test bit ordering
         opt_sol = -2
-        qaoa = QAOA(optimizer=COBYLA(), quantum_instance=quantum_instance, reps=2)
+        qaoa = QAOA(sampler=sampler, optimizer=COBYLA(), reps=2)
         min_eigen_optimizer = MinimumEigenOptimizer(qaoa)
         result = min_eigen_optimizer.solve(self.op_ordering)
         self.assertEqual(result.fval, opt_sol)
@@ -319,17 +301,16 @@ def test_samples_qaoa(self, simulator):
         self.assertAlmostEqual(result.raw_samples[0].fval, opt_sol)
         self.assertEqual(result.raw_samples[0].status, success)
 
-    @data("sv", "qasm")
-    def test_samples_vqe(self, simulator):
+    def test_samples_vqe(self):
         """Test samples for VQE"""
         # test minimize
         algorithm_globals.random_seed = 1
-        quantum_instance = self.sv_simulator if simulator == "sv" else self.qasm_simulator
         opt_sol = -2
         success = OptimizationResultStatus.SUCCESS
         optimizer = SPSA(maxiter=100)
         ry_ansatz = TwoLocal(5, "ry", "cz", reps=3, entanglement="full")
-        vqe_mes = VQE(ry_ansatz, optimizer=optimizer, quantum_instance=quantum_instance)
+        sampler = Sampler()
+        vqe_mes = SamplingVQE(sampler, ry_ansatz, optimizer=optimizer)
         vqe = MinimumEigenOptimizer(vqe_mes)
         results = vqe.solve(self.op_ordering)
         self.assertEqual(results.fval, opt_sol)
@@ -360,7 +341,7 @@ def test_runtime(self, subroutine):
         if subroutine == "vqe":
             ry_ansatz = TwoLocal(5, "ry", "cz", reps=3, entanglement="full")
             initial_point = np.random.default_rng(42).random(ry_ansatz.num_parameters)
-            solver = VQEProgram(
+            solver = VQEClient(
                 ansatz=ry_ansatz,
                 optimizer=optimizer,
                 initial_point=initial_point,
@@ -370,7 +351,7 @@ def test_runtime(self, subroutine):
         else:
             reps = 2
             initial_point = np.random.default_rng(42).random(2 * reps)
-            solver = QAOAProgram(
+            solver = QAOAClient(
                 optimizer=optimizer,
                 reps=reps,
                 initial_point=initial_point,
diff --git a/test/algorithms/test_recursive_optimization.py b/test/algorithms/test_recursive_optimization.py
old mode 100755
new mode 100644
index 53e206b7b..ec8c6d9cc
--- a/test/algorithms/test_recursive_optimization.py
+++ b/test/algorithms/test_recursive_optimization.py
@@ -10,36 +10,33 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Test Recursive Min Eigen Optimizer."""
+"""Test Recursive Min Eigen Optimizer with the primitive-based minimum eigensolver."""
 
 import unittest
 from test import QiskitOptimizationTestCase
 
 import numpy as np
-
-from qiskit import BasicAer
-from qiskit.utils import algorithm_globals, QuantumInstance
-
-from qiskit.algorithms import NumPyMinimumEigensolver, QAOA
+from qiskit.algorithms.minimum_eigensolvers import QAOA, NumPyMinimumEigensolver
+from qiskit.algorithms.optimizers import SLSQP
+from qiskit.primitives import Sampler
+from qiskit.utils import algorithm_globals
 
 import qiskit_optimization.optionals as _optionals
 from qiskit_optimization.algorithms import (
-    MinimumEigenOptimizer,
     CplexOptimizer,
+    MinimumEigenOptimizer,
     RecursiveMinimumEigenOptimizer,
-    WarmStartQAOAOptimizer,
     SlsqpOptimizer,
+    WarmStartQAOAOptimizer,
 )
-from qiskit_optimization.algorithms.recursive_minimum_eigen_optimizer import (
-    IntermediateResult,
-)
-from qiskit_optimization.problems import QuadraticProgram
+from qiskit_optimization.algorithms.recursive_minimum_eigen_optimizer import IntermediateResult
 from qiskit_optimization.converters import (
-    IntegerToBinary,
     InequalityToEquality,
+    IntegerToBinary,
     LinearEqualityToPenalty,
     QuadraticProgramToQubo,
 )
+from qiskit_optimization.problems import QuadraticProgram
 
 
 class TestRecursiveMinEigenOptimizer(QiskitOptimizationTestCase):
@@ -133,11 +130,9 @@ def test_recursive_warm_qaoa(self):
         """Test the recursive optimizer with warm start qaoa."""
         seed = 1234
         algorithm_globals.random_seed = seed
-        backend = BasicAer.get_backend("statevector_simulator")
         qaoa = QAOA(
-            quantum_instance=QuantumInstance(
-                backend=backend, seed_simulator=seed, seed_transpiler=seed
-            ),
+            sampler=Sampler(),
+            optimizer=SLSQP(),
             reps=1,
         )
         warm_qaoa = WarmStartQAOAOptimizer(
diff --git a/test/algorithms/test_warm_start_qaoa.py b/test/algorithms/test_warm_start_qaoa.py
index 39f0dba0c..be83a392b 100644
--- a/test/algorithms/test_warm_start_qaoa.py
+++ b/test/algorithms/test_warm_start_qaoa.py
@@ -10,22 +10,20 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-""" Test warm start QAOA optimizer. """
+""" Test warm start QAOA optimizer with the primitive-based minimum eigensolver. """
 
 import unittest
 from test import QiskitOptimizationTestCase
 
 import numpy as np
-
 from docplex.mp.model import Model
-from qiskit import BasicAer
-from qiskit.algorithms import QAOA
+from qiskit.algorithms.minimum_eigensolvers import QAOA
+from qiskit.algorithms.optimizers import SLSQP
+from qiskit.primitives.sampler import Sampler
 
 import qiskit_optimization.optionals as _optionals
 from qiskit_optimization.algorithms import SlsqpOptimizer
-from qiskit_optimization.algorithms.goemans_williamson_optimizer import (
-    GoemansWilliamsonOptimizer,
-)
+from qiskit_optimization.algorithms.goemans_williamson_optimizer import GoemansWilliamsonOptimizer
 from qiskit_optimization.algorithms.warm_start_qaoa_optimizer import (
     MeanAggregator,
     WarmStartQAOAOptimizer,
@@ -52,8 +50,7 @@ def test_max_cut(self):
         presolver = GoemansWilliamsonOptimizer(num_cuts=10)
         problem = Maxcut(graph).to_quadratic_program()
 
-        backend = BasicAer.get_backend("statevector_simulator")
-        qaoa = QAOA(quantum_instance=backend, reps=1)
+        qaoa = QAOA(sampler=Sampler(), optimizer=SLSQP(), reps=1)
         aggregator = MeanAggregator()
         optimizer = WarmStartQAOAOptimizer(
             pre_solver=presolver,
@@ -85,8 +82,7 @@ def test_constrained_binary(self):
 
         problem = from_docplex_mp(model)
 
-        backend = BasicAer.get_backend("statevector_simulator")
-        qaoa = QAOA(quantum_instance=backend, reps=1)
+        qaoa = QAOA(sampler=Sampler(), optimizer=SLSQP(), reps=1)
         aggregator = MeanAggregator()
         optimizer = WarmStartQAOAOptimizer(
             pre_solver=SlsqpOptimizer(),
@@ -113,8 +109,7 @@ def test_simple_qubo(self):
         model.minimize((u - v + 2) ** 2)
         problem = from_docplex_mp(model)
 
-        backend = BasicAer.get_backend("statevector_simulator")
-        qaoa = QAOA(quantum_instance=backend, reps=1)
+        qaoa = QAOA(sampler=Sampler(), optimizer=SLSQP(), reps=1)
         optimizer = WarmStartQAOAOptimizer(
             pre_solver=SlsqpOptimizer(),
             relax_for_pre_solver=True,
diff --git a/test/applications/test_optimization_application.py b/test/applications/test_optimization_application.py
new file mode 100644
index 000000000..7c5e5f020
--- /dev/null
+++ b/test/applications/test_optimization_application.py
@@ -0,0 +1,44 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2022.
+#
+# 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 OptimizationApplication class"""
+
+import unittest
+from test.optimization_test_case import QiskitOptimizationTestCase
+
+import numpy as np
+from ddt import data, ddt
+from qiskit.opflow import StateFn
+from qiskit.result import QuasiDistribution
+
+from qiskit_optimization.applications import OptimizationApplication
+
+
+@ddt
+class TestOptimizationApplication(QiskitOptimizationTestCase):
+    """Test OptimizationApplication class"""
+
+    @data(
+        np.array([0, 0, 1, 0]),
+        StateFn([0, 0, 1, 0]),
+        {"10": 0.8, "01": 0.2},
+        QuasiDistribution({"10": 0.8, "01": 0.2}),
+    )
+    def test_sample_most_likely(self, state_vector):
+        """Test sample_most_likely"""
+
+        result = OptimizationApplication.sample_most_likely(state_vector)
+        np.testing.assert_allclose(result, [0, 1])
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/applications/test_sk_model.py b/test/applications/test_sk_model.py
index d4027f091..921d17a08 100644
--- a/test/applications/test_sk_model.py
+++ b/test/applications/test_sk_model.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2021.
+# (C) Copyright IBM 2021, 2022.
 #
 # 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
diff --git a/test/converters/test_converters.py b/test/converters/test_converters.py
index 7290f3877..8086166b0 100644
--- a/test/converters/test_converters.py
+++ b/test/converters/test_converters.py
@@ -17,15 +17,12 @@
 
 import numpy as np
 from docplex.mp.model import Model
-from qiskit.algorithms import NumPyMinimumEigensolver
-from qiskit.opflow import Z, I
+from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver
+from qiskit.opflow import I, Z
+
 import qiskit_optimization.optionals as _optionals
-from qiskit_optimization import QuadraticProgram, QiskitOptimizationError
-from qiskit_optimization.algorithms import (
-    MinimumEigenOptimizer,
-    CplexOptimizer,
-    ADMMOptimizer,
-)
+from qiskit_optimization import QiskitOptimizationError, QuadraticProgram
+from qiskit_optimization.algorithms import ADMMOptimizer, CplexOptimizer, MinimumEigenOptimizer
 from qiskit_optimization.algorithms.admm_optimizer import ADMMParameters
 from qiskit_optimization.converters import (
     InequalityToEquality,
@@ -36,7 +33,6 @@
 from qiskit_optimization.problems import Constraint, Variable
 from qiskit_optimization.translators import from_docplex_mp
 
-
 QUBIT_OP_MAXIMIZE_SAMPLE = (
     -199999.5 * (I ^ I ^ I ^ Z)
     + -399999.5 * (I ^ I ^ Z ^ I)