Deflation for semismooth equations

TL;DR

This paper introduces a deflation technique combined with semismooth Newton methods to find multiple solutions of variational inequalities without changing initial guesses, proven effective in diverse applications.

Contribution

It presents a novel deflation approach that enables the discovery of multiple solutions in variational inequalities using a single initial guess.

Findings

Deflation effectively eliminates known solutions from consideration.

The method guarantees local superlinear convergence.

Demonstrated success on problems from various fields.

Abstract

Variational inequalities can in general support distinct solutions. In this paper we study an algorithm for computing distinct solutions of a variational inequality, without varying the initial guess supplied to the solver. The central idea is the combination of a semismooth Newton method with a deflation operator that eliminates known solutions from consideration. Given one root of a semismooth residual, deflation constructs a new problem for which a semismooth Newton method will not converge to the known root, even from the same initial guess. This enables the discovery of other roots. We prove the effectiveness of the deflation technique under the same assumptions that guarantee locally superlinear convergence of a semismooth Newton method. We demonstrate its utility on various finite- and infinite-dimensional examples drawn from constrained optimization, game theory, economics and solid mechanics.

Figures

Tables

Table 1. Table 1 : The values of ε 𝜀 \varepsilon at which new solutions were discovered via deflation, along with the number of semismooth Newton iterations required.

Solution	Discovered at	Semismooth Newton iterations
1	$\epsilon =0$	3
2	$\epsilon =4.0\times {10}^{-2}$	13
3	$\epsilon =4.1\times {10}^{-2}$	9
4	$\epsilon =7.0\times {10}^{-2}$	7
5	$\epsilon =7.1\times {10}^{-2}$	13
6	$\epsilon =1.44\times {10}^{-1}$	16
7	$\epsilon =1.45\times {10}^{-1}$	17

Table 2. Table 2 : Number of mesh refinements, degrees of freedom, nonlinear and linear iteration counts required for the continuation in γ 𝛾 \gamma along ε = 0.15 𝜀 0.15 \varepsilon=0.15 . The semismooth Newton solver exhibits γ 𝛾 \gamma - and mesh-independence, while the number of Krylov iterations per semismooth Newton step grows very slowly as γ 𝛾 \gamma and h ℎ h are refined.

$\gamma$	# Refs.	Dofs	Avg SSN its	Avg GMRES/MG its per SSN step
$1.00\times {10}^2$	0	$3.36\times {10}^3$	3.28	-
$1.33\times {10}^2$	0	$3.36\times {10}^3$	4.00	-
$2.23\times {10}^2$	0	$3.36\times {10}^3$	4.50	-
$4.47\times {10}^2$	0	$3.36\times {10}^3$	5.00	-
$1.05\times {10}^3$	0	$3.36\times {10}^3$	4.25	-
$2.80\times {10}^3$	1	$1.31\times {10}^4$	4.14	19.03
$8.42\times {10}^3$	2	$5.18\times {10}^4$	3.85	19.63
$2.81\times {10}^4$	3	$2.06\times {10}^5$	3.85	22.30
$1.03\times {10}^5$	4	$8.21\times {10}^5$	4.14	24.62
$4.13\times {10}^5$	5	$3.28\times {10}^6$	5.00	26.17
$1.00\times {10}^6$	5	$3.28\times {10}^6$	4.00	25.75

Equations

$$⟨Q(u),v-u⟩\ge 0\text{for all }v\in K.$$

$$F(z)=0,$$

$$z_{k+1}=z_k-F^{\prime }(z_k{)}^{-1}F(z_k)$$

$$(Mx+q{)}^T(y-x)\ge 0\text{for all }y\in {\mathbb{R}}_{+}^n,$$

$$G(z)=M(z;r)F(z).$$

$$\underset{z\to r}{\mathrm{liminf}}\Vert M(z;r)F(z){\Vert }_V>0.$$

$$M(z;r)=\frac{{\mathcal{I}}_V}{\Vert z-r{\Vert }_Z^p},$$

$$M(z;r)=\left(\frac{1}{\Vert z-r{\Vert }_Z^p}+1\right){\mathcal{I}}_V.$$

$$F(z+h)-F(z)-H(z+h)h=o(h)$$

$$z_{k+1}=z_k-H(z_k{)}^{-1}F(z_k).$$

$$M(z;r)=\left(\frac{1}{\Vert z-r{\Vert }_Z^p}+\sigma \right){\mathcal{I}}_V,$$

$$T(z)=F(z)-F(r)-H(z)(z-r).$$

$$\Vert M(z;r)F(z){\Vert }_V$$

$$\ge \frac{|\Vert H(z)(z-r){\Vert }_V-\Vert T(z){\Vert }_V|}{\Vert z-r{\Vert }_Z^p}-\sigma \Vert F(z){\Vert }_V$$

$$\ge \frac{\Vert H(z)(z-r){\Vert }_V-\Vert T(z){\Vert }_V}{\Vert z-r{\Vert }_Z^p}-\sigma \Vert F(z){\Vert }_V$$

$$\ge \frac{\gamma \Vert z-r{\Vert }_Z-\Vert T(z){\Vert }_V}{\Vert z-r{\Vert }_Z^p}-\sigma \Vert F(z){\Vert }_V,$$

$$\underset{z\to r}{\mathrm{liminf}}\Vert M(z;r)F(z){\Vert }_V\ge \underset{z\to r}{\mathrm{liminf}}\gamma \Vert z-r{\Vert }_Z^{1-p}>0.$$

$$M(z;r)=\left(\frac{1}{\Vert z-r{\Vert }_Z^p}+\sigma \right){\mathcal{I}}_V,$$

$$H_G(z)h=\left(\sigma +\frac{1}{\Vert z-r{\Vert }_Z^p}\right)H_F(z)h-p\frac{⟨z^{\ast },h{⟩}_{Z^{\ast },Z}}{\Vert z-r{\Vert }_Z^{p+1}}F(z),$$

$$\Vert z^{\ast }{\Vert }_{Z^{\ast }}\le 1,⟨z^{\ast },z-r{⟩}_{Z^{\ast },Z}=\Vert z-r{\Vert }_Z.$$

$$F_{\gamma }(z):=A(z)+\gamma \Phi (z)-f,$$

