---
abstract: 'We analyze backward step control globalization for finding zeros of Gâteaux-differentiable functions that map from a Banach space to a Hilbert space. The results include global convergence to a distinctive solution characterized by propagating the initial guess by a generalized Newton flow with guaranteed bounds on the discrete nonlinear residual norm decrease and an (also numerically) easily controllable asymptotic linear residual convergence rate. The convergence theory can be exploited to construct efficient numerical methods, which we demonstrate for the case of a Krylov–Newton method and an approximation-by-discretization multilevel framework. Both approaches optimize the asymptotic linear residual convergence rate, either over the Krylov subspace or through adaptive discretization, which in turn yields practical and efficient stopping criteria and refinement strategies that balance the nonlinear residuals with the relative residuals of the linear systems. We apply these methods to the class of nonlinear elliptic boundary value problems and present numerical results for the Carrier equation and the minimum surface equation.'
author:
- Andreas Potschka
title: '[Backward step control for Hilbert space problems]{}'
---

Introduction {#sec:introduction}
============

Let $U$ be a Banach space with norm ${\left\lVert.\right\rVert}_U$ and $V$ be a Hilbert space (we discuss generalizations to Banach spaces in section \[sec:Banach\]) with inner product ${\left(.,.\right)}_V$ and norm $\smash[t]{{\left\lVertv\right\rVert}_V = \smash[b]{\sqrt{{\left(v,v\right)}_V}}}$ for $v \in V$. For some open subset $D \subseteq U$, let $F: D \to V$ be continuously Gâteaux differentiable [@Hamilton1982 3.1.1. Def.] with derivative $F': D
\times U \to V$, i.e., $$F'(u, \delta u) = \lim_{h \to 0} \frac{1}{h} \left( F(u+h \delta u) - F(u)
  \right)$$ exists for all $u \in D$, $\delta u \in U$ and $F'$ is continuous as a function from the product space $D \times U$ to $V$ (which is a weaker requirement than continuity of $F'$ considered as a map from $D$ to $\mathcal{L}(U, V)$, the Banach space of all bounded linear operators from $U$ to $V$). Note that continuous Gâteaux differentiability implies linearity of $F'$ in the second argument [@Hamilton1982 3.2.5 Thm.]. In the following, we write $A(u) \delta u$ shorthand for $A(u, \delta u)$ whenever an operator $A$ is linear in the second argument. We consider the problem of finding an unknown $u \in D$ such that $$\label{eqn:FOfxIsZero}
  F(u) = 0_V$$ with a Newton-type iteration: Given $u_0 \in D$, find a suitable approximation $M: D \times V \to U$ (linear in the second argument) of the inverse of $F'(u)$ and a step size sequence $(t_k)_{k \in \mathbb{N}}$ satisfying $t_k \in [0, 1]$ such that the iteration $$\label{eqn:Newton}
  u_{k+1} = u_k + t_k \delta u_k \quad 
  \text{with } \delta u_k = -M(u_k) F(u_k)$$ converges to a solution $u^\ast \in D$ of . The first and main part of this article is devoted to finding a suitable step size sequence $(t_k)_{k \in \mathbb{N}}$ in section \[sec:convergence\]. Out of the many ways to construct $M(u)$, we elaborate on two choices in section \[sec:designM\]. For convenience, we define the negative generalized Newton flow $f: D \to U$ as $$f(u) = M(u)F(u).$$ As in [@Potschka2016], our convergence analysis will be based on generalized Newton paths $u^k: [0, \infty) \to U$, which are defined as the solutions of the initial value problems $$\label{eqn:genDavidenko}
  \frac{{{\mathrm d}}u^k}{{{\mathrm d}}t}(t) = -f(u^k(t)) \quad 
  \text{for } t \in [0, \infty) \quad 
  \text{with } u^k(0) = u_k.$$ We shall prove existence and uniqueness of solutions to in our setting in Theorem \[thm:arclength\]. Note that is an explicit Euler discretization of with step sizes $(t_k)_{k \in \mathbb{N}}$. If $M$ is chosen as the inverse of $F'$ then is a damped Newton method and is the Davidenko differential equation [@Davidenko1953]. The equivalence of the Newton method with explicit Euler on the Davidenko differential equation has been exploited by various authors (see, e.g., [@Deuflhard1974; @Brezinski1975; @Ascher1987; @Bock2000b; @Deuflhard2006]).

In the theory and the numerics below, the operator $M$ does not appear explicitly anymore and only the function $f$ will be required, which implicitly defines $M(u)$ in the direction $F(u)$. As it turns out, all other directions of $M(u)$ are not important. This observation leads to crucial improvements in the assumptions stated in [@Potschka2016].

The convergence theory below lends itself immediately to the construction of numerical algorithms for the approximate solution of via . In particular, it can be used to construct $M(u)$ from $F'(u)$ by finite-dimensional approximation, which can then be exploited to construct adaptive discretization schemes that optimize the contraction rate of the algorithm in $V$. In the case of Finite Element analysis, our approach delivers a multi-level Newton adaptive mesh refinement algorithm [@Hohmann1994] that can be used straight-forward as another tool complementing refinement strategies based on a posteriori error estimation (see, e.g., [@Graetsch2005; @Ainsworth2000; @Becker2001]).

Historically, Newton-type methods come in two flavors: Either, we do not solve the linearized systems with operator $F'(u_k)$ and right-hand side $-F(u_k)$ exactly, which is usually known under the name *inexact Newton* method [@Dembo1982; @Ypma1984]. Alternatively, we apply an approximation $M(u)$ of the inverse of $F'(u)$ directly, an approach which is sometimes called *approximate Newton* method and which is the classical form of Quasi-Newton methods [@Dennis1977]. As pointed out by Bank and Rose [@Bank1981], the two flavors are in fact different formulations of the same class of methods. As mentioned above, our analysis is based on the approximate Newton formulation in order to define the generalized Newton path , but the control of the linearized residual in the sense of inexact Newton methods emerges as the $\kappa$-condition A\[ass:kappa\], which is identical to a choice of $\eta_k = \kappa < 1$ for the residual forcing sequence $(\eta_k)$ first proposed in [@Dembo1982]. In turn, the classical analysis of the local rates of convergence in [@Dembo1982] carries over to our setting if $\kappa$ is allowed to be reduced from iteration to iteration in the sense of a forcing sequence once we are close to a solution. This shall, however, not be the focus of this paper, where we focus on the preasymptotic global convergence and are content with locally linear convergence rates.

