diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index f4060b4d..4fc1e90d 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -26,8 +26,14 @@ of the maintenance releases, please take a look at
 
 * :py:func:`cashocs.import_mesh` can now also directly import a Gmsh mesh file. Internally, the mesh is directly converted to xdmf and then read. At the moment, this only supports the conversion mode `physical`.
 
+* Add a new parameter :python:`min_inner_iter` to :py:func:`cashocs.newton_solve`. This parameter is used to specify the minimum amount of iterations that the inner (i.e. linear) solver should perform, regardless of termination criteria.
+
 * New configuration file parameters:
 
+  * Section StateSystem
+
+    * :ini:`min_inner_iter` is a parameter controlling the minimum number of iterations that a inner (i.e. linear) solver should perform when using Newton's method. Mainly intended to be used for inexact Newton's method.
+
   * Section LineSearch
 
     * :ini:`fail_if_not_converged` determines, whether the line search is cancelled once the state system cannot be solved or if a new iterate is tried instead.
diff --git a/cashocs/_pde_problems/state_problem.py b/cashocs/_pde_problems/state_problem.py
index 223d47c2..c640c7a7 100755
--- a/cashocs/_pde_problems/state_problem.py
+++ b/cashocs/_pde_problems/state_problem.py
@@ -66,6 +66,7 @@ def __init__(
         self.newton_atol = self.config.getfloat("StateSystem", "newton_atol")
         self.newton_damped = self.config.getboolean("StateSystem", "newton_damped")
         self.newton_inexact = self.config.getboolean("StateSystem", "newton_inexact")
+        self.min_inner_iter = self.config.getint("StateSystem", "min_inner_iter")
         self.newton_verbose = self.config.getboolean("StateSystem", "newton_verbose")
         self.newton_iter = self.config.getint("StateSystem", "newton_iter")
 
@@ -149,6 +150,7 @@ def solve(self) -> List[fenics.Function]:
                             A_tensor=self.A_tensors[i],
                             b_tensor=self.b_tensors[i],
                             preconditioner_form=self.db.form_db.preconditioner_forms[i],
+                            min_inner_iter=self.min_inner_iter,
                         )
 
             else:
diff --git a/cashocs/_utils/linalg.py b/cashocs/_utils/linalg.py
index 89ca05ec..e2c24c7e 100644
--- a/cashocs/_utils/linalg.py
+++ b/cashocs/_utils/linalg.py
@@ -24,10 +24,13 @@
 
 import fenics
 import numpy as np
+from packaging import version
+import petsc4py
 from petsc4py import PETSc
 import ufl
 
 from cashocs import _exceptions
+from cashocs import _loggers
 from cashocs._utils import forms as forms_module
 
 if TYPE_CHECKING:
@@ -323,6 +326,7 @@ def solve_linear_problem(
     atol: Optional[float] = None,
     comm: Optional[MPI.Comm] = None,
     P: Optional[PETSc.Mat] = None,  # pylint: disable=invalid-name
+    min_iter: int = 0,
 ) -> PETSc.Vec:
     """Solves a finite dimensional linear problem.
 
@@ -344,6 +348,8 @@ def solve_linear_problem(
             solving the linear problem. Overrides the specification in the ksp object
             and ksp_options.
         comm: The MPI communicator for the problem.
+        min_iter: The minimum number of iterations to use for the linear solver,
+            regardless of tolerances.
 
     Returns:
         The solution vector.
@@ -365,6 +371,14 @@ def solve_linear_problem(
     options = define_ksp_options(ksp_options)
 
     setup_fieldsplit_preconditioner(fun, ksp, options)
+    if min_iter > 0:
+        if version.parse(petsc4py.__version__) < version.parse("3.20"):
+            _loggers.warning(
+                "The ability to specify a minimum number of iterations for a linear "
+                "solver is only available with PETSc >= 3.20. The specified minimum "
+                "iteration number will be ignored."
+            )
+        options["ksp_min_it"] = min_iter
     setup_petsc_options([ksp], [options])
 
     if rtol is not None:
diff --git a/cashocs/io/config.py b/cashocs/io/config.py
index 2dd5928c..b6f1a8af 100644
--- a/cashocs/io/config.py
+++ b/cashocs/io/config.py
@@ -129,6 +129,9 @@ def __init__(self, config_file: Optional[str] = None) -> None:
                 "newton_inexact": {
                     "type": "bool",
                 },
+                "min_inner_iter": {
+                    "type": "int",
+                },
                 "newton_verbose": {
                     "type": "bool",
                 },
@@ -542,6 +545,7 @@ def __init__(self, config_file: Optional[str] = None) -> None:
 newton_iter = 50
 newton_damped = False
 newton_inexact = False
+min_inner_iter = 0
 newton_verbose = False
 picard_iteration = False
 picard_rtol = 1e-10
diff --git a/cashocs/nonlinear_solvers/newton_solver.py b/cashocs/nonlinear_solvers/newton_solver.py
index 17a7cd16..656cb7f3 100644
--- a/cashocs/nonlinear_solvers/newton_solver.py
+++ b/cashocs/nonlinear_solvers/newton_solver.py
@@ -57,6 +57,7 @@ def __init__(
         b_tensor: Optional[fenics.PETScVector] = None,
         is_linear: bool = False,
         preconditioner_form: ufl.Form = None,
+        min_inner_iter: int = 0,
     ) -> None:
         r"""Initializes self.
 
