diff --git a/docs/source/api/quad/solvers.acquisition_functions.rst b/docs/source/api/quad/solvers.acquisition_functions.rst
new file mode 100644
index 000000000..802557d47
--- /dev/null
+++ b/docs/source/api/quad/solvers.acquisition_functions.rst
@@ -0,0 +1,5 @@
+Acquisition Functions
+---------------------
+.. automodapi:: probnum.quad.solvers.acquisition_functions
+    :no-heading:
+    :headings: "*"
diff --git a/docs/source/api/quad/solvers.rst b/docs/source/api/quad/solvers.rst
index 955111967..501174120 100644
--- a/docs/source/api/quad/solvers.rst
+++ b/docs/source/api/quad/solvers.rst
@@ -15,6 +15,11 @@ probnum.quad.solvers
 
     solvers.policies
 
+.. toctree::
+    :hidden:
+
+    solvers.acquisition_functions
+
 .. toctree::
     :hidden:
 
diff --git a/src/probnum/quad/_bayesquad.py b/src/probnum/quad/_bayesquad.py
index 6dd889989..393a41bb7 100644
--- a/src/probnum/quad/_bayesquad.py
+++ b/src/probnum/quad/_bayesquad.py
@@ -6,6 +6,7 @@
 methods return a random variable, specifying the belief about the true value of the
 integral.
 """
+
 from __future__ import annotations
 
 from typing import Callable, Optional, Tuple
@@ -70,10 +71,11 @@ def bayesquad(
     policy
         Type of acquisition strategy to use. Defaults to 'bmc'. Options are
 
-        ==========================  =======
-         Bayesian Monte Carlo [2]_  ``bmc``
-         van Der Corput points      ``vdc``
-        ==========================  =======
+        ============================================  ===========
+         Bayesian Monte Carlo [2]_                    ``bmc``
+         Van Der Corput points                        ``vdc``
+         Uncertainty Sampling with random candidates  ``us_rand``
+        ============================================  ===========
 
     initial_design
         The type of initial design to use. If ``None`` is given, no initial design is
@@ -112,6 +114,9 @@ def bayesquad(
             num_initial_design_nodes : Optional[IntLike]
                 The number of nodes created by the initial design. Defaults to
                 ``input_dim * 5`` if an initial design is given.
+            us_rand_num_candidates : Optional[IntLike]
+                The number of candidate nodes used by the policy 'us_rand'. Defaults
+                to 1e2.
 
     Returns
     -------
diff --git a/src/probnum/quad/solvers/__init__.py b/src/probnum/quad/solvers/__init__.py
index 02ec29ec4..a2c996044 100644
--- a/src/probnum/quad/solvers/__init__.py
+++ b/src/probnum/quad/solvers/__init__.py
@@ -1,6 +1,12 @@
 """Bayesian quadrature methods and their components."""
 
-from . import belief_updates, initial_designs, policies, stopping_criteria
+from . import (
+    acquisition_functions,
+    belief_updates,
+    initial_designs,
+    policies,
+    stopping_criteria,
+)
 from ._bayesian_quadrature import BayesianQuadrature
 from ._bq_state import BQIterInfo, BQState
 
diff --git a/src/probnum/quad/solvers/_bayesian_quadrature.py b/src/probnum/quad/solvers/_bayesian_quadrature.py
index bdeef6257..7e06cfbcb 100644
--- a/src/probnum/quad/solvers/_bayesian_quadrature.py
+++ b/src/probnum/quad/solvers/_bayesian_quadrature.py
@@ -10,9 +10,15 @@
 from probnum.quad.integration_measures import IntegrationMeasure, LebesgueMeasure
 from probnum.quad.kernel_embeddings import KernelEmbedding
 from probnum.quad.solvers._bq_state import BQIterInfo, BQState
+from probnum.quad.solvers.acquisition_functions import WeightedPredictiveVariance
 from probnum.quad.solvers.belief_updates import BQBeliefUpdate, BQStandardBeliefUpdate
 from probnum.quad.solvers.initial_designs import InitialDesign, LatinDesign, MCDesign
-from probnum.quad.solvers.policies import Policy, RandomPolicy, VanDerCorputPolicy
+from probnum.quad.solvers.policies import (
+    Policy,
+    RandomMaxAcquisitionPolicy,
+    RandomPolicy,
+    VanDerCorputPolicy,
+)
 from probnum.quad.solvers.stopping_criteria import (
     BQStoppingCriterion,
     ImmediateStop,
@@ -38,7 +44,7 @@ class BayesianQuadrature:
     Parameters
     ----------
     kernel
-        The kernel used for the GP model.
+        The kernel used for the Gaussian process model.
     measure
         The integration measure.
     policy
@@ -139,6 +145,9 @@ def from_problem(
                 num_initial_design_nodes : Optional[IntLike]
                     The number of nodes created by the initial design. Defaults to
                     ``input_dim * 5`` if an initial design is given.
+                us_rand_num_candidates : Optional[IntLike]
+                    The number of candidate nodes used by the policy 'us_rand'. Defaults
+                    to 1e2.
 
         Returns
         -------
@@ -175,6 +184,7 @@ def from_problem(
         num_initial_design_nodes = options.get(
             "num_initial_design_nodes", int(5 * input_dim)
         )
+        us_rand_num_candidates = options.get("us_rand_num_candidates", int(1e2))
 
         # Set up integration measure
         if domain is None and measure is None:
@@ -198,6 +208,12 @@ def from_problem(
             policy = RandomPolicy(batch_size, measure.sample)
         elif policy == "vdc":
             policy = VanDerCorputPolicy(batch_size, measure)
+        elif policy == "us_rand":
+            policy = RandomMaxAcquisitionPolicy(
+                batch_size=1,
+                acquisition_func=WeightedPredictiveVariance(),
+                n_candidates=us_rand_num_candidates,
+            )
         else:
             raise NotImplementedError(f"The given policy ({policy}) is unknown.")
 
diff --git a/src/probnum/quad/solvers/_bq_state.py b/src/probnum/quad/solvers/_bq_state.py
index 337d5c571..095fa3d34 100644
--- a/src/probnum/quad/solvers/_bq_state.py
+++ b/src/probnum/quad/solvers/_bq_state.py
@@ -37,6 +37,8 @@ class BQState:
         Function evaluations at nodes.
     gram
         The kernel Gram matrix.
+    gram_cho_factor
+        The output of BQBeliefUpdate.compute_gram_cho_factor.
     kernel_means
         All kernel mean evaluations at ``nodes``.
 
@@ -55,6 +57,7 @@ def __init__(
         nodes: Optional[np.ndarray] = None,
         fun_evals: Optional[np.ndarray] = None,
         gram: np.ndarray = np.array([[]]),
+        gram_cho_factor: Tuple[np.ndarray, bool] = (np.array([[]]), False),
         kernel_means: np.ndarray = np.array([]),
     ):
         self.measure = measure
@@ -73,6 +76,7 @@ def __init__(
             self.fun_evals = fun_evals
 
         self.gram = gram
+        self.gram_cho_factor = gram_cho_factor
         self.kernel_means = kernel_means
 
     @classmethod
@@ -85,6 +89,7 @@ def from_new_data(
         integral_belief: Normal,
         prev_state: "BQState",
         gram: np.ndarray,
+        gram_cho_factor: Tuple[np.ndarray, bool],
         kernel_means: np.ndarray,
     ) -> "BQState":
         r"""Initialize state from updated data.