Similar to [@Deuflhard1991], we base our globalization approach on a continuous curve, but we substitute here the exact Newton path (obtained with the choice $M(u) = \smash[t]{F'(u)}^{-1}$ from ) by a generalized Newton path, which is allowed to have $M(u) \neq
\smash[t]{F'(u)}^{-1}.$ Following ideas of [@Hohmann1994], we can use the (contravariant) $\kappa$-condition A\[ass:kappa\] to design a multilevel Newton-type method for the construction of $M$ based on adaptive discretization. This leads to adaptive discretizations solely based on balancing the discretization residual with the nonlinear residual of the Newton-type method. We remark that this convenient black-box approach might be inferior to more involved schemes that balance the discretization errors and nonlinear errors (instead of residuals) [@Rannacher2013] or exploit underlying (energy) minimization properties [@Deuflhard1998] for particular problem classes. Nonetheless, our method is the first for which convergence to the closest solution in the sense of the generalized Newton flow can be proven.

#### Contributions

In this article, we extend the convergence analysis of backward step control for from the finite-dimensional to the Hilbert space setting. We provide reasonable assumptions and convergence results for with backward step control. The main result is convergence to a distinct solution characterized by the propagation of the initial guess by the generalized Newton flow provided that no singularity of the problem interferes. In addition, we prove an a priori bound on the nonlinear reduction of the residual norm. The convergence theory can be exploited to construct efficient numerical algorithms, which we discuss for the case of a Krylov–Newton method and a Finite Element approximation. Both are based on the optimization of the residual contraction constant, which yields in the latter case an efficient adaptive mesh refinement strategy. The results are demonstrated for the numerical solution of the Carrier equation and the minimal surface equation.

#### Overview

In section \[sec:convergence\], we discuss the general assumptions, derive and motivate the method of backward step control, provide reasons why implicit and higher-order time stepping methods for are not advisable, provide step size bounds, and establish the notion of generalized Newton paths, which are the central ingredient for the analysis of local and global convergence of backward step control. We then discuss extensions to Banach spaces and present an algorithmic realization with a minimal working Matlab example code for the solution of $\arctan(u) = 0$. In section \[sec:designM\], we exploit the convergence analysis to construct two Newton-type methods, a Krylov–Newton method and a method based on approximation-by-discretization. We apply these methods in section \[sec:ellipticPDE\] to the class of nonlinear elliptic boundary value problems and provide numerical results for the Carrier equation and the minimum surface equation.

#### Notation

We denote the open ball of radius $r>0$ around $u\in U$ by $B(u,r)$ and the Laplace operator by ${\Delta}= \nabla \cdot \nabla$. As usual, we write $C^0$ for the space of continuous functions, $H^1_0(\Omega)$ for the Sobolev space of square integrable functions on a bounded domain $\Omega \subset \mathbb{R}^n$ that vanish at the boundary and admit square integrable derivatives, and $H^{-1}(\Omega)$ for its dual space. The Euler number is denoted by $e =
\sum_{k=0}^{\infty} \frac{1}{k!}$.

Convergence analysis {#sec:convergence}
====================

The overall structure of the backward step control convergence analysis in Hilbert spaces is similar to the finite-dimensional case [@Potschka2016]. The intricate interplay of the changes in the details, however, advises us to present the convergence analysis in a self-contained fashion.

Discussion of assumptions
-------------------------

We start with the following definitions.

The *level function* $T: D \to \mathbb{R}$ is $T(u) = \frac{1}{2} {\left\lVertF(u)\right\rVert}_V^2.$

The *level set of $u \in D$* is $\smash[t]{\widetilde{\mathcal{T}}}(u) = \left\{ \bar{u} \in D \mid
  T(\bar{u}) \le T(u) \right\}.$

The *path connected level set of $u \in D$* is $$\mathcal{T}(u) = \left\{ \bar{u} \in \widetilde{\mathcal{T}}(u)
  \mid \exists 
  c \in C^0\left([0, 1], \widetilde{\mathcal{T}}(u)\right)
      \text{with } c(0) = u, c(1) = \bar{u} \right\}.$$

For $r \in (1, \infty)$ the *set of $r$-regular points* is $$\mathcal{R}_r = \left\{ u \in D \mid r^{-1} {\left\lVertF(u)\right\rVert}_V <
  {\left\lVertf(u)\right\rVert}_U < r {\left\lVertF(u)\right\rVert}_V \right\}.$$

The *set of $\infty$-regular points* is $\mathcal{R}_\infty = \bigcup_{r \in (1, \infty)} \mathcal{R}_r$.

We remark that if $u \in D \setminus \mathcal{R}_\infty$, which means that $u
\not\in \mathcal{R}_r$ for all $r
\in (1, \infty)$, then $M(u)$ is either not bounded or does not admit a bounded inverse [@Yosida1978 §I.6, Cor. 2, 3]. The contrary is, however, not true: $M(u)$ may be unbounded or not admit a bounded inverse although $u \in
\mathcal{R}_r$ for some $r \in (1, \infty)$, because in the definition of $\mathcal{R}_r$ only the action of $M(u)$ in direction $F(u)$ is of interest.

We require the following assumptions to hold true:

\[ass:validIni\] There exists an $r \in (1, \infty)$ such that $u_0 \in \mathcal{R}_r$, and ${\left\lVertF(u_0)\right\rVert}_V > 0$.

\[ass:kappa\] There exists a $\kappa < 1$ such that for all $u \in \mathcal{R}_r \cap
  \mathcal{T}(u_0)$ $${\left\lVertF(u) - F'(u) f(u)\right\rVert}_V \le \kappa {\left\lVertF(u)\right\rVert}_V.$$

\[ass:omega\] There exists an $\omega < \infty$ such that for all $u \in \mathcal{T}(u_0), t
  \in [0, 1]$ $${\left\lVert\left[ F'(u) - F'(u-tf(u)) \right] f(u)\right\rVert}_V \le \omega t {\left\lVertf(u)\right\rVert}_U
  {\left\lVertF(u)\right\rVert}_V.$$

\[ass:fLipschitz\] There exists an $L < \infty$ such that for all $u, \bar{u} \in
  \mathcal{T}(u_0)$ $${\left\lVertf(u) - f(\bar{u})\right\rVert}_U \le L {\left\lVertu - \bar{u}\right\rVert}_U.$$

\[ass:gamma\] For all $\eta > 0$ there exist constants $\gamma$, $t_\gamma > 0$ such that for all $t \in [0, t_\gamma]$, $u \in \mathcal{R}_r \cap \mathcal{T}(u_0)$ with ${\left\lVertf(u)\right\rVert}_U > \eta$ $${\left\lVertf(u-tf(u)) - f(u)\right\rVert}_U \ge \gamma t.$$

The main difference in the assumptions here compared to the finite-dimensional setting in [@Potschka2016] is the weakening of A\[ass:kappa\] and A\[ass:omega\] from a formulation with matrices to a formulation which requires the properties to hold only in the direction of the residual $F(u)$. Thus, all requirements can be postulated without using norms for operators that map between $U$ and $V$. Apart from the avoidance of operator norms, we had to replace all arguments based on compactness of bounded sets by other means for the proofs in the Hilbert space case. The discussion of the assumptions in [@Potschka2016] still applies to a large extent here: We require in A\[ass:validIni\] that $u_0$ is an $r$-regular point but not a solution. The central $\kappa$-condition A\[ass:kappa\] is a contravariant version of Bock’s covariant $\kappa$-condition [@Bock1987] in the sense that it quantifies on the one hand the deviation of the inexact increment $\delta u = -M(u) F(u)$ from the Newton increment $\delta
u^\mathrm{Newton} = -F(u)^{-1} F(u)$ in the $V$-norm $${\left\lVertF'(u) \left[ \delta u^\mathrm{Newton} - \delta u \right]\right\rVert}_V
  = {\left\lVertF(u) - F'(u)f(u)\right\rVert}_V
  \le \kappa {\left\lVertF(u)\right\rVert}_V$$ in comparison to Bock’s covariant $\kappa^\mathrm{cov} < 1$ in the $U$-norm (see also [@Potschka2013 section 5.2]) $$\label{eqn:kappa_Bock}
  {\left\lVertM(u - f(u)) \left[ F(u) - F'(u)f(u) \right]\right\rVert}_U 
  \le \kappa^\mathrm{cov} {\left\lVertf(u)\right\rVert}_U.$$ On the other hand, $\kappa$ in A\[ass:kappa\] characterizes the asymptotic Q-linear convergence rate of the residual norms ${\left\lVertF(u_k)\right\rVert}_V$, whereas $\kappa^\mathrm{cov}$ in characterizes the asymptotic R-linear convergence rate of the error ${\left\lVertu_k - u_\ast\right\rVert}_U$ if $u_\ast =
\lim_{k \to \infty} u_k$ (for a discussion of different affine invariances see [@Deuflhard2006]). The $\omega$-condition A\[ass:omega\] measures a combination of the nonlinearity and the well-posedness of because if $F'$ is Lipschitz continuous with Lipschitz constant $L'$, then we obtain $${\left\lVert\left[ F'(u) - F'(u - tf(u)) \right] f(u)\right\rVert}_V \le L' t {\left\lVertf(u)\right\rVert}_U^2$$ and boundedness of $M(u)$ in direction $F(u)$ with constant $C$ implies A\[ass:omega\] with $\omega = CL'$. The Lipschitz condition A\[ass:fLipschitz\] is classical. The nonstandard assumption A\[ass:gamma\] follows, for instance, if $f$ is bi-Lipschitz with constant $\ell$ $${\left\lVertf(u - tf(u)) - f(u)\right\rVert}_U \ge \ell t {\left\lVertf(u)\right\rVert}_U$$ with $\gamma = \eta \ell$ and $t_\gamma$ arbitrary.

Backward step control
---------------------

Newton-type methods are explicit Euler discretizations with step sizes $t_k$ of the generalized Newton flow . Thus, the convergence of Newton-type methods is strongly connected to the stability problem of the explicit Euler method. Implict Euler, in contrast, has ideal stability properties: It is an L-stable method (see, e.g., [@Hairer1996]). Hence, in order to determine $t_k$, we consider the backward iterate $$\bar{u}_k(t_k) := u_{k+1} + t_k f(u_{k+1}) = u_k + t_k g(u_k, t_k) \quad
\text{with } g(u, t) := f(u - tf(u)) - f(u).$$ The point $\bar{u}_k(t_k)$ is the starting point of a (stable) implicit Euler step for that arrives exactly at $u_{k+1}$, the result of a possibly unstable explicit Euler step starting from $u_k$. The idea of backward step control is based on a backward error argument: If a small perturbation of the starting point $u_k$ can be found from which a stable implicit Euler step arrives exactly at $u_{k+1}$, we can accept the step size. We thus require that the distance between $u_k$ and $\bar{u}_k(t_k)$ is bounded by some fixed constant $H>0$ through the choice $$\label{eqn:BSC} \tag{BSC}
  t_k = \min \mathcal{B}_H(u_k) \quad \text{where }
  \mathcal{B}_H(u) = \left\{ t \in [0, 1] \mid H = t {\left\lVertg(u,t)\right\rVert}_U
  \right\} \cup \{1\},$$ which implies ${\left\lVert\bar{u}_k(t_k) - u_k\right\rVert}_U \le H$ (with equality for $t_k <
1$) by continuity of $g$.

Implicit and higher-order time stepping methods
-----------------------------------------------

The question whether explicit Euler is really the best method to solve arises naturally. We can answer this question affirmatively for two reasons: First, all implicit methods have the drawback that an approximated inverse of an operator involving derivatives of the approximated inverse $M(u)$ would be required, e.g., in the case of the implicit Euler method $$\begin{aligned}
  0 &= u_{k+1} + t_k f(u_{k+1}) - u_k,
  \intertext{with a local Newton corrector}
  u_{k+1}^{i+1} &= u_{k+1}^i - \left[ {{\mathrm I}}_U + t_k f'(u_{k+1}^i) \right]^{-1}
  \left( u_{k+1}^i + t_k f(u_{k+1}^i) - u_k \right),\end{aligned}$$ which is not readily available and would require higher regularity of $M$ than guaranteed by the assumptions above. Second, higher order methods would destroy the well-known locally quadratic convergence of the Newton method, where $M(u) = (F'(u))^{-1}$. This can be seen from the homotopy formulation $$\label{eqn:uHomotopy}
  F(u(t)) - e^{-t} F(u_0) = 0,$$ which we can differentiate with respect to $t$ to arrive exactly at provided that $F'(u(t))$ stays invertible. Thus, the second order truncation error of explicit Euler is required to obtain locally quadratic convergence, because higher consistency orders would result in the locally linear convergence dictated by .

We have not explored multi-step methods of order one further. The possible outcome of this line of research is unfortunately unclear at present and exceeds the scope of this paper.

Step size bounds
----------------

\[lem:localFullSteps\] If A\[ass:validIni\] and A\[ass:fLipschitz\] hold, then delivers full steps $t_k = 1$ in the vicinity of a solution $u^\ast \in
  \mathcal{R}_r \cap \mathcal{T}(u_0)$.

Let $u_k \in B(u^\ast, L^{-2} H)$. Hence, it holds for all $t \in [0, 1]$ that $$t {\left\lVertg(u_k, t)\right\rVert}_U 
  \overset{\text{A\ref{ass:fLipschitz}}}{\le} L t^2 {\left\lVertf(u_k)\right\rVert}_U
  = L t^2 {\left\lVertf(u_k) - f(u^\ast)\right\rVert}_U 
  \overset{\text{A\ref{ass:fLipschitz}}}{\le} L^2 t^2 {\left\lVertu_k - u^\ast\right\rVert}_U < H.$$ Thus, $\mathcal{B}_H(u_k) = \{1\}$ and $t_k = 1$ by virtue of .

\[lem:lowerStepsizeBound\] If A\[ass:validIni\] and A\[ass:fLipschitz\] hold, then generates for all $u_k \in \mathcal{R}_r \cap \mathcal{T}(u_0)$ step sizes that are either $t_k = 1$ or have the lower bounds $$t_k \ge \frac{\sqrt{H}}{\sqrt{L {\left\lVertf(u_k)\right\rVert}_U}}
  > \frac{\sqrt{H}}{\sqrt{rL{\left\lVertF(u_k)\right\rVert}_V}}
  \ge \frac{\sqrt{H}}{\sqrt{rL{\left\lVertF(u_0)\right\rVert}_V}}.$$

*Proof * If $t_k < 1$, then $$t_k^2 \overset{\eqref{eqn:BSC}}{=} \frac{t_k H}{{\left\lVertg(u_k, t_k)\right\rVert}_U}
  \overset{\text{A\ref{ass:fLipschitz}}}{\ge}
  \frac{H}{L {\left\lVertf(u_k)\right\rVert}_U}
  > \frac{H}{rL{\left\lVertF(u_k)\right\rVert}_V}
  \ge \frac{H}{rL{\left\lVertF(u_0)\right\rVert}_V}. \qedhere$$

\[lem:upperStepsizeBound\] Let A\[ass:validIni\] and A\[ass:gamma\] hold and let $\bar{t}
  \in (0, 1)$ and $\eta > 0$. Then there exists an ${\overline{H}} > 0$ such that for all $H \in (0, {\overline{H}}]$ and $u \in \mathcal{R}_r \cap
  \mathcal{T}(u_0)$ with ${\left\lVertf(u)\right\rVert}_U \ge \eta$ it holds that $\min
  \mathcal{B}_H(u) \le \bar{t}$.

*Proof by contradiction * We assume to the contrary that for all ${\overline{H}} > 0$ there exists an $H
  \in (0, {\overline{H}}]$ and a $u \in \mathcal{R}_r \cap \mathcal{T}(u_0)$ satisfying ${\left\lVertf(u)\right\rVert}_U \ge \eta$ and $\min \mathcal{B}_H(u) > \bar{t}$. Then, A\[ass:gamma\] guarantees the existence of $\gamma$, $t_\gamma > 0$ such that for $t := \min \{ t_\gamma, \bar{t} \} < \min \mathcal{B}_H(u)$ we obtain from that $${\overline{H}} \ge H \ge t {\left\lVertg(u, t)\right\rVert}_U \ge \gamma t^2 > 0.$$ Because $\eta$ and thus $\gamma$ and $t$ are independent of ${\overline{H}}$, we obtain a contradiction for ${\overline{H}} \to 0$.

\[lem:dampedStepBound\] Let A\[ass:validIni\], A\[ass:fLipschitz\], and A\[ass:gamma\] hold and let $\eta > 0$. Then there exists an ${\overline{H}} > 0$ such that for all $H
  \in (0, {\overline{H}}]$ and $u \in \mathcal{R}_r \cap \mathcal{T}(u_0)$ it holds that $${\left\lVertf(u)\right\rVert}_U \min \mathcal{B}_H(u) \le \eta.$$

We choose $\bar{t} \in (0, 1)$ sufficiently small so that it satisfies $r
  {\left\lVertF(u_0)\right\rVert}_V \bar{t} \le \eta.$ Then, Lemma \[lem:upperStepsizeBound\] yields the existence of an ${\overline{H}} > 0$ such that for all $H \in (0, {\overline{H}}]$ and all $u \in
  \mathcal{R}_r \cap \mathcal{T}(u_0)$ with ${\left\lVertf(u)\right\rVert}_U \ge \eta$ it holds that $\min \mathcal{B}_H(u) \le \bar{t}$. Hence, $${\left\lVertf(u)\right\rVert}_U \min \mathcal{B}_H(u)
  \overset{\text{A\ref{ass:validIni}}}{<} 
  r {\left\lVertF(u)\right\rVert}_V \bar{t} \le r {\left\lVertF(u_0)\right\rVert} \bar{t} \le \eta.$$ For the remaining $u \in \mathcal{R}_r \cap \mathcal{T}(u_0)$ the assertion holds by virtue of ${\left\lVertf(u)\right\rVert}_U < \eta$.

Finite arclength of generalized Newton paths
--------------------------------------------

In the next step, we study the generalized Newton paths given by .

\[lem:Fdescent\] If A\[ass:validIni\], A\[ass:kappa\], A\[ass:fLipschitz\] and $u_k \in
  \mathcal{R}_r \cap \mathcal{T}(u_0)$ hold, then there exists $\bar{t} > 0$ such that has a unique local solution $u^k(t) \in
  \mathcal{R}_r \cap \mathcal{T}(u_k)$ for $t \in [0, \bar{t})$ which satisfies $${\left\lVertF(u^k(t))\right\rVert}_V \le e^{-(1-\kappa)t} {\left\lVertF(u_k)\right\rVert}_V \quad 
  \text{for all } t \in [0, \bar{t}).$$

The Picard–Lindelöf theorem [@Amann1990 II.7, exercise 3] yields with A\[ass:fLipschitz\] the existence of a unique local solution $u^k(t)$ to in some neighborhood $(-\bar{t}, \bar{t})$ of $t =
  0$. Without loss of generality, $\bar{t} > 0$ is small enough to ensure $u^k(t) \in \mathcal{R}_r$ for $t \in [0, \bar{t})$ because $\mathcal{R}_r$ is open. For ease of notation, we abbreviate $u^k(t)$ by $u$. The Cauchy–Schwarz inequality and A\[ass:kappa\] show that the level function is nonincreasing along this solution because $$\begin{aligned}
    \frac{{{\mathrm d}}}{{{\mathrm d}}t} T(u) &= {\left(F(u), F'(u) \frac{{{\mathrm d}}u}{{{\mathrm d}}t}\right)}_V
    = -{\left(F(u), F'(u) f(u)\right)}_V\\
    &= -{\left\lVertF(u)\right\rVert}_V^2 + {\left(F(u), F(u) - F'(u) f(u)\right)}_V\\
    &\le -{\left\lVertF(u)\right\rVert}_V^2 + \kappa {\left\lVertF(u)\right\rVert}_V^2
    = -2(1-\kappa) T(u) \le 0.
  \end{aligned}$$ Gronwall’s inequality (see, e.g., [@Amann1990]) yields $$T(u^k(t)) \le e^{-2 (1-\kappa) t} T(u_k)$$ and thus $u^k(t) \in \mathcal{T}(u_k)$ for $t \in [0, \bar{t})$. The assertion follows after multiplication by two and taking square roots.

We show in Theorem \[thm:arclength\] below that the quantities in the following definition are well-defined under suitable assumptions.

\[def:regularpart\] For $r \in (1,\infty)$, we define the *$r$-regular part $u_r^k$* of the generalized Newton path $u^k$ as the solution to the initial value problem $$\frac{{{\mathrm d}}u_r^k}{{{\mathrm d}}t}(t) = 
  \begin{cases}
    -f(u_r^k(t)) & \text{for } u_r^k(t) \in \mathcal{R}_r,\\
    0 & \text{otherwise},
  \end{cases}
  \quad \text{for } t \in [0, \infty),
  \quad \text{with } u_r^k(0) = u_k.$$ We denote its limit by $u_k^\ast = \lim_{t \to \infty} u^k_r(t)$ and define $t_k^\ast = \inf \{t \in [0, \infty) \mid u_r^k(t) \not\in \mathcal{R}_r \}$ with the usual convention that $\inf \varnothing = \infty$.

\[thm:arclength\] Let A\[ass:validIni\], A\[ass:kappa\], and A\[ass:fLipschitz\] hold. If $u_k \in \mathcal{R}_r \cap \mathcal{T}(u_0)$, then the $r$-regular part of the generalized Newton path exists uniquely and has a finite arclength $\ell(u_r^k)$ satisfying $${\left\lVertu_k - u_k^\ast\right\rVert}_U \le
  \ell(u_r^k) < \frac{r}{1 - \kappa} {\left\lVertF(u_k)\right\rVert}_V 
  < \frac{r^2}{1 - \kappa} {\left\lVertf(u_k)\right\rVert}_U.$$ If $u_k^\ast \in \mathcal{R}_r$, then $F(u_k^\ast) = 0$.

*Proof * The unique local solution of Lemma \[lem:Fdescent\] can be extended uniquely in $\mathcal{T}(u_k)$ by repeated application of the Picard–Lindelöf theorem either until $u^k(t) \not \in \mathcal{R}_r$ for some $t = t_k^\ast$ or to the whole interval $t \in [0, \infty)$. In the first case, the $r$-regular part is uniquely determined by $u_r^k(t) = u_r^k(t_k^\ast)$ for all $t \ge t_k^\ast$. We can now use the definition of $\mathcal{R}_r$ and Lemma \[lem:Fdescent\] in order to show $$\begin{aligned}
    \ell(u_r^k) &= \int_{0}^{\infty} {\left\lVert\frac{{{\mathrm d}}u_r^k}{{{\mathrm d}}t}(t)\right\rVert}_U {{\mathrm d}}t
    = \int_{0}^{t_k^\ast} {\left\lVertf(u^k(t))\right\rVert}_U {{\mathrm d}}t
    < r \int_{0}^{t_k^\ast} {\left\lVertF(u^k(t))\right\rVert}_V {{\mathrm d}}t \\
    &\le r \int_{0}^{\infty} e^{-(1-\kappa)t} {{\mathrm d}}t {\left\lVertF(u_k)\right\rVert}_V
    = \frac{r}{1-\kappa} {\left\lVertF(u_k)\right\rVert}_V < \frac{r^2}{1-\kappa}
    {\left\lVertf(u_k)\right\rVert}_U. 
  \end{aligned}$$ We obtain the lower arclength bound by noting that the shortest path between $u_k$ and $u_k^\ast$ has arclength ${\left\lVertu_k - u_k^\ast\right\rVert}_U$. If $u_k^\ast \in
  \mathcal{R}_r$, then $t_k^\ast = \infty$ and Lemma \[lem:Fdescent\] reveals $${\left\lVertF(u_k^\ast)\right\rVert}_V \le \lim_{t \to \infty} e^{-(1-\kappa)t}{\left\lVertF(u_k)\right\rVert}_V =
  0.
  \qedhere$$

Local convergence
-----------------

We use the next lemma to prove discrete descent of the residual norm.

\[lem:kappaScaled\] A\[ass:kappa\] holds if and only if for all $u \in \mathcal{R}_r \cap
  \mathcal{T}(u_0)$ and $t \in [0, 1]$ $${\left\lVertF(u) - t F'(u) f(u)\right\rVert}_V 
  \le \left[ 1 - (1-\kappa) t \right] {\left\lVertF(u)\right\rVert}_V.$$

*Proof * As in [@Potschka2016], the nontrivial direction of the proof follows from the convexity of the functional $\varphi(t) = {\left\lVertF(u) - t F'(u) f(u)\right\rVert}_V$ and A\[ass:kappa\] according to $$\varphi(t) \le (1-t)\varphi(0) + t \varphi(1) 
  \le \left[ (1-t) + \kappa t \right] {\left\lVertF(u)\right\rVert}_V.
  \qedhere$$

\[lem:Fdecrease\] Let A\[ass:validIni\], A\[ass:kappa\] and A\[ass:omega\] hold. If $u_k
  \in \mathcal{R}_r \cap \mathcal{T}(u_0)$, then $${\left\lVertF(u_{k+1})\right\rVert}_V \le \left[ 1 - (1-\kappa) t_k + \frac{\omega}{2}
  {\left\lVertf(u_k)\right\rVert}_U t_k^2 \right] {\left\lVertF(u_k)\right\rVert}_V.$$ Furthermore, if there exists a $\theta < 1$ such that the step size sequence satisfies $$\omega t_k {\left\lVertf(u_k)\right\rVert}_U \le 2 \theta (1 - \kappa),$$ then $${\left\lVertF(u_{k+1})\right\rVert}_V \le \left[ 1 - (1-\theta) (1-\kappa) t_k \right]
  {\left\lVertF(u_k)\right\rVert}_V.$$

Because $F'$ is continuously Gâteaux differentiable, we can apply [@Hamilton1982 3.2.2. Thm. and 2.1.4. Thm.] to show that for $u \in D$ and $\delta u \in U$ $$F(u + t_k \delta u) - F(u) =
    \int_{0}^{1} F'(u + \tau t_k \delta u) t_k \delta u {{\mathrm d}}\tau
    \overset{\tau = \frac{t}{t_k}}{=} \int_{0}^{t_k} F'(u + t \delta u) \delta
    u {{\mathrm d}}t.$$ Using Lemma \[lem:kappaScaled\] and A\[ass:omega\] we obtain the first assertion from $$\begin{aligned}
    & \quad\, {\left\lVertF(u_{k+1})\right\rVert}
    = {\left\lVertF(u_k) - \int_{0}^{t_k} F'(u_k - \tau f(u_k)) f(u_k) {{\mathrm d}}\tau\right\rVert}_V\\
    &= {\left\lVertF(u_k) - t_k F'(u_k) f(u_k) + \int_{0}^{t_k} \left[ F'(u_k) -
    F'(u_k - \tau f(u_k)) \right] f(u_k) {{\mathrm d}}\tau\right\rVert}_V\\
    &\le {\left\lVertF(u_k) - t_k F'(u_k) f(u_k)\right\rVert}_V
    + \int_{0}^{t_k} {\left\lVert\left[ F'(u_k) - F'(u_k - \tau f(u_k)) \right]
    f(u_k)\right\rVert}_V {{\mathrm d}}\tau\\
    &\le \left[ 1 - (1 - \kappa) t_k + \frac{\omega}{2} {\left\lVertf(u_k)\right\rVert}_U t_k^2
    \right] {\left\lVertF(u_k)\right\rVert}_V.
  \end{aligned}$$ The second assertion follows immediately.

We can now state a local convergence theorem.

\[thm:localConvergence\] Let A\[ass:validIni\], A\[ass:kappa\], and A\[ass:omega\] hold. If there exists a $\bar{t} \in (0,1)$ such that $t_k \ge \bar{t}$ for all $k \in \mathbb{N}$ and if there exists a $\theta \in (0, 1)$ such that for some $k \in
  \mathbb{N}$ the iterate $u_k$ satisfies $$\mathcal{T}(u_k) \subseteq \mathcal{R}_r
  \quad \text{and} \quad
  \omega r {\left\lVertF(u_k)\right\rVert}_V \le 2 \theta (1 - \kappa),$$ then $(u_k)_{k \in \mathbb{N}}$ converges to some point $u^\ast \in
  \mathcal{T}(u_k)$ with $F(u^\ast) = 0$.

Because $t_k \in [0, 1]$, we have $$\omega t_k {\left\lVertf(u_k)\right\rVert}_U \le \omega {\left\lVertf(u_k)\right\rVert}_U
  \overset{\text{A\ref{ass:validIni}}}{\le} \omega r {\left\lVertF(u_k)\right\rVert}_V
  \le 2 \theta (1 - \kappa).$$ Hence, repeated application of Lemma \[lem:Fdecrease\] yields for all $j \in
  \mathbb{N}$ $$\label{eqn:geomFConv}
    {\left\lVertF(u_{k+j})\right\rVert}_V 
    \le q^j {\left\lVertF(u_k)\right\rVert}_V \quad \text{with }
    q := 1 - (1-\theta)(1-\kappa) \bar{t}.$$ Because $q < 1$, ${\left\lVertF(u_k)\right\rVert}_V$ converges geometrically. In addition, we obtain that $(u_k)_{k \in
  \mathbb{N}}$ is a Cauchy sequence by virtue of $$\begin{aligned}
    {\left\lVertu_k - u_{k+j}\right\rVert}_U &\le \sum_{i=0}^{j-1} {\left\lVertu_{k+i} - u_{k+i-1}\right\rVert}_U
    = \sum_{i=0}^{j-1} t_{k+i} {\left\lVertf(u_{k+i})\right\rVert}_U\\
    &\le r \sum_{i=0}^{j-1} {\left\lVertF(u_{k+i})\right\rVert}_V 
    \le r {\left\lVertF(u_{k})\right\rVert}_V \sum_{i=0}^{\infty} q^i 
    = \frac{r}{1-q} {\left\lVertF(u_k)\right\rVert} \overset{k \to \infty}{\longrightarrow} 0.
  \end{aligned}$$ Thus, $(u_{k})_{k \in \mathbb{N}}$ converges to some $u^\ast \in \mathcal{T}(u_k)
  \subseteq \mathcal{R}_r$ and implies $F(u^\ast) = 0$.

For the rate of convergence, we obtain the following result:

\[lem:linearconvergence\] Under the assumptions of Theorem \[thm:localConvergence\], ${\left\lVertF(u_k)\right\rVert}_V$ converges linearly with asymptotic linear convergence rate $\kappa < 1$.

Because $u_k \in \mathcal{R}_r$, it follows that ${\left\lVertf(u_k)\right\rVert}_U \le r
  {\left\lVertF(u_k)\right\rVert}_V \to 0$. Hence, there is a sequence $(\theta_K)_{K \in
  \mathbb{N}}$ with $\theta_K \in [0, 1)$ and $\theta_K \to 0$ such that $$\omega t_k {\left\lVertf(u_k)\right\rVert}_U \le 2 \theta_K (1 - \kappa) \quad \text{for all } k
  \ge K.$$ Lemma \[lem:localFullSteps\] and repeated application of Lemma \[lem:Fdecrease\] then deliver $${\left\lVertF(u_{K+1})\right\rVert}_V \le \left[ 1 - (1-\theta_K)(1-\kappa) \right]
  {\left\lVertF(u_K)\right\rVert}_V,$$ where $1 - (1-\theta_K)(1-\kappa) \to \kappa$ as $K \to \infty$.

In order to obtain methods with guaranteed superlinear or quadratic local convergence, it is necessary to appropriately drive $\kappa$ to zero as $F(u_k)
\to 0$ as exhaustively described by the means of forcing sequences $\eta_k =
\kappa$ in [@Dembo1982].

Global convergence
------------------

The following lemmas are required for the main theorem, which assures convergence of the iterates to $u_0^\ast$. As a prerequisite, we prove that every neighborhood of an isolated zero $u^\ast$ of $F$ contains a path connected level set that contains a neighborhood of $u^\ast$.

\[lem:boundedLevelSets\] Let A\[ass:validIni\], A\[ass:kappa\], and A\[ass:fLipschitz\] hold. If there exist $\varepsilon > 0$ and $u^\ast \in D$ such that $u^\ast$ is the only zero of $F$ on $B(u^\ast, \varepsilon) \subseteq \mathcal{R}_r \cap
  \mathcal{T}(u_0),$ then there exists an $\tilde{\varepsilon} > 0$ with $$\bigcup_{u \in B(u^\ast, \tilde{\varepsilon})} \mathcal{T}(u) 
  \subseteq B(u^\ast, \varepsilon).$$

*Proof by contradiction * We assume to the contrary that there exists a sequence $(u_n)_{n \in
  \mathbb{N}}$ with ${\left\lVertu_n - u^*\right\rVert}_U < \frac{\varepsilon}{2n}$ and $\mathcal{T}(u_n) \not\subseteq B(u^*, \varepsilon)$. Hence, there exists a sequence $(\tilde{v}_n)_{n \in \mathbb{N}}$ with $\tilde{v}_n \in \mathcal{T}(u_n)$ and ${\left\lVert\tilde{v}_n - u^\ast\right\rVert}_U \ge
  \varepsilon$. Because $\mathcal{T}(u_n)$ is path connected, there exist continuous functions $c_n: [0, 1] \to \mathcal{T}(u_n)$ with $c_n(0) = u_n$ and $c_n(1) = \tilde{v}_n$. Because ${\left\lVertc_n(0) - u^*\right\rVert}_U <
  \frac{\varepsilon}{2}$ and ${\left\lVertc_n(1) - u^*\right\rVert}_U \ge \varepsilon$, the intermediate value theorem yields the existence of $v_n = c_n(\tau_n) \in
  \mathcal{T}(u_n)$ for some $\tau_n \in [0, 1]$ satisfying $$\label{eqn:annulus}
    {\left\lVertv_n - u^*\right\rVert}_U = \frac{\varepsilon}{2}.$$ By Theorem \[thm:arclength\], we obtain for the distance to the limit $v_n^\ast$ of the $r$-regular part of the generalized Newton path emanating from $v_n$ that $$\begin{aligned}
    {\left\lVertv_n - v_n^\ast\right\rVert}_U 
    &\le \ell(v_r^n) < \frac{r}{1-\kappa} {\left\lVertF(v_n)\right\rVert}_V
    \le \frac{r}{1-\kappa} {\left\lVertF(u_n)\right\rVert}_V
    < \frac{r^2}{1-\kappa} {\left\lVertf(u_n)\right\rVert}_U\\
    &= \frac{r^2}{1-\kappa} {\left\lVertf(u_n) - f(u^*)\right\rVert}_U
    \le \frac{r^2L}{1-\kappa} {\left\lVertu_n - u^*\right\rVert}_U <
    \frac{r^2L\varepsilon}{2(1-\kappa)n} \to 0,
  \end{aligned}$$ which implies for some sufficiently large $n$ that $v^*_n \in B(u^*,
  \varepsilon) \subseteq \mathcal{R}_r$ and thus $F(v^*_n) = 0$. By we get $v_n^* \neq u^*$ in contradiction to the uniqueness of $u^*$. We also need a bound on the deviation of two neighboring generalized Newton paths emanating from $u_k$ and $u_{k+1}$.

\[lem:oneStepError\] Let A\[ass:fLipschitz\] hold. If $u^k(t_k + \tau)$, $u^{k+1}(\tau) \in
  \mathcal{T}(u_k)$ for all $\tau \in [0, t]$, then $${\left\lVertu^k(t_k + t) - u^{k+1}(t)\right\rVert}_U 
  \le \frac{1}{2} {\left\lVertf(u_k)\right\rVert}_U L e^{L(t_k+t)} t_k^2.$$

We use the integral form of the Gronwall inequality as in [@Potschka2016 Lem. 8.5].

In order to prove an a priori bound on the decrease of the nonlinear residual for , we need the following Lemma.

\[lem:sqrtBound\] Let $h > 0$. If a sequence $(a_k)_{k \in \mathbb{N}}$ of nonnegative numbers satisfies $a_{k+1}^2 \le a_k^2 - 2 h a_k$ for all $k \in \mathbb{N},$ then $a_k \le \max \{a_0 - k h, 0\}$ for all $k \in \mathbb{N}.$

If $a_k \le 2 h$, it follows that $a_{k+1} = 0$ by virtue of $0 \le a_{k+1}^2 \le a_k (a_k - 2h) \le 0.$ If $a_k > 2 h$, we immediately obtain $a_{k+1}^2 \le a_k^2 - 2 h a_k + h^2 = (a_k - h)^2.$ Taking the square root completes the proof.

We can now prove the main theorem of this article.

(0,0) circle (1); (-0.9,0) circle (0.7); (0,-2) arc (-90:90:2); (0,2) arc (90:270:2); (0,2) – (-1,2); (-1,2) arc (90:150:5); (0,-2) – (-1,-2); (-1,-2) arc (90:150:1); (0,0) – (-1,0); (-1,0) arc (90:150:3); (-1,-0.1) – (-1,0.1); (-0.8,-0.1) – (-0.8,0.1); (-0.8,-0.7) arc (-90:90:0.7); (-0.8,-0.7) – (-1,-0.7); (-0.8,0.7) – (-1,0.7); (-1,0.7) arc(90:150:3.7); (-1,-0.7) arc(90:150:2.3); (0,0) circle (2pt) node\[right\] [$u_0^\ast$]{}; (-0.9,-1.2) – (-0.9,0) node [(0,0) circle (1pt);]{}; at (0,-1.4) [$u^0\left( \sum_{i=0}^{k-1} t_i \right)$]{}; at (0,1) [$B(u_0^\ast, \frac{\varepsilon}{2})$]{}; at (-3.8,1.5) [$\mathcal{N}_\varepsilon$]{}; at (-3,0.5) [$\mathcal{N}_{\tilde{\varepsilon}}$]{}; at (-2.5,-0.1) [$u^0(t)$]{};

\[thm:convergenceBSC\] Let A\[ass:validIni\], A\[ass:kappa\], A\[ass:omega\], A\[ass:fLipschitz\], and A\[ass:gamma\] hold. If there exists an $\varepsilon > 0$ such that $$\mathcal{N}_\varepsilon := \bigcup_{t \in [0, \infty)} B(u^0(t), \varepsilon)
  \subseteq \mathcal{R}_r$$ and $u_0^\ast$ is an isolated zero of $F$, then there exist constants ${\overline{H}} >0$, $K < \infty$, and $c > 0$ such that the following statements hold true:

1.  For all $H \in (0, {\overline{H}}]$ iteration with step size selection converges to $u_0^\ast$.

2.  After $\smash[t]{\bar{k}} = \lceil K / \smash[t]{\sqrt{H}} \rceil$ steps, $u_{\bar{k}}$ lies in the region of local full step convergence.

3.  The sequence of residual norms ${\left\lVertF(u_k)\right\rVert}_V$ decreases geometrically and satisfies the additional a priori bound $$\sqrt{{\left\lVertF(u_k)\right\rVert}_V}
            \le \sqrt{{\left\lVertF(u_0)\right\rVert}_V} - k c \sqrt{H}$$ for all $k \in \{ i \in \mathbb{N} \mid t_j < 1, j = 0, \dotsc, i \}$, which includes $k = \smash[t]{\bar{k}}.$

Because $u_0^\ast$ is an isolated zero of $F$, we can assume without loss of generality that $\varepsilon > 0$ was chosen small enough such that $u_0^\ast$ is the only zero of $F$ in $\mathcal{N}_\varepsilon$. Let now $\theta \in (0, 1)$. We further assume without loss of generality $\varepsilon$ to be sufficiently small to satisfy $\omega r^2 L \varepsilon \le 2 \theta (1-\kappa)$. It follows with A\[ass:validIni\] and A\[ass:fLipschitz\] that for all $u
  \in B(u_0^\ast, \varepsilon)$ $$\begin{aligned}
    \omega r {\left\lVertF(u)\right\rVert}_V 
    &< \omega r^2 {\left\lVertf(u)\right\rVert}_U
    = \omega r^2 {\left\lVertf(u) - f(u_0^\ast)\right\rVert}_U\\
    &\le \omega r^2 L {\left\lVertu - u_0^\ast\right\rVert}_U
    < \omega r^2 L \varepsilon
    \le 2 \theta (1 - \kappa).
  \end{aligned}$$ Thus, Theorem \[thm:localConvergence\] establishes that $B(u_0^\ast,
  \varepsilon)$ is contained in the region of local full step convergence. By possibly reducing $\varepsilon$ further according to Lemma \[lem:boundedLevelSets\], we can guarantee that if $u_{k} \in
  B(u_0^\ast, \varepsilon)$ for some $k \in \mathbb{N}$, then $(u_k)_k$ converges to the unique zero $u_0^\ast$ in $\mathcal{N}_\varepsilon.$ For the first statement, it now remains to show that $u_k \in B(u_0^\ast,
  \varepsilon)$ for some $k \in \mathbb{N}.$

To this end, we choose $T_\ast < \infty$ such that $${\left\lVertu^0(t) - u_0^\ast\right\rVert}_U \le \frac{\varepsilon}{2} \quad \text{for all } t
  \ge T_\ast - 1.$$ Because $u_0^\ast$ is the only zero of $F$ on $\mathcal{N}_\varepsilon$, we can choose an $\tilde{\varepsilon} \in (0, \frac{\varepsilon}{2})$ such that there is an $\eta > 0$ satisfying $${\left\lVertf(u)\right\rVert}_U > \eta \quad \text{for all } 
  u \in \mathcal{T}(u_0) \cap \mathcal{N}_{\tilde{\varepsilon}}, \text{ where }
  \mathcal{N}_{\tilde{\varepsilon}} := \bigcup_{t \in [0, T_\ast]} B(u^0(t),
  \tilde{\varepsilon}).$$ From A\[ass:gamma\] we obtain the existence of constants $\gamma > 0$ and $t_\gamma \in (0, 1)$ such that $$\label{eqn:lbg}
    {\left\lVertg(u, t)\right\rVert}_U \ge \gamma t \quad \text{for all } t \in [0, t_\gamma]
    \text{ and } u \in \mathcal{N}_{\tilde{\varepsilon}}.$$ We can then use $\bar{t} = t_\gamma$ in Lemma \[lem:upperStepsizeBound\] to obtain a constant ${\overline{H}} > 0$ such that $$\label{eqn:ubtk}
    \min \mathcal{B}_H(u) \le t_\gamma < 1 \quad \text{for all } H \in (0,
    {\overline{H}}] \text{ and } u \in \mathcal{N}_{\tilde{\varepsilon}}.$$ In anticipation of a later argument in the proof, we can assume without loss of generality that ${\overline{H}}$ is sufficiently small to satisfy $$\label{eqn:smallHbar}
    T_\ast e^{L T_\ast}\left( rL {\left\lVertF(u_0)\right\rVert}_V \right)^{\frac{3}{2}}
    \smash[t]{{\overline{H}}}^{\frac{1}{2}} \le 2 \gamma \tilde{\varepsilon}.$$ Combining the inequalities and with , we see that $$\label{eqn:lbtk2}
    \gamma t_k^2 \le t_k {\left\lVertg(u_k, t_k)\right\rVert}_U \le H \quad \text{for all } H \in
    (0, {\overline{H}}] \text{ and } u_k \in \mathcal{N}_{\tilde{\varepsilon}}.$$ We now choose an $H$-dependent $\bar{k} \in \mathbb{N}$ that satisfies $$T_\ast - 1 \le \sum_{i=0}^{\bar{k}-1} t_i \le T_\ast$$ and show by induction that $u_k \in \mathcal{T}(u_0) \cap
  \mathcal{N}_{\tilde{\varepsilon}}$ for all $k \le \bar{k}$, which clearly holds true for $k=0$ because $u_0 = u^0(0) \in B(u_0, \tilde{\varepsilon}).$ We can now assume inductively that $u_i \in \mathcal{T}(u_0) \cap
  \mathcal{N}_{\tilde{\varepsilon}}$ for all $i \le k-1$ with $k \le \bar{k}$ in order to show $u_k \in  \mathcal{T}(u_0) \cap
  \mathcal{N}_{\tilde{\varepsilon}}.$ Because $\mathcal{N}_{\tilde{\varepsilon}} \subseteq \mathcal{R}_r$, it follows from that $t_i < 1$. Hence, Lemma \[lem:lowerStepsizeBound\] yields $$T_\ast \ge
  \sum_{i=0}^{k-1} t_i \ge k \frac{\sqrt{H}}{\sqrt{r L {\left\lVertF(u_0)\right\rVert}_V}},$$ which implies the bound $$\label{eqn:ubkbar}
    k \le \frac{T_\ast \sqrt{rL{\left\lVertF(u_0)\right\rVert}_V}}{\sqrt{H}}.$$ We then obtain by a telescope argument, Lemma \[lem:oneStepError\], A\[ass:validIni\], , , and that $$\begin{aligned}
    {\left\lVertu^0 \left( \sum_{i=0}^{k-1} t_i \right) - u_k\right\rVert}_U
    &\le \sum_{j=0}^{k-1} {\left\lVertu^j \left( \sum_{i=j}^{k-1} t_i \right) -
    u^{j+1} \left( \sum_{i=j+1}^{k-1} t_i \right)\right\rVert}_U\\
    &\le
    \sum_{j=0}^{k-1} \frac{1}{2} {\left\lVertf(u_{j})\right\rVert}_U L e^{L T_\ast} t_j^2
    < \frac{rL e^{LT_\ast} {\left\lVertF(u_0)\right\rVert}_V}{2 \gamma} H k\\
    &\le \frac{T_\ast e^{LT_\ast} \left( rL {\left\lVertF(u_0)\right\rVert}_V
    \right)^{\frac{3}{2}} }{2 \gamma} \sqrt{H} \le \tilde{\varepsilon}.
  \end{aligned}$$ Because $\sum_{i=0}^{k-1} t_i \le T_\ast$, we have established by induction that $u_k \in \mathcal{T}(u_0) \cap \mathcal{N}_{\tilde{\varepsilon}}$ for all $k \le \bar{k}$. Finally, $u_{\bar{k}} \in B(u_0^\ast, \varepsilon)$ by virtue of $${\left\lVertu_0^\ast - u_{\bar{k}}\right\rVert}_U
  \le {\left\lVert u_0^\ast - u^0 \left( \sum_{i=0}^{\bar{k}-1} t_i \right)\right\rVert}_U
  + {\left\lVert u^0 \left( \sum_{i=0}^{\bar{k}-1} t_i \right) - u_{\bar{k}}\right\rVert}_U
  < \frac{\varepsilon}{2} + \tilde{\varepsilon} < \varepsilon.$$ This proves the first statement. The second assertion follows from by choosing $K = T_\ast
  \sqrt{r L {\left\lVertF(u_0)\right\rVert}_V}$ for $k = \bar{k}$. For the third statement, we obtain from the lower step size bound of Lemma \[lem:lowerStepsizeBound\] and Lemma \[lem:Fdecrease\] that for all $k$ with $t_k < 1$ it holds that $$\begin{gathered}
    {\left\lVertF(u_{k+1})\right\rVert}_V \le \left[ 1 - (1-\theta)(1-\kappa)t_k \right]
    {\left\lVertF(u_k)\right\rVert}_V
    \le {\left\lVertF(u_k)\right\rVert}_V - 2 c \sqrt{H {\left\lVertF(u_k)\right\rVert}_V},\\
    \text{with } c = \frac{(1-\theta)(1-\kappa)}{2 \sqrt{rL}} > 0.
  \end{gathered}$$ Hence, the residual norm sequence converges geometrically and Lemma \[lem:sqrtBound\] yields the a priori bound with $a_k =
  \sqrt{{\left\lVertF(u_k)\right\rVert}_V}$ and $h = c \sqrt{H}$. Because we established $u_k \in
  \mathcal{N}_{\tilde{\varepsilon}}$ for all $k \le \smash[t]{\bar{k}}$ in the proof of the first statement, we obtain $t_k < 1$ for all $k \le
  \smash[t]{\bar{k}}$ from , which completes the proof.

A note on generalizations to Banach spaces {#sec:Banach}
------------------------------------------

The only step in the convergence proof that exploits the Hilbert space structure of $V$ is Lemma \[lem:Fdescent\]. All remaining steps can be carried out even if $V$ is only a Banach space. The general approach here is to modify the used level function and to require an additional assumption akin to condition A\[ass:kappa\].

In this section, we shortly present the necessary modifications for the special case of the Lebesgue space $V =
L^p(\Omega)$ of real-valued $p$-integrable functions with $ 2 < p <
\infty$ and $\Omega \subset \mathbb{R}^n$. First, we need to consider a different level function $$T_p(u) = \frac{1}{p} {\left\lVertF(u)\right\rVert}_V^p 
= \frac{1}{p} \int_{\Omega} {\left\lvertF(u)(s)\right\rvert}^p {{\mathrm d}}s,$$ with analogously defined level sets $\mathcal{T}_p(u)$. In addition to A\[ass:kappa\], we require pointwise that $$\left[F(u) - F'(u)f(u)\right](s) \le \kappa {\left\lvertF(u)(s)\right\rvert} \quad \text{for almost all } s
\in \Omega \text{ and all } u \in \mathcal{R}_r \cap \mathcal{T}_p(u_0).$$ Standard arguments on the differentiability of $p$-norms (see, e.g., [@Lieb2001 Thm. 2.6]) then deliver with the abbreviation $u = u^k(t)$ that $$\begin{aligned}
  & \quad\, \frac{{{\mathrm d}}}{{{\mathrm d}}t} T_p(u^k(t))
  = -\int_{\Omega} {\left\lvertF(u)(s)\right\rvert}^{p-2} \left[ F(u)(s) \cdot ( F'(u) f(u)) (s)
  \right] {{\mathrm d}}s\\
  &= -\int_{\Omega} {\left\lvertF(u)(s)\right\rvert}^p {{\mathrm d}}s 
  + \int_{\Omega} {\left\lvertF(u)(s)\right\rvert}^{p-2} \left[ F(u)(s) \cdot (F(u) - F'(u) f(u))
  (s) \right] {{\mathrm d}}s\\
  &\le -p T_p(u)
  + \int_{\Omega} {\left\lvertF(u)(s)\right\rvert}^{p-1} (F(u) - F'(u) f(u)) (s) {{\mathrm d}}s\\
  &\le -p(1-\kappa) T_p(u) \le 0.\end{aligned}$$ After application of Gronwall’s inequality, multiplication by $p$ and taking the $p$-th root, we establish the result of Lemma \[lem:Fdescent\] $${\left\lVertF(u^k(t))\right\rVert}_V \le e^{-(1-\kappa)t} {\left\lVertF(u_k)\right\rVert}_V.$$

Algorithmic realization {#sec:realization}
-----------------------

The algorithmic realization of can be carried over verbatim from the finite-dimensional setting laid out in [@Potschka2016 section 10] with the use of ${\left\lVert.\right\rVert}_U$ for all occuring norms. As in [@Potschka2016], we do not use monotone iterations [@Galantai2014] for numerical computations here but use the bisection procedure with exponentially smoothed step size prediction. For convenience, we sketch it again: In order to determine $t_k$ from , we approximately compute a zero of the Lipschitz continuous scalar function $t \mapsto t {\left\lVertg(u_k, t)\right\rVert}_U - H$ by a bracketing procedure. Numerically, we are content with a $t_k$ that satisfies $$\label{eqn:HlHu}
  H'_k := t_k {\left\lVertg(u_k, t_k)\right\rVert}_U \in [H^\mathrm{l}, H^\mathrm{u}] 
  \quad \text{or} \quad
  t_k \approx 1 \text{ and } H'_k < H^\mathrm{l},$$ where $H^\mathrm{l} < H$ and $H^\mathrm{u} > H$ are close to $H$.

For illustration purposes, we provide in Table \[alg:example\_atan\] example code in Matlab, which computes for the real-valued function $F(u) =
\arctan(u)$ the iterates $u_1, \dotsc, u_5$ of iteration with $M(u) = F'(u)^{-1} = u^2 + 1$ and step sizes $t_k$ satisfying , starting from $u_0 = 2$ for $H = 0.8$. Full step Newton diverges for this choice of $u_0$. For the sake of brevity, the remaining algorithmic parameters suggested in [@Potschka2016 section 10.2] are used as explict values in the code and termination and error checks are omitted.

      k       t         u        du       dup    Hprime
      0  1.0000   2.0e+00  -5.5e+00   1.7e+01   2.3e+01  decrease t
      0  0.5000   2.0e+00  -5.5e+00   1.0e+00   3.3e+00  decrease t
      0  0.2500   2.0e+00  -5.5e+00  -7.6e-01   1.2e+00  accept t
      1  0.2335   6.2e-01  -7.6e-01  -4.9e-01   6.3e-02  increase t
      1  0.6168   6.2e-01  -7.6e-01  -1.5e-01   3.8e-01  accept t
      2  0.7543   1.5e-01  -1.5e-01  -3.4e-02   8.6e-02  accept t
      3  1.0000   3.4e-02  -3.4e-02   2.7e-05   3.4e-02  accept t
      4  1.0000  -2.7e-05   2.7e-05  -1.3e-14   2.7e-05  accept t
      5  1.0000   1.3e-14  -1.3e-14  -0.0e+00   1.3e-14  accept t
      

We display the output in Table \[tab:output\]. Each line in the output corresponds to one evaluation of $f(u) = M(u) F(u)$, not counting the initial evaluation of $f(u_0)$. We observe that for $k = 0$, two bisection steps are required to reduce $t_0$ to $\frac{1}{4}$. For $k = 1$, one bisection step is required to increase $t_1$ to above $0.6$. For $k \ge 2$, the predicted step sizes are accepted without further bisection steps. The iterate $u_5$, which is known already at the end of iteration 4, is an acceptable solution candidate, because ${\left\lVert\delta u_5\right\rVert}_2 \approx 1.3 \cdot 10^{-14}$.

The smoothed step size prediction in line 5 is vital for reducing the number of bracketing steps. The predicted step size often already satisfies in all but a few iterations and thus almost no extra computational effort in terms of increment evaluations $f(u_k)$, which typically comprise setting up and solving one linear system, is required for the globalization procedure in most iterations.

Discussion of drawbacks
-----------------------

Backward step control appears to suffer from two substantial drawbacks: First, it depends on a problem specific parameter $H > 0$ that needs to be chosen sufficiently small for convergence, but not too small to cause unecessarily many iterations and thus inefficiency of the method. Second, it may not converge in one step on affine linear problems. We argue here that these two seemingly detrimental properties are actually *necessary* for the class of methods that converge to the closest solution (in the sense of the Newton flow): Let us assume we want to solve a one-dimensional nonlinear equation $F(u) = 0$, $F:
\mathbb{R} \to \mathbb{R}$, with a certain *parameter-free* nonlinear method. Assume it produces a sequence of iterates $(u^k)_{k \in \mathbb{N}}
\subset \mathbb{R}$ that converges to some $u^\ast \neq u^0$. Because $u^0$ cannot be an accumulation point of $(u^k)_{k \in \mathbb{N}}$, we can find the closest iterate $u^j$ to $u^0$ for some $j \in \mathbb{N}$. Any function that differs from $F$ only on the open interval $I$ between $u^0$ and $u^j$ will inevitably lead to the same iterates $(u^k)_{k \in \mathbb{N}}$ if we apply the parameter-free nonlinear method to it. Thus, we can modify $F$ smoothly on $I$ to introduce zeros of $F$ that are closer to $u^0$ than $u^\ast$. We illustrate this construction for an affine linear $F$ in Fig. \[fig:drawback\]. We conclude that any method that provably converges to the closest zero must depend on a parameter such as $H$ to account for problem specific quantities that are virtually impossible to estimate numerically.

Moreover, existing globalization methods [@Deuflhard1991; @Hohmann1994; @Deuflhard1998; @Deuflhard2006] are not parameter free either, because they require an initial step size guess $t_0$ that needs to be sufficiently small.

Design of Newton-type methods {#sec:designM}
=============================

The convergence analysis of backward step control lends itself immediately to the design of globally convergent Newton-type methods. The two required steps are:

1.  Define $M$ (or directly $f$) respecting the $\kappa$-condition A\[ass:kappa\].

2.  Use to determine the step size sequence $(t_k)$.

The second step is generic. For the first step, we give two important examples in the following two sections. Both methods find approximations $\delta u_k$ of the Newton increment $\delta
u_k^\mathrm{Newton}$ determined by the (infinite-dimensional) linear system $$\label{eqn:HilbertSpaceNewton}
  F'(u_k) \delta u_k^\mathrm{Newton} = -F(u_k).$$ We emphasize that the operator $M(u_k)$ is implicitly determined by requiring $\delta u_k = -M(u_k) F(u_k)$. It is not required to actually compute $M(u_k)$, as long as we have $\delta u_k$. However, we need to make sure that $f$ is Lipschitz continuous with respect to $u$.

Krylov–Newton methods {#sec:KrylovNewton}
---------------------

Krylov subspace methods like CG, MINRES, or GMRES [@Hestenes1952; @Paige1975; @Saad1986] for the iterative solution of linear systems, originally developed for large but finite-dimensional sparse systems, can also be stated for infinite-dimensional linear operators. The convergence theory is more complicated in the infinite-dimensional case (see, e.g., [@Nevanlinna1993; @Gasparo2008]). The structure of the linear mapping $F'(u_k): U \to V$ does usually not immediately admit the application of Krylov subspace methods, unless a left preconditioner $P_k: V \to U$ is available, such that $A_k := P_k F'(u_k)$ is an endomorphism on $U$. Then, the $m$-th iterate $\delta u_k^{m}$ of a Krylov subspace method applied to is an approximate solution restricted to the (at most $m$-dimensional) $m$-th Krylov subspace $$\mathcal{K}^m(A_k, P_k F(u_k)) = \left\{ q(A_k) P_k F(u_k) \mid
  q \text{ is a polynomial of degree less than } m \right\}.$$ We now focus on residual minimizing Krylov subspace like GMRES and MINRES, which are constructed on the basis of the additional optimality requirement $$\label{eqn:residual_minimizing}
  \delta u_k^m = \operatorname*{arg~min}_{\delta u \in \mathcal{K}^m(A_k, P_k F(u_k))}
  {\left\lVertP_k F(u_k) + A_k \delta u_k)\right\rVert}_U.$$ As in [@Guennel2014], we investigate the important special case where $V = U^\ast$. In this setting, we can choose $P_k$ as the Riesz isomorphism and immediately obtain $$\label{eqn:Krylov}
  \begin{aligned}
    {\left\lVertP_k F(u_k) + A_k \delta u_k)\right\rVert}_U &=
    {\left\lVertP_k \left[ F(u_k) + F'(u_k) \delta u_k \right]\right\rVert}_U\\
    &= {\left\lVertF(u_k) + F'(u_k) \delta u_k)\right\rVert}_V.
  \end{aligned}$$ In the light of and , the $\kappa$-condition A\[ass:kappa\] is nothing but the classical relative termination condition used by Krylov subspace methods with given relative tolerance $\kappa$. In other words, GMRES and MINRES choose from the Krylov subspace the increment that achieves the smallest left-hand side in the condition of A\[ass:kappa\] and thus yield (with Lemma \[lem:linearconvergence\]) an asymptotic linear convergence rate bounded above by $\kappa$ for the nonlinear Krylov–Newton method. Also for other Krylov space methods like preconditioned CG, the most commonly used termination criterion is that of the relative residual, even though the relative residual is not guaranteed to decrease monotonically. Thus, the backward step convergence theory of section \[sec:convergence\] delivers suitable termination criteria for the inner linear iterations. In particular, a rather loose relative stopping criterion of, say, $\kappa =
\frac{1}{10}$ delivers asymptotically already one decimal digit of accuracy per nonlinear iteration.

With this approach, there is one theoretic gap we need to close: The increment $-f(u_k) = -M(u_k) F(u_k) = \delta u_k^m$ depends on the number of Krylov subspace iterations $m$, which is determined adaptively. As a concatenation of a finite number of Lipschitz continuous operations, the $m$-th Krylov iterate depends Lipschitz continuously on $u_k$, but changes in $m$ from one nonlinear $k$-iteration to another can lead to discontinuities in the operator $M$. However, a small modification of the above approach can make sure that $M(u)$ and thus $f(u)$ are Lipschitz continuous with respect to $u$, as required for A\[ass:fLipschitz\]: Instead of using the final iterate $\delta u_k^m$, we could use a linear combination $f(u_k) = (1-\alpha_k) \delta
u_k^{m-1} + \alpha_k \delta u_k^m$ of the two last iterates such that instead of the inequality A\[ass:kappa\] the equality $$\nu(\alpha_k) := {\left\lVertF(u_k) + F'(u_k) \left[ (1-\alpha_k) \delta u_k^{m-1} +
\alpha_k \delta u_k^m \right]\right\rVert}_V = \kappa {\left\lVertF(u_k)\right\rVert}_V$$ holds for some $\alpha_k \in [0, 1]$. This is always possible by virtue of the intermediate value theorem applied to the continuous function $\nu(\alpha)$, which satisfies $\nu(0) > \kappa {\left\lVertF(u_k)\right\rVert}$ and $\nu(1) \le \kappa
{\left\lVertF(u_k)\right\rVert}$ because the Krylov subspace method has terminated in step $m$ but not yet in step $m-1$. If we now assume that there is an upper bound on the number of Krylov iterations required to satisfy the relative termination criterion A\[ass:kappa\] on the level set of $u^0$, we can establish Lipschitz continuity of $M$.

Based on our experience in practical computations, however, it is more efficient to always use $\alpha_k = 1$ and robustify the bisection procedure for the approximate solution to against discontinuities of $g$ by relaxing the lower bound $H^\mathrm{l}$ in closer to zero. This might give rise to smaller than necessary step sizes $t_k$, which has not been observed to be problematic in practical computations, but usually delivers faster local residual contraction once $t_k = 1$.

We report numerical results of a Krylov–Newton method for the Carrier equation in section \[sec:carrier\].

Approximation by discretization {#sec:M_by_discretization}
-------------------------------

Following the multilevel Newton approach of [@Hohmann1994], we can also compute an approximate solution $\delta u_k$ by first discretizing and then (approximately) solving the discretized system. The $\kappa$-condition A\[ass:kappa\] yields a computable criterion for checking if the approximation is accurate enough to ensure convergence in $U$. If not, we need to improve the discretization (and possibly the accuracy of the approximate solution to the resulting finite-dimensional linear system).

In this conceptually simple approach, challenges can arise in the evaluation of the $V$-norms in A\[ass:kappa\]. We address these issues for the case of an elliptic partial differential equation in section \[sec:minsurf\]. Moreover, the evaluation of the $V$-norm in A\[ass:kappa\] can provide a means to adaptively discretize and in turn also the original nonlinear problem.

As in section \[sec:KrylovNewton\], the procedure for approximating a solution to usually involves some discrete decisions, for instance the marking and refinement of certain discretization cells as $k$ increases. Thus, the so constructed $f(u)$ is not Lipschitz continuous. A smoothed formulation akin to the interpolation construction in section \[sec:KrylovNewton\] exceeds the scope of this article and shall be investigated in future work.

Numerical examples from nonlinear elliptic boundary value problems {#sec:ellipticPDE}
==================================================================

In this section, we illustrate the general paradigms presented in section \[sec:M\_by\_discretization\] for the class of elliptic boundary value problems on a bounded domain $\Omega \subset {{\mathbb R}}^n$ with continuously differentiable boundary $\partial \Omega$ that can be cast as nonlinear root-finding problems with the Sobolev spaces $U = H^1_0(\Omega)$ and $V =
H^{-1}(\Omega)$.

Based on the Poincaré inequality (see, e.g., [@Evans2010]), we can use the inner product $${\left(u,v\right)}_U = \int_{\Omega} \nabla u \cdot \nabla v$$ for the Hilbert space $U$. For the Hilbert space $V$, we can then compute norms via the Riesz representation theorem [@Yosida1978 §III.3]: For every $v
\in V$, there exists a uniquely determined $r_v \in U$ such that $$\label{eqn:RieszRepresentation}
  {\left(u, r_v\right)}_U = \int_{\Omega} v u \text{ for all } u \in U
  \quad \text{and} \quad
  {\left\lVertv\right\rVert}_V = {\left\lVertr_v\right\rVert}_U.$$

We investigate the numerical performance of backward step control on two nonlinear elliptic boundary value problems on bounded domains $\Omega$ with continuously differentiable boundary.

All algorithmic parameters of backward step control are chosen as in [@Potschka2016] unless otherwise stated.

Preconditioned Krylov subspace methods {#sec:PGMRES}
--------------------------------------

We first illustrate the general Krylov subspace method approach presented in section \[sec:KrylovNewton\]. As our focus in the case of Krylov–Newton methods here lies on a concise algorithmic statement rather than ultimate computational speed, we use Chebfun [@Battles2004; @Driscoll2014] as an algorithmic tool for the numerical results in section \[sec:carrier\], because it allows to compute numerically with functions (represented as adaptively truncated Chebyshev expansions) instead of numbers [@Trefethen2007] and supports the automatic computation of Fréchet derivatives by the use of automatic differentiation in function space [@Birkisson2012].

Because the linear operators in Chebfun are implemented in strong form, we also use the strong form of the inner product in $U$ $${\left(u,v\right)}_U = -\int_{-1}^{1} u {\Delta}v = -\int_{-1}^{1} v {\Delta}u,$$ (which requires $u$ or $v$ to have square integrable second derivatives). It follows from that $$-\int_{-1}^{1} u {\Delta}r_v = {\left(u, r_v\right)}_U = \int_{-1}^{1} v u 
\quad \text{for all } u \in U,$$ and, thus, $r_v = -{\Delta}^{-1} v$ and $${\left\lVertv\right\rVert}_V^2 = {\left\lVertr_v\right\rVert}_U^2 = {\left\lVert-{\Delta}^{-1} v\right\rVert}_U^2.
$$ In Chebfun, ${\Delta}^{-1}$ can be evaluated efficiently using ultraspherical spectral collocation [@Olver2013]. Based on these prerequisites, we modified Chebfun’s builtin GMRES to use the inner product and norm of $U$ (instead of $L^2(\Omega)$) in combination with ${\Delta}^{-1}$ as a preconditioner, which yields the correct residual norm ${\left\lVert{\Delta}^{-1}v\right\rVert}_U = {\left\lVertv\right\rVert}_V$. We note that we could have used MINRES or even CG because the resulting left-preconditioned linear system is self-adjoint and positive definite in the case at hand [@Guennel2014]. This does not affect our results dramatically, because GMRES produces the same iterates as MINRES for self-adjoint systems. However, the current version of Chebfun does not ship a MINRES implementation. The numerical results for the Carrier equation in section \[sec:carrier\] indicate that using the Riesz isomorphism as a preconditioner works satisfactorily in our example.

### Application to the Carrier equation {#sec:carrier}

For $\varepsilon > 0$, we want to determine a real-valued function $u(x)$ on $x \in [-1, 1]$ that solves the nonlinear second order boundary value problem $$\label{eqn:carrier}
  \varepsilon {\Delta}u + 2 (1-x^2) u + u^2 = 1, \quad u(\pm 1) = 0,$$ which—according to [@Bender1999 §9.7]—is due to Carrier. For small $\varepsilon$, it becomes challenging to solve because of the existence of many local solutions (compare Fig. \[fig:carrier\_solution\]).

We apply the GMRES–Newton method described in section \[sec:PGMRES\] to and illustrate the convergence of the method in Fig. \[fig:carrier\_convergence\] for $\varepsilon
= 10^{-3}$, $\kappa = 10^{-2}$, and varying values of $H = H_\mathrm{rel}
{\left\lVert\delta u_0\right\rVert}_U$, where we choose $H_\mathrm{rel} \in \{ 0.5, 0.1, 0.05,
0.01 \}$. The initial guess is $u_0 = 0$ and the termination criterion is ${\left\lVertF(x_k)\right\rVert}_V \le 10^{-11}$. The relative GMRES termination tolerance $\kappa$ was chosen rather large but at the same time small enough to ensure sufficiently fast local convergence with two decimal digits per iteration. We first observe that no convergence can be obtained for $H_\mathrm{rel} = 0.5$, even though ${\left\lVertF(u_{41})\right\rVert}_V$ drops below $2.6 \cdot 10^{-5}$ and ${\left\lVert\delta u_7\right\rVert}_U \approx 0.28$. From these numbers we can estimate the nonlinearity of the problem in terms of $\omega$ and its well-posedness in terms of $r$ based on Lemma \[lem:Fdecrease\], which yields $$\omega \ge 2 (1 - \kappa) / (t_7 {\left\lVertf(u_7)\right\rVert}_U) \approx 7.1
\quad \text{and} \quad
\omega r \ge 2 (1 - \kappa) / (t_{41} {\left\lVertF(u_{41})\right\rVert}_V) \approx 7.6 \cdot
10^{4}$$ and shows that the problem is highly nonlinear.

For the remaining choices of $H_\mathrm{rel}$ we obtain convergence, albeit a different local solution is found for $H_\mathrm{rel} = 0.1$ than for the others (compare Fig. \[fig:carrier\_solution\]), which nicely confirms the theory of Theorem \[thm:convergenceBSC\]. As guaranteed by Lemma \[lem:localFullSteps\], full steps $t_k = 1$ are taken in the vicinity of a solution and we can clearly observe the asymptotic linear convergence rate of $\kappa = 10^{-2}$ for the residual norm predicted by Lemma \[lem:linearconvergence\]. The final increment norms seem rather large because we do not compute $\delta u_k$ if ${\left\lVertF(u_k)\right\rVert}_V$ is already below $10^{-11}$. Thus, the last increment norm lags behind by one iteration and would be much smaller if we computed it again for the final iterate.

In Fig. \[fig:carrier\_gmres\], we see that the number of GMRES iterations needed in each nonlinear iteration stays moderately small. In each GMRES iteration, one operator-vector-multiplication must be carried out, which we compute via Chebfun as a directional derivative of $F$. In total, 1255 ($H_\mathrm{rel} = 0.1$), 1455 ($H_\mathrm{rel} = 0.05$), and 2471 ($H_\mathrm{rel} = 0.01$) directional derivatives of $F$ are required to solve the problem.

In addition, we can observe from Fig. \[fig:carrier\_gmres\] that the predicted step size $t_k$ needs to be corrected only in few iterations $k$ by the backward step control bisection procedure outlined in section \[sec:realization\]. In the first iteration, four ($H_\mathrm{rel} = 0.1, 0.05$) and five ($H_\mathrm{rel} = 0.01$) bisection steps are required to reduce the initial step size guess $t_0 = 1$. In all other iterations marked by $\bullet$ in Fig. \[fig:carrier\_gmres\], only one additional bisection step is necessary, except in iteration $k=8$ for $H_\mathrm{rel} = 0.1$ and $k=21$ for $H_\mathrm{rel} = 0.01$, where two bisection steps need to be taken. This backs up our claim that the computational overhead of backward step control is small.

Approximation by Finite Elements {#sec:FE}
--------------------------------

In contrast to section \[sec:PGMRES\], we now explicitly discretize the increment $\delta u_k$ and the step determination equation by Finite Elements: To this end, let $\mathcal{C}$ be a partition of $\Omega$ into cells $C \in
\mathcal{C}$. We can then construct the finite-dimensional Finite Element subspace $$U_{\mathcal{C}}^p = \left\{ u \in C^0(\Omega) \mid u \text{ is a polynomial
of degree $p$ on each } C \in \mathcal{C} \right\} \subset U.$$ The increment is then determined by finding $\delta u_k \in U_{\mathcal{C}}^p$ such that $$\label{eqn:discretizedincrement}
  (F'(u_k) \delta u_k) \varphi = -F(u_k) \varphi \quad 
  \text{for all } \varphi \in U_{\mathcal{C}}^p.$$ By fixing a nodal basis of $U_{\mathcal{C}}^p$, we obtain a linear system with a large but sparse (typically symmetric positive definite) ${\left\lvertU_{\mathcal{C}}^p\right\rvert}$-by-${\left\lvertU_{\mathcal{C}}^p\right\rvert}$ matrix.

### Computation of norms in $V$

The discretization in the previous paragraph is completely standard. We now describe the non-standard part, which comprises the computation of $$\label{eqn:kappak}
  {\left\lVertF(u_k)\right\rVert}_V \quad \text{and} \quad 
  \kappa_k := \frac{{\left\lVertF(u_k) + F'(u_k) \delta u_k\right\rVert}_V}{{\left\lVertF(u_k)\right\rVert}_V}$$ in the space $V = H^{-1}(\Omega)$. To this end, we use again the Riesz representation . However, using the same Finite Element subspace $U_{\mathcal{C}}^p$ for the discretization of would yield the wrong value $\kappa_k = 0$ because the numerator vanishes. Moreover, using $U_{\mathcal{C}}^p$ would also give wrong results for ${\left\lVertF(u_k)\right\rVert}_V$ because the residual projected on $U_{\mathcal{C}}^p$ converges locally quadratic (as a Newton method on a finite-dimensional space) if we solve exactly, but does not see the discretization error. Thus, we need to solve on richer Finite Element spaces. The numerical results in section \[sec:minsurf\] indicate that choosing $U_{\mathcal{C}}^{p+1}$ seems to be sufficient for good estimates of the required $V$-norms.

### Adaptive mesh refinement to minimize the contraction rate $\kappa$

Using $p$-refinement for the solution of instead of refinement of $\mathcal{C}$ has two advantages: First, the increase in the degrees of freedom ${\left\lvert\smash[b]{\smash[t]{U_{\mathcal{C}}^{p+1}}}\right\rvert}$ with respect to ${\left\lvertU_{\mathcal{C}}^p\right\rvert}$ is only moderate if $p$ is moderately large, e.g., $p = 3$. Second, the squared norm $\smash[t]{{\left\lVertr_v\right\rVert}}_U^2$ can then be written as a sum of contributions from each cell $C \in \mathcal{C}$, which indicate which cells should ideally be refined if $\kappa_k$ is larger than a desired residual contraction rate $\kappa < 1$ prescribed by the user. The cellwise contributions $\kappa_k$ can be treated in the same fashion as existing cellwise error indicators.

### Application to the minimum surface equation {#sec:minsurf}

In this section, we consider the classical minimum surface problem in the following special form: Let ${\left\lvert.\right\rvert}$ denote the Euclidean norm in $\mathbb{R}^2$ and let $\Omega = \left\{ x \in \mathbb{R}^2 \mid
{\left\lvertx\right\rvert}_2 < 1 \right\}$ and $u^\partial(x) = \sin(2 \pi (x_1 + x_2))$. We seek a function $u$ on $\Omega$ that equals $u^\partial$ on $\partial
\Omega$ and minimizes the area of its graph $$\min I(u) = \int_{\Omega} \sqrt{1 + \smash[t]{{\left\lvert\nabla u\right\rvert}}^2} \quad
  \text{s.t.} \quad u\big|_{\partial\Omega} = u^\partial\big|_{\partial\Omega}.$$ With the spaces $U = H^1_0(\Omega)$ and $V = H^{-1}(\Omega)$ as before, the minimum is described as the solution $u \in u^\partial + U$ to the variational problem $$F(u) \varphi := \int_{\Omega} \nabla \varphi \cdot {\mathbf{g}}(\nabla u)
 = 0 \quad \text{for all } \varphi \in U,
 \quad \text{where } {\mathbf{g}}({\mathbf{v}}) := \left(1 +
 \smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2\right)^{-\frac{1}{2}} {\mathbf{v}}.$$ Thus, $F$ maps $u^\partial + U$ to $V$. Its Gâteaux derivative $F':
(u^\partial + U) \times U \to V$ can then be expressed as $$\left( F'(u) \delta u \right) \varphi
= \int_{\Omega} \nabla \varphi \cdot {\mathbf{g}}'(\nabla u) \nabla \delta u,$$ where the Jacobian of ${\mathbf{g}}$ is given by $${\mathbf{g}}'({\mathbf{v}}) = \frac{1}{\sqrt{1 + \smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2}} \left( {{\mathrm I}}_2 -
  \frac{1}{1 + \smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2} {\mathbf{v}} {\mathbf{v}}^T \right).$$ We recall that $F$ is continuously Gâteaux differentiable from $U = H^1_0(\Omega)$ (with norm $\smash[t]{{\left\lVertu\right\rVert}}_U^2 = \int_{\Omega}
\smash[t]{{\left\lvert\nabla u\right\rvert}}^2$) to $V = H^{-1}(\Omega)$. To see this, we use the chain rule [@Hamilton1982 3.3.4. Thm.] on the continuously Gâteaux differentiable $L^2(\Omega, \mathbb{R}^2)$ inner product and the Nemytskii operator defined by ${\mathbf{g}}$. The Nemytskii operator defined by ${\mathbf{g}}$ mapping from $L^2(\Omega,
\mathbb{R}^2)$ to itself is continuously Gâteaux differentiable by virtue of [@Goldberg1992 Thm. 8, Rem. 6], because ${\mathbf{g}}$ and its Jacobian are uniformly bounded ${\left\lvert{\mathbf{g}}({\mathbf{v}})\right\rvert} \le 1$ and ${\left\lvert{\mathbf{g}}'({\mathbf{v}})\right\rvert}_{2 \times 2}
\le 1$ (where ${\left\lvert.\right\rvert}_{2 \times 2}$ denotes the spectral norm of 2-by-2 matrices), implied by the eigenvalues $$\frac{1}{ \sqrt{1 + \smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2}}
  \left( 1 - \frac{\smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2}{1 + \smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2} \right)
  \quad \text{and} \quad
  \frac{1}{\sqrt{1 + \smash[t]{{\left\lvert{\mathbf{v}}\right\rvert}}^2}},$$ of ${\mathbf{g}}'({\mathbf{v}})$ corresponding to the eigenspace spanned by ${\mathbf{v}}$ and its complement. We remark here that $F$ does not satisfy the stronger property of being continuously Fréchet differentiable as a mapping from $H^1_0(\Omega)$ to $H^{-1}(\Omega)$ as noted in [@Wachsmuth2014].

We solve the resulting system only approximately with a preconditioned Conjugate Gradient (PCG) method [@Hestenes1952]. The resulting inexactness also contributes to the computations of $\kappa_k$ in .

All computations were obtained with the software package deal.II [@Bangerth2016; @Bangerth2007].

### Efficient computation of Riesz representations

As described in section \[sec:FE\], the Riesz representations $r_v$ need to be computed from in order to compute the $V$-norms entering $\kappa_k$. The following algorithmical and computational approaches are important to prevent the computation times for the solution of on the richer Finite Element subspace $U_{\mathcal{C}}^{p+1}$ from dominating the overall computational effort:

1.  We found that PCG with a symmetric Gauss–Seidel smoother as preconditioner delivers good results. In the computations reported below, we employed a Chebyshev smoother of degree four because it yields a slightly better performance than Gauss–Seidel when running the code on several processors in parallel. The use of multigrid does not pay off because it usually involves a rather expensive setup machinery on unstructured meshes (see, e.g., [@Janssen2011]) and the right-hand sides of consist mainly of high-frequency residuals in $U_{\mathcal{C}}^{p+1} \setminus U_{\mathcal{C}}^p$, the low frequency residuals having been mostly eliminated on $U_{\mathcal{C}}^p$ already.

2.  In order to avoid a possible memory bottleneck caused by storing the stiffness matrix discretized on the high-dimensional space $U_{\mathcal{C}}^{p+1}$, we use a matrix-free realization of the Laplacian [@Kronbichler2012].

3.  Because the resulting computation times are then dominated by the bandwidth of the access to main memory, we perform all computations involved in the (matrix-free) matrix-vector products for the solution of only with single instead of double precision floating point arithmetic. The numerical results below indicate that this approach is still sufficiently accurate, while being considerably faster.

4.  Because $\kappa_k$ only steers the algorithm but does not affect the quality of the iterates $x_k$ directly, it can be computed with rather low accuracy requirements in the PCG method. We use relative stopping criteria of $0.1$ and $0.05$ for the numerator and the denominator of $\kappa_k$ in .

### Details about the numerical setup

![Left: The initial mesh for the solution of the minimum surface equation has 425 cells. With $p=3$, the resulting Finite Element space has 2929 degrees of freedom. Right: The accepted and discarded values of $\kappa_k$ in the numerical solution of the minimum surface equation. Discarded values are marked with $\circ$ at fractional iteration numbers, the values of $\kappa_k$ on the finally successful mesh are marked with $\times$ at integer iteration numbers. The discarded values of $\kappa_k$ converge to 1, while the accepted ones approach approximately 0.08.[]{data-label="fig:initialmesh"}](initialmesh){height="30mm"}

We discretize and by nodal Finite Elements of order $p = 3$ and $p
+ 1 = 4$, respectively, on quadrilaterals with tensor product polynomials using Gauss–Lobatto nodes. For the elements on the curved boundary, we employ polynomial tranformations of degree seven from the reference cell to the physical cells. The relative tolerance of PCG for the solution of is $0.001$ (in the Euclidean norm on the discretized vectors).

As a starting guess, we let $u_0$ be the Finite Element interpolation of $u^\partial$ on the coarsest mesh depicted in Fig. \[fig:initialmesh\]. In a first phase, we iterate until the increment norm in $U$ is below 0.01 without computing any $V$-norms. This first phase is a finite-dimensional Newton method ($\kappa = 0$) globalized with backward step control ($H_{\mathrm{rel}} = 0.05$, $t_0 = 1$) on the Finite Element subspace $U_{\mathcal{C}}^p$ belonging to the initial mesh.

For the successive phase of nonlinear adaptive mesh refinement described in section \[sec:FE\], we choose $\kappa = 0.5$. If $\kappa_k > \kappa$, we mark all cells for refinement that have a contribution of more than $2^{-p}$ times the maximum cell contribution to $\kappa_k$, up to a given maximum number of 200,000 cells. The resulting number of cells might be slightly higher in the final mesh due to mesh smoothing in deal.II. In the first eight iterations, the step size is gradually increased from $t_0 =
0.0625$ to $t_7 = 0.9738$. All other iterations are performed with full steps $t_k = 1$. The mesh is refined for the first time in iteration 11, kept for iteration 12, and then successively refined in each further step until the maximum number of cells is reached in iteration 21. We can furthermore observe from Fig. \[fig:initialmesh\] that from iteration 16 on, the first trial value of $\kappa_k$ before the refinement is always one (up to three decimal digits), which shows that no further improvement can be achieved by performing more nonlinear iterations on the current discretization, which is automatically detected correctly by the algorithm.

We compare in Fig. \[fig:minsurfConvergence\] the convergence of the residual norm ${\left\lVertF(u_k)\right\rVert}_V$ (computed afterwards with high accuracy on a refined mesh, which is generated by one additional global refinement step of the triangulation of the final mesh) of backward step control $\kappa$-optimizing adaptive mesh refinement with the convergence when using mesh refinement with the deal.II builtin Kelly error indicator [@Kelly1983; @Gago1983]. In contrast to the theory of backward step control, there is no theoretical guideline for the Kelly indicator on how many nonlinear iterations to run before another round of refinement is triggered. We choose to refine as soon as the increment norm ${\left\lVert\delta
u_k\right\rVert}_U$ becomes less than or equal to a factor $\rho > 0$ of the error estimate returned on the last mesh by the Kelly indicator. Fig. \[fig:minsurfConvergence\] shows that $\kappa$-optimizing adaptive mesh refinement delivers the best ratio of residual norm versus CPU time and versus the number of degrees of freedom compared to mesh refinement based on the Kelly indicator for varying values of $\rho = 0.5, 0.1, 0.01$.

The computations for the solution of the minimum surface equation with a final number of 1.5 million degrees of freedom using $\kappa$-optimizing mesh refinement took 112s wall clock time on the four cores of a mid 2012 MacBook Pro, 2.3 GHz Intel Core i7, 8 GB. Out of this grand total, the computations necessary for estimating $\kappa_k$ took only 18s, even though they need to be performed on the high-dimensional Finite Element space $U_{\mathcal{C}}^{p+1}$.

Conclusions
===========

We presented a comprehensive convergence analysis for , a method that globalizes the convergence of Newton-type methods for the solution of in a Hilbert space setting. We proved that under the reasonable assumptions A\[ass:validIni\]–A\[ass:gamma\] the iterates $u_k$ either leave the region of $r$-regular points $\mathcal{R}_r$ (in which case we need to adjust $M$ or embed $F$ in a suitable homotopy in order to prevent attraction to singularities) or converge to the distinctive solution $u_0^\ast$ (the initial guess $u_0$ propagated by the generalized Newton flow ) provided that $H > 0$ is chosen sufficiently small. Moreover, we provided an $H$-dependent a priori bound on the decrease of ${\left\lVertF(u_k)\right\rVert}_V$ and characterized the asymptotic linear residual convergence rate by $\kappa$. We provided efficient numerical methods based on the blueprint of bounding and optimizing $\kappa$ in each iteration, either over a finite-dimensional subspace in a Krylov–Newton method or through an adaptive Finite Element discretization, in order to balance the nonlinear residual norm with the residual norm of the linear systems. We applied these methods to the class of nonlinear elliptic boundary value problems and presented numerical results for the Carrier equation in a Chebfun implementation and for the minimum surface equation in deal.II. The challenge to efficiently compute norms in $V = H^{-1}(\Omega)$ via the Riesz preconditioner can be addressed by suitable numerical methods and techniques.

The author is grateful to Felix Lenders and to Gerd Wachsmuth for comments on an earlier draft of this manuscript and to the anonymous reviewers for their fruitful comments. This work was funded by the European Research Council through S. Engell’s and H.G. Bock’s ERC Advanced Investigator Grant MOBOCON (291 458) and by the German Federal Ministry of Education and Research under grants 05M2013-GOSSIP and 05M2016-MOPhaPro.

[10]{} \[1\][[\#1]{}]{} urlstyle \[1\][DOI \#1]{}

Ainsworth, M., Oden, J.: A posteriori error estimation in finite element analysis. Pure and Applied Mathematics (New York). Wiley-Interscience \[John Wiley & Sons\], New York (2000)

Amann, H.: Ordinary differential equations, *de Gruyter Studies in Mathematics*, vol. 13. Walter de Gruyter & Co., Berlin (1990)

Ascher, U., Osborne, M.: [A]{} note on solving nonlinear equations and the natural criterion function. ournal of [O]{}ptimization [T]{}heory and [A]{}pplications **55**(1), 147–152 (1987)

Bangerth, W., Davydov, D., Heister, T., Heltai, L., Kanschat, G., Kronbichler, M., Maier, M., Turcksin, B., Wells, D.: The `deal.II` library, version 8.4. Journal of Numerical Mathematics **24** (2016)

Bangerth, W., Hartmann, R., Kanschat, G.: [deal.II]{} – a general purpose object oriented finite element library. ACM Trans. Math. Softw. **33**(4), 24/1–24/27 (2007)

Bank, R., Rose, D.: Global approximate [N]{}ewton methods. Numer. Math. **37**(2), 279–295 (1981)

Battles, Z., Trefethen, L.: An extension of [MATLAB]{} to continuous functions and operators. SIAM J. Sci. Comput. **25**(5), 1743–1770 (2004)

Becker, R., Rannacher, R.: An optimal control approach to a posteriori error estimation in finite element methods. Acta Numer. **10**, 1–102 (2001)

Bender, C., Orszag, S.: Advanced mathematical methods for scientists and engineers. [I]{}. Springer-Verlag, New York (1999). Asymptotic methods and perturbation theory, Reprint of the 1978 original

Birkisson, A., Driscoll, T.: Automatic [F]{}réchet differentiation for the numerical solution of boundary-value problems. ACM Trans. Math. Software **38**(4) (2012)

Bock, H.: [R]{}andwertproblemmethoden zur [P]{}arameteridentifizierung in [S]{}ystemen nichtlinearer [D]{}ifferentialgleichungen, *Bonner Mathematische Schriften*, vol. 183. Universität Bonn, Bonn (1987)

Bock, H., Kostina, E., Schlöder, J.: [O]{}n the [R]{}ole of [N]{}atural [L]{}evel [F]{}unctions to [A]{}chieve [G]{}lobal [C]{}onvergence for [D]{}amped [N]{}ewton [M]{}ethods. In: M. Powell, S. Scholtes (eds.) [S]{}ystem [M]{}odelling and [O]{}ptimization. [M]{}ethods, [T]{}heory and [A]{}pplications, pp. 51–74. Kluwer (2000)

Brezinski, C.: Numerical stability of a quadratic method for solving systems of nonlinear equations. Computing **14**(3), 205–211 (1975)

Davidenko, D.: [O]{}n a new method of numerical solution of systems of nonlinear equations. oklady [A]{}kad. [N]{}auk [SSSR]{} ([N]{}.[S]{}.) **88**, 601–602 (1953)

Dembo, R., Eisenstat, S., Steihaug, T.: Inexact [N]{}ewton methods. SIAM J. Numer. Anal. **19**(2), 400–408 (1982)

Dennis Jr., J., Moré, J.: Quasi-[N]{}ewton methods, motivation and theory. SIAM Rev. **19**(1), 46–89 (1977)

Deuflhard, P.: [A]{} [M]{}odified [N]{}ewton [M]{}ethod for the [S]{}olution of [I]{}ll-conditioned [S]{}ystems of [N]{}onlinear [E]{}quations with [A]{}pplications to [M]{}ultiple [S]{}hooting. umer. [M]{}ath. **22**, 289–311 (1974)

Deuflhard, P.: Global inexact [N]{}ewton methods for very large scale nonlinear problems. Impact Comput. Sci. Engrg. **3**(4), 366–393 (1991)

Deuflhard, P.: [N]{}ewton methods for nonlinear problems, *Springer Series in Computational Mathematics*, vol. 35. Springer, Heidelberg (2011). Affine invariance and adaptive algorithms

Deuflhard, P., Weiser, M.: Global inexact [N]{}ewton multilevel [FEM]{} for nonlinear elliptic problems. In: Multigrid methods [V]{} ([S]{}tuttgart, 1996), *Lect. Notes Comput. Sci. Eng.*, vol. 3, pp. 71–89. Springer, Berlin (1998)

Driscoll, T., Hale, N., Trefethen, L. (eds.): Chebfun guide. Pafnuty Publications, Oxford (2014)

Evans, L.: Partial differential equations, *Graduate Studies in Mathematics*, vol. 19, second edn. American Mathematical Society, Providence, RI (2010)

Gago, J.d.S., Kelly, D., Zienkiewicz, O., Babu[š]{}ka, I.: A posteriori error analysis and adaptive processes in the finite element method. [II]{}. [A]{}daptive mesh refinement. Internat. J. Numer. Methods Engrg. **19**(11), 1621–1656 (1983)

Galántai, A., Abaffy, J.: Always convergent iteration methods for nonlinear equations of [L]{}ipschitz functions. Numer. Algor. pp. 1–11 (2014)

Gasparo, M., Papini, A., Pasquali, A.: Some properties of [GMRES]{} in [H]{}ilbert spaces. Numer. Func. Anal. Opt. **29**(11–12), 1276–1285 (2008)

Goldberg, H., Kampowsky, W., Tröltzsch, F.: On [N]{}emytskij operators in [$L_p$]{}-spaces of abstract functions. Math. Nachr. **155**, 127–140 (1992)

Gr[ä]{}tsch, T., Bathe, K.J.: A posteriori error estimation techniques in practical finite element analysis. Comput. & Structures **83**(4-5), 235–265 (2005)

Günnel, A., Herzog, R., Sachs, E.: A note on preconditioners and scalar products in [K]{}rylov subspace methods for self-adjoint problems in [H]{}ilbert space. Electron. Trans. Numer. Anal. **41**, 13–20 (2014)

Hairer, E., Wanner, G.: Solving Ordinary Differential Equations [I]{}[I]{}, *Springer Series in Computational Mathematics*, vol. 14, second edn. Springer, Berlin (1996)

Hamilton, R.: The inverse function theorem of [N]{}ash and [M]{}oser. Bull. Amer. Math. Soc. (N.S.) **7**(1), 65–222 (1982)

Hestenes, M., Stiefel, E.: Methods of conjugate gradients for solving linear systems. J. Res. Nat. Bur. Standards **49**, 409–436 (1952)

Hohmann, A.: Inexact [G]{}auss [N]{}ewton methods for parameter dependent nonlinear problems. Ph.D. thesis, Freie Universität Berlin (1994)

Janssen, B., Kanschat, G.: Adaptive multilevel methods with local smoothing for [$H^1$]{}- and [$H^{\rm curl}$]{}-conforming high order finite element methods. SIAM J. Sci. Comput. **33**(4), 2095–2114 (2011)

Kelly, D., Gago, J.d.S., Zienkiewicz, O., Babu[š]{}ka, I.: A posteriori error analysis and adaptive processes in the finite element method. [I]{}. [E]{}rror analysis. Internat. J. Numer. Methods Engrg. **19**(11), 1593–1619 (1983)

Kronbichler, M., Kormann, K.: A generic interface for parallel cell-based finite element operator application. Comput. & Fluids **63**, 135–147 (2012)

Lieb, E., Loss, M.: Analysis, *Graduate Studies in Mathematics*, vol. 14, second edn. American Mathematical Society, Providence, RI (2001)

Nevanlinna, O.: Convergence of iterations for linear equations. Lectures in mathematics. Birkh[ä]{}user, Basel, Boston, Berlin (1993)

Olver, S., Townsend, A.: A fast and well-conditioned spectral method. SIAM Rev. **55**(3), 462–489 (2013)

Paige, C., Saunders, M.: [S]{}olutions of sparse indefinite systems of linear equations. . [N]{}umer. [A]{}nal. **12**(4), 617–629 (1975)

Potschka, A.: A direct method for parabolic [PDE]{} constrained optimization problems. Advances in [N]{}umerical [M]{}athematics. Springer (2013)

Potschka, A.: Backward step control for global [N]{}ewton-type methods. SIAM J. Numer. Anal. **54**(1), 361–387 (2016)

Rannacher, R., Vihharev, J.: Adaptive finite element analysis of nonlinear problems: balancing of discretization and iteration errors. J. Numer. Math. **21**(1), 23–61 (2013)

Saad, Y., Schultz, M.: [GMRES]{}: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. **7**(3), 856–869 (1986)

Trefethen, L.: Computing numerically with functions instead of numbers. Math. Comput. Sci. **1**(1), 9–19 (2007)

Wachsmuth, G.: Differentiability of implicit functions: [B]{}eyond the implicit function theorem. J. Math. Anal. Appl. **414**(1), 259–272 (2014)

Yosida, K.: Functional analysis, fifth edn. Springer-[V]{}erlag (1978)

Ypma, T.: Local convergence of inexact [N]{}ewton methods. SIAM J. Numer. Anal. **21**(3), 583–590 (1984)