@@ -95,6 +96,10 @@ def __init__(
             is_linear: A boolean flag, which indicates whether the problem is actually
                 linear.
             preconditioner_form: A UFL form which defines the preconditioner matrix.
+            min_inner_iter: Minimum number of iterations for the inner (linear) solver
+                to perform. This parameter is only intended for being used with an
+                inexact Newton's method (there, in the first steps, only a small handful
+                of inner iterations are needed to satisfy the termination criterion).
 
         """
         self.nonlinear_form = nonlinear_form
@@ -116,6 +121,8 @@ def __init__(
         else:
             self.preconditioner_form = preconditioner_form
 
+        self.min_inner_iter = min_inner_iter
+
         self.verbose = verbose if not self.is_linear else False
 
         temp_derivative = derivative or fenics.derivative(self.nonlinear_form, self.u)
@@ -293,6 +300,7 @@ def solve(self) -> fenics.Function:
                 atol=self.atol / 10.0,
                 comm=self.comm,
                 P=self.P_matrix,
+                min_iter=self.min_inner_iter,
             )
 
             if self.is_linear:
@@ -440,6 +448,7 @@ def newton_solve(
     b_tensor: Optional[fenics.PETScVector] = None,
     is_linear: bool = False,
     preconditioner_form: ufl.Form = None,
+    min_inner_iter: int = 0,
 ) -> fenics.Function:
     r"""Solves a nonlinear problem with Newton\'s method.
 
@@ -475,6 +484,11 @@ def newton_solve(
             sub-problem.
         is_linear: A boolean flag, which indicates whether the problem is actually
             linear.
+        preconditioner_form: A UFL form which defines the preconditioner matrix.
+        min_inner_iter: Minimum number of iterations for the inner (linear) solver
+            to perform. This parameter is only intended for being used with an
+            inexact Newton's method (there, in the first steps, only a small handful
+            of inner iterations are needed to satisfy the termination criterion).
 
     Returns:
         The solution of the nonlinear variational problem, if converged. This overwrites
@@ -523,6 +537,7 @@ def newton_solve(
         b_tensor=b_tensor,
         is_linear=is_linear,
         preconditioner_form=preconditioner_form,
+        min_inner_iter=min_inner_iter,
     )
 
     solution = solver.solve()
diff --git a/cashocs/nonlinear_solvers/picard_solver.py b/cashocs/nonlinear_solvers/picard_solver.py
index ae617f98..9c1aa839 100644
--- a/cashocs/nonlinear_solvers/picard_solver.py
+++ b/cashocs/nonlinear_solvers/picard_solver.py
@@ -94,6 +94,7 @@ def picard_iteration(
     b_tensors: Optional[List[fenics.PETScVector]] = None,
     inner_is_linear: bool = False,
     preconditioner_forms: Optional[Union[List[ufl.Form], ufl.Form]] = None,
+    min_inner_iter: int = 0,
 ) -> None:
     """Solves a system of coupled PDEs via a Picard iteration.
 
@@ -122,8 +123,10 @@ def picard_iteration(
         inner_is_linear: Boolean flag, if this is ``True``, all problems are actually
             linear ones, and only a linear solver is used.
         preconditioner_forms: The list of forms for the preconditioner. The default
-                is `None`, so that the preconditioner matrix is the same as the system
-                matrix.
+            is `None`, so that the preconditioner matrix is the same as the system
+            matrix.
+        min_inner_iter: Minimum iterations the inner (linear) solver should perform
+            regardless of termination criteria.
 
     """
     is_printing = verbose and fenics.MPI.rank(fenics.MPI.comm_world) == 0
