Repository: pymor/pymor

Issue:
Default tolerances for newton algorithm
The current defaults for `pymor.algorithms.newton.newton` force the algorithm to never converge. Previously, we had `stagnation_threshold` set to `0.9` which seemed toopermissive for a default convergence criterion. @pymor/pymor-devs, @TiKeil, @lbalicki, do you have any opinions on this?

Also, note that all current stopping criteria are based on the Euclidean norm of the residual, whereas the norm of the update is also frequently used. Should we implement both (this is what FEniCS does)? I have also seen the dual norm of the residual being used in an NGSolve tutorial. I am not sure if it makes sense to implement everything one might possibly think of ..

Discussion:

Patch:
diff --git a/src/pymor/algorithms/newton.py b/src/pymor/algorithms/newton.py
index b751fb310..7efd3580b 100644
--- a/src/pymor/algorithms/newton.py
+++ b/src/pymor/algorithms/newton.py
@@ -14,11 +14,11 @@ from pymor.core.logger import getLogger
 
 
 @defaults('miniter', 'maxiter', 'rtol', 'atol', 'relax', 'stagnation_window', 'stagnation_threshold')
-def newton(operator, rhs, initial_guess=None, mu=None, error_product=None, least_squares=False,
-           miniter=0, maxiter=100, rtol=0., atol=0., relax=1., line_search_params=None,
-           stagnation_window=3, stagnation_threshold=np.inf,
+def newton(operator, rhs, initial_guess=None, mu=None, range_product=None, source_product=None, least_squares=False,
+           miniter=0, maxiter=100, atol=0., rtol=1e-7, relax='armijo', line_search_params=None,
+           stagnation_window=3, stagnation_threshold=np.inf, error_measure='update',
            return_stages=False, return_residuals=False):
-    """Basic Newton algorithm.
+    """Newton algorithm.
 
     This method solves the nonlinear equation ::
 
@@ -38,31 +38,39 @@ def newton(operator, rhs, initial_guess=None, mu=None, error_product=None, least
         solution `U`.
     mu
         The |parameter values| for which to solve the equation.
-    error_norm
-        The product with which the norm of the residual is computed. If `None`, the
-        Euclidean product is used.
+    range_product
+        The inner product `Operator` on `operator.range` with which the norm
+        of the resiudal is computed. If `None`, the Euclidean inner product
+        is used.
+    source_product
+        The inner product `Operator` on `operator.source` with which the norm
+        of the solution and update vectors is computed. If `None`, the Euclidean inner
+        product is used.
     least_squares
         If `True`, use a least squares linear solver (e.g. for residual minimization).
     miniter
         Minimum amount of iterations to perform.
     maxiter
         Fail if the iteration count reaches this value without converging.
-    rtol
-        Finish when the residual norm has been reduced by this factor relative to the
-        norm of the initial residual.
     atol
-        Finish when the residual norm is below this threshold.
+        Finish when the error measure is below this threshold.
+    rtol
+        Finish when the error measure has been reduced by this factor
+        relative to the norm of the initial residual resp. the norm of the current solution.
     relax
-        If real valued, relaxation factor for Newton updates; otherwise 'armijo' to
-        indicate that the :func:~pymor.algorithms.line_search.armijo line search algorithm
+        If real valued, relaxation factor for Newton updates; otherwise `'armijo'` to
+        indicate that the :func:`~pymor.algorithms.line_search.armijo` line search algorithm
         shall be used.
     line_search_params
         Dictionary of additional parameters passed to the line search method.
     stagnation_window
-        Finish when the residual norm has not been reduced by a factor of
+        Finish when the error measure has not been reduced by a factor of
         `stagnation_threshold` during the last `stagnation_window` iterations.
     stagnation_threshold
         See `stagnation_window`.
+    error_measure
+        If `'residual'`, convergence depends on the norm of the residual. If
+        `'update'`, convergence depends on the norm of the update vector.
     return_stages
         If `True`, return a |VectorArray| of the intermediate approximations of `U`
         after each iteration.
@@ -77,15 +85,19 @@ def newton(operator, rhs, initial_guess=None, mu=None, error_product=None, least
     data
         Dict containing the following fields:
 
-            :error_sequence:  |NumPy array| containing the residual norms after each iteration.
+            :solution_norms:  |NumPy array| of the solution norms after each iteration.
+            :update_norms:    |NumPy array| of the norms of the update vectors for each iteration.
+            :residual_norms:  |NumPy array| of the residual norms after each iteration.
             :stages:          See `return_stages`.
             :residuals:       See `return_residuals`.
 
     Raises
     ------
     NewtonError
-        Raised if the Netwon algorithm failed to converge.
+        Raised if the Newton algorithm failed to converge.
     """
+    assert error_measure in ('residual', 'update')
+
     logger = getLogger('pymor.algorithms.newton')
 
     data = {}
@@ -102,60 +114,107 @@ def newton(operator, rhs, initial_guess=None, mu=None, error_product=None, least
     U = initial_guess.copy()
     residual = rhs - operator.apply(U, mu=mu)
 
-    err = residual.norm(error_product)[0]
-    logger.info(f'      Initial Residual: {err:5e}')
+    # compute norms
+    solution_norm = U.norm(source_product)[0]
+    solution_norms = [solution_norm]
+    update_norms = []
+    residual_norm = residual.norm(range_product)[0]
+    residual_norms = [residual_norm]
+
+    # select error measure for convergence criteria
+    err = residual_norm if error_measure == 'residual' else np.inf
+    err_scale_factor = err
+    errs = residual_norms if error_measure == 'residual' else update_norms
+
+    logger.info(f'     norm:{solution_norm:.3e}                                 res:{residual_norm:.3e}')
 
     iteration = 0
-    error_sequence = [err]
     while True:
+        # check for convergence / failure of convergence
         if iteration >= miniter:
-            if err <= atol:
-                logger.info(f'Absolute limit of {atol} reached. Stopping.')
+            if residual_norm == 0:
+                # handle the corner case where error_norm == update, U is the exact solution
+                # and the jacobian of operator is not invertible at the exact solution
+                logger.info(f'Norm of residual exactly zero. Converged.')
+                break
+            if err < atol:
+                logger.info(f'Absolute tolerance of {atol} for norm of {error_measure} reached. Converged.')
                 break
-            if err/error_sequence[0] <= rtol:
-                logger.info(f'Prescribed total reduction of {rtol} reached. Stopping.')
+            if err < rtol * err_scale_factor:
+                logger.info(f'Relative tolerance of {rtol} for norm of {error_measure} reached. Converged.')
                 break
-            if (len(error_sequence) >= stagnation_window + 1
-                    and err/max(error_sequence[-stagnation_window - 1:]) >= stagnation_threshold):
-                logger.info(f'Error is stagnating (threshold: {stagnation_threshold:5e}, window: {stagnation_window}). '
-                            f'Stopping.')
+            if (len(errs) >= stagnation_window + 1
+                    and err > stagnation_threshold * max(errs[-stagnation_window - 1:])):
+                logger.info(f'Norm of {error_measure} is stagnating (threshold: {stagnation_threshold:5e}, '
+                            f'window: {stagnation_window}). Converged.')
                 break
             if iteration >= maxiter:
-                raise NewtonError('Failed to converge')
+                raise NewtonError('Failed to converge after {iteration} iterations.')
+
+        iteration += 1
+
+        # store convergence history
         if iteration > 0 and return_stages:
             data['stages'].append(U)
         if return_residuals:
             data['residuals'].append(residual)
-        iteration += 1
+
+        # compute update
         jacobian = operator.jacobian(U, mu=mu)
         try:
-            correction = jacobian.apply_inverse(residual, least_squares=least_squares)
+            update = jacobian.apply_inverse(residual, least_squares=least_squares)
         except InversionError:
             raise NewtonError('Could not invert jacobian')
+
+        # compute step size
         if isinstance(relax, Number):
             step_size = relax
         elif relax == 'armijo':
-            logger.info(f'Using Armijo as line search method')
             def res(x):
                 residual_vec = rhs - operator.apply(x, mu=mu)
-                return residual_vec.norm(error_product)[0]
-            if error_product is None:
+                return residual_vec.norm(range_product)[0]
+
+            if range_product is None:
                 grad = - (jacobian.apply(residual) + jacobian.apply_adjoint(residual))
             else:
-                grad = - (jacobian.apply_adjoint(error_product.apply(residual)) + jacobian.apply(error_product.apply_adjoint(residual)))
-            step_size = armijo(res, U, correction, grad=grad, initial_value=err, **(line_search_params or {}))
+                grad = - (jacobian.apply_adjoint(range_product.apply(residual))
+                          + jacobian.apply(range_product.apply_adjoint(residual)))
+            step_size = armijo(res, U, update, grad=grad, initial_value=residual_norm, **(line_search_params or {}))
         else:
             raise ValueError('Unknown line search method')
-        U.axpy(step_size, correction)
+
+        # update solution and residual
+        U.axpy(step_size, update)
         residual = rhs - operator.apply(U, mu=mu)
 
-        err = residual.norm(error_product)[0]
-        logger.info(f'Iteration {iteration:2}: Residual: {err:5e},  '
-                    f'Reduction: {err / error_sequence[-1]:5e}, Total Reduction: {err / error_sequence[0]:5e}')
-        error_sequence.append(err)
-        if not np.isfinite(err):
+        # compute norms
+        solution_norm = U.norm(source_product)[0]
+        solution_norms.append(solution_norm)
+        update_norm = update.norm(source_product)[0] * step_size
+        update_norms.append(update_norm)
+        residual_norm = residual.norm(range_product)[0]
+        residual_norms.append(residual_norm)
+
+        # select error measure for next iteration
+        err = residual_norm if error_measure == 'residual' else update_norm
+        if error_measure == 'update':
+            err_scale_factor = solution_norm
+
+        logger.info(f'it:{iteration} '
+                    f'norm:{solution_norm:.3e} '
+                    f'upd:{update_norm:.3e} '
+                    f'rel_upd:{update_norm / solution_norm:.3e} '
+                    f'res:{residual_norm:.3e} '
+                    f'red:{residual_norm / residual_norms[-2]:.3e} '
+                    f'tot_red:{residual_norm / residual_norms[0]:.3e}')
+
+        if not np.isfinite(residual_norm) or not np.isfinite(solution_norm):
             raise NewtonError('Failed to converge')
 
-    data['error_sequence'] = np.array(error_sequence)
+    logger.info('')
+
+    data['solution_norms'] = np.array(solution_norms)
+    data['update_norms']   = np.array(update_norms)
+    data['residual_norms'] = np.array(residual_norms)
 
     return U, data
diff --git a/src/pymor/operators/interface.py b/src/pymor/operators/interface.py
index ce383e43a..0dceb5024 100644
--- a/src/pymor/operators/interface.py
+++ b/src/pymor/operators/interface.py
@@ -233,13 +233,16 @@ class Operator(ParametricObject):
                 options = {}
             options['least_squares'] = least_squares
 
-            R = V.empty(reserve=len(V))
-            for i in range(len(V)):
-                try:
-                    R.append(newton(self, V[i], initial_guess=initial_guess[i] if initial_guess is not None else None,
-                                    mu=mu, **options)[0])
-                except NewtonError as e:
-                    raise InversionError(e)
+            with self.logger.block('Solving nonlinear problem using newton algorithm ...'):
+                R = V.empty(reserve=len(V))
+                for i in range(len(V)):
+                    try:
+                        R.append(newton(self, V[i],
+                                        initial_guess=initial_guess[i] if initial_guess is not None else None,
+                                        mu=mu,
+                                        **options)[0])
+                    except NewtonError as e:
+                        raise InversionError(e)
             return R
 
     def apply_inverse_adjoint(self, U, mu=None, initial_guess=None, least_squares=False):


Test patch:
diff --git a/src/pymortests/algorithms/stuff.py b/src/pymortests/algorithms/stuff.py
index be73f2b3b..4c0f72308 100644
--- a/src/pymortests/algorithms/stuff.py
+++ b/src/pymortests/algorithms/stuff.py
@@ -19,9 +19,13 @@ def _newton(mop, initial_value=1.0, **kwargs):
     return newton(mop, rhs, initial_guess=guess, **kwargs)
 
 @pytest.mark.parametrize("order", list(range(1, 8)))
-def test_newton(order):
+@pytest.mark.parametrize("error_measure", ['update', 'residual'])
+def test_newton(order, error_measure):
     mop = MonomOperator(order)
-    U, _ = _newton(mop, atol=1e-15)
+    U, _ = _newton(mop,
+                   atol=1e-15 if error_measure == 'residual' else 1e-7,
+                   rtol=0.,
+                   error_measure=error_measure)
     assert float_cmp(mop.apply(U).to_numpy(), 0.0)
 
 def test_newton_fail():
@@ -37,7 +41,7 @@ def test_newton_with_line_search():
 def test_newton_fail_without_line_search():
     mop = MonomOperator(3) - 2 * MonomOperator(1) + 2 * MonomOperator(0)
     with pytest.raises(NewtonError):
-        _ = _newton(mop, initial_value=0.0, atol=1e-15)
+        _ = _newton(mop, initial_value=0.0, atol=1e-15, relax=1.)
 
 def test_newton_unknown_line_search():
     mop = MonomOperator(1)