$$(x)^{\varepsilon}_{+}:=\left\{\begin{array}[]{ll}x-\frac{\varepsilon}{2},&x\geq\varepsilon,\\ \frac{x^{2}}{2\varepsilon},&x\in(0,\varepsilon),\\ 0,&x\leq 0.\end{array}\right.$$

$$\Vert F_{\gamma ,\epsilon }^{\prime }(z_{i,\epsilon }{)}^{-1}\left(F_{\gamma ,\epsilon }^{\prime }(u)-F_{\gamma ,\epsilon }^{\prime }(v)\right)\Vert \le {\omega }_{i,\epsilon }\Vert u-v\Vert \forall u,v\in E_{i,\epsilon }.$$

$$\Vert [(\cdot {)}_{+}^{\epsilon }{]}^{\prime }(u)-[(\cdot {)}_{+}^{\epsilon }{]}^{\prime }(v){\Vert }_{(H^1{)}^{\ast }}$$

$$|[(\cdot {)}_{+}^{\epsilon }{]}^{\prime }(u(x))-[(\cdot {)}_{+}^{\epsilon }{]}^{\prime }(v(x))|$$

$$z\ge 0,F(z)\ge 0,z\perp F(z).$$

$$\text{find }z\in K\text{ s.t. }F(z{)}^T(y-z)\ge 0\text{ for all }y\in K,$$

$$K=\{z\in {\mathbb{R}}^n:z\ge 0\},$$

$${\varphi }_{FB}$$

$${\varphi }_{FB}(a,b)$$

$${\Phi }_i(z)={\varphi }_{FB}(z_i,F_i(z)).$$

$$G(z)=M(z;r_1)M(z;r_2)\cdots M(z;r_n)\Phi (z),$$

$$\displaystyle F(z)=\left[\begin{array}[]{l}3z_{1}^{2}+2z_{1}z_{2}+2z_{2}^{2}+z_{3}+3z_{4}-6\\ 2z_{1}^{2}+z_{2}^{2}+z_{1}+10z_{3}+2z_{4}-2\\ 3z_{1}^{2}+z_{1}z_{2}+2z_{2}^{2}+2z_{3}+9z_{4}-9\\ z_{1}^{2}+3z_{2}^{2}+2z_{3}+3z_{4}-3\end{array}\right].$$

Full text

{amscode}

65K15 Numerical methods for variational inequalities and related problems, 65P30 Bifurcation problems, 65H10 Systems of equations, 35M86 Nonlinear unilateral problems and nonlinear variational inequalities of mixed type, 90C33 Complementarity and equilibrium problems and variational inequalities (finite dimensions).

1 Introduction

Variational inequalities are a fundamental class of problem that arise in many branches of applied mathematics. The problem to be solved is: given a real reflexive Banach space $U$, a closed convex subset $K\subset U$, and an operator $Q:K\to U^{*}$ mapping to the dual space $U^{*}$ of $U$, find $u\in K$ such that

[TABLE]

This is denoted by VI$(Q,K)$. As an elementary example, consider the problem of minimizing a differentiable function $f:\mathbb{R}\to\mathbb{R}$ over the closed interval $I\subset\mathbb{R}$. The necessary condition for $z\in\mathbb{R}$ to be a (local) minimum is that $z$ satisfies VI($f^{\prime},I$). More generally, the minimizers of a general smooth nonlinear program satisfy a variational inequality, which is related to the familiar Karush–Kuhn–Tucker conditions under a suitable constraint qualification, see e.g. [nocedal2006]. Variational inequalities also arise naturally in problems of solid mechanics involving contact [fichera1973, kikuchi1988], in game theory for the calculation of Nash equilibria [harker1984], in phase separation with nonsmooth free energy [blowey1991], and other fields. For more details on variational inequalities, see [facchinei2003, ferris1997, glowinski1984, glowinski1981, harker1990, kinderlehrer2000, lions1967] and the references therein.

A very popular and successful strategy for computing a solution of a variational inequality is to reformulate it as a semismooth rootfinding problem [chen2000, deluca1996, hintermueller2002, kummer1988, klatte2002, qi1993, ulbrich2011], where it is possible to do so. That is, (1) is equivalently reformulated as the task of finding $z\in Z$ such that

[TABLE]

for a residual function $F:Z\to V$, where $Z$ and $V$ are real Banach spaces and $Z$ is reflexive. The space $Z$ is typically constructed via $Z=U\times\Lambda$, where $\Lambda$ is a suitable space of Lagrange multipliers. The residual $F$ may not be differentiable in the classical Fréchet sense but enjoys a weaker property called semismoothness (defined later in Definition 2.3). The problem (2) is constructed in such a way that there is a bijection between solutions of (1) and roots of (2). While the standard Newton–Kantorovich iteration

[TABLE]

requires the existence of the Fréchet derivative $F^{\prime}$ of $F$, it is possible to define a semismooth Newton iteration for semismooth residuals (defined later in section 2.3). This method exhibits locally superlinear convergence under certain regularity conditions on the solution $z$.

Variational inequalities often admit multiple solutions, and these are typically significant for the application at hand. For example, a nonconvex optimization problem may permit several local minima, while a game may permit multiple Nash equilibria. Identifying these distinct solutions is important for understanding the system as a whole. The question of calculating distinct roots of semismooth residuals such as (2) naturally follows. In this paper, we analyze a numerical technique called deflation that can successfully identify multiple solutions of variational inequalities, provided they exist and are isolated from each other.

The central idea of deflation is to compute distinct roots of the semismooth residual $F$ in the following manner. Let us suppose we are given a semismooth residual $F:Z\to V$ and a single known root $r\in Z$ that satisfies some regularity conditions to be made precise later. Deflation constructs a modified residual $G:Z\to V$ with the following properties:

1. (1)

$G$ preserves roots: $G({z})=0\iff F({z})=0$, for ${z}\in Z\setminus\{r\}$, 2. (2)

A semismooth Newton method applied to $G$ from any initial guess in $Z\setminus\{r\}$ will not converge to $r$.

This latter property holds even if semismooth Newton on $G$ is initialized from the same initial guess that led to the convergence to $r$ in the first instance. That is, if semismooth Newton is applied to $G$, and it converges, it will converge to a distinct solution. By enforcing nonconvergence of the semismooth Newton method to known solutions, deflation enables the discovery of unknown ones. The deflated problem is constructed via the application of a deflation operator to the underlying problem $F$.

The idea of deflation was first investigated in the context of differentiable maps $F:\mathbb{R}^{n}\to\mathbb{R}^{n}$ by Brown and Gearhart [brown1971], and was subsequently reinvented in the context of optimization as the tunneling method of Levy and Gómez [levy1985]. Birkisson [birkisson2013] and Farrell et al. [farrell2014] analyzed it in the context of Fréchet differentiable maps between Banach spaces, allowing for its application to smooth partial differential equations. The main contribution of this work is to extend the theory to the case where $F$ is semismooth, but not Fréchet differentiable.

The importance of multiple solutions of variational inequalities has motivated other authors to develop various approaches for computing them. A simple strategy is to vary the initial guess given to the solver [vonstengel2002], but this is heuristic and labour-intensive [tinloi2003]. By contrast, the deflation technique does not involve modifying the initial guess; the problem itself is modified. Judice and Mitra [judice1988] develop an algorithm for enumerating the solutions of linear complementarity problems, variational inequalities of the form: find $x\in\mathbb{R}^{n}_{+}$ such that

[TABLE]

where $\mathbb{R}^{n}_{+}$ is the non-negative orthant of Euclidean space, $M\in\mathbb{R}^{n\times n}$, and $q\in\mathbb{R}$. Their algorithm requires exhaustive exploration of a binary tree whose size is exponential in terms of the size of the problem, and is thus impractical for large problems. Tin-Loi and Tseng [tinloi2003] develop an algorithm for finding multiple solutions of linear complementarity problems by augmenting the problem with constraints that eliminate known solutions; while very successful on the problems considered, the size of each problem increases with each solution eliminated. By contrast, the deflation technique does not increase the size of the problems to be solved after each solution found. Another strategy is to extend classical path-following algorithms to parameter-dependent variational inequalities, and trace out the bifurcation diagram as the parameter is varied [conrad1988, mittelmann1983]. By continuing around turning points, distinct solutions for the same parameter values may be identified; however, this strategy will only identify distinct solutions that happen to lie on a connected branch. The deflation technique enables the discovery of solutions on disconnected branches.

Deflation was first applied in the context of finite-dimensional mixed complementarity problems (a particular kind of variational inequality) by Kanzow [kanzow2000]. Kanzow reports some success on rather difficult test problems, but remarks that deflation is “not …very reliable for larger problems”. This impression is more widely shared: Allgower & Georg [allgower2003] remark (in the context of nonlinear equations, not variational inequalities) that “it is often a matter of seeming chance whether one obtains an additional solution”. We hypothesize that these negative experiences are a consequence of using a poor deflation operator, the norm deflation operator proposed by Brown & Gearhart. We will demonstrate that deflation is much more robust and effective for semismooth problems with an elementary modification to the deflation operator that recovers the correct behaviour of the deflated problem at infinity.

This paper is laid out as follows. In section 2 we define a deflation operator, and prove its effectiveness in the semismooth case. The regularity conditions required on the known solution $z$ are exactly those used to prove locally superlinear convergence of the semismooth Newton method itself in [chen2000, hintermueller2002]; no additional assumptions are required. In section 3, we apply the technique to calculating distinct solutions of several illustrative finite-dimensional variational inequalities, while in section 4 we apply the technique to infinite-dimensional problems with a mesh-independent function-space-based algorithm. These examples demonstrate that the shifted deflation operator applied in this work is effective in numerical practice. We conclude with some remarks and open questions in section 5.

2 Deflation

2.1 Deflation operators

Given a known solution $r$, the deflated problem $G(z)=0$ is constructed by the application of a deflation operator to the original problem $F(z)=0$:

[TABLE]

The requirements for $M$ to enable the discovery of new solutions are captured in the definition of a deflation operator. The definition has been slightly modified from [farrell2014] to allow it to apply to semismooth problems in the sequel.

Definition 2.1** (Deflation operator for isolated solutions).**

Let $F:D\subset Z\to V$, with its domain of definition $D$ an open subset of $Z$. Let $r\in D$ be an isolated root of $F$, i.e. $F(r)=0$ and there exists an open ball around $r$ with no other roots of $F$. We say that $M(\cdot,r):Z\setminus\{r\}\to L(V,V)$ is a deflation operator for $F$ at $r$ if

$M(z;r)\in L(V,V)$ is invertible for all $z\neq r$ in a neighborhood of $r$. 2. 2.

The deflated residual does not converge to zero as $z\to r$:

[TABLE]

Property (6) is referred to as the deflation property. The fundamental example of a deflation operator, proposed by Brown & Gearhart in 1971 [brown1971], is

[TABLE]

where $\mathcal{I}_{V}$ is the identity map on $V$, and $p\geq 1$. The power $p$ controls the rate of blowup as $z\to r$. This is known to be a deflation operator in the case where $F$ is continuously Fréchet differentiable [farrell2014]. This operator was one of two considered by Kanzow, and the operator considered by Allgower & Georg. However, this operator has a major drawback: as $\|z-r\|\to\infty$ in any direction, $M(z;r)\to 0$. This often leads to $G(z)\to 0$ as well, depending on the behaviour of $F(z)$ at infinity. Farrell et al. [farrell2014] suggested a simple modification to the deflation operator to recover the behaviour $M(z;r)\to 1$ and hence $G(z)\to F(z)$ at infinity: the addition of a shift. The shifted deflation operator is

[TABLE]

This is much more effective in numerical practice than (7); the incorrect behaviour at infinity likely accounts for the unsatisfactory performance reported by Kanzow and Allgower & Georg. This will be investigated further in the examples in section 3.

Remark 1*.*

After the deflated problem (5) is constructed, the Newton–Kantorovich or semismooth Newton algorithms will be applied to it, and therefore the differentiability of the deflation operator should be established. These deflation operators (7) and (8) are differentiable away from $z=r$ if the norm used on the Banach space $Z$ is differentiable away from zero, i.e. if the Banach space $Z$ is Fréchet smooth. Note that a reflexive Banach space always admits an equivalent Fréchet smooth norm [fry2002], and hence this requirement is satisfied after possibly renorming. We therefore assume this property henceforth.

Remark 2*.*

In practice only an approximation $\tilde{r}\approx r$ is available for use in the deflation operator. The question then arises of how this approximation affects the computation of the roots of $G$ (e.g. Wilkinson [wilkinson1963, pp. 55] considered this issue in the context of unshifted deflation for polynomial rootfinding). Since the shifted deflation operator (8) satisfies $M(z;\tilde{r})\approx 1$ away from $\tilde{r}$ (and $r$), no difficulties are encountered for solutions that are sufficiently far apart. If two solutions are very close together, a simple remedy discussed by Wilkinson is to calculate a root of $G(z)$, then use that as initial guess for further Newton iterations on $F(z)=0$.

2.2 The Fréchet-differentiable case

For completeness, we state the result arguing that (7) and (8) are deflation operators in the Fréchet-differentiable case. Incidentally, this result can be proven analogously to the semismooth case as in Theorem 2.4 below, which provides an alternative proof to the one found in [farrell2014].

Theorem 2.2** (Deflation for Fréchet differentiable problems [farrell2014]).**

Let $F:D\to V$ be a continuously Fréchet differentiable operator with derivative $F^{\prime}:D\to L(Z,V)$, and let $M$ be given by (7) or (8). Let $r\in D$ be an isolated solution of $F$, i.e. satisfy $F(r)=0$ with $F^{\prime}(r)$ invertible. Then $M$ is a deflation operator for $F$ at $r$.

Remark 3*.*

It may be more convenient to use another norm $\|\cdot\|_{X}$ in the deflation operator, provided $Z\hookrightarrow X$. It is also possible to use a seminorm $|\cdot|_{X}$, provided $\lim_{z\to r}\|T(z)\|_{V}/|z-r|_{X}=0$, where $T(z)$ is the Taylor remainder associated with $F(z)=F(r)+F^{\prime}(r)(z-r)+T(z)$.

2.3 The semismooth case

We now consider the semismooth case.

Definition 2.3** (Semismoothness [chen2000, hintermueller2002, mifflin1977]).**

Let $Z$ and $V$ be Banach spaces. Let $F:D\subset Z\to V$, where $D$ is an open subset of $Z$. $F$ is semismooth at $z\in D$ if it is locally Lipschitz continuous at $z$ and there exists an open neighbourhood $N\subset D$ containing $z$ with a Newton derivative, i.e. a mapping $H:D\to L(Z,V)$ with the property that

[TABLE]

for all $z$ in $N$.

With this Newton derivative, the semismooth Newton iteration is given by

[TABLE]

We now state the main result of this work. The following theorem is novel.

Theorem 2.4** (Deflation for semismooth problems).**

Let $Z$ and $V$ be Banach spaces, and let $F:D\subset Z\to V$. Let $r$ be a root of $F$. Suppose $F$ is semismooth at $r$ with Newton derivative $H:D\to L(Z,V)$ in a neighbourhood $N\subset Z$ around $r$. Further assume that $H(z)$ is invertible for all $z\in N$ and $\{\|H(z)^{-1}\|:z\in N\}$ has a finite upper bound $\Gamma$. Then the operators (7) and (8) are deflation operators for $p\geq 1$.

Remark 4*.*

These are the same assumptions used to prove the locally superlinear convergence of the semismooth Newton method in [chen2000, hintermueller2002].

Proof.

For brevity, define

[TABLE]

with $\sigma=0$ corresponding to (7) and $\sigma=1$ corresponding to (8). Invertibility of $M(z;r)$ for all $z\neq r$ is obvious. Consider $z\in N\setminus\{r\}$. Let $\gamma=\Gamma^{-1}$, and define

[TABLE]

As before, $T(z)=o(z-r)$ from the definition of semismoothness. We then have

[TABLE]

Since $\|F(z)\|_{V}\to 0$ and $\|T(z)\|_{V}/\|z-r\|_{Z}\to 0$ as $z\to r$, we have

[TABLE]

as required. ∎

It remains to show that the deflated problem (5) is in fact semismooth. This property is verified in the following result.

Lemma 2.5**.**

Let $F:D\subset Z\to V$ be semismooth at $z\in D,z\neq r$ with Newton derivative $H_{F}$. Let $M$ be given by

[TABLE]

for some $\sigma\geq 0$. Then the product $G(z)=M(z;r)F(z)$ is also semismooth at $D$ with Newton derivative action

[TABLE]

where $z^{*}\in Z^{*}$ is the derivative of the norm $\|\cdot\|_{Z}$ at $z-r$, i.e. satisfies

[TABLE]

Remark 5*.*

If $Z$ is a Hilbert space, then the Riesz representation of $z^{*}$ is ${(z-r)}/{\|z-r\|}$.

Proof.

This follows from the well-known calculus rules for semismooth and continuously Fréchet differentiable mappings, see e.g., [hintermueller2002, ulbrich2011] ∎

Remark 6*.*

Since the deflated problem is semismooth, the usual sufficient conditions guaranteeing local superlinear convergence may be applied [hintermueller2002] to the deflated residual.

Remark 7*.*

Since the deflated problem is also semismooth, any devices developed for globalizing convergence may be applied, such as line search techniques and continuation e.g. [kanzow2000, ralph1994].

Remark 8*.*

It would be of significant interest to derive sufficient conditions guaranteeing the convergence of the same initial guess to two distinct solutions via deflation. Some initial results in this vein in the smooth case are discussed in [farrell2015d], where it is shown that repeated applications of the well-known Rall-Rheinboldt global convergence theorem [rall1974, rheinboldt1978] can assure convergence to two distinct solutions starting from the same initial guess. A Rall-Rheinboldt-type result, as opposed to Newton-Kantorovich, would be essential to prove convergence to multiple solutions in the context of deflation. This is because the Rall-Rheinboldt theorem places conditions on the radii of convergence of the balls centered at the solutions, rather than guaranteeing the existence of a unique solution in a ball around the initial guess. We briefly investigate the limitations of the classical theory in the framework of semismooth equations below, and in doing so we explain the need for a result of Rall-Rheinboldt-type that is native to the semismooth case.

Many infinite-dimensional semismooth equations of interest share a common structure. In particular, due to low multiplier regularity for bound constrained variational problems, one often resorts to a Moreau-Yosida-type approximation and considers a sequence of (semismooth) equations with residuals taking the form

[TABLE]

where $A$ is a continuously Fréchet differentiable operator, $\Phi$ is a semismooth superposition operator, $\gamma>0$ is a penalty parameter, and $f$ is constant, cf. [hintermueller2006]. For the sake of argument, assume that $\Omega\subset\mathbb{R}^{n}$ is a nonempty, open, and bounded set; $Z=H^{1}(\Omega)$ the usual Sobolev space of $L^{2}$-functions with weak derivatives in $L^{2}$, and $\Phi$ is generated by the function $\phi(x):=\max\{0,x\}$, i.e., $\Phi(z)(x):=\max\{0,z(x)\}$. Since $\phi$ and thus $\Phi$ are nonsmooth, we cannot directly employ the arguments in [farrell2015d].

However, by smoothing the $\max$-function, we can obtain a further approximation of the original problem that is regular enough to exploit the Rall-Rheinboldt theory. The remaining question is whether the convergence guarantees for the smooth problem are stable as $\varepsilon\downarrow 0$. Suppose we replace $\max\{0,x\}$ by

[TABLE]

and define $\Phi_{\varepsilon}(z)(x):=(z(x))^{\varepsilon}_{+}$ and $F_{\gamma,\varepsilon}(z):=A(z)+\gamma\Phi_{\varepsilon}(z)-f$. One of the essential ingredients of the sufficient conditions in the Rall-Rheinboldt theorem are the (local) Lipschitz properties of the derivative $F^{\prime}_{\gamma,\varepsilon}$ at the distinct solutions. In particular, we require an open neighborhood $E_{i,\varepsilon}$ of each solution $z_{i,\varepsilon}$ along with a constant $\omega_{i,\varepsilon}>0$ such that

[TABLE]

Whereas the operator $A$ can be assumed to be unproblematic, one readily derives the estimates

[TABLE]

for any two functions $u,v\in H^{1}(\Omega)$. Although the difference of the smoothed operators is uniformly bounded in $u,v,$ and $\varepsilon\geq 0$, the pointwise relation, which holds as an equality on the set $\left\{x\in\Omega\left|u(x),v(x)\in(0,\varepsilon)\right.\right\}$, would indicate that any form of affine-covariant Lipschitz constant $\omega_{i,\varepsilon}$ would unfavorably depend on $\varepsilon$. As a result, the $\omega_{i,\varepsilon}$-dependent radii associated with the balls of convergence for the distinct roots would converge to zero. Additional assumptions on the structure of $E_{i,\varepsilon}$ that would avoid these issues are unrealistic, e.g., suppose $\Omega\subset\mathbb{R}^{1}$ so that $H^{1}(\Omega)\hookrightarrow C(\overline{\Omega})$ and assume that there exists $\eta>0$ (independent of $\varepsilon$) such that $z_{i,\varepsilon}\leq-\eta<0$.

3 Finite-dimensional examples

We investigate the effectiveness of the deflation approach by applying it to various semismooth problems in the literature that exhibit distinct solutions.

3.1 Complementarity problems

We consider the nonlinear complementarity problem NCP$(F)$: given $F:\mathbb{R}^{n}\to\mathbb{R}^{n}$, find $z\in\mathbb{R}^{n}$ such that

[TABLE]

This is equivalent to the variational inequality

[TABLE]

where

[TABLE]

with all inequalities understood componentwise. We apply a standard semismooth reformulation of the problem using the Fischer–Burmeister NCP function

[TABLE]

which has the property that $\phi_{FB}(a,b)=0\iff a\geq 0,b\geq 0,ab=0$ [fischer1992]. The nonlinear complementarity problem (21) is equivalent to finding roots of the semismooth residual $\Phi:\mathbb{R}^{n}\to\mathbb{R}^{n}$ defined by

[TABLE]

This is then solved with a semismooth Newton method [kanzow2000, kummer1992, qi1993]. As roots $\{r_{1},\dots,r_{n}\}$ are discovered, semismooth Newton is applied to

[TABLE]

where $M$ is given by (8). That is, the deflation operators for each solution are concatenated to deflate all known solutions. Unless noted otherwise, the parameter choice $p=2$ was used.

3.2 Kojima and Shindoh (1986)

This problem was first proposed by Kojima and Shindoh [kojima1986] and is an NCP with $F:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ given by

[TABLE]

It admits two solutions,

[TABLE]

This problem was used again by Dirkse and Ferris [dirske1995] as an example of a problem in which classical Newton solvers struggle to find a solution. This is because one of the two solutions, ${z}^{(2)}$, has a degenerate third component, i.e. ${z}^{(2)}_{3}=F_{3}({z}^{(2)})=0$, and hence does not satisfy strict complementarity. Another feature of this problem is that the linear complementarity problem formed through linearization of the residual $F$ around zero has no solution, causing difficulties for the Josephy–Newton method there [josephy1979].

This is a relatively easy problem to solve and deflation with shifting ($\sigma=1$) successfully finds both solutions from many initial guesses. We chose initial guess $[7/10,\dots,7/10]^{T}$. With no line search, semismooth Newton converged to ${z}^{(1)}$ in 7 iterations; after deflation, semismooth Newton converged to ${z}^{(2)}$ in 12 iterations. By contrast, without shifting ($\sigma=0$) deflation did not identify any additional solutions.

3.3 Gould (2001)

This is a nonconvex quadratic programming problem with linear constraints suggested by N. I. M. Gould in an invited lecture to the 19th biennial conference on numerical analysis [gould2002]. It is a quadratic minimization problem with an indefinite Hessian of the form

[TABLE]

The first order Karush–Kuhn–Tucker optimality conditions yield an NCP with residual

[TABLE]

where $z=[x,\lambda]$, with $\lambda=[\lambda_{1},\lambda_{2}]$ the vector of the Lagrange multipliers associated with $F_{3}(z)\geq 0$ and $F_{4}(z)\geq 0$ respectively. Note that in this case it is not necessary to use Lagrange multipliers to enforce $x\geq 0$ as this is implicit in the NCP formulation. The nonconvexity of the function $f$ makes this problem difficult; it attains two minima with similar functional values and has a saddle point at $x=[1/4,1/2]^{T}$. The central path to be followed by an interior point method is pathological, with different paths converging to the different minima.

We directly solve the arising NCP with the semismooth Newton method with deflation with shifting and without line search. The initial guess was $[2/10,2/10,0,0]^{T}$. In order, the three solutions found were

[TABLE]

These are the saddle point, the global minimum and the local minimum respectively. The KKT conditions make no distinction between minima and saddle points, and hence the solver finds both kinds of stationary points. The number of iterations required was 5, 7 and 10 respectively. As before, without shifting deflation did not successfully identify any additional solutions.

3.4 Aggarwal (1973)

This is a Nash bimatrix equilibrium problem arising in game theory. This kind of problem was first introduced by von Neumann and Morgenstern [vonneumann1945] and the existence of its solutions was further studied by Nash [nash1951] and Lemke and Howson [lemke1964]. In the same paper, Lemke and Howson also presented a numerical algorithm for computing solutions to these kinds of problems. This example was introduced by Aggarwal [aggarwal1973] to prove that it is impossible to find all solutions of such problems using a modification of the Lemke–Howson method that had been conjectured to compute all solutions.

The problem consists of finding the equilibrium points of a bimatrix (non-zero sum, two person) game. Let $A$ and $B$ be the $n\times n$ payoff matrices of players $1$ and $2$ respectively. Let us assume that player $1$ plays the $i^{th}$ pure strategy and player $2$ selects the $j^{th}$ pure strategy amongst the $n$ strategies available to each. The entries of $A$ and $B$, $a_{i,j}$ and $b_{i,j}$ respectively, correspond to the payoff received by each player. It is then possible to define a mixed strategy for a player which consists of a $n\times 1$ vector $x$ such that $x_{i}\geq 0$ and $x_{1}+...+x_{n}=1$. Denote by $x$ and $y$ the mixed strategies for player $1$ and $2$ respectively. The entries of these vectors stand for the probability of the player adopting the corresponding pure strategy. The expected payoffs of the two players are then $x^{T}Ay$ and $x^{T}By$ respectively. An equilibrium point $(x^{*},y^{*})$ is reached when, for all $x$, $y$,

[TABLE]

i.e. neither player can unilaterally improve their payoff.

Aggarwal’s counterexample admits three Nash equilibria. These equilibria are related to the solutions of the NCP with residual

[TABLE]

where $z=[x,y]^{T}$ and $e=[1,1,\dots,1]^{T}$, $\overline{A}$ and $\overline{B}$ are positive-valued loss matrices related to $A$ and $B$ respectively, and $x$ and $y$ relate to the mixed strategy adopted by each player [murty1988, §1.4]. The data for this problem is

[TABLE]

This problem is quite difficult, and we therefore turned to continuation to aid convergence, as described below. The problem was modified to introduce an artificial parameter $\mu$

[TABLE]

with the original problem given by $\mu=1$. With deflation with shifting, three solutions were found for $\mu=1/1000$ from the initial guess $[0,\dots,0]^{T}$, in 5, 24 and 26 iterations of semismooth Newton respectively. (As in the previous examples, deflation without shifting did not identify any additional solutions.) All three branches were then successfully continued to $\mu=1$ using 50 equispaced continuation steps and simple zero-order continuation, i.e. the solution for the previous value $\mu_{-}$ is used as initial guess for the solution of the next value $\mu_{+}$. The three solutions found were

[TABLE]

Aggarwal observed that the conjectured scheme mentioned above could compute ${z}^{(1)}$ and ${z}^{(3)}$, but could not compute ${z}^{(2)}$.

3.5 Gérard, Leclère and Philpott (2017)

Gérard et al. describe a stochastic market where the agents are risk-averse, i.e. estimate their welfare using a coherent risk measure [gerard2017]. They give an example of an incomplete market with three different equilibria, two stable and one unstable. The authors examine the convergence of the well-known PATH solver [dirske1995b, ferris2000], and discover that PATH always yields the unstable equilibrium, even when initialized from many distinct initial guesses. (An alternative tâtonnement algorithm does discover all three equilibria when initialized from different initial guesses.) We therefore investigate whether deflation can assist a semismooth Newton method in discovering all three solutions from a single initial guess.

Mathematically, the problem is a mixed complementarity problem, a generalization of nonlinear complementarity problems. Let $\mathbb{R}_{\infty}:=\mathbb{R}\cup\{-\infty,+\infty\}$. Given $F:\mathbb{R}^{N}\to\mathbb{R}^{N}$, a lower bound $l\in\mathbb{R}_{\infty}^{N}$ and an upper bound $u\in\mathbb{R}_{\infty}^{N}$, the task is to find $z\in\mathbb{R}^{N}$ such that exactly one of the following holds for each $i=1,\dots,N$:

1. (a)

$l_{i}\leq z_{i}\leq u_{i}$ and $F_{i}(z)=0$; 2. (b)

$l_{i}=z_{i}$ and $F_{i}(z)>0$; 3. (c)

$z_{i}=u_{i}$ and $F_{i}(z)<0$.

This is referred to as MCP($F,l,u$). NCP($F$) is a special case with the particular choice $l=[0,\dots,0]^{T}$ and $u=[\infty,\dots,\infty]^{T}$.

The problem at hand is given by $F:\mathbb{R}^{10}\to\mathbb{R}^{10}$, where

[TABLE]

with $z=(x_{0},x_{11},x_{12},y_{1},y_{2},\pi_{1},\pi_{2},u_{4},u_{5},\theta_{P})$. The bounds are given by $l=[0,\dots 0,-\infty]^{T}$ and $u=[\infty,\dots,\infty]^{T}$, i.e. all variables except $\theta_{P}$ have lower bound 0, and $\theta_{P}$ is unconstrained.

To demonstrate that the deflation concept is not confined to a particular semismooth reformulation, in this example we use an alternative NCP function. Define

[TABLE]

which again has the property that $\phi_{\text{MP}}(a,b)=0\iff a\geq 0,b\geq 0,ab=0$. The semismooth reformulation of the MCP employed is

[TABLE]

Deflation with shifting was applied from the initial guess $z_{0}=[0,\dots,0]^{T}$ with $p=1$, and a line search algorithm was used to aid convergence (Alg. 2 of Brune et al. [brune2015]). With these parameters, the procedure identified all three solutions. With $p=2$ or without the line search, only two solutions were found. In order, the three solutions found were

[TABLE]

where only the equilibrium prices are shown for brevity. The solutions were found with 15, 9 and 17 semismooth Newton iterations respectively. The solution found by PATH is the latter. As with the previous examples, deflation without shifting did not identify any additional solutions.

4 Infinite-dimensional examples

4.1 Solving infinite-dimensional variational inequalities

When solving inequality-constrained infinite-dimensional problems, additional care must be taken. The main issue here is a general lack of regularity of the Lagrange multipliers. As a result, it is often impossible to derive a complementarity system (analogous to KKT-conditions for nonlinear programs) that can be reformulated as a single semismooth equation. Even in situations where the associated multiplier is regular enough to allow such a reformulation, we encounter insurmountable issues in the derivation of a function-space-based generalized Newton method. Ignoring these issues and taking a first-discretize-then-optimize approach will generally lead to mesh-dependent convergence, i.e. the number of iterations required to converge increases significantly as the mesh is refined.

This is illustrated in the following example. Suppose $\Omega\subset\mathbb{R}^{n}$ is a nonempty, open, and bounded subset and let $\mathcal{J}\!:\!H^{1}_{0}(\Omega)\to\mathbb{R}$ be Gâteaux differentiable. Consider the model problem

[TABLE]

We denote the feasible set by $K$. If (46) admits a solution $\overline{u}$, then we have

[TABLE]

with $\mathcal{J}^{\prime}(\overline{u})\in H^{-1}(\Omega)$. Since $K$ is a cone, an equivalent formulation holds:

[TABLE]

where $K^{\circ}$ is the polar cone to $K$ given by

[TABLE]

According to the Radon–Riesz theorem, $\lambda\in K^{\circ}$ is in fact a locally finite Radon measure on $\Omega$ [bonnans2000, pg. 564]. Therefore, $\lambda$ cannot in general be evaluated pointwise, in which case (47) cannot be reformulated as a semismooth system of equations.

Suppose further that $\Omega$ is a convex polyhedron and that $\mathcal{J}^{\prime}({u})=A{u}-f$, with $A$ a second-order linear elliptic operator with smooth coefficients and $f\in L^{2}(\Omega)$. Then $\overline{u}\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ [kinderlehrer2000, chap. IV] and thus $\lambda\in L^{2}(\Omega)$. We may then rewrite (47) as

[TABLE]

where $(x)_{+}:=\max\{0,x\}$. Even in this ideal case, the nonsmooth superposition operator $\Phi(\lambda,u):=\lambda-(\lambda-{u})_{+}$ must be defined from $L^{2}(\Omega)\times H^{1}_{0}(\Omega$) into $L^{2}(\Omega)$. The natural choice for the generalized derivative of $\Phi$ is given by

[TABLE]

In order for this to be a Newton derivative, the following approximation property must hold:

[TABLE]

However, this only holds true if $\Phi$ is defined from $L^{2+\varepsilon}(\Omega)\times H^{1}_{0}(\Omega)\to L^{2}(\Omega)$ for $\varepsilon>0$ [hintermueller2002, ulbrich2011], which is not available even in this ideal case. This so-called “missing norm gap” persists for all other known NCP-functions. As a result, the infinite-dimensional problem (48) is not semismooth, which manifests itself as mesh-dependence on the discrete level [hintermueller2004].

An alternative mesh-independent scheme can be constructed from the Moreau–Yosida regularization of the indicator functional for the constraints with respect to the $L^{2}(\Omega)$ topology. The key property of this scheme is that the plus function $(\cdot)_{+}$ is only applied to $u$, and since $u\in H^{1}_{0}(\Omega)\hookrightarrow L^{2+\varepsilon}(\Omega)$, a norm gap holds and this operator is semismooth.

We sketch the approach taken in our implementation; for more details, see [hintermueller2006, MHintermueller_KKunisch_2006]. Let

[TABLE]

where $i_{\mathbb{R}_{+}}(x)$ is the indicator function for $\mathbb{R}_{+}$ (i.e. 0 if $x\leq 0$ and $+\infty$ otherwise). The minimization problem (46) is equivalent to the unconstrained problem

[TABLE]

Approximating $\mathcal{I}(u)$ by its $L^{2}(\Omega)$ Moreau–Yosida regularization

[TABLE]

we obtain a sequence of $\gamma$-dependent problems of the form

[TABLE]

with penalty parameter $\gamma\to\infty$. The associated first-order necessary condition is

[TABLE]

which is semismooth for the reasons outlined above. An initial $\gamma$ is chosen and $\overline{u}_{\gamma}$ computed. Once this is found, the solver continues in $\gamma$, using an analytical path-following scheme to drive the penalty parameter $\gamma\to\infty$ efficiently [hintermueller2006, MHintermueller_KKunisch_2006, LAdam_HHintermueller_TMSurowiec_2018]. The mesh and $\gamma$ are linked; as $\gamma\to\infty$, the mesh is uniformly refined to ensure balanced error estimates (cf. [hintermueller2009]). At every update of $\gamma$, the mesh is refined zero or more times until

[TABLE]

is satisfied, where $h$ is the characteristic mesh size. The process terminates once $\gamma$ reaches a target value $\gamma_{\max}$, which in this work is taken to be $\gamma_{\max}=10^{6}$.

4.2 Computing multiple solutions of infinite-dimensional variational inequalities with deflation

Deflation can be combined with the Moreau–Yosida solver as follows. For the initial value of $\gamma$, a set of initial guesses is supplied. Each guess is used as a starting point for semismooth Newton applied to (53); if the guess is successful, the solution found is deflated, and the guess is attempted again. When all guesses have been exhausted, the analytical path-following strategy is applied to all solutions found, and the next value of $\gamma$ is taken to be the minimum of these (the most conservative update of all solutions found). If necessary, the mesh is refined and all solutions prolonged. The solutions found for the previous step are then used as initial guesses for the next, until the process terminates with $\gamma\geq\gamma_{\max}$.

4.3 Zeidler (1988)

Zeidler [zeidler1988, pg. 320] and Phú [phu1987] study a long thin elastic rod under the action of a compressive load constrained to lie in a channel of fixed width. Let $s\in[0,L]$ denote the arclength of the rod, $y\in H^{1}_{0}(0,L)$ denote its vertical displacement from centerline of the channel, and $\theta\in H^{1}(0,L)$ denote the angle between the rod and the centerline of the channel. The potential energy of the system is given by

[TABLE]

where $B\in\mathbb{R}$ is the bending stiffness of the rod, $P\in\mathbb{R}$ is the compressive load applied to the right end-point, $\rho\in\mathbb{R}$ is the mass per unit length of the rod, and $g\in\mathbb{R}$ is the acceleration due to gravity. The rod is placed in a channel of width $2\alpha$ such that

[TABLE]

is satisfied almost everywhere. Equilibria of the system are therefore given by the local minimizers of

[TABLE]

Zeidler and Phú consider small deformations and linearize the system using the Taylor expansions

[TABLE]

yielding the system

[TABLE]

where irrelevant constant terms in the functional have been neglected. Substituting the constraint $\theta=y^{\prime}$, we arrive at the final system

[TABLE]

In the absence of the inequality constraint on $y$, the optimality conditions for this linearized problem comprise a linear beam equation, with either a unique solution or a one-dimensional nullspace. However, in the presence of the inequality constraint $|y|\leq\alpha$, the problem remains nonlinear and can support distinct isolated solutions.

The system (60) was discretized using $H^{2}(0,L)$-conforming cubic Hermite finite elements [kirby2018a, kirby2018b] using Firedrake [rathgeber2016]. The arising linear systems were solved using the sparse LU factorization of PETSc [balay2017].

In the absence of gravity ($g=0$), Phú proved that the first bifurcation of the system (59) occurs at $P=B\pi^{2}/L^{2}$. In the presence of small gravity, the reflective symmetry of the system is broken and the zero solution is no longer a trivial solution, but the bifurcation point will be nearby. We therefore consider the system for parameter values $B=1,g=1,\rho=1,L=1,\alpha=0.4$ and $P=10.4$. This choice of $P$ is sufficiently greater than $B\pi^{2}/L^{2}\approx 9.87$ that it is reasonable to expect the system to support distinct solutions.

The mesh-independent Moreau–Yosida solver with deflation was applied with initial guess $y=0$ for $\gamma=10$. This converged in one iteration to the first solution, Figure 1(a). This is expected as the inequality constraints are inactive at this solution and the problem is therefore equivalent to the linear beam equation. Deflation with shifting was then applied with deflation operator

[TABLE]

where $r_{y}$ denotes the solution already known. The solver was reinitialized from the zero initial guess and converged after 6 semismooth Newton iterations to a second solution that violates the lower bound. (For this low value of $\gamma$, the bound constraints are only weakly enforced.) This solution was then deflated using the same operator (61) and the solver was re-initialized from the zero initial guess. The procedure then converged after 14 semismooth Newton iterations to the third solution that violates the upper bound111We also experimented with a mesh-dependent semismooth Newton method applied to this problem. The second solution was found after 42 iterations, while the third was found after 45, and both required a line search. In this case the convergence of the mesh-independent scheme is much more robust.. These three solutions were then continued to $\gamma=\gamma_{\max}$ in 9 continuation steps, with no further solutions found. The three solutions found for $\gamma=\gamma_{\max}$ are shown in Figures 1(a)–1(c).

This experiment demonstrates an important property of the deflation strategy: deflation is capable of computing distinct solutions of infinite-dimensional variational inequalities whose solutions exhibit both nontrivial active sets and no activity whatsoever, from the same initial guess.

4.4 A two-dimensional beam under axial compression with obstacle constraints

In this example we consider a two-dimensional analogue of the previous problem, and compute several equilibrium configurations of a hyperelastic beam under axial compression with obstacle constraints.

The physical model employed is compressible neo-Hookean hyperelasticity. Let $\Omega=(0,1)\times(0,1/10)$ denote the undeformed reference configuration, with boundary $\partial\Omega=\partial\Omega_{\textrm{left}}\cup\partial\Omega_{\textrm{bottom}}\cup\partial\Omega_{\textrm{right}}\cup\partial\Omega_{\textrm{top}}$. Homogeneous Dirichlet conditions are imposed on $\partial\Omega_{\textrm{left}}$, axial compression Dirichlet conditions are imposed on $\partial\Omega_{\textrm{right}}$, and natural boundary conditions are imposed on $\partial\Omega_{\textrm{top}}$ and $\partial\Omega_{\textrm{bottom}}$.

In addition, box constraints are imposed on the vertical component $u_{2}$ of the displacement vector field $u:\Omega\to\mathbb{R}^{2}$ on $\partial\Omega_{\textrm{top}}$ and $\partial\Omega_{\textrm{bottom}}$. Let

[TABLE]

where $\mathrm{tr}_{\Gamma}:H^{1}(\Omega)\to H^{1/2}(\Gamma)$ is the standard trace operator, and let $\tau_{\mathrm{bottom}}$ be defined analogously. Since $H^{1/2}(\Gamma)\hookrightarrow L^{2}(\Gamma)$, we may impose pointwise bound constraints of the type $\tau_{\mathrm{top}}(u)\leq\alpha$ and $\tau_{\mathrm{bottom}}(u)\geq\alpha$.

For a given axial compression $\varepsilon$, we seek a displacement vector

[TABLE]

that satisfies

[TABLE]

where $\psi(u)$ is the isotropic compressible neo–Hookean strain energy density with Young’s modulus $E=10^{6}$ and Poisson ratio $\nu=0.3$, $B=(0,-1000)$ is the body force density due to gravity, and $\alpha$ is the value of the bound constraint enforced. The problem is discretized using piecewise linear finite elements; the coarsest grid employed has 3200 triangular elements.

The goal is to solve this problem for $\varepsilon=0.15$ with $\alpha=8\times 10^{-2}$. To do this, continuation is employed. The Moreau–Yosida regularization of (64) is solved with fixed $\gamma=100$ from $\varepsilon=10^{-3}$ to $\varepsilon=0.15$ with steps of $\Delta\varepsilon=10^{-3}$, with deflation employed at each continuation step. This initial continuation process yields seven solutions at $\varepsilon=0.15$, three inactive solutions and four active solutions. The continuation procedure with mesh refinement described in section 4.2 is then applied to continue these solutions from $\gamma=100$ to $\gamma=10^{6}$. On the coarsest grid, the LU algorithm of MUMPS [amestoy2001] is used to solve the linear systems arising in semismooth Newton; on finer meshes GMRES-accelerated geometric multigrid is employed, using all levels in the hierarchy, with three iterations of Chebyshev-accelerated point-block SOR as a smoother (see Ulbrich et al. [ulbrich2017] for rigorous analysis of multigrid in a Moreau–Yosida regularization context). For all solves, the full undamped semismooth Newton step is used, i.e. no line search is found to be necessary for this problem.

This procedure yields seven solutions for $\varepsilon=0.15$ and $\gamma=10^{6}$. To give a sense of the work involved, we report the number of semismooth Newton iterations required for the initial discovery of each solution as $\varepsilon$ is continued in Table 1, and the average number of semismooth Newton iterations per solve and GMRES iterations per semismooth Newton step as $\gamma$ is continued in Table 2. The absolute and relative tolerances of the nonlinear solver were both set to $10^{-8}$, while the absolute and relative tolerances of the linear solver were set to 0 and $10^{-8}$. In all cases iteration counts are modest. In particular, the results of Table 2 show that the number of semismooth Newton iterations required does not increase as $\gamma$ and $h$ are refined, while the number of GMRES-accelerated multigrid V-cycles grows very slowly.

The active solutions are shown in Figure 2 and the inactive solutions are shown in Figure 3. For each active solution, we solve (64) without the obstacle constraints; the corresponding solutions are also plotted to indicate the extent to which the obstacle constraints influence the solutions. As can be seen, the bound constraints significantly change the solutions, and are active on a set of positive measure on the boundary.

These results are encouraging. The function-space-based semismooth Newton method combined with analytical path-following, parameter continuation, multigrid and deflation appears very promising for constrained non-convex variational problems with multiple solutions.

5 Conclusion

Deflation is a useful technique for identifying distinct solutions of variational inequalities with semismooth Newton methods. In particular, employing shifted deflation operators significantly improves the robustness of the approach. The main strengths of the deflation method are that it is effective, straightforward to implement and that it does not significantly increase the cost per Newton iteration.

While the method is found to be effective in numerical experiments, at present no sufficient conditions are known that guarantee convergence of the method to additional solutions. While such conditions are unlikely to be necessary, and may be difficult to verify a priori in computational practice, their availability would establish the foundations of the method and give insight into the design of appropriate deflation operators. The identification of such sufficient conditions forms an important open question and a direction for future research.