@@ -188,6 +191,7 @@ def picard_iteration(
                 b_tensor=b_tensor,
                 is_linear=inner_is_linear,
                 preconditioner_form=preconditioner_form_list[j],
+                min_inner_iter=min_inner_iter,
             )
 
     if is_printing:
diff --git a/docs/source/user/demos/optimal_control/doc_config.rst b/docs/source/user/demos/optimal_control/doc_config.rst
index a6bd6d3a..96cb08dd 100755
--- a/docs/source/user/demos/optimal_control/doc_config.rst
+++ b/docs/source/user/demos/optimal_control/doc_config.rst
@@ -110,6 +110,14 @@ Additionally, we have the boolean parameter :ini:`newton_inexact`, defined via
 
 which sets up an inexact Newton method for solving nonlinear problems in case this is :ini:`newton_inexact = True`. The default is :ini:`newton_inexact = False`.
 
+For the inexact Newton's method, it can be beneficial to set a minimum number of iterations that have to be performed regardless of termination criteria. This is sensible, as the tolerance for the early iterations is usually very low, meaning that the linear solver converges in a handful of iterations. This leads to unwanted behavior, as the Jacobian is re-assembled very often initially. The parameter
+
+.. code-block:: ini
+
+    min_inner_iter = 0
+
+can be used to set the minimum amount of iterations that a linear solver should perform (when used in Newton's method). If this is :ini:`0` (which is the default), there is no minimum number of iterations prescribed.
+
 The parameter 
 
 .. code-block:: ini
@@ -633,6 +641,8 @@ in the following.
       - if :ini:`newton_damped = True`, damping is enabled
     * - :ini:`newton_inexact = False`
       - if :ini:`newton_inexact = True`, an inexact Newton's method is used
+    * - :ini:`min_inner_iter = 0`
+      - Specifies the minimum amount of iterations the inner (linear) solver has to perform in a Newton method
     * - :ini:`newton_verbose = False`
       - :ini:`newton_verbose = True` enables verbose output of Newton's method
     * - :ini:`picard_iteration = False`
diff --git a/docs/source/user/demos/shape_optimization/doc_config.rst b/docs/source/user/demos/shape_optimization/doc_config.rst
index 130d0474..8d301ad6 100755
--- a/docs/source/user/demos/shape_optimization/doc_config.rst
+++ b/docs/source/user/demos/shape_optimization/doc_config.rst
@@ -161,6 +161,14 @@ Additionally, we have the boolean parameter :ini:`newton_inexact`, defined via
 
 which sets up an inexact Newton method for solving nonlinear problems in case this is :ini:`newton_inexact = True`. The default is :ini:`newton_inexact = False`.
 
+For the inexact Newton's method, it can be beneficial to set a minimum number of iterations that have to be performed regardless of termination criteria. This is sensible, as the tolerance for the early iterations is usually very low, meaning that the linear solver converges in a handful of iterations. This leads to unwanted behavior, as the Jacobian is re-assembled very often initially. The parameter
+
+.. code-block:: ini
+
+    min_inner_iter = 0
+
+can be used to set the minimum amount of iterations that a linear solver should perform (when used in Newton's method). If this is :ini:`0` (which is the default), there is no minimum number of iterations prescribed.
+
 Next, we have the parameter
 
 .. code-block:: ini
@@ -1129,6 +1137,8 @@ in the following.
       - if :ini:`newton_damped = True`, damping is enabled
     * - :ini:`newton_inexact = False`
       - if :ini:`newton_inexact = True`, an inexact Newton's method is used
+    * - :ini:`min_inner_iter = 0`
+      - Specifies the minimum amount of iterations the inner (linear) solver has to perform in a Newton method
     * - :ini:`newton_verbose = False`
       - :ini:`newton_verbose = True` enables verbose output of Newton's method
     * - :ini:`picard_iteration = False`
diff --git a/pyproject.toml b/pyproject.toml
index 6ede151b..d2c548ce 100755
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -47,6 +47,7 @@ dependencies = [
 	"numpy",
 	"typing_extensions",
     "matplotlib",
+    "packaging",
 ]
 
 [project.urls]
@@ -136,5 +137,6 @@ module = [
 	"meshio",
 	"h5py",
     "matplotlib.*",
+    "packaging.*",
 ]
 ignore_missing_imports = true