diff --git a/benchmarks/random_variables.py b/benchmarks/random_variables.py
index 4e00242214..1e09ea2ef0 100644
--- a/benchmarks/random_variables.py
+++ b/benchmarks/random_variables.py
@@ -13,7 +13,7 @@
     "multivar_normal",
     "matrixvar_normal",
     "symmatrixvar_normal",
-    # "operatorvar_normal"
+    "operatorvar_normal",
 ]
 
 
@@ -59,15 +59,15 @@ class Functions:
     Benchmark various functions of distributions.
     """
 
-    param_names = ["rv_name", "property"]
+    param_names = ["randvar", "property"]
     params = [rv_names, ["pdf", "logpdf", "cdf", "logcdf"]]
 
-    def setup(self, rv_name, property):
-        self.randvar = get_randvar(rv_name=rv_name)
-        self.eval_point = np.random.uniform(self.randvar.shape)
-        self.quantile = np.random.uniform(self.randvar.shape)
+    def setup(self, randvar, property):
+        self.randvar = get_randvar(rv_name=randvar)
+        self.eval_point = np.random.uniform(size=self.randvar.shape)
+        self.quantile = np.random.uniform(size=self.randvar.shape)
 
-    def time_distr_functions(self, rv_name, property):
+    def time_distr_functions(self, randvar, property):
         """Times evaluation of the pdf, logpdf, cdf and logcdf."""
         try:
             if property == "pdf":
@@ -87,18 +87,18 @@ class Sampling:
     Benchmark sampling routines for various distributions.
     """
 
-    param_names = ["rv_name"]
+    param_names = ["randvar"]
     params = [rv_names]
 
-    def setup(self, rv_name):
+    def setup(self, randvar):
         np.random.seed(42)
         self.n_samples = 1000
-        self.randvar = get_randvar(rv_name=rv_name)
+        self.randvar = get_randvar(rv_name=randvar)
 
-    def time_sample(self, rv_name):
+    def time_sample(self, randvar):
         """Times sampling from this distribution."""
         self.randvar.sample(self.n_samples)
 
-    def peakmem_sample(self, rv_name):
+    def peakmem_sample(self, randvar):
         """Peak memory of sampling process."""
         self.randvar.sample(self.n_samples)
diff --git a/src/probnum/random_variables/_normal.py b/src/probnum/random_variables/_normal.py
index 776bbc2181..81f0baa839 100644
--- a/src/probnum/random_variables/_normal.py
+++ b/src/probnum/random_variables/_normal.py
@@ -3,6 +3,7 @@
 from typing import Callable, Optional, Union
 
 import numpy as np
+import scipy.linalg
 import scipy.stats
 
 from probnum import utils as _utils
@@ -464,8 +465,9 @@ def _univariate_entropy(self: _ValueType) -> np.float_:
     def _dense_cov_cholesky(self) -> np.ndarray:
         eps = 10 ** -12  # damping needed to avoid negative definite covariances
 
-        return np.linalg.cholesky(
-            self._dense_cov + eps * np.eye(self._dense_cov.shape[0], dtype=self.dtype)
+        return scipy.linalg.cholesky(
+            self._dense_cov + eps * np.eye(self._dense_cov.shape[0], dtype=self.dtype),
+            lower=True,
         )
 
     def _dense_sample(self, size: ShapeType = ()) -> np.ndarray:
@@ -538,8 +540,8 @@ def _kronecker_cov_cholesky(self) -> linops.Kronecker:
         eps = 10 ** -12  # damping needed to avoid negative definite covariances
 
         return linops.Kronecker(
-            A=np.linalg.cholesky(A + eps * np.eye(A.shape[0])),
-            B=np.linalg.cholesky(B + eps * np.eye(B.shape[0])),
+            A=scipy.linalg.cholesky(A + eps * np.eye(A.shape[0]), lower=True),
+            B=scipy.linalg.cholesky(B + eps * np.eye(B.shape[0]), lower=True),
             dtype=self.dtype,
         )
 
@@ -551,7 +553,9 @@ def _symmetric_kronecker_identical_factors_cov_cholesky(self) -> linops.Kronecke
         A = self._cov.A.todense()
 
         eps = 10 ** -12  # damping needed to avoid negative definite covariances
-        A_cholesky = np.linalg.cholesky(A + eps * np.eye(A.shape[0], dtype=self.dtype))
+        A_cholesky = scipy.linalg.cholesky(
+            A + eps * np.eye(A.shape[0], dtype=self.dtype), lower=True
+        )
 
         return linops.SymmetricKronecker(A=A_cholesky, dtype=self.dtype)