diff --git a/superdsm/dsm.py b/superdsm/dsm.py
index 3269a92c..b6856b5b 100644
--- a/superdsm/dsm.py
+++ b/superdsm/dsm.py
@@ -259,11 +259,9 @@ class Energy:
     :param epsilon: Corresponds to the constant :math:`\\epsilon` which is used for the smooth approximation of the regularization term :math:`\\|\\xi\\|_1 \\approx \\mathbb 1^\\top_\\Omega \\sqrt{\\xi^2 + \\epsilon} - \\sqrt{\\epsilon} \\cdot \\#\\Omega` (see Supplemental Material 2 of :ref:`Kostrykin and Rohr, TPAMI 2023 <references>`).
     :param alpha: Governs the regularization of the deformations and corresponds to :math:`\\alpha` described in :ref:`pipeline_theory_cvxprog`. Increasing this value leads to a smoother segmentation result.
     :param smooth_matrix_factory: An object with a ``get`` method which yields the matrix :math:`\\tilde G_\\omega` for any image region :math:`\\omega` (represented as a binary mask and passed as a parameter).
-    :param sparsity_tol: Absolute values below this threshold will be treated as zeros for computation of the gradient.
-    :param hessian_sparsity_tol: Absolute values below this threshold will be treated as zeros for computation of the Hessian.
     """
 
-    def __init__(self, roi, epsilon, alpha, smooth_matrix_factory, sparsity_tol=0, hessian_sparsity_tol=0):
+    def __init__(self, roi, epsilon, alpha, smooth_matrix_factory):
         self.roi = roi
         self.p   = None
 
@@ -279,12 +277,6 @@ def __init__(self, roi, epsilon, alpha, smooth_matrix_factory, sparsity_tol=0, h
         assert alpha >= 0, 'alpha must be positive'
         self.alpha = alpha
 
-        assert sparsity_tol >= 0, 'sparsity_tol must be positive'
-        self.sparsity_tol = sparsity_tol
-
-        assert hessian_sparsity_tol >= 0, 'hessian_sparsity_tol must be positive'
-        self.hessian_sparsity_tol = hessian_sparsity_tol
-
         # pre-compute common terms occuring in the computation of the derivatives
         self.q = _compute_polynomial_derivatives(self.x)
     
@@ -343,7 +335,6 @@ def grad(self, params):
         self._update_theta()
         term1 = -self.y * self.theta
         grad = np.asarray([term1 * q for q in self.q]) @ self.w
-        term1[abs(term1) < self.sparsity_tol] = 0
         term1_sparse = coo_matrix(term1).transpose(copy=False)
         if self.smooth_mat.shape[1] > 0:
             grad2  = (self.w.reshape(-1)[None, :] @ self.smooth_mat.multiply(term1_sparse)).reshape(-1)
@@ -360,7 +351,6 @@ def hessian(self, params):
         self._update_maps(params)
         self._update_theta()
         kappa = self.theta - np.square(self.theta)
-        kappa[kappa < self.sparsity_tol] = 0
         pixelmask = (kappa != 0)
         term4 = np.sqrt(kappa[pixelmask] * self.w[pixelmask])[None, :]
         D1 = np.asarray([-self.y * qi for qi in self.q])[:, pixelmask] * term4
@@ -373,13 +363,6 @@ def hessian(self, params):
             assert np.allclose(0, g[g < 0])
             g[g < 0] = 0
             H += sparse_diag(np.concatenate([np.zeros(6), g]))
-            if self.hessian_sparsity_tol > 0:
-                H = H.tocoo()
-                H_mask = (np.abs(H.data) >= self.hessian_sparsity_tol)
-                H_mask = np.logical_or(H_mask, H.row == H.col)
-                H.data = H.data[H_mask]
-                H.row  = H.row [H_mask]
-                H.col  = H.col [H_mask]
         else:
             H = D1 @ D1.T
         return H
diff --git a/superdsm/dsmcfg.py b/superdsm/dsmcfg.py
index 1efeb3be..e9d32314 100644
--- a/superdsm/dsmcfg.py
+++ b/superdsm/dsmcfg.py
@@ -6,7 +6,6 @@
 DSM_CONFIG_DEFAULTS = {
     'cachesize': 1,
     'cachetest': None,
-    'sparsity_tol': 0,
     'init': 'elliptical',
     'smooth_amount': 10,
     'epsilon': 1.0,
@@ -37,9 +36,6 @@ class DSM_Config(Stage):
     ``dsm/cachetest``
         The test function to be used for cache testing. If ``None``, then ``numpy.array_equal`` will be used. Using other functions like ``numpy.allclose`` has shown to introduce numerical instabilities. Defaults to ``None``.
 
-    ``dsm/sparsity_tol``
-        Absolute values below this threshold will be treated as zeros during optimization. Defaults to 0.
-
     ``dsm/init``
         Either a function or a string. If this is function, then it will be called to determine the initialization, and the dimension of the vector :math:`\\xi` will be passed as a parameter. If this is a string, then the initialization corresponds to the result of convex programming using elliptical models (if set to ``elliptical``, see Supplemental Material 6 of :ref:`Kostrykin and Rohr, TPAMI 2023 <references>`) or a zeros vector of is used (otherwise). Defaults to ``elliptical``.
 
diff --git a/superdsm/objects.py b/superdsm/objects.py
index 113d0577..0726ac27 100644
--- a/superdsm/objects.py
+++ b/superdsm/objects.py
@@ -358,7 +358,7 @@ def _compute_elliptical_solution(J_elliptical, CP_params):
     return solution_array
 
 
-def cvxprog(region, scale, epsilon, alpha, smooth_amount, smooth_subsample, gaussian_shape_multiplier, smooth_mat_allocation_lock, smooth_mat_dtype, sparsity_tol=0, hessian_sparsity_tol=0, init=None, cachesize=0, cachetest=None, cp_timeout=None):
+def cvxprog(region, scale, epsilon, alpha, smooth_amount, smooth_subsample, gaussian_shape_multiplier, smooth_mat_allocation_lock, smooth_mat_dtype, init=None, cachesize=0, cachetest=None, cp_timeout=None):
     """Fits a deformable shape model to the intensities of an image region.
     
     Performs convex programming in an image region :math:`X` to determine the value of the set energy function :math:`\\nu(X)` and the optimal parameters :math:`\\theta` and :math:`\\xi` (see :ref:`pipeline_theory_cvxprog` and :ref:`pipeline_theory_jointsegandclustersplit`).