@@ -105,6 +110,8 @@ def from_new_data(
             Previous state of the BQ loop.
         gram
             The Gram matrix of the given nodes.
+        gram_cho_factor
+            The output of BQBeliefUpdate.compute_gram_cho_factor for ``gram``.
         kernel_means
             The kernel means at the given nodes.
 
@@ -123,6 +130,7 @@ def from_new_data(
             nodes=nodes,
             fun_evals=fun_evals,
             gram=gram,
+            gram_cho_factor=gram_cho_factor,
             kernel_means=kernel_means,
         )
 
diff --git a/src/probnum/quad/solvers/acquisition_functions/__init__.py b/src/probnum/quad/solvers/acquisition_functions/__init__.py
new file mode 100644
index 000000000..377256f6f
--- /dev/null
+++ b/src/probnum/quad/solvers/acquisition_functions/__init__.py
@@ -0,0 +1,13 @@
+"""Acquisition functions for Bayesian quadrature."""
+
+from ._acquisition_function import AcquisitionFunction
+from ._predictive_variance import WeightedPredictiveVariance
+
+# Public classes and functions. Order is reflected in documentation.
+__all__ = [
+    "AcquisitionFunction",
+    "WeightedPredictiveVariance",
+]
+
+# Set correct module paths. Corrects links and module paths in documentation.
+WeightedPredictiveVariance.__module__ = "probnum.quad.solvers.acquisition_functions"
diff --git a/src/probnum/quad/solvers/acquisition_functions/_acquisition_function.py b/src/probnum/quad/solvers/acquisition_functions/_acquisition_function.py
new file mode 100644
index 000000000..6f689234d
--- /dev/null
+++ b/src/probnum/quad/solvers/acquisition_functions/_acquisition_function.py
@@ -0,0 +1,43 @@
+"""Abstract base class for BQ acquisition functions."""
+
+from __future__ import annotations
+
+import abc
+from typing import Optional, Tuple
+
+import numpy as np
+
+from probnum.quad.solvers._bq_state import BQState
+
+
+class AcquisitionFunction(abc.ABC):
+    """An abstract class for an acquisition function for Bayesian quadrature."""
+
+    @property
+    @abc.abstractmethod
+    def has_gradients(self) -> bool:
+        """Whether the acquisition function exposes gradients."""
+        raise NotImplementedError
+
+    @abc.abstractmethod
+    def __call__(
+        self, x: np.ndarray, bq_state: BQState
+    ) -> Tuple[np.ndarray, Optional[np.ndarray]]:
+        """Evaluates the acquisition function and optionally its gradients.
+
+        Parameters
+        ----------
+        x
+            *shape=(batch_size, input_dim)* -- The nodes where the acquisition function
+            is being evaluated.
+        bq_state
+            State of the BQ belief.
+
+        Returns
+        -------
+        acquisition_values :
+            *shape=(batch_size, )* -- The acquisition values at nodes ``x``.
+        acquisition_gradients :
+            *shape=(batch_size, input_dim)* -- The corresponding gradients (optional).
+        """
+        raise NotImplementedError
diff --git a/src/probnum/quad/solvers/acquisition_functions/_predictive_variance.py b/src/probnum/quad/solvers/acquisition_functions/_predictive_variance.py
new file mode 100644
index 000000000..64c09a350
--- /dev/null
+++ b/src/probnum/quad/solvers/acquisition_functions/_predictive_variance.py
@@ -0,0 +1,48 @@
+"""Uncertainty sampling for Bayesian Monte Carlo."""
+
+from __future__ import annotations
+
+from typing import Optional, Tuple
+
+import numpy as np
+
+from probnum.quad.solvers._bq_state import BQState
+from probnum.quad.solvers.belief_updates import BQStandardBeliefUpdate
+
+from ._acquisition_function import AcquisitionFunction
+
+# pylint: disable=too-few-public-methods, fixme
+
+
+class WeightedPredictiveVariance(AcquisitionFunction):
+    r"""The predictive variance acquisition function that yields uncertainty sampling.
+
+    The acquisition function is
+
+    .. math::
+        a(x) = \operatorname{Var}(f(x)) p(x)^2
+
+    where :math:`\operatorname{Var}(f(x))` is the predictive variance of the model and
+    :math:`p(x)` is the density of the integration measure :math:`\mu`.
+
+    """
+
+    @property
+    def has_gradients(self) -> bool:
+        # Todo (#581): this needs to return True, once gradients are available
+        return False
+
+    def __call__(
+        self,
+        x: np.ndarray,
+        bq_state: BQState,
+    ) -> Tuple[np.ndarray, Optional[np.ndarray]]:
+        predictive_variance = bq_state.kernel(x, x)
+        if bq_state.fun_evals.shape != (0,):
+            kXx = bq_state.kernel.matrix(bq_state.nodes, x)
+            regression_weights = BQStandardBeliefUpdate.gram_cho_solve(
+                bq_state.gram_cho_factor, kXx
+            )
+            predictive_variance -= np.sum(regression_weights * kXx, axis=0)
+        values = bq_state.scale_sq * predictive_variance * bq_state.measure(x) ** 2
+        return values, None
diff --git a/src/probnum/quad/solvers/belief_updates/_belief_update.py b/src/probnum/quad/solvers/belief_updates/_belief_update.py
index 4f9d3dc42..20cc28582 100644
--- a/src/probnum/quad/solvers/belief_updates/_belief_update.py
+++ b/src/probnum/quad/solvers/belief_updates/_belief_update.py
@@ -58,7 +58,7 @@ def __call__(
         """
         raise NotImplementedError
 
-    def _compute_gram_cho_factor(self, gram: np.ndarray) -> np.ndarray:
+    def compute_gram_cho_factor(self, gram: np.ndarray) -> Tuple[np.ndarray, bool]:
         """Compute the Cholesky decomposition of a positive-definite Gram matrix for use
         in scipy.linalg.cho_solve
 
@@ -68,19 +68,36 @@ def _compute_gram_cho_factor(self, gram: np.ndarray) -> np.ndarray:
 
         Parameters
         ----------
-        gram :
+        gram
             symmetric pos. def. kernel Gram matrix :math:`K`, shape (nevals, nevals)
 
         Returns
         -------
         gram_cho_factor :
             The upper triangular Cholesky decomposition of the Gram matrix. Other
-            parts of the matrix contain random data.
+            parts of the matrix contain random data. A boolean that indicates whether
+            the matrix is lower triangular (always False but needed for scipy).
         """
         return cho_factor(gram + self.jitter * np.eye(gram.shape[0]))
 
-    # pylint: disable=no-self-use
-    def _gram_cho_solve(self, gram_cho_factor: np.ndarray, z: np.ndarray) -> np.ndarray:
+    @staticmethod
+    def gram_cho_solve(
+        gram_cho_factor: Tuple[np.ndarray, bool], z: np.ndarray
+    ) -> np.ndarray:
         """Wrapper for scipy.linalg.cho_solve. Meant to be used for linear systems of
-        the gram matrix. Requires the solution of scipy.linalg.cho_factor as input."""
+        the gram matrix. Requires the solution of scipy.linalg.cho_factor as input.
+
+        Parameters
+        ----------
+        gram_cho_factor
+            The return object of compute_gram_cho_factor.
+        z
+            An array of appropriate shape.
+
+        Returns
+        -------
+        solution :
+            The solution ``x`` to the linear system ``gram x = z``.
+
+        """
         return cho_solve(gram_cho_factor, z)
diff --git a/src/probnum/quad/solvers/belief_updates/_standard_update.py b/src/probnum/quad/solvers/belief_updates/_standard_update.py
index 2991d0b48..8c95d3973 100644
--- a/src/probnum/quad/solvers/belief_updates/_standard_update.py
+++ b/src/probnum/quad/solvers/belief_updates/_standard_update.py
@@ -72,13 +72,13 @@ def __call__(
             )
 
         # Cholesky factorisation of the Gram matrix
-        gram_cho_factor = self._compute_gram_cho_factor(gram)
+        gram_cho_factor = self.compute_gram_cho_factor(gram)
 
         # Estimate scaling parameter
         new_scale_sq = self._estimate_scale(fun_evals, gram_cho_factor, bq_state)
 
         # Integral mean and variance
-        weights = self._gram_cho_solve(gram_cho_factor, kernel_means)
+        weights = self.gram_cho_solve(gram_cho_factor, kernel_means)
         integral_mean = weights @ fun_evals
         initial_integral_variance = new_kernel_embedding.kernel_variance()
         integral_variance = new_scale_sq * (
@@ -94,6 +94,7 @@ def __call__(
             integral_belief=new_belief,
             prev_state=bq_state,
             gram=gram,
+            gram_cho_factor=gram_cho_factor,
             kernel_means=kernel_means,
         )
 
@@ -108,7 +109,10 @@ def _estimate_kernel(self, kernel: Kernel) -> Tuple[Kernel, bool]:
         return new_kernel, kernel_was_updated
 
     def _estimate_scale(
-        self, fun_evals: np.ndarray, gram_cho_factor: np.ndarray, bq_state: BQState
+        self,
+        fun_evals: np.ndarray,
+        gram_cho_factor: Tuple[np.ndarray, bool],
+        bq_state: BQState,
     ) -> FloatLike:
         """Estimate the scale parameter."""
         if self.scale_estimation is None:
@@ -116,7 +120,7 @@ def _estimate_scale(
         elif self.scale_estimation == "mle":
             new_scale_sq = (
                 fun_evals
-                @ self._gram_cho_solve(gram_cho_factor, fun_evals)
+                @ self.gram_cho_solve(gram_cho_factor, fun_evals)
                 / fun_evals.shape[0]
             )
         else:
diff --git a/src/probnum/quad/solvers/policies/__init__.py b/src/probnum/quad/solvers/policies/__init__.py
index e53c1b18c..eb78873e1 100644
--- a/src/probnum/quad/solvers/policies/__init__.py
+++ b/src/probnum/quad/solvers/policies/__init__.py
@@ -1,5 +1,6 @@
 """Policies for Bayesian quadrature."""
 
+from ._max_aquisition_policy import RandomMaxAcquisitionPolicy
 from ._policy import Policy
 from ._random_policy import RandomPolicy
 from ._van_der_corput_policy import VanDerCorputPolicy
@@ -9,9 +10,11 @@
     "Policy",
     "RandomPolicy",
     "VanDerCorputPolicy",
+    "RandomMaxAcquisitionPolicy",
 ]
 
 # Set correct module paths. Corrects links and module paths in documentation.
 Policy.__module__ = "probnum.quad.solvers.policies"
 RandomPolicy.__module__ = "probnum.quad.solvers.policies"
 VanDerCorputPolicy.__module__ = "probnum.quad.solvers.policies"
+RandomMaxAcquisitionPolicy.__module__ = "probnum.quad.solvers.policies"
diff --git a/src/probnum/quad/solvers/policies/_max_aquisition_policy.py b/src/probnum/quad/solvers/policies/_max_aquisition_policy.py
new file mode 100644
index 000000000..b374237e6
--- /dev/null
+++ b/src/probnum/quad/solvers/policies/_max_aquisition_policy.py
@@ -0,0 +1,71 @@
+"""Max acquisition policy for Bayesian quadrature."""
+
+from __future__ import annotations
+
+from typing import Optional
+
+import numpy as np
+
+from probnum.quad.solvers._bq_state import BQState
+from probnum.quad.solvers.acquisition_functions import AcquisitionFunction
+from probnum.typing import IntLike
+
+from ._policy import Policy
+
+# pylint: disable=too-few-public-methods, fixme
+
+
+class RandomMaxAcquisitionPolicy(Policy):
+    """Policy that maximizes an acquisition function by sampling random candidate nodes.
+
+    The candidate nodes are random draws from the integration measure. The node with the
+    largest acquisition value is chosen.
+
+    Parameters
+    ----------
+    batch_size
+        Size of batch of nodes when calling the policy once (must be equal to 1).
+    acquisition_func
+        The acquisition function.
+    n_candidates
+        The number of candidate nodes.
+
+    Raises
+    ------
+    ValueError
+        If ``batch_size`` is not 1, or if ``n_candidates`` is too small.
+    """
+
+    def __init__(
+        self,
+        batch_size: IntLike,
+        acquisition_func: AcquisitionFunction,
+        n_candidates: IntLike,
+    ) -> None:
+
+        if batch_size != 1:
+            raise ValueError(
+                f"RandomMaxAcquisitionPolicy can only be used with batch "
+                f"size 1 ({batch_size})."
+            )
+        if n_candidates < 1:
+            raise ValueError(
+                f"The number of candidates ({n_candidates}) must be equal "
+                f"or larger than 1."
+            )
+
+        super().__init__(batch_size=batch_size)
+        self.acquisition_func = acquisition_func
+        self.n_candidates = int(n_candidates)
+
+    @property
+    def requires_rng(self) -> bool:
+        return True
+
+    def __call__(
+        self, bq_state: BQState, rng: Optional[np.random.Generator]
+    ) -> np.ndarray:
+        random_nodes = bq_state.measure.sample(n_sample=self.n_candidates, rng=rng)
+        values = self.acquisition_func(random_nodes, bq_state)[0]
+        idx_max = int(np.argmax(values))
+        return random_nodes[idx_max, :][None, :]
diff --git a/src/probnum/quad/solvers/policies/_van_der_corput_policy.py b/src/probnum/quad/solvers/policies/_van_der_corput_policy.py
index 63ab5a3b1..ec31bcee1 100644
--- a/src/probnum/quad/solvers/policies/_van_der_corput_policy.py
+++ b/src/probnum/quad/solvers/policies/_van_der_corput_policy.py
@@ -30,6 +30,13 @@ class VanDerCorputPolicy(Policy):
     measure
         The integration measure with finite domain.
 
+    Raises
+    ------
+    ValueError
+        If input dimension is not 1.
+    ValueError
+        If measure domain is not bounded.
+
     References
     ----------
     .. [1] https://en.wikipedia.org/wiki/Van_der_Corput_sequence
diff --git a/tests/test_quad/test_acquisition_functions.py b/tests/test_quad/test_acquisition_functions.py
new file mode 100644
index 000000000..020ed5540
--- /dev/null
+++ b/tests/test_quad/test_acquisition_functions.py
@@ -0,0 +1,62 @@
+"""Basic tests for BQ acquisition functions."""
+
+# New acquisition functions need to be added to the fixtures 'acquisition'.
+
+
+import numpy as np
+import pytest
+
+from probnum.quad.integration_measures import LebesgueMeasure
+from probnum.quad.solvers import BQState
+from probnum.quad.solvers.acquisition_functions import WeightedPredictiveVariance
+from probnum.randprocs.kernels import ExpQuad
+from probnum.randvars import Normal
+
+
+@pytest.fixture(params=[pytest.param(n, id=f"nevals{n}") for n in [1, 3, 5]])
+def nevals(request):
+    return request.param
+
+
+@pytest.fixture(params=[pytest.param(a) for a in [WeightedPredictiveVariance]])
+def acquisition(request):
+    return request.param()
+
+
+def test_acquisition_shapes(acquisition, input_dim, nevals, rng):
+
+    m = LebesgueMeasure(input_dim=input_dim, domain=(0, 1))
+    k = ExpQuad(input_shape=(input_dim,))
+    bq_state_no_data = BQState(
+        measure=m,
+        kernel=k,
+    )
+    bq_state = BQState(
+        measure=m,
+        kernel=k,
+        nodes=np.zeros([nevals, input_dim]),
+        fun_evals=np.ones(nevals),
+        gram=np.eye(nevals),
+        gram_cho_factor=(np.eye(nevals), False),
+        kernel_means=np.ones(nevals),
+        previous_integral_beliefs=nevals * (Normal(mean=0.0, cov=1.0),),
+        integral_belief=Normal(mean=0.0, cov=1.0),
+    )
+
+    n_nodes = 3
+    x = np.ones([n_nodes, input_dim])
+
+    # no data yet in bq_state
+    res = acquisition(x, bq_state_no_data)
+    assert res[0].shape == (n_nodes,)  # values
+    assert res[1] is None  # no gradients yet
+
+    # data has been collected previously
+    res = acquisition(x, bq_state)
+    assert res[0].shape == (n_nodes,)  # values
+    assert res[1] is None  # no gradients yet
+
+
+def test_acquisition_property_values(acquisition):
+    # no gradients yet. this may change (#581).
+    assert not acquisition.has_gradients
diff --git a/tests/test_quad/test_bayesian_quadrature.py b/tests/test_quad/test_bayesian_quadrature.py
index 838a3931c..29fdba289 100644
--- a/tests/test_quad/test_bayesian_quadrature.py
+++ b/tests/test_quad/test_bayesian_quadrature.py
@@ -8,7 +8,11 @@
 from probnum.quad.solvers import BayesianQuadrature
 from probnum.quad.solvers.belief_updates import BQStandardBeliefUpdate
 from probnum.quad.solvers.initial_designs import LatinDesign, MCDesign
-from probnum.quad.solvers.policies import RandomPolicy, VanDerCorputPolicy
+from probnum.quad.solvers.policies import (
+    RandomMaxAcquisitionPolicy,
+    RandomPolicy,
+    VanDerCorputPolicy,
+)
 from probnum.quad.solvers.stopping_criteria import (
     ImmediateStop,
     IntegralVarianceTolerance,
@@ -80,7 +84,12 @@ def test_bayesian_quadrature_wrong_input(input_dim):
 
 
 @pytest.mark.parametrize(
-    "policy, policy_type", [("bmc", RandomPolicy), ("vdc", VanDerCorputPolicy)]
+    "policy, policy_type",
+    [
+        ("bmc", RandomPolicy),
+        ("vdc", VanDerCorputPolicy),
+        ("us_rand", RandomMaxAcquisitionPolicy),
+    ],
 )
 def test_bq_from_problem_policy_assignment(policy, policy_type):
     """Test if correct policy is assigned from string identifier."""
@@ -121,7 +130,7 @@ def test_bq_from_problem_stopping_criterion_assignment(max_evals, var_tol, rel_t
     assert isinstance(bq.stopping_criterion, t)
 
 
-def test_bq_from_problem_defaults(bq, bq_no_policy):
+def test_bq_from_problem_default_attribute_types(bq, bq_no_policy):
 
     # defaults if policy is available
     assert isinstance(bq.measure, LebesgueMeasure)
@@ -138,6 +147,78 @@ def test_bq_from_problem_defaults(bq, bq_no_policy):
     assert bq_no_policy.initial_design is None
 
 
+def test_bq_from_problem_options_default_values():
+    """Check if default values of the options dictionary are set correctly."""
+
+    bq = BayesianQuadrature.from_problem(
+        input_dim=2,
+        domain=(0, 1),
+    )
+
+    # batch_size default
+    assert bq.policy.batch_size == 1
+
+    # jitter default
+    assert bq.belief_update.jitter == 1e-8
+
+    # n_candidates for policy 'us_rand'
+    bq = BayesianQuadrature.from_problem(input_dim=2, domain=(0, 1), policy="us_rand")
+    assert bq.policy.n_candidates == int(1e2)
+
+    # num_initial_design_nodes for initial design
+    input_dim = 5
+    bq = BayesianQuadrature.from_problem(
+        input_dim=input_dim,
+        domain=(0, 1),
+        initial_design="mc",
+    )
+    assert bq.initial_design.num_nodes == int(5 * input_dim)
+
+
+def test_bq_from_problem_options_custom_values(bq, bq_no_policy):
+    """Check if custom values of the options dictionary are set correctly."""
+
+    # batch_size manual value
+    batch_size = 3
+    bq = BayesianQuadrature.from_problem(
+        input_dim=2,
+        domain=(0, 1),
+        options=dict(batch_size=batch_size),
+    )
+    assert bq.policy.batch_size == batch_size
+
+    # jitter manual value
+    jitter = 1.3
+    bq = BayesianQuadrature.from_problem(
+        input_dim=2,
+        domain=(0, 1),
+        options=dict(jitter=jitter),
+    )
+    assert bq.belief_update.jitter == jitter
+
+    # n_candidates for policy 'us_rand'
+    us_rand_num_candidates = 5
+    bq = BayesianQuadrature.from_problem(
+        input_dim=2,
+        domain=(0, 1),
+        policy="us_rand",
+        options=dict(us_rand_num_candidates=us_rand_num_candidates),
+    )
+    assert bq.policy.n_candidates == us_rand_num_candidates
+
+    # num_initial_design_nodes for initial design
+    input_dim = 5
+    num_initial_design_nodes = 3
+    assert int(input_dim * 5) != num_initial_design_nodes
+    bq = BayesianQuadrature.from_problem(
+        input_dim=input_dim,
+        domain=(0, 1),
+        initial_design="mc",
+        options=dict(num_initial_design_nodes=num_initial_design_nodes),
+    )
+    assert bq.initial_design.num_nodes == num_initial_design_nodes
+
+
 # Tests for input checks and exception raises start here.
 
 
diff --git a/tests/test_quad/test_bq_state.py b/tests/test_quad/test_bq_state.py
index 85123f06c..1300bc8a0 100644
--- a/tests/test_quad/test_bq_state.py
+++ b/tests/test_quad/test_bq_state.py
@@ -105,6 +105,8 @@ def test_state_defaults_types(state, request):
     assert isinstance(s.nodes, np.ndarray)
     assert isinstance(s.fun_evals, np.ndarray)
     assert isinstance(s.gram, np.ndarray)
+    assert isinstance(s.gram_cho_factor[0], np.ndarray)
+    assert isinstance(s.gram_cho_factor[1], bool)
     assert isinstance(s.kernel_means, np.ndarray)
     assert isinstance(s.previous_integral_beliefs, tuple)
     assert isinstance(s.scale_sq, float)
@@ -122,6 +124,7 @@ def test_state_defaults_shapes(state, request):
     s = request.getfixturevalue(state)
     assert len(s.previous_integral_beliefs) == 0
     assert s.gram.shape == (1, 0)
+    assert s.gram_cho_factor[0].shape == (1, 0)
     assert s.kernel_means.shape == (0,)
 
 
@@ -151,6 +154,7 @@ def test_state_from_new_data(state, request):
     y = np.ones(new_nevals)
     integral = Normal(0, 1)
     gram = np.eye(new_nevals)
+    gram_cho_factor = (np.eye(new_nevals), False)
     kernel_means = np.ones(new_nevals)
     kernel = ExpQuad(input_shape=(old_state.input_dim,))
     scale_sq = 1.7
@@ -164,6 +168,7 @@ def test_state_from_new_data(state, request):
         integral_belief=integral,
         prev_state=old_state,
         gram=gram,
+        gram_cho_factor=gram_cho_factor,
         kernel_means=kernel_means,
     )
 
@@ -174,6 +179,8 @@ def test_state_from_new_data(state, request):
     assert isinstance(s.nodes, np.ndarray)
     assert isinstance(s.fun_evals, np.ndarray)
     assert isinstance(s.gram, np.ndarray)
+    assert isinstance(s.gram_cho_factor[0], np.ndarray)
+    assert isinstance(s.gram_cho_factor[1], bool)
     assert isinstance(s.kernel_means, np.ndarray)
     assert isinstance(s.integral_belief, Normal)
     assert isinstance(s.previous_integral_beliefs, tuple)
@@ -183,6 +190,7 @@ def test_state_from_new_data(state, request):
     assert s.fun_evals.shape == (new_nevals,)
     assert len(s.previous_integral_beliefs) == 1
     assert s.gram.shape == (new_nevals, new_nevals)
+    assert s.gram_cho_factor[0].shape == (new_nevals, new_nevals)
     assert s.kernel_means.shape == (new_nevals,)
 
     # values
diff --git a/tests/test_quad/test_policy.py b/tests/test_quad/test_policy.py
index 76e5a9e94..3b2dd0d81 100644
--- a/tests/test_quad/test_policy.py
+++ b/tests/test_quad/test_policy.py
@@ -1,7 +1,8 @@
 """Basic tests for BQ policies."""
 
 # New policies need to be added to the fixtures 'policy_name' and 'policy_params'
-# and 'policy'.
+# and 'policy'. Further, add the new policy to the assignment test in
+# test_bayesian_quadrature.py.
 
 
 import numpy as np
@@ -9,18 +10,24 @@
 
 from probnum.quad.integration_measures import GaussianMeasure, LebesgueMeasure
 from probnum.quad.solvers import BQState
-from probnum.quad.solvers.policies import RandomPolicy, VanDerCorputPolicy
+from probnum.quad.solvers.acquisition_functions import WeightedPredictiveVariance
+from probnum.quad.solvers.policies import (
+    RandomMaxAcquisitionPolicy,
+    RandomPolicy,
+    VanDerCorputPolicy,
+)
 from probnum.randprocs.kernels import ExpQuad
 
 
-@pytest.fixture
-def batch_size():
-    return 3
+@pytest.fixture(params=[pytest.param(s, id=f"bs{s}") for s in [1, 3, 5]])
+def batch_size(request):
+    return request.param
 
 
 @pytest.fixture(
     params=[
-        pytest.param(name, id=name) for name in ["RandomPolicy", "VanDerCorputPolicy"]
+        pytest.param(name, id=name)
+        for name in ["RandomPolicy", "VanDerCorputPolicy", "RandomMaxAcquisitionPolicy"]
     ]
 )
 def policy_name(request):
@@ -40,10 +47,12 @@ def _get_bq_states(ndim):
             kernel=ExpQuad(input_shape=(ndim,)),
             nodes=np.zeros([nevals, ndim]),
             fun_evals=np.ones(nevals),
+            gram=np.eye(nevals),
+            gram_cho_factor=(np.eye(nevals), False),
         )
         return bq_state, bq_state_no_data
 
-    params = dict(name=policy_name, ndim=input_dim)
+    params = dict(name=policy_name, input_dim=input_dim)
     params["bq_state"], params["bq_state_no_data"] = _get_bq_states(input_dim)
 
     if policy_name == "RandomPolicy":
@@ -61,8 +70,16 @@ def _get_bq_states(ndim):
             measure=LebesgueMeasure(input_dim=1, domain=(0, 1)),
         )
         params["bq_state"], params["bq_state_no_data"] = _get_bq_states(1)
-        params["ndim"] = 1
+        params["input_dim"] = 1
         params["requires_rng"] = False
+    elif policy_name == "RandomMaxAcquisitionPolicy":
+        # Since RandomMaxAcquisitionPolicy requires batch_size=1, we override it here.
+        input_params = dict(
+            batch_size=1,
+            acquisition_func=WeightedPredictiveVariance(),
+            n_candidates=10,
+        )
+        params["requires_rng"] = True
     else:
         raise NotImplementedError
 
@@ -74,29 +91,31 @@ def _get_bq_states(ndim):
 @pytest.fixture()
 def policy(policy_params):
     name = policy_params.pop("name")
-    input_params = policy_params.pop("input_params")
-
-    if name == "RandomPolicy":
-        return RandomPolicy(**input_params), policy_params
-    elif name == "VanDerCorputPolicy":
-        return VanDerCorputPolicy(**input_params), policy_params
-    else:
-        raise NotImplementedError
+    input_params = policy_params["input_params"]
+    # add new policy to this dict
+    policy = dict(
+        RandomPolicy=RandomPolicy,
+        VanDerCorputPolicy=VanDerCorputPolicy,
+        RandomMaxAcquisitionPolicy=RandomMaxAcquisitionPolicy,
+    )
+    return policy[name](**input_params), policy_params
 
 
 # Tests shared by all policies start here.
 
 
-def test_policy_shapes(policy, batch_size, rng):
+def test_policy_shapes(policy, rng):
     policy, params = policy
-    bq_state, bq_state_no_data = params["bq_state"], params["bq_state_no_data"]
-    ndim = params["ndim"]
+    bq_state = params["bq_state"]
+    bq_state_no_data = params["bq_state_no_data"]
+    input_dim = params["input_dim"]
+    batch_size = params["input_params"]["batch_size"]
 
     # bq state contains data
-    assert policy(bq_state, rng).shape == (batch_size, ndim)
+    assert policy(bq_state, rng).shape == (batch_size, input_dim)
 
     # bq state contains no data yet
-    assert policy(bq_state_no_data, rng).shape == (batch_size, ndim)
+    assert policy(bq_state_no_data, rng).shape == (batch_size, input_dim)
 
 
 def test_policy_property_values(policy):
@@ -104,6 +123,39 @@ def test_policy_property_values(policy):
     assert policy.requires_rng is params["requires_rng"]
 
 
+# Tests specific to RandomMaxAcquisitionPolicy start here
+
+
+def test_random_max_acquisition_attribute_values():
+    n_cadidates = 10
+    policy = RandomMaxAcquisitionPolicy(
+        batch_size=1,
+        acquisition_func=WeightedPredictiveVariance(),
+        n_candidates=n_cadidates,
+    )
+    assert policy.n_candidates == n_cadidates
+
+
+def test_random_max_acquisition_raises():
+    # batch size larger than 1
+    with pytest.raises(ValueError):
+        wrong_batch_size = 2
+        RandomMaxAcquisitionPolicy(
+            batch_size=wrong_batch_size,
+            acquisition_func=WeightedPredictiveVariance(),
+            n_candidates=10,
+        )
+
+    # n_candidates too small
+    with pytest.raises(ValueError):
+        wrong_n_candidates = 0
+        RandomMaxAcquisitionPolicy(
+            batch_size=1,
+            acquisition_func=WeightedPredictiveVariance(),
+            n_candidates=wrong_n_candidates,
+        )
+
+
 # Tests specific to VanDerCorputPolicy start here