diff --git a/.github/workflows/CI_testing.yml b/.github/workflows/CI_testing.yml
index 70f69762..3fc778f3 100644
--- a/.github/workflows/CI_testing.yml
+++ b/.github/workflows/CI_testing.yml
@@ -35,6 +35,21 @@ jobs:
           cd ./hippylibX/test &&
           mpirun -n 2 python3 test_eigendecomposition.py
 
+      - name: low rank Hessian testing
+        run: |
+          cd ./hippylibX/test &&
+          mpirun -n 2 python3 test_lowRankHessian.py
+
+      - name: low rank Hessian preconditioner testing
+        run: |
+          cd ./hippylibX/test &&
+          mpirun -n 2 python3 test_lowRankHessian_preconditioner.py
+
+      - name: sample testing
+        run: |
+          cd ./hippylibX/test &&
+          mpirun -n 2 python3 test_sampling.py
+
       - name: run serial check
         run: |
           cd ./hippylibX/test && 
diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py
index af19e6d2..b62d3d69 100644
--- a/example/poisson_dirichlet_example.py
+++ b/example/poisson_dirichlet_example.py
@@ -10,7 +10,6 @@
 from matplotlib import pyplot as plt
 from typing import Sequence, Dict
 
-
 sys.path.append(os.environ.get("HIPPYLIBX_BASE_DIR", "../"))
 import hippylibX as hpx
 
@@ -106,7 +105,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     misfit_form = PoissonMisfitForm(d, noise_variance)
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form, [bc0])
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
     prior_mean = prior_mean.x
 
     prior = hpx.BiLaplacianPrior(
@@ -184,7 +183,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
 
     with dlx.io.VTXWriter(
         msh.comm,
-        "poisson_Dirichlet_BiLaplacian_prior_np{0:d}_Prior.bp".format(nproc),
+        "poisson_Dirichlet_BiLaplacian_prior_np{0:d}.bp".format(nproc),
         [m_fun, m_true_fun, u_map_fun, u_true_fun, d_fun],
     ) as vtx:
         vtx.write(0.0)
@@ -215,17 +214,48 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
 
     hpx.parRandom.normal(1.0, Omega)
 
-    d, U = hpx.doublePassG(
-        Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
-    )
+    d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)
 
-    eigen_decomposition_results = {
-        "A": Hmisfit.mat,
-        "B": prior.R,
-        "k": k,
-        "d": d,
-        "U": U,
-    }
+    # generating prior and posterior samples
+    lap_aprx = hpx.LaplaceApproximator(prior, d, U)
+    lap_aprx.mean = prior.generate_parameter(0)
+    lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:]
+
+    noise = prior.generate_parameter("noise")
+
+    num_samples_generate = 5
+    use_vtx = False
+    prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample")
+    posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample")
+    if use_vtx:
+        with dlx.io.VTXWriter(
+            msh.comm,
+            "poisson_Dirichlet_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc),
+            [prior_sample, posterior_sample],
+        ) as vtx:
+            for i in range(num_samples_generate):
+                hpx.parRandom.normal(1.0, noise)
+                lap_aprx.sample(noise, prior_sample.x, posterior_sample.x)
+                prior_sample.x.scatter_forward()
+                posterior_sample.x.scatter_forward()
+                vtx.write(float(i))
+    else:
+        ############################################
+        with dlx.io.XDMFFile(
+            msh.comm,
+            "poisson_Dirichlet_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc),
+            "w",
+        ) as file:
+            file.write_mesh(msh)
+            for i in range(num_samples_generate):
+                hpx.parRandom.normal(1.0, noise)
+                lap_aprx.sample(noise, prior_sample.x, posterior_sample.x)
+                prior_sample.x.scatter_forward()
+                posterior_sample.x.scatter_forward()
+                file.write_function(prior_sample, float(i))
+                file.write_function(posterior_sample, float(i))
+
+    eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
     final_results = {
         "data_misfit_True": data_misfit_True,
@@ -235,7 +265,6 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     }
 
     return final_results
-    #######################################
 
 
 if __name__ == "__main__":
@@ -243,7 +272,9 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     ny = 64
     noise_variance = 1e-4
     prior_param = {"gamma": 0.03, "delta": 0.3}
+
     final_results = run_inversion(nx, ny, noise_variance, prior_param)
+
     k, d = (
         final_results["eigen_decomposition_results"]["k"],
         final_results["eigen_decomposition_results"]["d"],
diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py
index b11d7439..50a8649b 100644
--- a/example/poisson_dirichlet_example_reg.py
+++ b/example/poisson_dirichlet_example_reg.py
@@ -109,7 +109,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form, [bc0])
 
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
 
     prior_gamma = prior_param["gamma"]
     prior_delta = prior_param["delta"]
@@ -225,17 +225,9 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
 
     hpx.parRandom.normal(1.0, Omega)
 
-    d, U = hpx.doublePassG(
-        Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
-    )
+    d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)
 