@@ -376,7 +376,7 @@ def cvxprog(region, scale, epsilon, alpha, smooth_amount, smooth_subsample, gaus
     """
     _print_heading('initializing')
     smooth_matrix_factory = SmoothMatrixFactory(smooth_amount, gaussian_shape_multiplier, smooth_subsample, smooth_mat_allocation_lock, smooth_mat_dtype)
-    J = Energy(region, epsilon, alpha, smooth_matrix_factory, sparsity_tol, hessian_sparsity_tol)
+    J = Energy(region, epsilon, alpha, smooth_matrix_factory)
     CP_params = {'cachesize': cachesize, 'cachetest': cachetest, 'scale': scale / J.smooth_mat.shape[0], 'timeout': cp_timeout}
     print(f'scale: {CP_params["scale"]:g}')
     print(f'region: {str(region.model.shape)}, offset: {str(region.offset)}')
diff --git a/superdsm/version.py b/superdsm/version.py
index 02f26d92..f7abe76c 100644
--- a/superdsm/version.py
+++ b/superdsm/version.py
@@ -1,5 +1,5 @@
 VERSION_MAJOR = 0
-VERSION_MINOR = 3
+VERSION_MINOR = 4
 VERSION_PATCH = 0
 
 VERSION = '%d.%d%s' % (VERSION_MAJOR, VERSION_MINOR, '.%d' % VERSION_PATCH if VERSION_PATCH > 0 else '')