-    eigen_decomposition_results = {
-        "A": Hmisfit.mat,
-        "B": prior.R,
-        "k": k,
-        "d": d,
-        "U": U,
-    }
+    eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
     final_results = {
         "data_misfit_True": data_misfit_True,
diff --git a/example/poisson_example.py b/example/poisson_example.py
index 162d575d..01177d37 100644
--- a/example/poisson_example.py
+++ b/example/poisson_example.py
@@ -87,7 +87,7 @@ def run_inversion(
     misfit_form = PoissonMisfitForm(d, noise_variance)
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form)
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
     prior_mean = prior_mean.x
 
     prior = hpx.BiLaplacianPrior(
@@ -195,17 +195,48 @@ def run_inversion(
 
     hpx.parRandom.normal(1.0, Omega)
 
-    d, U = hpx.doublePassG(
-        Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
-    )
+    d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)
 
-    eigen_decomposition_results = {
-        "A": Hmisfit.mat,
-        "B": prior.R,
-        "k": k,
-        "d": d,
-        "U": U,
-    }
+    # generating prior and posterior samples
+    lap_aprx = hpx.LaplaceApproximator(prior, d, U)
+    lap_aprx.mean = prior.generate_parameter(0)
+    lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:]
+
+    noise = prior.generate_parameter("noise")
+
+    num_samples_generate = 5
+    use_vtx = False
+    prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample")
+    posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample")
+    if use_vtx:
+        with dlx.io.VTXWriter(
+            msh.comm,
+            "poisson_Robin_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc),
+            [prior_sample, posterior_sample],
+        ) as vtx:
+            for i in range(num_samples_generate):
+                hpx.parRandom.normal(1.0, noise)
+                lap_aprx.sample(noise, prior_sample.x, posterior_sample.x)
+                prior_sample.x.scatter_forward()
+                posterior_sample.x.scatter_forward()
+                vtx.write(float(i))
+    else:
+        ############################################
+        with dlx.io.XDMFFile(
+            msh.comm,
+            "poisson_Robin_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc),
+            "w",
+        ) as file:
+            file.write_mesh(msh)
+            for i in range(num_samples_generate):
+                hpx.parRandom.normal(1.0, noise)
+                lap_aprx.sample(noise, prior_sample.x, posterior_sample.x)
+                prior_sample.x.scatter_forward()
+                posterior_sample.x.scatter_forward()
+                file.write_function(prior_sample, float(i))
+                file.write_function(posterior_sample, float(i))
+
+    eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
     final_results = {
         "data_misfit_True": data_misfit_True,
@@ -222,9 +253,10 @@ def run_inversion(
 if __name__ == "__main__":
     nx = 64
     ny = 64
-    noise_variance = 1e-4
+    noise_variance = 1e-6
     prior_param = {"gamma": 0.07, "delta": 0.7}
     final_results = run_inversion(nx, ny, noise_variance, prior_param)
+
     k, d = (
         final_results["eigen_decomposition_results"]["k"],
         final_results["eigen_decomposition_results"]["d"],
diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py
index 0da1e975..a658b948 100644
--- a/example/poisson_example_reg.py
+++ b/example/poisson_example_reg.py
@@ -85,7 +85,7 @@ def run_inversion(
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form)
 
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
 
     prior_gamma = prior_param["gamma"]
     prior_delta = prior_param["delta"]
@@ -198,17 +198,9 @@ def run_inversion(
 
     hpx.parRandom.normal(1.0, Omega)
 
-    d, U = hpx.doublePassG(
-        Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
-    )
+    d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)
 
-    eigen_decomposition_results = {
-        "A": Hmisfit.mat,
-        "B": prior.R,
-        "k": k,
-        "d": d,
-        "U": U,
-    }
+    eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
     final_results = {
         "data_misfit_True": data_misfit_True,
@@ -224,7 +216,7 @@ def run_inversion(
 if __name__ == "__main__":
     nx = 64
     ny = 64
-    noise_variance = 1e-4
+    noise_variance = 1e-6
     prior_param = {"gamma": 0.02, "delta": 0.2}
     final_results = run_inversion(nx, ny, noise_variance, prior_param)
     k, d = (
diff --git a/example/sfsi_toy_gaussian.py b/example/sfsi_toy_gaussian.py
index 742dbbbd..e63bbb20 100644
--- a/example/sfsi_toy_gaussian.py
+++ b/example/sfsi_toy_gaussian.py
@@ -205,11 +205,6 @@ def run_inversion(
         optimizer_results["optimizer"] = True
     else:
         optimizer_results["optimizer"] = False
-    final_results = {
-        "data_misfit_True": data_misfit_True,
-        "data_misfit_False": data_misfit_False,
-        "optimizer_results": optimizer_results,
-    }
 
     Hmisfit = hpx.ReducedHessian(model, misfit_only=True)
 
@@ -228,17 +223,48 @@ def run_inversion(
 
     hpx.parRandom.normal(1.0, Omega)
 
-    d, U = hpx.doublePassG(
-        Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
-    )
+    d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)
 
-    eigen_decomposition_results = {
-        "A": Hmisfit.mat,
-        "B": prior.R,
-        "k": k,
-        "d": d,
-        "U": U,
-    }
+    # generating prior and posterior samples
+    lap_aprx = hpx.LaplaceApproximator(prior, d, U)
+    lap_aprx.mean = prior.generate_parameter(0)
+    lap_aprx.mean.array[:] = x[hpx.PARAMETER].array[:]
+
+    noise = prior.generate_parameter("noise")
+
+    num_samples_generate = 5
+    use_vtx = False
+    prior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="prior_sample")
+    posterior_sample = dlx.fem.Function(Vh[hpx.PARAMETER], name="posterior_sample")
+    if use_vtx:
+        with dlx.io.VTXWriter(
+            msh.comm,
+            "pact_prior_Bilaplacian_samples_np{0:d}.bp".format(nproc),
+            [prior_sample, posterior_sample],
+        ) as vtx:
+            for i in range(num_samples_generate):
+                hpx.parRandom.normal(1.0, noise)
+                lap_aprx.sample(noise, prior_sample.x, posterior_sample.x)
+                prior_sample.x.scatter_forward()
+                posterior_sample.x.scatter_forward()
+                vtx.write(float(i))
+    else:
+        ############################################
+        with dlx.io.XDMFFile(
+            msh.comm,
+            "pact_prior_Bilaplacian_samples_np{0:d}.xdmf".format(nproc),
+            "w",
+        ) as file:
+            file.write_mesh(msh)
+            for i in range(num_samples_generate):
+                hpx.parRandom.normal(1.0, noise)
+                lap_aprx.sample(noise, prior_sample.x, posterior_sample.x)
+                prior_sample.x.scatter_forward()
+                posterior_sample.x.scatter_forward()
+                file.write_function(prior_sample, float(i))
+                file.write_function(posterior_sample, float(i))
+
+    eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
     final_results = {
         "data_misfit_True": data_misfit_True,
@@ -248,7 +274,7 @@ def run_inversion(
     }
 
     return final_results
-    #######################################
+    ######################################
 
 
 if __name__ == "__main__":
diff --git a/example/sfsi_toy_gaussian_reg.py b/example/sfsi_toy_gaussian_reg.py
index 05e829d3..9a74503a 100644
--- a/example/sfsi_toy_gaussian_reg.py
+++ b/example/sfsi_toy_gaussian_reg.py
@@ -229,17 +229,9 @@ def run_inversion(
 
     hpx.parRandom.normal(1.0, Omega)
 
-    d, U = hpx.doublePassG(
-        Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1, check=False
-    )
+    d, U = hpx.doublePassG(Hmisfit.mat, prior.R, prior.Rsolver, Omega, k, s=1)
 
-    eigen_decomposition_results = {
-        "A": Hmisfit.mat,
-        "B": prior.R,
-        "k": k,
-        "d": d,
-        "U": U,
-    }
+    eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
     final_results = {
         "data_misfit_True": data_misfit_True,
diff --git a/hippylibX/algorithms/NewtonCG.py b/hippylibX/algorithms/NewtonCG.py
index be345a55..1406bde8 100644
--- a/hippylibX/algorithms/NewtonCG.py
+++ b/hippylibX/algorithms/NewtonCG.py
@@ -240,7 +240,7 @@ def _solve_ls(self, x: list) -> list:
             solver.parameters["print_level"] = print_level - 1
             mg_neg = self.model.generate_vector(PARAMETER)
             mg_neg.array[:] = -1 * mg.array[:]
-            solver.solve(mhat, mg_neg)
+            solver.solve(mg_neg, mhat)
             self.total_cg_iter += HessApply.ncalls
             alpha = 1.0
             descent = 0
@@ -374,7 +374,8 @@ def _solve_tr(self, x):
             solver.parameters["zero_initial_guess"] = True
             solver.parameters["print_level"] = print_level - 1
 
-            solver.solve(mhat, -mg)
+            # solver.solve(mhat, -mg)
+            solver.solve(-mg, mhat)
             self.total_cg_iter += HessApply.ncalls
             if self.it == 1:
                 self.model.prior.R.mult(mhat, R_mhat)
diff --git a/hippylibX/algorithms/__init__.py b/hippylibX/algorithms/__init__.py
index 7649993a..f198b04c 100644
--- a/hippylibX/algorithms/__init__.py
+++ b/hippylibX/algorithms/__init__.py
@@ -2,3 +2,4 @@
 from .NewtonCG import ReducedSpaceNewtonCG, ReducedSpaceNewtonCG_ParameterList  # noqa
 from .multivector import MultiVector, MatMvMult, MatMvTranspmult, MvDSmatMult  # noqa
 from .randomizedEigensolver import doublePassG  # noqa
+from .lowRankOperator import LowRankOperator  # noqa
diff --git a/hippylibX/algorithms/cgsolverSteihaug.py b/hippylibX/algorithms/cgsolverSteihaug.py
index 10720413..d378a50d 100644
--- a/hippylibX/algorithms/cgsolverSteihaug.py
+++ b/hippylibX/algorithms/cgsolverSteihaug.py
@@ -160,14 +160,14 @@ def update_x_with_TR(self, x, alpha, d):
             x.axpy(tau * alpha, d)
             return True
 
-    def solve(self, x: dlx.la.Vector, b: dlx.la.Vector) -> None:
+    def solve(self, b: dlx.la.Vector, x: dlx.la.Vector) -> None:
         temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
         temp_petsc_vec_b = dlx.la.create_petsc_vector_wrap(b)
-        self.solve_petsc(temp_petsc_vec_x, temp_petsc_vec_b)
+        self.solve_petsc(temp_petsc_vec_b, temp_petsc_vec_x)
         temp_petsc_vec_x.destroy()
         temp_petsc_vec_b.destroy()
 
-    def solve_petsc(self, x: petsc4py.PETSc.Vec, b: petsc4py.PETSc.Vec) -> None:
+    def solve_petsc(self, b: petsc4py.PETSc.Vec, x: petsc4py.PETSc.Vec) -> None:
         """
         Solve the linear system :math:`Ax = b`
         """
diff --git a/hippylibX/algorithms/lowRankOperator.py b/hippylibX/algorithms/lowRankOperator.py
new file mode 100644
index 00000000..803e7503
--- /dev/null
+++ b/hippylibX/algorithms/lowRankOperator.py
@@ -0,0 +1,91 @@
+# Copyright (c) 2016-2018, The University of Texas at Austin
+# & University of California--Merced.
+# Copyright (c) 2019-2020, The University of Texas at Austin
+# University of California--Merced, Washington University in St. Louis.
+#
+# All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the hIPPYlib library. For more information and source code
+# availability see https://hippylib.github.io.
+#
+# hIPPYlib is free software; you can redistribute it and/or modify it under the
+# terms of the GNU General Public License (as published by the Free
+# Software Foundation) version 2.0 dated June 1991.
+
+import petsc4py.PETSc
+from .multivector import MultiVector, MatMvMult
+import numpy as np
+import petsc4py
+
+
+class LowRankOperator:
+    """
+    This class model the action of a low rank operator :math:`A = U D U^T`.
+    Here :math:`D` is a diagonal matrix, and the columns of are orthonormal
+    in some weighted inner-product.
+    """
+
+    def __init__(
+        self, d: np.array, U: MultiVector, createVecLeft=None, createVecRight=None
+    ):
+        """
+        Construct the low rank operator given :code:`d` and :code:`U`.
+        """
+        self.d = d
+        self.U = U
+        self.createVecLeft = createVecLeft
+        self.createVecRight = createVecRight
+
+    def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None:
+        """
+        Compute :math:`y = Ax = U D U^T x`
+        """
+        Utx = self.U.dot(x)
+        dUtx = self.d * Utx  # elementwise mult
+        y.scale(0.0)
+        self.U.reduce(y, dUtx)
+
+    def solve(self, rhs: petsc4py.PETSc.Vec, sol: petsc4py.PETSc.Vec) -> None:
+        """
+        Compute :math:`\mbox{sol} = U D^-1 U^T x`
+        """
+        Utr = self.U.dot(rhs)
+        dinvUtr = Utr / self.d
+        sol.scale(0.0)
+        self.U.reduce(sol, dinvUtr)
+
+    def get_diagonal(self, diag: petsc4py.PETSc.Vec) -> None:
+        """
+        Compute the diagonal of :code:`A`.
+        """
+        diag.scale(0.0)
+        tmp = self.U[0].duplicate()
+        for i in range(self.U.nvec):
+            tmp.pointwiseMult(self.U[i], self.U[i])
+            diag.axpy(self.d[i], tmp)
+
+    def trace(self, W=None) -> float:
+        """
+        Compute the trace of :code:`A`.
+        If the weight :code:`W` is given, compute the trace of :math:`W^{1/2} A W^{1/2}`.
+        This is equivalent to :math:`\mbox{tr}_W(A) = \sum_i \lambda_i`,
+        where :math:`\lambda_i` are the generalized eigenvalues of
+        :math:`A x = \lambda W^{-1} x`.
+
+        .. note:: If :math:`U` is a :math:`W`-orthogonal matrix then :math:`\mbox{tr}_W(A) = \sum_i D(i,i)`.
+        """
+        if W is None:
+            tmp = self.U[0].duplicate()
+            tmp.scale(0.0)
+            self.U.reduce(tmp, np.sqrt(self.d))
+            tr = tmp.dot(tmp)
+        else:
+            WU = MultiVector.createFromVec(self.U[0], self.U.nvec)
+            MatMvMult(W, self.U, WU)
+            diagWUtU = np.zeros_like(self.d)
+            for i in range(self.d.shape[0]):
+                diagWUtU[i] = WU[i].dot(self.U[i])
+            tr = np.sum(self.d * diagWUtU)
+
+        return tr
diff --git a/hippylibX/algorithms/multivector.py b/hippylibX/algorithms/multivector.py
index 67828d1a..b732b150 100644
--- a/hippylibX/algorithms/multivector.py
+++ b/hippylibX/algorithms/multivector.py
@@ -56,18 +56,14 @@ def dot(self, v: Union[petsc4py.PETSc.Vec, Type["MultiVector"]]) -> np.array:
 
         return return_values
 
-    def reduce(self, alpha: np.array) -> petsc4py.PETSc.Vec:
-        return_vec = self[0].duplicate()
-        return_vec.scale(0.0)
+    def reduce(self, y: petsc4py.PETSc.Vec, alpha: np.array) -> None:
         for i in range(self.nvec):
-            return_vec.axpy(alpha[i], self[i])
-        return return_vec
+            y.axpy(alpha[i], self[i])
 
     def axpy(self, alpha: Union[float, np.array], Y: Type["MultiVector"]) -> None:
         if isinstance(alpha, float):
             for i in range(self.nvec):
                 self[i].axpy(alpha, Y[i])
-
         else:
             for i in range(self.nvec):
                 self[i].axpy(alpha[i], Y[i])
@@ -154,5 +150,4 @@ def MvDSmatMult(X: MultiVector, A: np.array, Y: MultiVector) -> None:
     ), "Y Number of vecs incompatible with number of cols in A"
     for j in range(Y.nvec):
         Y[j].scale(0.0)
-        reduced_vec = X.reduce(A[:, j].flatten())
-        Y.data[j] = X.data[0].duplicate(reduced_vec.getArray())
+        X.reduce(Y[j], A[:, j].flatten())
diff --git a/hippylibX/algorithms/randomizedEigensolver.py b/hippylibX/algorithms/randomizedEigensolver.py
index 5f4e8d24..4af6c539 100644
--- a/hippylibX/algorithms/randomizedEigensolver.py
+++ b/hippylibX/algorithms/randomizedEigensolver.py
@@ -13,7 +13,6 @@ def doublePassG(
     Omega: MultiVector,
     k: int,
     s=1,
-    check=False,
 ) -> tuple[np.array, MultiVector]:
     nvec = Omega.nvec
     assert nvec >= k
diff --git a/hippylibX/modeling/Regularization.py b/hippylibX/modeling/Regularization.py
index ffea59a2..37408d61 100644
--- a/hippylibX/modeling/Regularization.py
+++ b/hippylibX/modeling/Regularization.py
@@ -77,6 +77,9 @@ def __del__(self):
         self.Msolver.destroy()
         self.M.destroy()
 
+    def generate_parameter(self, dim: int) -> dlx.la.Vector:
+        return dlx.la.Vector(self.Vh.dofmap.index_map)
+
     def _createsolver(self, petsc_options: dict) -> petsc4py.PETSc.KSP:
         ksp = petsc4py.PETSc.KSP().create(self.Vh.mesh.comm)
         problem_prefix = f"dolfinx_solve_{id(self)}"
diff --git a/hippylibX/modeling/__init__.py b/hippylibX/modeling/__init__.py
index 4493c511..57a5493e 100644
--- a/hippylibX/modeling/__init__.py
+++ b/hippylibX/modeling/__init__.py
@@ -6,3 +6,8 @@
 from .Regularization import H1TikhonvFunctional, VariationalRegularization  # noqa
 from .variables import STATE, PARAMETER, ADJOINT, NVAR  # noqa
 from .reducedHessian import ReducedHessian  # noqa
+from .laplaceApproximation import (
+    LowRankHessian,  # noqa
+    LowRankPosteriorSampler,  # noqa
+    LaplaceApproximator,  # noqa
+)
diff --git a/hippylibX/modeling/laplaceApproximation.py b/hippylibX/modeling/laplaceApproximation.py
new file mode 100644
index 00000000..5e999c43
--- /dev/null
+++ b/hippylibX/modeling/laplaceApproximation.py
@@ -0,0 +1,262 @@
+# Copyright (c) 2016-2018, The University of Texas at Austin
+# & University of California--Merced.
+# Copyright (c) 2019-2020, The University of Texas at Austin
+# University of California--Merced, Washington University in St. Louis.
+#
+# All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the hIPPYlib library. For more information and source code
+# availability see https://hippylib.github.io.
+#
+# hIPPYlib is free software; you can redistribute it and/or modify it under the
+# terms of the GNU General Public License (as published by the Free
+# Software Foundation) version 2.0 dated June 1991.
+
+import petsc4py.PETSc
+from ..algorithms.lowRankOperator import LowRankOperator
+from ..algorithms.multivector import MultiVector
+import numpy as np
+import petsc4py
+import dolfinx as dlx
+
+
+def not_implemented(func):
+    def wrapper(*args, **kwargs):
+        raise NotImplementedError(f"{func.__name__} Function has not been implemented.")
+
+    return wrapper
+
+
+class LowRankHessian:
+    """
+    Operator that represents the action of the low rank approximation
+    of the Hessian and of its inverse.
+    """
+
+    def __init__(self, prior, d: np.array, U: MultiVector):
+        self.U = U
+        self.prior = prior
+        self.LowRankH = LowRankOperator(d, self.U)
+        dsolve = d / (np.ones(d.shape, dtype=d.dtype) + d)
+        self.LowRankHinv = LowRankOperator(dsolve, self.U)
+
+        self.help = self.prior.R.createVecRight()
+        self.help1 = self.prior.R.createVecLeft()
+
+    def __del__(self) -> None:
+        self.help.destroy()
+        self.help1.destroy()
+
+    def createVecRight(self) -> petsc4py.PETSc.Vec:
+        return self.prior.R.createVecRight()
+
+    def createVecLeft(self) -> petsc4py.PETSc.Vec:
+        return self.prior.R.createVecLeft()
+
+    def mult(self, x: petsc4py.PETSc.Vec, y: petsc4py.PETSc.Vec) -> None:
+        self.prior.R.mult(x, y)
+        self.LowRankH.mult(y, self.help)
+        self.prior.R.mult(self.help, self.help1)
+        y.axpy(1, self.help1)
+
+    def solve(self, rhs: petsc4py.PETSc.Vec, sol: petsc4py.PETSc.Vec) -> None:
+        self.prior.Rsolver.solve(rhs, sol)
+        self.LowRankHinv.mult(rhs, self.help)
+        sol.axpy(-1, self.help)
+
+
+class LowRankPosteriorSampler:
+    """
+    Object to sample from the low rank approximation
+    of the posterior.
+
+        .. math:: y = ( I - U S U^{T}) x,
+
+    where
+
+    :math:`S = I - (I + D)^{-1/2}, x \\sim \\mathcal{N}(0, R^{-1}).`
+    """
+
+    def __init__(
+        self,
+        prior,
+        d: np.array,
+        U: MultiVector,
+    ):
+        self.U = U
+        self.prior = prior
+
+        ones = np.ones(d.shape, dtype=d.dtype)
+        self.d = ones - np.power(ones + d, -0.5)
+        self.lrsqrt = LowRankOperator(self.d, self.U)
+        self.help = self.prior.R.createVecLeft()
+
+    def __del__(self) -> None:
+        self.help.destroy()
+
+    def createVecRight(self) -> petsc4py.PETSc.Vec:
+        return self.prior.R.createVecRight()
+
+    def createVecLeft(self) -> petsc4py.PETSc.Vec:
+        return self.prior.R.createVecLeft()
+
+    def mult(self, noise: dlx.la.Vector, s: dlx.la.Vector):
+        temp_petsc_vec_noise = dlx.la.create_petsc_vector_wrap(noise)
+        temp_petsc_vec_s = dlx.la.create_petsc_vector_wrap(s)
+
+        self.prior.R.mult(temp_petsc_vec_noise, self.help)
+        self.lrsqrt.mult(self.help, temp_petsc_vec_s)
+        temp_petsc_vec_s.axpby(-1.0, 1.0, temp_petsc_vec_noise)
+        temp_petsc_vec_noise.destroy()
+        temp_petsc_vec_s.destroy()
+
+
+class LaplaceApproximator:
+    """
+    Class for the low rank Gaussian Approximation of the Posterior.
+    This class provides functionality for approximate Hessian
+    apply, solve, and Gaussian sampling based on the low rank
+    factorization of the Hessian.
+
+    In particular if :math:`d` and :math:`U` are the dominant eigenpairs of
+    :math:`H_{\\mbox{misfit}} U[:,i] = d[i] R U[:,i]`
+    then we have:
+
+    - low rank Hessian apply: :math:`y = ( R + RU D U^{T}) x.`
+    - low rank Hessian solve: :math:`y = (R^-1 - U (I + D^{-1})^{-1} U^T) x.`
+    - low rank Hessian Gaussian sampling: :math:`y = ( I - U S U^{T}) x`, where :math:`S = I - (I + D)^{-1/2}` and :math:`x \\sim \\mathcal{N}(0, R^{-1}).`
+    """
+
+    def __init__(self, prior, d: np.array, U: MultiVector, mean=None):
+        """
+        Construct the Gaussian approximation of the posterior.
+        Input:
+        - :code:`prior`: the prior mode.
+        - :code:`d`:     the dominant generalized eigenvalues of the Hessian misfit.
+        - :code:`U`:     the dominant generalized eigenvector of the Hessian misfit :math:`U^T R U = I.`
+        - :code:`mean`:  the MAP point.
+        """
+        self.prior = prior
+        self.d = d
+        self.U = U
+        self.Hlr = LowRankHessian(prior, d, U)
+        self.sampler = LowRankPosteriorSampler(self.prior, self.d, self.U)
+        self.mean = mean
+
+        if self.mean is None:
+            self.mean = self.prior.generate_parameter(0)
+
+    def cost(self, m: dlx.la.Vector) -> float:
+        if self.mean is None:
+            dm = m
+        else:
+            dm = m - self.mean
+        temp_petsc_vec_dm = dlx.la.create_petsc_vector_wrap(dm)
+        self.Hlr.mult(temp_petsc_vec_dm, self.Hlr.help1)
+        return_value = 0.5 * self.Hlr.help1.dot(temp_petsc_vec_dm)
+        temp_petsc_vec_dm.destroy()
+        return return_value
+
+    def sample(self, *args, **kwargs):
+        """
+        possible calls:
+
+        1) :code:`sample(s_prior, s_post, add_mean=True)`
+
+           Given a prior sample  :code:`s_prior` compute a sample :code:`s_post` from the posterior.
+
+           - :code:`s_prior` is a sample from the prior centered at 0 (input).
+           - :code:`s_post` is a sample from the posterior (output).
+           - if :code:`add_mean=True` (default) then the samples will be centered at the map point.
+
+        2) :code:`sample(noise, s_prior, s_post, add_mean=True)`
+
+           Given :code:`noise` :math:`\\sim \\mathcal{N}(0, I)` compute a sample :code:`s_prior` from the prior and
+           :code:`s_post` from the posterior.
+
+           - :code:`noise` is a realization of white noise (input).
+           - :code:`s_prior` is a sample from the prior (output).
+           - :code:`s_post`  is a sample from the posterior.
+           - if :code:`add_mean=True` (default) then the prior and posterior samples will be centered at the respective means.
+
+        """
+        add_mean = True
+        for name, value in kwargs.items():
+            if name == "add_mean":
+                add_mean = value
+            else:
+                raise NameError(name)
+
+        if len(args) == 2:
+            self._sample_given_prior(args[0], args[1])
+            if add_mean:
+                args[1].array[:] += self.mean.array
+
+        elif len(args) == 3:
+            self._sample_given_white_noise(args[0], args[1], args[2])
+            if add_mean:
+                args[1].array[:] += self.prior.mean.array
+                args[2].array[:] += self.mean.array
+        else:
+            raise NameError("Invalid number of parameters in Posterior::sample")
+
+    def _sample_given_white_noise(
+        self, noise: dlx.la.Vector, s_prior: dlx.la.Vector, s_post: dlx.la.Vector
+    ):
+        self.prior.sample(noise, s_prior, add_mean=False)
+        self.sampler.mult(s_prior, s_post)
+
+    def _sample_given_prior(self, s_prior: dlx.la.Vector, s_post: dlx.la.Vector):
+        self.sampler.mult(s_prior, s_post)
+
+    @not_implemented
+    def trace(self, **kwargs):
+        """
+        Compute/estimate the trace of the posterior, prior distribution
+        and the trace of the data informed correction.
+
+        See :code:`_Prior.trace` for more details.
+        """
+        pr_trace = self.prior.trace(**kwargs)
+        corr_trace = self.trace_update()
+        post_trace = pr_trace - corr_trace
+        return post_trace, pr_trace, corr_trace
+
+    @not_implemented
+    def trace_update(self):
+        return self.Hlr.LowRankHinv.trace(self.prior.M)
+
+    @not_implemented
+    def pointwise_variance(self, **kwargs):
+        """
+        Compute/estimate the pointwise variance of the posterior, prior distribution
+        and the pointwise variance reduction informed by the data.
+
+        See :code:`_Prior.pointwise_variance` for more details.
+        """
+        pr_pointwise_variance = self.prior.pointwise_variance(**kwargs)
+        # correction_pointwise_variance = Vector(self.prior.R.mpi_comm())
+        correction_pointwise_variance = None
+        self.init_vector(correction_pointwise_variance, 0)
+        self.Hlr.LowRankHinv.get_diagonal(correction_pointwise_variance)
+        post_pointwise_variance = pr_pointwise_variance - correction_pointwise_variance
+        return (
+            post_pointwise_variance,
+            pr_pointwise_variance,
+            correction_pointwise_variance,
+        )
+
+    def klDistanceFromPrior(self, sub_comp=False):
+        dplus1 = self.d + np.ones_like(self.d)
+
+        c_logdet = 0.5 * np.sum(np.log(dplus1))
+        c_trace = -0.5 * np.sum(self.d / dplus1)
+        c_shift = self.prior.cost(self.mean)
+
+        kld = c_logdet + c_trace + c_shift
+
+        if sub_comp:
+            return kld, c_logdet, c_trace, c_shift
+        else:
+            return kld
diff --git a/hippylibX/test/test_eigendecomposition.py b/hippylibX/test/test_eigendecomposition.py
index c05dbd4a..be7553c5 100644
--- a/hippylibX/test/test_eigendecomposition.py
+++ b/hippylibX/test/test_eigendecomposition.py
@@ -83,7 +83,7 @@ def test_eigendecomposition(self):
             out["eigen_decomposition_results"]["d"],
             out["eigen_decomposition_results"]["U"],
         )
-        result = check_g(A, B, U, d)
+        result = check_g(A.mat, B.R, U, d)
         check_output(self, result)
 
 
diff --git a/hippylibX/test/test_lowRankHessian.py b/hippylibX/test/test_lowRankHessian.py
new file mode 100644
index 00000000..47b7c91c
--- /dev/null
+++ b/hippylibX/test/test_lowRankHessian.py
@@ -0,0 +1,96 @@
+# Test to check if Hlr^(-1).Hlr == Hlr.Hlr(-1) == Identity
+# Hlr -> hpx.LowRankHessian
+
+# x ,y z as petsc4py Vectors:
+# y = Hlr(x)
+# z = Hlr^(-1)y => z == x
+
+# y = Hlr(-1)x
+# z = Hlr(y) => z == x
+
+import unittest
+import sys
+import os
+import numpy as np
+import dolfinx as dlx
+import dolfinx.fem.petsc
+import petsc4py
+from typing import Any
+
+sys.path.append(os.path.abspath("../.."))
+
+import hippylibX as hpx
+
+sys.path.append(os.path.abspath("../../example"))
+from example import poisson_dirichlet_example
+
+
+def low_Rank_Hessian_mult_solve(
+    prior: Any, d: np.array, U: hpx.MultiVector
+) -> tuple[float, float]:
+    Hlr = hpx.LowRankHessian(prior, d, U)
+    vec1 = dlx.la.vector(prior.Vh.dofmap.index_map)
+    vec2 = dlx.la.vector(prior.Vh.dofmap.index_map)
+    vec3 = dlx.la.vector(prior.Vh.dofmap.index_map)
+
+    hpx.parRandom.normal(1.0, vec1)
+    temp_petsc_vec1 = dlx.la.create_petsc_vector_wrap(vec1)
+    temp_petsc_vec2 = dlx.la.create_petsc_vector_wrap(vec2)
+    temp_petsc_vec3 = dlx.la.create_petsc_vector_wrap(vec3)
+
+    Hlr.mult(temp_petsc_vec1, temp_petsc_vec2)
+    Hlr.solve(temp_petsc_vec2, temp_petsc_vec3)
+    temp_petsc_vec3.axpy(-1.0, temp_petsc_vec1)
+
+    value1 = temp_petsc_vec3.norm(petsc4py.PETSc.NormType.N2)
+
+    hpx.parRandom.normal(1.0, vec1)
+
+    Hlr.solve(temp_petsc_vec1, temp_petsc_vec2)
+    Hlr.mult(temp_petsc_vec2, temp_petsc_vec3)
+    temp_petsc_vec3.axpy(-1.0, temp_petsc_vec1)
+
+    value2 = temp_petsc_vec3.norm(petsc4py.PETSc.NormType.N2)
+
+    temp_petsc_vec1.destroy()
+    temp_petsc_vec2.destroy()
+    temp_petsc_vec3.destroy()
+
+    return value1, value2
+
+
+def check_output(self, result: tuple[float, float]):
+    self.assertLessEqual(
+        np.abs(result[0]),
+        1e-6,
+        "Hlr(-1).Hlr not equal to identity",
+    )
+
+    self.assertLessEqual(
+        np.abs(result[1]),
+        1e-6,
+        "Hlr.Hlr(-1) not equal to identity",
+    )
+
+
+class Testing_Execution(unittest.TestCase):
+    def test_lowRankHessian(self):
+        hpx.parRandom.replay()
+        nx = 64
+        ny = 64
+        noise_variance = 1e-4
+        prior_param = {"gamma": 0.1, "delta": 1.0}
+        out = poisson_dirichlet_example.run_inversion(
+            nx, ny, noise_variance, prior_param
+        )
+        prior, d, U = (
+            out["eigen_decomposition_results"]["B"],
+            out["eigen_decomposition_results"]["d"],
+            out["eigen_decomposition_results"]["U"],
+        )
+        result = low_Rank_Hessian_mult_solve(prior, d, U)
+        check_output(self, result)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/hippylibX/test/test_lowRankHessian_preconditioner.py b/hippylibX/test/test_lowRankHessian_preconditioner.py
new file mode 100644
index 00000000..2974ca1c
--- /dev/null
+++ b/hippylibX/test/test_lowRankHessian_preconditioner.py
@@ -0,0 +1,64 @@
+# Low rank Hessian should be a very good preconditioner for the data
+# misfit Hessian.
+
+import unittest
+import sys
+import os
+import dolfinx as dlx
+
+sys.path.append(os.path.abspath("../.."))
+
+import hippylibX as hpx
+
+sys.path.append(os.path.abspath("../../example"))
+from example import poisson_dirichlet_example
+
+
+class Testing_Execution(unittest.TestCase):
+    def test_lowRankHessian_precondtioner(self):
+        hpx.parRandom.replay()
+        nx = 64
+        ny = 64
+        noise_variance = 1e-4
+        prior_param = {"gamma": 0.1, "delta": 1.0}
+        out = poisson_dirichlet_example.run_inversion(
+            nx, ny, noise_variance, prior_param
+        )
+        Hmisfit, prior, d, U = (
+            out["eigen_decomposition_results"]["A"],
+            out["eigen_decomposition_results"]["B"],
+            out["eigen_decomposition_results"]["d"],
+            out["eigen_decomposition_results"]["U"],
+        )
+        lowRank_Hessian = hpx.LowRankHessian(prior, d, U)
+
+        Hmisfit.misfit_only = False
+        matr = Hmisfit.mat
+
+        vec1 = dlx.la.vector(prior.Vh.dofmap.index_map)
+        hpx.parRandom.normal(1.0, vec1)
+
+        vec2 = dlx.la.vector(prior.Vh.dofmap.index_map)
+        parameters = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug_ParameterList()
+        solver = hpx.algorithms.cgsolverSteihaug.CGSolverSteihaug(
+            parameters, prior.Vh.mesh.comm
+        )
+        solver.set_operator(matr)
+
+        solver.set_preconditioner(prior.Rsolver)
+        solver.solve(vec1, vec2)
+        Rsolver_iter_count = solver.iter
+
+        solver.set_preconditioner(lowRank_Hessian)
+        solver.solve(vec1, vec2)
+        Hlr_iter_count = solver.iter
+
+        self.assertLess(
+            Hlr_iter_count,
+            Rsolver_iter_count,
+            "lowRank Hessian not better pc than prior.Rsolver",
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/hippylibX/test/test_sampling.py b/hippylibX/test/test_sampling.py
new file mode 100644
index 00000000..fdf427f7
--- /dev/null
+++ b/hippylibX/test/test_sampling.py
@@ -0,0 +1,82 @@
+# testing the lowRankPosteriorSampler.sample method from modeling/laplaceApproximation.py
+# x : random vector
+# prior.Rsolver.solve(x, x1)
+# Hlrsampler.sample(x1, x2)
+# Hlrsampler.sample(x2, x3)
+# Hlr.solve(x, y)
+# assert x3 == y
+
+import unittest
+import sys
+import os
+import dolfinx as dlx
+import petsc4py
+import dolfinx.fem.petsc
+import numpy as np
+
+sys.path.append(os.path.abspath("../.."))
+
+import hippylibX as hpx
+
+sys.path.append(os.path.abspath("../../example"))
+from example import poisson_dirichlet_example
+
+
+class Testing_Execution(unittest.TestCase):
+    def test_sampling(self):
+        hpx.parRandom.replay()
+
+        nx = 64
+        ny = 64
+        noise_variance = 1e-4
+        prior_param = {"gamma": 0.1, "delta": 1.0}
+        out = poisson_dirichlet_example.run_inversion(
+            nx, ny, noise_variance, prior_param
+        )
+        prior, d, U = (
+            out["eigen_decomposition_results"]["B"],
+            out["eigen_decomposition_results"]["d"],
+            out["eigen_decomposition_results"]["U"],
+        )
+
+        x = dlx.la.vector(prior.Vh.dofmap.index_map)
+        x1 = dlx.la.vector(prior.Vh.dofmap.index_map)
+        x2 = dlx.la.vector(prior.Vh.dofmap.index_map)
+        x3 = dlx.la.vector(prior.Vh.dofmap.index_map)
+        y = dlx.la.vector(prior.Vh.dofmap.index_map)
+
+        hpx.parRandom.normal(1.0, x)
+
+        temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
+        temp_petsc_vec_x1 = dlx.la.create_petsc_vector_wrap(x1)
+        prior.Rsolver.solve(temp_petsc_vec_x, temp_petsc_vec_x1)
+
+        PS_lr = hpx.LowRankPosteriorSampler(prior, d, U)
+        PS_lr.mult(x1, x2)
+        PS_lr.mult(x2, x3)
+
+        Hlr = hpx.LowRankHessian(prior, d, U)
+
+        temp_petsc_vec_y = dlx.la.create_petsc_vector_wrap(y)
+        Hlr.solve(temp_petsc_vec_x, temp_petsc_vec_y)
+
+        # assert x3 == y
+        temp_petsc_vec_x3 = dlx.la.create_petsc_vector_wrap(x3)
+        temp_petsc_vec_y.axpy(-1.0, temp_petsc_vec_x3)
+
+        value = temp_petsc_vec_y.norm(petsc4py.PETSc.NormType.N2)
+
+        temp_petsc_vec_x.destroy()
+        temp_petsc_vec_x1.destroy()
+        temp_petsc_vec_x3.destroy()
+        temp_petsc_vec_y.destroy()
+
+        self.assertLessEqual(
+            np.abs(value),
+            1e-6,
+            "lowRankPosteriorSampler.sample failed",
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()