# Semismooth and smoothing Newton methods for nonlinear systems with complementarity constraints: Adaptivity and inexact resolution

Ibtihel Ben Gharbia, Joëlle Ferzly, Martin Vohralík, Soleiman Yousef

► To cite this version:

Ibtihel Ben Gharbia, Joëlle Ferzly, Martin Vohralík, Soleiman Yousef. Semismooth and smoothing Newton methods for nonlinear systems with complementarity constraints: Adaptivity and inexact resolution. Journal of Computational and Applied Mathematics, 2023, 420, pp.114765. 10.1016/j.cam.2022.114765. hal-03355116v2

HAL Id: hal-03355116

https://inria.hal.science/hal-03355116v2

Submitted on 17 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Semismooth and smoothing Newton methods for nonlinear systems with complementarity constraints: adaptivity and inexact resolution*

Ibtihel Ben Gharbia† Joëlle Ferzly†‡§ Martin Vohralík‡§ Soleiman Yousef†

June 17, 2022

## Abstract

We consider nonlinear algebraic systems with complementarity constraints stemming from numerical discretizations of nonlinear complementarity problems. The particularity is that they are non-differentiable, so that classical linearization schemes like the Newton method cannot be applied directly. To approximate the solution of such nonlinear systems, an iterative linearization algorithm like the semismooth Newton-min or an interior-point algorithm can be used. Alternatively, the non-differentiable nonlinearity can be smoothed, which allows a direct application of the Newton method. Corresponding linear systems can be solved only approximately using an iterative linear algebraic solver, leading to inexact approaches. In this work, we design a general framework to systematically steer these different ingredients. We first derive an a posteriori error estimate given by the norm of the considered system's residual. We then, relying on smoothing, design a simple strategy of tightening the smoothing parameter. We finally distinguish the smoothing, linearization, and algebraic error components, which enables us to formulate an adaptive algorithm where the linear and nonlinear solvers are stopped when the corresponding error components do not affect significantly the overall error. Numerical experiments indicate that the proposed algorithm, possibly in combination with the GMRES algebraic solver, ensures important savings in terms of the number of iterations and execution time. It appears rather promising in comparison with the other methods, namely since its performance seems remarkably stable over a range of academic and industrial problems.

**Key words:** nonlinear complementarity constraints, semismooth smoothing Newton methods, interior-point method, a posteriori error estimate, adaptivity, stopping criteria

## 1 Introduction

Consider a system of algebraic equations with complementarity constraints written in the following form: Find a vector $X \in \mathbb{R}^n$ such that

$$ \mathbb{E}X = F, \qquad (1.1a) $$

$$ K(X) \geq 0, \quad G(X) \geq 0, \quad K(X) \cdot G(X) = 0, \qquad (1.1b) $$

where for two integers $n > 1$ and $0 < m < n$, $\mathbb{E} \in \mathbb{R}^{n-m,n}$ is a matrix, $K: \mathbb{R}^n \to \mathbb{R}^m$ and $G: \mathbb{R}^n \to \mathbb{R}^m$ are (linear) operators, and $F \in \mathbb{R}^{n-m}$ is a given vector. The first line (1.1a) typically represents the discretization of a linear partial differential equation. The second line (1.1b) then represents the complementarity constraints. It states that the vectors $K(X)$ and $G(X)$ have nonnegative components and are orthogonal.

Complementarity problems have important applications in many fields: economics, engineering, operations research, nonlinear analysis... In the literature, many theoretical results and numerical methods have been proposed to solve problem (1.1), see for example the books of Facchini and Pang [23, 24], Ito and Kunisch [30], Ulbrich [43], Bonnans et al. [14], and the study of Aganagić [1].

<footnote>*This project has partly received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR).</footnote>

<footnote>†IFP Énergies nouvelles, 1 et 4 avenue du Bois Préau, 92852 Rueil-Malmaison, France (ibtihel.ben-gharbia@ifpen.fr, joelle.ferzly@ifpen.fr, soleiman.yousef@ifpen.fr).</footnote>

<footnote>‡Inria, 2 Rue Simone Iff, 75589 Paris, France (martin.vohralik@inria.fr).</footnote>

<footnote>§CERMICS, École des Ponts, 77455 Marne-la-Vallée, France</footnote>
By means of so-called C-functions (C for complementarity), see [23, 24], the complementarity constraints (1.1b) can be rewritten as a system of equations $C(X) = 0$, where $C: \mathbb{R}^n \to \mathbb{R}^n$ is nonlinear and semismooth. We then obtain the following equivalent formulation of problem (1.1): Find $X \in \mathbb{R}^n$ such that

$$ \begin{aligned} \mathbb{E}X &= F, \\ C(X) &= 0. \end{aligned} \tag{1.2} $$

A direct application of the standard Newton method to (1.2) is, however, impeded by the fact that $C(X)$ is not differentiable. An introduction of the Clarke differential [16] allows to give a weaker differentiability meaning and leads to the class of semismooth Newton methods, with reputedly good convergence properties [33, 24, 10, 11, 12, 19, 20]. These methods are in certain cases equivalent to primal-dual active set strategies, see Hintermüller et al. [27]. Moreover, in [47], a regularized semismooth Newton method combined with a hyperplane projection technique was proposed.

Augmented Lagrangian method is one of the commonly used algorithms for constrained optimization, see, e.g., [29] and the references therein. It seeks a solution by replacing the original constrained problem by a series of unconstrained problems and add to the objective function a penalty term, and another term designed to mimic a Lagrange multiplier.

An additional technique, often used in a function space setting, consists in introducing a proper regularization, motivated by the augmented Lagrangian method. It allows to apply an infinite-dimensional semismooth Newton method for the solution of the regularized problem, see, e.g., [43]. In the present context, this leads to replacing the complementarity conditions (1.1b) by

$$ K(X) + \min\{0, -K(X) + \gamma G(X)\} = 0, $$

for a parameter $\gamma > 0$. This method can be combined with a path-following strategy to update the regularization parameter $\gamma$, see for instance [41, 28, 40].

Another important class of methods for constrained optimization problems of the form (1.1) is formed by interior-point methods. These methods consist in generating a sequence in the feasible region $K(X) \ge 0$ and $G(X) \ge 0$, under the assumption of knowing a feasible initial point. We refer to the work of Wright [46], Bellavia et al. [4], and the references therein for a review.

Lastly, an additional notable method is the smoothing Newton method. The main idea of this approach is to approximate the semismooth (non-differentiable) function $C$ from (1.2) by a smooth (differentiable) function that depends on a smoothing parameter. The problem is reformulated as a sequence of regularized smooth equations that can be solved by applying the standard Newton method, and where one drives the smoothing parameter down to zero, cf. [39, 36, 35] and the references therein.

In this work, we design a general framework to systematically steer the above different ingredients. Our main philosophy is adaptive smoothing (regularization). For $\mu^j > 0$, let a smoothed function $C_{\mu^j}(\cdot)$, satisfy $\|C_{\mu^j}(X) - C(X)\| \to 0$ as $\mu^j \to 0$, for $X \in \mathbb{R}^n$. The smoothing parameter $\mu^j$ is reduced at each smoothing iteration $j \ge 1$. Thus, problem (1.1), or equivalently (1.2), can be reformulated as a system of smooth (differentiable) equations written in the form: Find $X^j \in \mathbb{R}^n$ such that

$$ \begin{aligned} \mathbb{E}X^j &= F, \\ C_{\mu^j}(X^j) &= 0. \end{aligned} \tag{1.3} $$

Hence, Newton-type methods can be applied to solve system (1.3), yielding, at each linearization step $k \ge 1$, a linear system

$$ A_{\mu^j}^{j,k-1} X^{j,k} = B_{\mu^j}^{j,k-1}, \tag{1.4} $$

where $A_{\mu^j}^{j,k-1} \in \mathbb{R}^{n,n}$ is a matrix and $B_{\mu^j}^{j,k-1} \in \mathbb{R}^n$ is a vector.

Solving (1.4) with a direct method may be very expensive. A popular approach is to solve it approximately by applying only a few steps of an iterative algebraic solver. Such inexact approaches can be found in [22, 32] for semismooth Newton methods, in [39, 26] for smoothing Newton methods, in [49] for augmented Lagrangian methods, and in [3] for interior-point methods. In the algorithms introduced therein, the iterations of different solvers are stopped according to a fixed maximal number of iterations, the Euclidean norm of the residual vector, or other parameters-dependant stopping criteria. In this work, the a posteriori estimate constitute a distinctive element at the heart of the proposed smoothing method. Importantly, it ensures the desired balance between each source of error at any resolution step, unlike existing approaches based on classical stopping criteria.

At each linear algebraic step $i \ge 1$ for (1.4), one in particular obtains $X^{j,k,i} \in \mathbb{R}^n$ such that

$$ A_{\mu^j}^{j,k-1} X^{j,k,i} = B_{\mu^j}^{j,k-1} - R_{\text{alg}}^{j,k,i}, $$
where $R_{alg}^{j,k,i} \in \mathbb{R}^n$ is the algebraic residual vector of (1.4).

Our principal aim is to reduce the computational cost of the numerical resolution of (1.1) by employing an adaptive strategy based on a posteriori error estimates. There is a well-developed literature on a posteriori error estimates and mesh adaptivity for partial differential equations, see for instance the books of Ainsworth and Oden [2], Repin [37], and Nochetto et al. [34]. For variational inequalities, we can mention the contributions of Repin [38], Ben Belgacem et al. [5], Bürg and Schröder [15], and Dabaghi et al. [17]. Although smoothing Newton approaches have been widely studied, to the best of our knowledge, almost no work has been done to this day on a posteriori error estimates and adaptivity for solvers applied to discrete problems of the form (1.1).

We first derive an upper bound on the norm of the residual of system (1.2), given by

$$R(X^{j,k,i}) := \begin{bmatrix} F - \mathbb{E}X^{j,k,i} \\ -C(X^{j,k,i}) \end{bmatrix}.$$

Then, decomposing $R(X^{j,k,i})$, we distinguish the different error components. This leads to an a posteriori control of the form

$$\|R(X^{j,k,i})\|_{r} \leq \eta^{j,k,i} = \eta_{\text{sm}}^{j,k,i} + \eta_{\text{lin}}^{j,k,i} + \eta_{\text{alg}}^{j,k,i}. \quad (1.5)$$

Here, $\eta^{j,k,i}$ is a fully computable upper bound that holds true at any smoothing (regularization) step $j$, linearization step $k$, and algebraic solver step $i$, whereas the role of the estimators $\eta_{\text{sm}}^{j,k,i}$, $\eta_{\text{lin}}^{j,k,i}$, and $\eta_{\text{alg}}^{j,k,i}$ is to identify the smoothing, linearization, and algebraic components of the error. This error bound allows to define adaptive stopping criteria for the nonlinear and linear algebraic solvers, in the spirit of [21, 17], and the references therein. These criteria, as well as a simple way to tighten the smoothing parameter $\mu^j$, are incorporated in a three-level adaptive algorithm. In contrast to common approaches, where the termination requires reaching a fixed threshold, the particularity of this adaptive algorithm is that the iterations are stopped when the error component of the concerned solver is smaller than the total error, up to a desired fraction. The efficiency of the proposed adaptive algorithm for (inexact) smoothing Newton methods and (inexact) interior-point methods is showcased numerically on practical problems.

It is relevant to mention that this work is extented in [8], where the present approach is applied to a system of PDEs with complementarity constraints in infinite-dimensional space. In particular, taking into account the discretization error allows to adaptively steer the smoothing in system (1.3). Although we do not address mesh adaptivity in our work, we underline that a posteriori estimators are an important tool for adaptive mesh refinement strategies, see, e.g., [18] and the references therein. Consequently, algorithms based on the previous criteria ensure significant computational gains in terms of total number of iterations and mesh cells.

Our manuscript is organized as follows. In Section 2, we recall a semismooth Newton method based on an equivalent reformulation of the complementarity constraints in the form (1.2). Section 3 is devoted to introduce our adaptive inexact smoothing Newton method based on the reformulation as a system of smooth equations as in (1.3). We establish here the a posteriori error estimates (1.5) and propose an adaptive algorithm with a posteriori stopping criteria. We survey a nonparametric interior-point method in Section 4, and introduce its adaptive version in Section 5. Finally, a detailed numerical study is presented in Sections 6 and 7.

## 2 Semismooth Newton method

The purpose of this section is to briefly recall the semismooth Newton method to approximate the solution of the nonlinear system of equations (1.1), see, e.g., [33, 23, 17]. The complementarity constraints represented by (1.1b) as algebraic inequalities are here rewritten as non-differentiable algebraic equalities, using a complementarity function (C-function). A function $\tilde{C}: \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}^m$, $m \ge 1$, is called a C-function if

$$\tilde{C}(x, y) = 0 \iff x \ge 0, y \ge 0, x \cdot y = 0 \quad \forall (x, y) \in \mathbb{R}^m \times \mathbb{R}^m.$$

A variety of C-functions can be found in the literature, see, e.g., [42, 25]. We give as examples the minimum (min) function and the Fischer-Burmeister (F-B) function: for $l = 1, \dots, m$,

$$(\tilde{C}_{\min}(x, y))_l := (\min(x, y))_l = (x_l + y_l)/2 - |x_l - y_l|/2, \quad (2.1)$$

$$(\tilde{C}_{\text{FB}}(x, y))_l := \sqrt{x_l^2 + y_l^2} - (x_l + y_l). \quad (2.2)$$

In general, the C-functions are not Fréchet differentiable. The min and the Fischer-Burmeister functions are, for example, differentiable everywhere except in $x = y$ and $(0,0)$, respectively. Let us introduce a function
C : $\mathbb{R}^n \to \mathbb{R}^m$ defined as $C(X) := \tilde{C}(K(X), G(X))$, where $\tilde{C} : \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}^m$ is any C-function. By using this reformulation in (1.1b), it is obvious that problem (1.1) can be equivalently rewritten as: Find a vector $X \in \mathbb{R}^n$, such that

115

116

117

118

$$ \begin{aligned} \mathbb{E}X &= F, \\ C(X) &= 0. \end{aligned} \tag{2.3a} $$

119

Next, we detail the semismooth Newton linearization. Let an initial vector $X^0 \in \mathbb{R}^n$ be given. At the step $k \ge 1$, one looks for $X^k \in \mathbb{R}^n$ such that

121

122

123

$$ A^{k-1}X^k = B^{k-1}, \tag{2.4} $$

where the square matrix $A^{k-1} \in \mathbb{R}^{n,n}$ and the right-hand side vector $B^{k-1} \in \mathbb{R}^n$ are given by

124

125

$$ A^{k-1} := \left[ \mathbb{E} \mathbf{J}_C(X^{k-1}) \right], \quad B^{k-1} := \left[ \mathbf{F} \mathbf{J}_C(X^{k-1})X^{k-1} - \mathbf{C}(X^{k-1}) \right]. \tag{2.5} $$

Note that the Jacobian corresponding to (2.3a) is constant and equal to $\mathbb{E}$ since it is linear. The semismooth nonlinearity occurs in the second line (2.3b): the notation $\mathbf{J}_C$ in (2.5) stands for the Jacobian matrix in the sense of Clarke of the function $C$, cf. [23, 24]. To give an example, consider the semismooth min function (2.1) and define the matrices $\mathbb{K}$ and $\mathbb{G} \in \mathbb{R}^{m,n}$ respectively by $\mathbb{K} := [\nabla \mathbf{K}(X)]$ and $\mathbb{G} := [\nabla \mathbf{G}(X)]$. Then the $l^{th}$ row of the Jacobian matrix in the sense of Clarke $\mathbf{J}_C$ is either given by the $l^{th}$ row of $\mathbb{K}$, if $(\mathbf{K}(X^{k-1}))_l \le (\mathbf{G}(X^{k-1}))_l$, or by the $l^{th}$ row of $\mathbb{G}$, if $(\mathbf{G}(X^{k-1}))_l < (\mathbf{K}(X^{k-1}))_l$.

126

127

128

129

131

We will need below the total residual vector of problem (2.3), defined by

132

133

$$ \mathbf{R}(V) := \left[ \begin{array}{c} \mathbf{F} - \mathbb{E}V \\ -\mathbf{C}(V) \end{array} \right], \quad V \in \mathbb{R}^n. \tag{2.6} $$

In this context, the relative norm of a vector $V \in \mathbb{R}^n$ is given by $||V||_r := ||V||/||\mathbf{R}(X^0)||$, where $|| \cdot ||$ is the $L_2$-norm.

134

135

# 3 Adaptive inexact smoothing Newton method

136

In this section we introduce our adaptive inexact smoothing Newton method. Based on a posteriori error estimators, adaptive stopping criteria are formulated to conceive an adaptive iterative algorithm.

137

138

## 3.1 Smoothing of the C-functions

139

The key of our developments is to smooth the non-differentiable equation formulation (2.3b) of the complementarity constraints (1.1b) with the help of a smooth (i.e. continuously differentiable) function. This smoothing allows us to approximately transform the nonsmooth nonlinear system (2.3) to a smooth system of nonlinear equations to be solved by using the standard Newton method.

140

141

142

143

Let $\mu > 0$ be a (small) smoothing parameter. We construct an approximation function $\tilde{C}_\mu : \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}^m$ of a C-function $\tilde{C}$ such that $\tilde{C}_\mu(\cdot, \cdot)$ is of class $C^1$ on $\mathbb{R}^m \times \mathbb{R}^m$ and satisfies

144

145

146

$$ \| \tilde{C}(x, y) - \tilde{C}_{\mu}(x, y) \| \to 0 \text{ as } \mu \to 0 \text{ for all } (x, y) \in \mathbb{R}^m \times \mathbb{R}^m. $$

For example, for $l = 1, \dots, m$, a possible smoothing of the min and the Fischer-Burmeister functions (2.1) and (2.2) can be

147

148

$$ (\tilde{C}_{\min_{\mu}}(x, y))_l = \frac{x_l + y_l}{2} - \frac{(|x - y|_{\mu})_l}{2}, \quad \text{with } (|z|_{\mu})_l = \sqrt{z_l^2 + \mu^2}, \tag{3.1} $$

149

$$ (\tilde{C}_{\text{FB}_\mu}(x, y))_l = \sqrt{\mu^2 + x_l^2 + y_l^2 - (x_l + y_l)}, \tag{3.2} $$

150

151

where the $\mu$-smooth absolute value function $| \cdot |_\mu : \mathbb{R}^m \to \mathbb{R}_+^m$, $m \ge 0$, replaces the absolute value function (not differentiable at 0), see Figure 1. Note that both functions $| \cdot |_\mu$ and $\tilde{C}_{\text{FB},\mu}$ are of class $C^\infty$.

152

153
<figure><img src="image_2.png" /><figcaption>Figure 1: Left: Absolute value function $|\cdot|$ and smoothed absolute value function $|\cdot|_{\mu}$. Right: Fischer-Burmeister function $\tilde{C}_{\text{FB}}(\cdot)$ and smoothed Fischer-Burmeister function $\tilde{C}_{\text{FB}_{\mu}}(\cdot)$, for different values of the smoothing parameter $\mu$.</figcaption></figure>

We define the function $C_{\mu}: \mathbb{R}^n \to \mathbb{R}^m$ as $C_{\mu}(X) := \tilde{C}_{\mu}(K(X), G(X))$, where $\tilde{C}_{\mu}: \mathbb{R}^m \times \mathbb{R}^m \to \mathbb{R}^m$ is any smoothed C-function of at least class $C^1$. This allows to approximate problem (1.1) or (1.2) by a system of smooth equations: Find a vector $X \in \mathbb{R}^n$, such that

$$ \begin{aligned} \mathbb{E}X &= F, \\ C_{\mu}(X) &= 0. \end{aligned} \qquad (3.3) $$

Thus, Newton-type methods can be applied to solve the system of nonlinear algebraic equations (3.3).

Fixing $\mu^1 > 0$, we now describe an iterative method for solving problem (2.3). At the beginning of each smoothing iteration (outer iteration) denoted hereafter by $j \ge 1$, an initial guess $X^j \in \mathbb{R}^n$ is given, and a smoothing parameter $\mu^j$ is determined; $\mu^j$ will be driven down to zero. Then some iterative nonlinear solver like the Newton method is employed to solve the smoothed problem written in the form: Find $X^j \in \mathbb{R}^n$ such that

$$ \begin{aligned} \mathbb{E}X^j &= F, \\ C_{\mu^j}(X^j) &= 0. \end{aligned} \qquad (3.4) $$

## 3.2 Newton linearization of the nonlinear algebraic system

In what follows, we detail the Newton method employed to solve problem (3.4) at a fixed outer smoothing step $j \ge 1$. Given an initial vector $X^{j,0}$ (typically $X^{j,0} = X^{j-1}$), Newton's algorithm generates a sequence $(X^{j,k})_{k \ge 1}$ with $X^{j,k} \in \mathbb{R}^n$ given by the following system of linear algebraic equations

$$ A_{\mu^j}^{j,k-1} X^{j,k} = B_{\mu^j}^{j,k-1}, \qquad (3.5) $$

where the Jacobian matrix $A_{\mu^j}^{j,k-1} \in \mathbb{R}^{n,n}$ and the right-hand side vector $B_{\mu^j}^{j,k-1} \in \mathbb{R}^n$ are defined by

$$ A_{\mu^j}^{j,k-1} := \left[ \mathbb{E} \mathbf{J}_{C_{\mu^j}}(\mathbf{X}^{j,k-1}) \right], \quad B_{\mu^j}^{j,k-1} := \left[ \mathbf{F} \mathbf{J}_{C_{\mu^j}}(\mathbf{X}^{j,k-1}) \mathbf{X}^{j,k-1} - C_{\mu^j}(\mathbf{X}^{j,k-1}) \right], \qquad (3.6) $$

with $\mathbf{J}_{C_{\mu^j}}(\mathbf{X}^{j,k-1})$ the Jacobian matrix of the smooth function $C_{\mu^j}$ at $\mathbf{X}^{j,k-1}$.

## 3.3 Inexact solution of the linear algebraic system

The linearized system (3.5) may not be solved exactly, since the use of a direct method may be expensive. For this reason, we consider in this work also an inexact resolution. For a fixed smoothing step $j \ge 1$, a fixed Newton step $k \ge 1$, and an initial guess $X^{j,k,0}$ (typically $X^{j,k,0} = X^{j,k-1}$), only a few steps of an iterative linear algebraic solver can be applied to find an approximate solution to (3.5), yielding, on step $i \ge 1$, an approximation $X^{j,k,i}$ to $X^{j,k}$. This satisfies (3.5) up to the residual vector given by

$$ B_{\mu^j}^{j,k-1} - A_{\mu^j}^{j,k-1} X^{j,k,i}. \qquad (3.7) $$
Define now the linearization function $C_{\mu j}^{j,k-1} : \mathbb{R}^n \to \mathbb{R}^m$ of $C_{\mu j}$ at smoothing step $j$ and Newton step $k$ as

$$C_{\mu j}^{j,k-1}(\mathbf{V}) := C_{\mu j}(\mathbf{X}^{j,k-1}) + \mathbf{J}_{C_{\mu j}}(\mathbf{X}^{j,k-1})(\mathbf{V} - \mathbf{X}^{j,k-1}) \quad \forall \mathbf{V} \in \mathbb{R}^n. \quad (3.8)$$

This allows us to write the algebraic residual vector for $\mathbf{V} \in \mathbb{R}^n$ as

$$\mathbf{R}_{\text{alg}}^{\text{AISN}}(\mathbf{V}) := \mathbf{B}_{\mu j}^{j,k-1} - \mathbf{A}_{\mu j}^{j,k-1} \mathbf{V} = \begin{bmatrix} \mathbf{F} - \mathbb{E}\mathbf{V} \\ -\mathbf{C}_{\mu j}^{j,k-1}(\mathbf{V}) \end{bmatrix}. \quad (3.9)$$

## 3.4 An upper bound for the norm of the residual

We consider the total residual vector of problem (2.3) given in (2.6). By adding and subtracting $C_{\mu j}(\mathbf{X}^{j,k,i})$ and its linearization $C_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i})$ given by (3.8), the total residual vector can be decomposed as follows:

$$\begin{aligned} \mathbf{R}(\mathbf{X}^{j,k,i}) &= \left[ \begin{array}{c} \mathbf{F} - \mathbb{E}\mathbf{X}^{j,k,i} \\ -\mathbf{C}(\mathbf{X}^{j,k,i}) \pm \mathbf{C}_{\mu j}(\mathbf{X}^{j,k,i}) \pm \mathbf{C}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) \end{array} \right] \\ &= \underbrace{\begin{bmatrix} \mathbf{0} \\ \mathbf{C}_{\mu j}(\mathbf{X}^{j,k,i}) - \mathbf{C}(\mathbf{X}^{j,k,i}) \end{bmatrix}}_{\text{smoothing}} + \underbrace{\begin{bmatrix} \mathbf{0} \\ \mathbf{C}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) - \mathbf{C}_{\mu j}(\mathbf{X}^{j,k,i}) \end{bmatrix}}_{\text{linearization}} + \underbrace{\begin{bmatrix} \mathbf{F} - \mathbb{E}\mathbf{X}^{j,k,i} \\ -\mathbf{C}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) \end{bmatrix}}_{\text{algebraic}}. \end{aligned}$$

It is reasonable to get these three terms. Indeed, the first one reflects the error due to the approximation of the semismooth function $C$ by the smoothed function $C_{\mu j}$. The second term is related to the linearization of the nonlinear smooth problem (3.4). Taking into account that the resolution of the smooth linearized problem (3.5) is possibly done “inexactly”, the remaining term represents the error of the inexact algebraic resolution. By the triangle inequality, the relative norm of $\mathbf{R}(\mathbf{X}^{j,k,i})$ is thus bounded by the smoothing, linearization, and algebraic estimators respectively defined as

$$\eta_{\text{sm,AISN}}^{j,k,i} := \|C_{\mu j}(\mathbf{X}^{j,k,i}) - \mathbf{C}(\mathbf{X}^{j,k,i})\|_{\mathrm{r}}, \quad (3.10a)$$

$$\eta_{\text{lin,AISN}}^{j,k,i} := \|C_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) - \mathbf{C}_{\mu j}(\mathbf{X}^{j,k,i})\|_{\mathrm{r}}, \quad (3.10b)$$

$$\eta_{\text{alg,AISN}}^{j,k,i} := \left( \|\mathbf{F} - \mathbb{E}\mathbf{X}^{j,k,i}\|_{\mathrm{r}}^2 + \|\mathbf{C}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i})\|_{\mathrm{r}}^2 \right)^{\frac{1}{2}}. \quad (3.10c)$$

Note that $\eta_{\text{alg,AISN}}^{j,k,i}$ is exactly equal to the relative norm of $\mathbf{R}_{\text{alg}}^{\text{AISN}}(\mathbf{X}^{j,k,i})$ given by (3.9). From these developments we conclude:

**Theorem 3.1** Let $\mathbf{X}^{j,k,i} \in \mathbb{R}^n$ arise from an inexact solve of (3.5). We have

$$\|\mathbf{R}(\mathbf{X}^{j,k,i})\|_{\mathrm{r}} \le \eta_{\text{AISN}}^{j,k,i} := \eta_{\text{sm,AISN}}^{j,k,i} + \eta_{\text{lin,AISN}}^{j,k,i} + \eta_{\text{alg,AISN}}^{j,k,i}.$$

## 3.5 Adaptive inexact smoothing Newton algorithm

Theorem 3.1 motivates the following. Let two real parameters $\alpha_{\text{lin}}$ and $\alpha_{\text{alg}}$ be given in $[0, 1]$, representing the desired relative size of the algebraic and linearization errors, and let $\varepsilon > 0$ be a given desired tolerance for the total error. The stopping criteria for the linearization, algebraic, and smoothing steps, with the bars denoting the stopping indices, are respectively set as

$$\eta_{\text{alg,AISN}}^{j,k,\bar{i}} < \alpha_{\text{alg}} \eta_{\text{lin,AISN}}^{j,k,\bar{i}}, \quad (3.11a)$$

$$\eta_{\text{lin,AISN}}^{j,\bar{k},\bar{i}} < \alpha_{\text{lin}} \eta_{\text{sm,AISN}}^{j,\bar{k},\bar{i}}, \quad (3.11b)$$

$$\|\mathbf{R}(\mathbf{X}^{\bar{j},\bar{k},\bar{i}})\|_{\mathrm{r}} < \varepsilon. \quad (3.11c)$$

The first criterion (3.11a) for the algebraic iterative solver expresses that there is no need to continue with the algebraic steps when the linearization error becomes dominant. Similarly, the second one (3.11b) aims at stopping the linearization iterations when the linearization error does not substantially contribute to the smoothing error. Finally, the termination criterion for the smoothing steps (3.11c) is of the standard type, that is when we stop the entire procedure, when the relative norm of the total residual vector lies below the desired tolerance $\varepsilon$.

The entire method is described by the following adaptive algorithm, which drives the smoothing parameter $\mu^j$ to zero as $\mu^j := \alpha \mu^{j-1}$ at each smoothing iteration. Other common empirical ways to progressively reduce $\mu^j$ can be found, e.g., in [45]. The adaptive inexact smoothing Newton algorithm is the following:
**Algorithm 1** Adaptive inexact smoothing Newton algorithm

**1. Initialization**

Choose a tolerance $\epsilon > 0$ and parameters $\alpha \in ]0, 1[$ and $\alpha_{\text{lin}}, \alpha_{\text{alg}} \in ]0, 1]$. Fix $\mu^1 > 0$ and an initial approximation $\mathbf{X}^0 \in \mathbb{R}^n$. Set $j := 1$.

**2. Smoothing loop**

2.1 Set $\mathbf{X}^{j,0} := \mathbf{X}^0$ as an initial guess for the nonlinear solver. Set $k := 1$.

**2.2 Newton linearization loop**

2.2.1 From $\mathbf{X}^{j,k-1}$ define $A_{\mu^j}^{j,k-1} \in \mathbb{R}^{n,n}$ and $B_{\mu^j}^{j,k-1} \in \mathbb{R}^n$ by (3.6).

2.2.2 Consider the problem of finding a solution $\mathbf{X}^{j,k}$ to

$$A_{\mu^j}^{j,k-1} \mathbf{X}^{j,k} = B_{\mu^j}^{j,k-1}. \qquad (3.12)$$

2.2.3 Set $\mathbf{X}^{j,k,0} := \mathbf{X}^{j,k-1}$ as initial guess for the iterative algebraic solver. Set $i := 1$.

**2.2.4 Algebraic solver loop**

i) Starting from $\mathbf{X}^{j,k-1}$, perform a step of the iterative algebraic solver for the solution of (3.12), yielding, on step $i$ an approximation $\mathbf{X}^{j,k,i}$ to $\mathbf{X}^{j,k}$ satisfying

$$A_{\mu^j}^{j,k-1} \mathbf{X}^{j,k,i} = B_{\mu^j}^{j,k-1} - R_{\text{alg}}^{\text{AISN}}(\mathbf{X}^{j,k,i}).$$

ii) Compute the estimators given in (3.10).

iii) If $\eta_{\text{alg,AISN}}^{j,k,i} < \alpha_{\text{alg}} \eta_{\text{lin,AISN}}^{j,k,i}$, set $\bar{i} := i$ and stop. If not, set $i := i + 1$ and go to i).

2.2.5 If $\eta_{\text{lin,AISN}}^{j,k,i} < \alpha_{\text{lin}} \eta_{\text{sm,AISN}}^{j,k,i}$, set $\bar{k} := k$ and stop. If not, set $k := k + 1$ and go to 2.2.1.

2.3 If $\|R(\mathbf{X}^{j,\bar{k},\bar{i}})\|_{\text{r}} < \epsilon$, set $\bar{j} := j$ and stop.

If not, set $j := j + 1$, $\mathbf{X}^{j,0} := \mathbf{X}^{j-1,\bar{k},\bar{i}}$, and $\mu^j := \alpha \mu^{j-1}$. Then set $k := 1$ and go to 2.2.1.
## 4 Nonparametric interior-point method

Now we employ a nonparametric interior-point method to problem (1.1). More precisely, we consider the method introduced in [44] where a systematic strategy is used to steer the sequence of smoothing parameters towards zero.

We introduce a vector $\mu = \mu 1 \in \mathbb{R}^m$, where $\mu > 0$ is the smoothing parameter and $1 \in \mathbb{R}^m$ is the vector with all components equal to 1. The original nonsmooth problem (1.1) is replaced by a smoothed problem written in the form: Find $\mathbf{X} \in \mathbb{R}^n$ such that

$$ \mathbb{E}\mathbf{X} = \mathbf{F}, \qquad (4.1a) $$

$$ \mathbf{K}(\mathbf{X}) \geq 0, \mathbf{G}(\mathbf{X}) \geq 0, \mathbf{K}(\mathbf{X})\mathbf{G}(\mathbf{X}) = \mu, \qquad (4.1b) $$

where $[(\mathbf{K}(\mathbf{X})\mathbf{G}(\mathbf{X}))]_m = [\mathbf{K}(\mathbf{X})]_m[\mathbf{G}(\mathbf{X})]_m$. In order to properly adjust the sequence of smoothing parameters, the smoothing parameter $\mu$ is treated as an unknown, by introducing the following new equation into system (4.1)

$$ \theta\mu + \mu^2 = 0, \qquad (4.2) $$

where $\theta$ is a small positive real parameter, chosen once and for all. This equation prevents $\mu$ from rushing to zero in just one iteration, and ensures quadratic convergence, see [44]. The unknown of system (4.1) is now the enlarged vector $\mathbf{X} = (\mathbf{X}, \mu)^T \in \mathbb{R}^{n+1}$. We are thus brought back to applying the standard Newton method to a smooth problem.

Let $\mathbf{X}^0 \in \mathbb{R}^n$ such that $\mathbf{K}(\mathbf{X}^0) \geq 0$ and $\mathbf{G}(\mathbf{X}^0) \geq 0$ be given. To update the iterate $\mathbf{X}^{k-1}$, we compute a search direction denoted by $\mathbf{d}^k = [\mathbf{d}_{\mathbf{X}}^k, d_{\mu}^k] \in \mathbb{R}^{n+1}$, where $\mathbf{d}_{\mathbf{X}}^k \in \mathbb{R}^n$ and $d_{\mu}^k \in \mathbb{R}$. Then, to preserve positivity of $\mathbf{K}(\mathbf{X}^k)$ and $\mathbf{G}(\mathbf{X}^k)$ at each step of the nonlinear solver, a truncation of the Newton direction $\mathbf{d}^k$ is performed so that the corresponding update satisfies $\mathbf{K}(\mathbf{X}^{k-1} + \kappa^k \mathbf{d}_{\mathbf{X}}^k) \geq 0$ and $\mathbf{G}(\mathbf{X}^{k-1} + \kappa^k \mathbf{d}_{\mathbf{X}}^k) \geq 0$ for some $\kappa^k \in ]0, 1]$, as close to 1 as possible. After this, we can set

$$ \mathbf{X}^k := \mathbf{X}^{k-1} + \kappa^k \mathbf{d}^k. $$

Recall that our goal is to make $\mu$ equal to 0 in the limit while ensuring the positivity of the updated iterate. Another choice for the additional equation (4.2) added to system (4.1) was developed and introduced in a recent work, see [45, Section 3]. The proposed equation does not require to truncate the Newton direction, and couples $\mu$ and $\mathbf{X}$ in a tighter way.

We rewrite system (4.1) as an enlarged system of $n+1$ equations

$$ \begin{aligned} \mathbb{E}\mathbf{X} &= \mathbf{F}, \\ \mathbf{K}(\mathbf{X})\mathbf{G}(\mathbf{X}) - \mu &= \mathbf{0}, \\ \theta\mu + \mu^2 &= \mathbf{0}. \end{aligned} \qquad (4.3) $$

## 5 Adaptive inexact interior-point method

We present in this section our adaptive inexact version of the nonparametric interior point method of Section 4. In contrast to Section 4, we consider, however, $\mu > 0$ as a parameter, and not as an unknown. At each smoothing step $j \geq 1$, we may solve the system of smoothing equations written as: Find $\mathbf{X}^j \in \mathbb{R}^n$ such that $\mathbf{K}(\mathbf{X}^j) \geq 0, \mathbf{G}(\mathbf{X}^j) \geq 0$, and

$$ \mathbb{E}\mathbf{X}^j = \mathbf{F}, \qquad (5.1a) $$

$$ \mathbf{H}_{\mu^j}(\mathbf{X}^j) := \mathbf{K}(\mathbf{X}^j)\mathbf{G}(\mathbf{X}^j) - \mu^j = \mathbf{0}. \qquad (5.1b) $$

The values of $\mu^j$ are gradually decreased at each smoothing iteration, creating a sequence of suitable $\mu^j$ converging to zero.

### 5.1 Newton linearization of the nonlinear algebraic system

Let $\mathbf{X}^0 \in \mathbb{R}^n$ such that $\mathbf{K}(\mathbf{X}^0) \geq 0$ and $\mathbf{G}(\mathbf{X}^0) \geq 0$ be given. At each smoothing iteration $j \geq 1$ and each linearization step $k \geq 1$, starting with an initial approximation $\mathbf{X}^{j,0}$ such that $\mathbf{K}(\mathbf{X}^{j,0}) \geq 0$ and $\mathbf{G}(\mathbf{X}^{j,0}) \geq 0$ (typically $\mathbf{X}^{j,0} = \mathbf{X}^{j-1}$), we try to approach the solution of problem (5.1) by finding $\mathbf{X}^{j,k} \in \mathbb{R}^n$ such that

$$ \mathbb{A}_{\mu^j}^{j,k-1} \mathbf{X}^{j,k} = \mathbf{B}_{\mu^j}^{j,k-1}, \qquad (5.2) $$
where the Jacobian matrix $A_{\mu j}^{j,k-1} \in \mathbb{R}^{n,n}$ and the right-hand side vector $B_{\mu j}^{j,k-1} \in \mathbb{R}^n$ are defined by

$$A_{\mu j}^{j,k-1} := \left[ \mathbf{J}_{H_{\mu j}}(\mathbf{X}^{j,k-1}) \right], \quad B_{\mu j}^{j,k-1} := \left[ \mathbf{J}_{H_{\mu j}}(\mathbf{X}^{j,k-1})\mathbf{X}^{j,k-1} - \mathbf{H}_{\mu j}(\mathbf{X}^{j,k-1}) \right], \quad (5.3)$$

with $\mathbf{J}_{H_{\mu j}}$ the Jacobian matrix of $\mathbf{H}_{\mu j}$. To ensure the positivity of the complementarity constraints, we then define the direction $\mathbf{d}^{j,k} := \mathbf{X}^{j,k} - \mathbf{X}^{j,k-1} \in \mathbb{R}^n$ and find $\kappa^{j,k} \in ]0, 1]$ such that

$$K(\mathbf{X}^{j,k-1} + \kappa^{j,k} \mathbf{d}^{j,k}) \geq 0 \quad \text{and} \quad G(\mathbf{X}^{j,k-1} + \kappa^{j,k} \mathbf{d}^{j,k}) \geq 0.$$

## 5.2 Inexact solution of the linear algebraic system

For a fixed smoothing iteration $j \ge 1$, a fixed Newton step $k \ge 1$, and an initial guess $\mathbf{X}^{j,k,0}$ (typically $\mathbf{X}^{j,k,0} = \mathbf{X}^{j,k-1}$), an iterative algebraic solver can be applied to approach the solution of (5.2), yielding, on step $i \ge 1$, an approximation $\mathbf{X}^{j,k,i}$ to $\mathbf{X}^{j,k}$. This satisfies (5.2) up to a residual vector defined by

$$B_{\mu j}^{j,k-1} - A_{\mu j}^{j,k-1} \mathbf{X}^{j,k,i}. \quad (5.4)$$

Introduce the linearization $\mathbf{H}_{\mu j}^{j,k-1} : \mathbb{R}^n \to \mathbb{R}^m$ of $\mathbf{H}_{\mu j}(\cdot)$ such that for $\mathbf{V} \in \mathbb{R}^n$,

$$\mathbf{H}_{\mu j}^{j,k-1}(\mathbf{V}) := \mathbf{H}_{\mu j}(\mathbf{X}^{j,k-1}) + \mathbf{J}_{H_{\mu j}}(\mathbf{X}^{j,k-1})(\mathbf{V} - \mathbf{X}^{j,k-1}). \quad (5.5)$$

Using (5.5), the algebraic residual vector can be written as follows

$$\mathbf{R}_{\text{alg}}^{\text{AIIP}}(\mathbf{V}) := \mathbf{B}_{\mu j}^{j,k-1} - A_{\mu j}^{j,k-1} \mathbf{V} = \left[ \begin{array}{c} \mathbf{F} - \mathbf{E}\mathbf{V} \\ -\mathbf{H}_{\mu j}^{j,k-1}(\mathbf{V}) \end{array} \right], \quad \mathbf{V} \in \mathbb{R}^n. \quad (5.6)$$

We now define the function $\mathbf{H} : \mathbb{R}^n \to \mathbb{R}^m$ by

$$\mathbf{H}(\mathbf{V}) := \mathbf{K}(\mathbf{V})\mathbf{G}(\mathbf{V}), \quad \mathbf{V} \in \mathbb{R}^n. \quad (5.7)$$

and the total residual vector associated to the adaptive inexact interior-point method by

$$\mathbf{R}^{\text{AIIP}}(\mathbf{V}) := \begin{bmatrix} \mathbf{F} - \mathbf{E}\mathbf{V} \\ -\mathbf{H}(\mathbf{V}) \end{bmatrix}, \quad \mathbf{V} \in \mathbb{R}^n. \quad (5.8)$$

Here again, the relative norm of a given vector $\mathbf{V} \in \mathbb{R}^n$ is given by $||\mathbf{V}||_r := ||\mathbf{V}|| / ||\mathbf{R}^{\text{AIIP}}(\mathbf{X}^0)||_r$.

## 5.3 An upper bound for the norm of the residual

In the same spirit as in Section 3.4, we decompose at each smoothing step $j \ge 1$, each linearization step $k \ge 1$, and each algebraic step $i \ge 1$ the total residual vector given by (5.8)

$$\mathbf{R}^{\text{AIIP}}(\mathbf{X}^{j,k,i}) = \underbrace{\begin{bmatrix} \mathbf{0} \\ \mathbf{H}_{\mu j}(\mathbf{X}^{j,k,i}) - \mathbf{H}(\mathbf{X}^{j,k,i}) \end{bmatrix}}_{\text{smoothing}} + \underbrace{\begin{bmatrix} \mathbf{0} \\ \mathbf{H}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) - \mathbf{H}_{\mu j}(\mathbf{X}^{j,k,i}) \end{bmatrix}}_{\text{linearization}} + \underbrace{\begin{bmatrix} \mathbf{F} - \mathbf{E}\mathbf{X}^{j,k,i} \\ -\mathbf{H}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) \end{bmatrix}}_{\text{algebraic}}.$$

We then define the smoothing, linearization, and algebraic estimators by

$$\eta_{\text{sm,AIIP}}^{j,k,i} := ||\mathbf{H}_{\mu j}(\mathbf{X}^{j,k,i}) - \mathbf{H}(\mathbf{X}^{j,k,i})||_r = ||\boldsymbol{\mu}^j||_r, \quad (5.9a)$$

$$\eta_{\text{lin,AIIP}}^{j,k,i} := ||\mathbf{H}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i}) - \mathbf{H}_{\mu j}(\mathbf{X}^{j,k,i})||_r, \quad (5.9b)$$

$$\eta_{\text{alg,AIIP}}^{j,k,i} := \left( \|\mathbf{F} - \mathbf{E}\mathbf{X}^{j,k,i}\|_r^2 + \|\mathbf{H}_{\mu j}^{j,k-1}(\mathbf{X}^{j,k,i})\|_r^2 \right)^{\frac{1}{2}}. \quad (5.9c)$$

Then we have an upper bound for the norm $||\mathbf{R}^{\text{AIIP}}(\mathbf{X}^{j,k,i})||_r$:

**Theorem 5.1** Let $\mathbf{X}^{j,k,i} \in \mathbb{R}^n$ be the approximation of $\mathbf{X}$ given by an iterative algebraic solver. Then we have

$$||\mathbf{R}^{\text{AIIP}}(\mathbf{X}^{j,k,i})||_r \le \eta_{\text{AIIP}}^{j,k,i} := \eta_{\text{sm,AIIP}}^{j,k,i} + \eta_{\text{lin,AIIP}}^{j,k,i} + \eta_{\text{alg,AIIP}}^{j,k,i}.$$
## 5.4 Adaptive inexact interior-point algorithm

Our proposed adaptive inexact interior-point algorithm implements adaptive stopping criteria formulated using the error component estimators given by (5.9) is as follows:

<caption>**Algorithm 2** Adaptive inexact interior-point algorithm</caption>

**1. Initialization**

Choose a tolerance $\varepsilon > 0$ and parameters $\alpha \in ]0, 1[$ and $\alpha_{\text{lin}}, \alpha_{\text{alg}} \in ]0, 1]$.

Fix $\mu^1 > 0$ and an initial vector $\mathbf{X}^0 \in \mathbb{R}^n$ such that $\mathbf{K}(\mathbf{X}^0) \ge \mathbf{0}$ and $\mathbf{G}(\mathbf{X}^0) \ge \mathbf{0}$. Set $j := 1$.

**2. Smoothing loop**

2.1 Set $\mathbf{X}^{j,0} := \mathbf{X}^0$ as an initial guess for the linearization loop and $k := 1$.

**2.2 Interior-point linearization loop**

2.2.1 From $\mathbf{X}^{j,k-1}$ define $\mathbf{A}_{\mu^j}^{j,k-1} \in \mathbb{R}^{n,n}$ and $\mathbf{B}_{\mu^j}^{j,k-1} \in \mathbb{R}^n$ by (5.3).

2.2.2 Consider the problem of finding $\mathbf{X}^{j,k} \in \mathbb{R}^n$ such that

$$ \mathbf{A}_{\mu^j}^{j,k-1} \mathbf{X}^{j,k} = \mathbf{B}_{\mu^j}^{j,k-1}. \qquad (5.10) $$

2.2.3 Set $\mathbf{X}^{j,k,0} := \mathbf{X}^{j,k-1}$ as initial guess for the iterative algebraic solver. Set $i := 1$.

**2.2.4 Algebraic solver loop**

i) Starting from $\mathbf{X}^{j,k-1}$ perform a step of the iterative algebraic solver for (5.10), yielding, at step $i \ge 1$, a vector $\mathbf{X}^{j,k,i} \in \mathbb{R}^n$ such that

$$ \mathbf{A}_{\mu^j}^{j,k-1} \mathbf{X}^{j,k} = \mathbf{B}_{\mu^j}^{j,k-1} - \mathbf{R}_{\text{alg}}^{\text{AIIP}}(\mathbf{X}^{j,k,i}). $$

ii) Set $\mathbf{d}^{j,k,i} := \mathbf{X}^{j,k} - \mathbf{X}^{j,k-1}$ and compute $\kappa^{j,k,i} \in ]0, 1]$ such that

$$ \mathbf{K}(\mathbf{X}^{j,k-1} + \kappa^{j,k,i} \mathbf{d}^{j,k,i}) \ge \mathbf{0} \text{ and } \mathbf{G}(\mathbf{X}^{j,k-1} + \kappa^{j,k,i} \mathbf{d}^{j,k,i}) \ge \mathbf{0}. $$

Then set $\mathbf{X}^{j,k,i} := \mathbf{X}^{j,k-1} + \kappa^{j,k,i} \mathbf{d}^{j,k,i}$.

iii) Compute the estimators given by (5.9).

iv) If $\eta_{\text{alg,AIIP}}^{j,k,i} < \alpha_{\text{alg}} \eta_{\text{lin,AIIP}}^{j,k,i}$, set $\bar{i} := i$ and stop. If not, set $i := i + 1$ and go to i).

2.2.5 If $\eta_{\text{lin,AIIP}}^{j,k,\bar{i}} < \alpha_{\text{lin}} \eta_{\text{sm,AIIP}}^{j,k,\bar{i}}$, set $\bar{k} := k$ and stop. If not, set $k := k + 1$ and go to 2.2.1.

2.3 If $\| \mathbf{R}^{\text{AIIP}}(\mathbf{X}^{j,\bar{k},\bar{i}}) \|_{\text{r}} < \varepsilon$, set $\bar{j} := j$ and stop. If not, set $j := j + 1$, $\mathbf{X}^{j,0} := \mathbf{X}^{j-1,\bar{k},\bar{i}}$, and $\mu^j := \alpha \mu^{j-1}$. Then set $k := 1$ and go to 2.2.1.

# 6 Numerical experiments: Problem of contact between two membranes

This section reports some numerical illustrations obtained using the algorithms previously presented. We consider here the model problem of contact between two membranes.

## 6.1 Problem statement

Let $\Omega = (a, b)$ be a one-dimensional domain. The problem reads: Find $u_1, u_2$, and $\lambda$ such that

$$ \left\{ \begin{array}{lll} -\mu_1 \Delta u_1 - \lambda = f_1 & \text{in } \Omega, \\ -\mu_2 \Delta u_2 + \lambda = f_2 & \text{in } \Omega, \\ (u_1 - u_2)\lambda = 0, & u_1 - u_2 \ge 0, & \lambda \ge 0 \text{ in } \Omega, \\ u_1 = g & \text{on } \partial\Omega, \\ u_2 = 0 & \text{on } \partial\Omega, \end{array} \right. \qquad (6.1) $$

306

307
where $u_1$ and $u_2$ represent the vertical displacements of the two membranes and $\lambda$ is a Lagrange multiplier characterizing the action of the second membrane on the first one, $-\lambda$ being the reaction. The constant parameters $\mu_1, \mu_2 > 0$ correspond to the tension of each membrane, whereas $f_1, f_2 \in L^2(\Omega)$ are given external forces. The boundary condition prescribed by a constant $g > 0$ ensures that, on the boundary $\partial\Omega$, the first membrane is above the second one. The third line of (6.1) represents the linear complementarity conditions which serve to distinguish two different physical situations: either the membranes are separated ($u_1 > u_2$ and $\lambda = 0$), or they are in contact ($u_1 = u_2$ and $\lambda > 0$). We discretize this problem by the finite volume method. The corresponding discretization can be written under the form of problem (1.1).

## 6.2 Test problem setting

Following [5], we set $\Omega = (-1, 1)$ and consider the following analytical solution for $x \in \Omega$

$$u_1(x) := g(2x^2 - 1), \quad u_2(x) := \begin{cases} 2g(1 - x^2)(2x^2 - 1) & \text{if } x < \frac{-1}{\sqrt{2}} \text{ or } x > \frac{1}{\sqrt{2}}, \\ g(2x^2 - 1) & \text{otherwise,} \end{cases}$$

$$\lambda(x) := \begin{cases} 0 & \text{if } x < \frac{-1}{\sqrt{2}} \text{ or } x > \frac{1}{\sqrt{2}}, \\ 2g & \text{otherwise.} \end{cases}$$

This triple is the solution of (6.1) for the data $f_1$ and $f_2$ given by

$$f_1(x) := \begin{cases} -4g & \text{if } x < \frac{-1}{\sqrt{2}} \text{ or } x > \frac{1}{\sqrt{2}}, \\ -6g & \text{otherwise,} \end{cases} \quad \text{and} \quad f_2(x) := \begin{cases} -12g(1 - 4x^2) & \text{if } x < \frac{-1}{\sqrt{2}} \text{ or } x > \frac{1}{\sqrt{2}}, \\ -2g & \text{otherwise.} \end{cases}$$

Throughout the computational experiments, the parameters $\mu_1$ and $\mu_2$ are set to 1 and the boundary condition $g$ for the first membrane is taken equal to 0.1. Let $N$ be the number of mesh elements. The initial guess $\mathbf{X}^0 \in \mathbb{R}^{3N}$ has its first $N$ components equal to $g$ and its other components equal to zero for the semismooth and smoothing Newton methods. For the nonparametric interior-point method (resp. the adaptive interior-point method), the initialization is given by $\mathbf{X}^0 = [0.1 \ 0.5 \ 0.05]^T \in \mathbb{R}^{3N+1}$ (resp. $\mathbf{X}^0 = [0.1 \ 0.5]^T \in \mathbb{R}^{3N}$). All the simulations are performed in MATLAB. We consider $N = 25000$ elements, leading to the matrix $\mathbf{A}$ of size $n = 75000$.

## 6.3 Semismooth Newton method

We start by presenting the numerical results of the semismooth Newton method described in Section 2, using the F-B function (2.2). The stopping criterion is on the total residual vector (2.6)

$$\| \mathbf{R}(\mathbf{X}^k) \|_r < 10^{-8}. \qquad (6.2)$$

To achieve this stopping criterion, 527 semismooth Newton-F-B iterations (CPU time: 68.9s) and 2232 Newton-min iterations (CPU time: 338.9s) are needed. Figure 2 represents the evolution of $\| \mathbf{R}(\mathbf{X}^k) \|_r$ as a function of the semismooth Newton-F-B iterations. We can see that it decreases slowly during iterations, then the convergence gets extremely fast at the end.
<figure><img src="image_3.png" /><figcaption>Figure 2: [Semismooth Newton method, F-B function (2.2), stopping criterion (6.2)] Relative norm of the total residual vector (2.6) as a function of semismooth Newton iterations.</figcaption></figure>

## 6.4 Adaptive smoothing Newton method

We now test the adaptive smoothing Newton method, denoted by ASN, with the smoothed F-B function (3.2). This consists in employing the method presented in Section 3, summarized in Algorithm 1, but with an exact resolution of the nonlinear system (3.5). The linearization and smoothing estimators are respectively defined by

$$ \eta_{\text{lin,ASN}}^{j,k} := \|C_{\mu^j}(\mathbf{X}^{j,k})\|_{\text{r}}, \qquad (6.3a) $$

$$ \eta_{\text{sm,ASN}}^{j,k} := \|C_{\mu^j}(\mathbf{X}^{j,k}) - \mathbf{C}(\mathbf{X}^{j,k})\|_{\text{r}}, \qquad (6.3b) $$

and the total estimator by $\eta_{\text{ASN}}^{j,k} := \eta_{\text{sm,ASN}}^{j,k} + \eta_{\text{lin,ASN}}^{j,k}$.

First, we analyze the performance of the adaptive stopping criterion based on the estimators for stopping the linearization steps. We compare it with the classical approach in where the linearization is continued until the relative norm of the linearization estimator becomes smaller than a threshold taken as $10^{-8}$, i.e.,

$$ \text{Classical stopping criterion: } \eta_{\text{lin,ASN}}^{j,k} < 10^{-8}, \qquad (6.4) $$

$$ \text{Adaptive stopping criterion: } \eta_{\text{lin,ASN}}^{j,k} < \alpha_{\text{lin}} \eta_{\text{sm,ASN}}^{j,k} \qquad (6.5) $$

We set $\mu^1 = 1$, $\varepsilon = 10^{-8}$, $\alpha_{\text{lin}} = 1$, and $\alpha = 0.1$ in Algorithm 1. Figure 3 depicts the evolution of the estimators and the relative norm of the total residual vector $\mathbf{R}(\mathbf{X}^{j,k})$ given in (2.6) as a function of the smoothing Newton-F-B iterations, at a specific smoothing iteration $j = 1$ ($\mu^1 = 1$), left, and $j = 3$ ($\mu^3 = 10^{-2}$), right. We can observe from Figure 3, left, that, as expected, the smoothing estimator and $\|\mathbf{R}(\mathbf{X}^{j,k})\|_{\text{r}}$ stagnates after few steps, since here the smoothing parameter $\mu^1$ is equal to 1, whereas the linearization estimator steadily decreases. If we consider the classical stopping criterion (6.4), the linearization will only be stopped at step $k = 8$. On the other hand, with our adaptive stopping criterion (6.5), only one iteration is necessary. Clearly after a few linearization steps, the linearization estimator no longer affects significantly the smoothing estimator, and we can economize many useless iterations.

Next, we provide in Table 1 the results obtained using the adaptive stopping criterion (6.5) to stop the nonlinear solver. We terminate the smoothing iterations using the relative norm of the total residual vector (2.6)

$$ \|\mathbf{R}(\mathbf{X}^{j,\bar{k}})\|_{\text{r}} < 10^{-8}. \qquad (6.6) $$

We present the cumulated number of Newton iterations Niter, the estimators (6.3), and the relative norm of the total residual vector (2.6) at each smoothing step $j$. In terms of numbers, 10 smoothing iterations and 36 cumulated Newton iterations (CPU time: 6.9s) are needed to achieve the stopping criterion (6.6). From Table 1, one can see that for each value of $\mu^j$, the Newton iterations are stopped according to (6.5). $\|\mathbf{R}(\mathbf{X}^{j,\bar{k}})\|_{\text{r}}$ decreases until lying below $10^{-8}$. Figure 4 displays the curve of the estimators as a function of cumulated Newton iterations and smoothing iterations, as well as the relative norm of the total residual vector as a function of smoothing iterations. The improvement of the performance with respect to the semismooth Newton-F-B method of Section 6.3 is spectacular.
<figure><img src="image_4.png" /><figcaption>Figure 3: [Adaptive smoothing Newton method, smoothed F-B function (3.2), classical and adaptive stopping criteria (6.4) and (6.5)] Relative norm of the total residual vector (2.6) and estimators (6.3) as a function of Newton iterations $k$, at a specific smoothing iteration $j = 1$ ($\mu^1 = 1$), left, and at $j = 3$ ($\mu^3 = 10^{-2}$), right.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: [Adaptive smoothing Newton method, smoothed F-B function (3.2), adaptive stopping criterion (6.5)] Estimators (6.3) as a function of cumulated Newton iterations (left). Estimators (6.3) (middle) and relative norm of the total residual vector (2.6) (right) as a function of smoothing iterations $j$ at convergence of the linearization solver.</figcaption></figure>

<table><caption>Table 1: [Adaptive smoothing Newton method, smoothed F-B function (3.2), adaptive stopping criterion (6.5)] Number of Newton iterations Niter, estimators (6.3), and relative norm of the total residual vector (2.6) at each smoothing iteration $j$, at convergence of the linearization solver.</caption><thead><tr><th>μ<sup>j</sup></th><th>Niter</th><th>η<sup>j,bar</sup><sub>lin,ASN</sub></th><th>η<sup>j,bar</sup><sub>sm,ASN</sub></th><th>||R(X<sup>j,bar</sup>)||<sub>r</sub></th></tr></thead><tbody><tr><td>1e+00</td><td>1</td><td>2.17e+03</td><td>4.24e+03</td><td>2.17e+03</td></tr><tr><td>1e-01</td><td>3</td><td>6.00e+01</td><td>2.37e+02</td><td>2.03e+02</td></tr><tr><td>1e-02</td><td>4</td><td>9.73e+00</td><td>1.53e+01</td><td>1.01e+01</td></tr><tr><td>1e-03</td><td>5</td><td>3.18e-01</td><td>6.84e-01</td><td>6.00e-01</td></tr><tr><td>1e-04</td><td>7</td><td>9.87e-03</td><td>3.58e-02</td><td>3.43e-02</td></tr><tr><td>1e-05</td><td>4</td><td>1.06e-03</td><td>2.33e-03</td><td>1.87e-03</td></tr><tr><td>1e-06</td><td>3</td><td>1.14e-04</td><td>1.50e-04</td><td>7.45e-05</td></tr><tr><td>1e-07</td><td>3</td><td>4.85e-06</td><td>8.04e-06</td><td>3.84e-06</td></tr><tr><td>1e-08</td><td>3</td><td>3.23e-07</td><td>4.72e-07</td><td>1.83e-07</td></tr><tr><td>1e-09</td><td>3</td><td>1.43e-08</td><td>2.15e-08</td><td>8.04e-09</td></tr></tbody></table>
With the intention to compare the proposed method to existing methods, we complete the semismooth Newton method by a path-following strategy to solve problem (1.1), following [48]. For the sake of brevity, we shall not detail this here. The following test compares the semismooth Newton method (SSN) and the semismooth Newton method with path-following (SSN-pf) in which the linearization is stopped when the criterion (6.2) is satisfied, to the adaptive smoothing Newton method, using the smoothed min and F-B functions (3.1) and (3.2) and the stopping criteria (6.5) and (6.2) respectively for the linearization and smoothing iterations. We compare the number of cumulated linearization iterations and the global CPU time of the simulation for the different strategies. The results are displayed in Figure 5. They confirm the expected reduction of the computational cost of the numerical resolution with our adaptive approaches. Actually, we notice that the semismooth Newton method with path-following (red curve) and the adaptive smoothing Newton method (purple and dark blue curves) require significantly fewer cumulated Newton iterations and time to converge, in comparison with the semismooth Newton method (green and orange curves). Therefore, employing the path-following strategy or the adaptive strategy based on a posteriori error estimates enables to save many unnecessary additional iterations, and yield much better results than the pure semismooth Newton method. We note that, using the adaptive smoothing Newton method, one obtains similar computational results using both the smoothed F-B or the smoothed min function.

<figure><img src="image_6.png" /><figcaption>Figure 5: [Semismooth Newton method (with and without a path-following strategy) and adaptive smoothing method] Cumulated number of Newton iterations (left) and CPU time (right) as a function of the number of mesh elements.</figcaption></figure>

## 6.5 Adaptive inexact smoothing Newton method

We focus in this section on the adaptive inexact Newton method introduced in Section 3 and investigate the performance of Algorithm 1 using the smoothed F-B function (3.2) together with the restarted GMRES method. Typically, we use a fixed restart parameter equal to 300. The behavior of the adaptive smoothing solvers can be improved dramatically by using good preconditioners. Here, we merely use an ILU preconditioner to speed-up the GMRES solver. For other possibilities for preconditioners, we refer to, e.g., [31] and the references therein. To point out the efficiency of the adaptivity, we test two approaches. First, we stop the algebraic iterations using the standard GMRES stopping criterion on the relative residual given by

$$ R_{\text{rel}} := \frac{\|\mathbb{M}_2 \setminus (\mathbb{M}_1 \setminus (B_{\mu j}^{j,k-1} - A_{\mu j}^{j,k-1} X^{j,k,i}))\|}{\|\mathbb{M}_2 \setminus (\mathbb{M}_1 \setminus (B_{\mu j}^{j,k-1} - A_{\mu j}^{j,k-1} X^{j,k-1}))\|} \le 10^{-10}, \quad (6.7) $$

where $\mathbb{M}_1$ and $\mathbb{M}_2$ are the preconditioner matrices. Second, we incorporate the adaptive stopping criteria (3.11a) for the algebraic solver in Algorithm 1. We set the parameters $\mu^1 = 1, \varepsilon = 10^{-5}, \alpha_{\text{alg}} = 10^{-3}, \alpha_{\text{lin}} = 1$, and $\alpha = 0.1$. Figure 6 depicts the evolution of the algebraic and linearization estimators and the GMRES relative residual during the algebraic resolution, for specific smoothing step $j$ and linearization step $k$. For $j = 2$ and $k = 2$, we see that 22 GMRES iterations are needed to achieve the standard stopping criterion (6.7), whereas in the adaptive resolution case, only 10 GMRES iterations are required to satisfy the adaptive stopping criterion (3.11a). In this case, we can avoid many unnecessary iterations. One can also see from the right part of Figure 6, for $j = 3$ and $k = 1$, that the overall gain in terms of algebraic iterations obtained using our stopping criteria is quite significant.

Figure 7, left, shows the evolution of the estimators during smoothing iterations, at convergence of the nonlinear and linear solvers. As expected, the estimators decrease when $\mu$ decreases at each smoothing step. In the middle
<figure><img src="image_7.png" /><figcaption>Figure 6: [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Algebraic and linearization estimators (3.10) and GMRES relative residual as a function of the GMRES iterations $i$, for a fixed smoothing and linearization iterations, $j = 2, k = 2, i$ varies, left, and $j = 3, k = 1, i$ varies, right, using the classical stopping criterion (6.7) and the adaptive one (3.11a).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Estimators (3.10) as a function of smoothing iterations $j$ at convergence of the algebraic and linearization solvers, left. Estimators as a function of cumulated Newton iterations at convergence of the algebraic solver, middle. Estimators as a function of cumulated GMRES iterations during the first two smoothing iterations ($j = 1$ and $j = 2$), right.</figcaption></figure>

405 part of Figure 7, we can observe the behavior of the estimators at the end of the algebraic iterations, during the 406 linearization iterations. We present 8 curves, each one corresponding to a specific value of $\mu^j$. We can see that at 407 each smoothing iteration $j$, the smoothing estimator $\eta_{\text{sm,AISN}}^{j,k,\tilde{i}}$ stagnates after about two iterations. The linearization 408 estimator $\eta_{\text{lin,AISN}}^{j,k,\tilde{i}}$ decreases until becoming smaller than the smoothing estimator, satisfying the stopping criterion 409 (3.11b). Finally, the detected behavior in terms of all smoothing iterations $j$, linearization iterations $k$, and algebraic 410 solver iterations $i$ is presented in Figure 7, right, for $j \le 2$. The overall results are collected in Table 2. We present 411 in particular the number of linearization and cumulated algebraic iterations per smoothing step $j$, Niter and Giter 412 respectively, as well as the estimators (3.10) and the relative norm of the total residual vector (2.6) at the end of 413 each smoothing step $j$. Using the adaptive stopping criteria (3.11), 8 smoothing iterations, 39 cumulated Newton 414 iterations, and 5999 cumulated GMRES iterations are needed to ensure convergence. Figure 8 illustrates the 415 performance of the adaptive inexact smoothing Newton method. It represents the ratio between: 1) the number 416 of algebraic iterations (left) and the CPU time (right) using the classical GMRES stopping criterion (6.7) and 2) 417 the number of algebraic iterations and the CPU time using the adaptive stopping criterion (3.11a) for GMRES, as 418 a function of the number of elements. For larger systems, 20-times fewer iterations and 18-times faster execution 419 time are achieved.
<table><caption>Table 2: [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Number of Newton iterations and cumulated GMRES iterations, estimators (3.10), and relative norm of the total residual vector (2.6) at each smoothing iteration j, at convergence of the algebraic and linearization solvers.</caption><thead><tr><th>μ<sup>j</sup></th><th>Niter</th><th>Giter</th><th>η<sup>j,k,̄i</sup><sub>lin,AISN</sub></th><th>η<sup>j,k,̄i</sup><sub>sm,AISN</sub></th><th>η<sup>j,k,̄i</sup><sub>alg,AISN</sub></th><th>||R(X<sup>j,k,̄i</sup>)||<sub>r</sub></th></tr></thead><tbody><tr><td>1e+00</td><td>1</td><td>8</td><td>2.16e+03</td><td>4.24e+03</td><td>1.80e+00</td><td>2.19e+03</td></tr><tr><td>1e-01</td><td>4</td><td>34</td><td>5.95e+01</td><td>2.31e+02</td><td>1.89e-02</td><td>1.80e+02</td></tr><tr><td>1e-02</td><td>3</td><td>391</td><td>1.54e+01</td><td>1.73e+01</td><td>1.41e-02</td><td>6.75e+00</td></tr><tr><td>1e-03</td><td>4</td><td>198</td><td>5.04e-01</td><td>8.16e-01</td><td>4.60e-04</td><td>5.95e-01</td></tr><tr><td>1e-04</td><td>10</td><td>796</td><td>7.99e-03</td><td>3.53e-02</td><td>5.58e-06</td><td>3.43e-02</td></tr><tr><td>1e-05</td><td>10</td><td>684</td><td>8.54e-04</td><td>2.12e-03</td><td>7.61e-07</td><td>1.94e-03</td></tr><tr><td>1e-06</td><td>4</td><td>513</td><td>9.03e-05</td><td>1.48e-04</td><td>7.42e-08</td><td>1.05e-04</td></tr><tr><td>1e-07</td><td>3</td><td>3375</td><td>6.04e-06</td><td>8.14e-06</td><td>4.27e-09</td><td>4.26e-06</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>Figure 8: [Adaptive inexact smoothing Newton method, smoothed F-B function (3.2), Algorithm 1] Ratio between: the number of algebraic iterations (left) and CPU time (right) needed by the classical stopping criterion (6.7) to converge to the number and time needed by the adaptive stopping criterion (3.11a), as a function of the number of mesh elements.</figcaption></figure>

## 6.6 Nonparametric interior-point method

We consider here the nonparametric interior-point approach of Section 4, where the dimension of the corresponding problem is $n = 3N + 1$. The value of the constant $\theta$ in the additional equation (4.2) is $10^{-1}$. The stopping criterion is on the relative norm of the linearization residual vector

$$ \|\mathbf{R}^{\text{IP}}(\mathbf{X}^k)\|_r := \|\mathbf{R}^{\text{IP}}(\mathbf{X}^k)\|/\|\mathbf{R}^{\text{IP}}(\mathbf{X}^0)\| < 10^{-8}, \quad (6.8) $$

with

$$ \mathbf{R}^{\text{IP}}(\mathbf{X}^k) := \left[ \begin{array}{c} \mathbf{F} - \mathbb{E}\mathbf{X}^k \\ \mu - \mathbf{K}(\mathbf{X}^k)\mathbf{G}(\mathbf{X}^k) \\ -\theta\mu^k - (\mu^k)^2 \end{array} \right]. $$

Using this method, 19 Newton iterations (CPU time: 6.7s) are needed to reach the end of the simulation. Figure 9 shows that $||\mathbf{R}^{\text{IP}}(\mathbf{X}^k)||_r$ decreases during the Newton interior-point iterations until satisfying the stopping criterion (6.8).

## 6.7 Adaptive interior-point method

Next, we consider the adaptive interior-point method, which is the method presented in Section 5, Algorithm 2 without applying an algebraic iterative solver to approximate the solution of the linear system (5.2). In this case, we can define the linearization and smoothing estimators respectively by

$$ \eta_{\text{lin,AIP}}^{j,k} := ||\mathbf{H}_{\mu^j}(\mathbf{X}^{j,k})||_r, \quad (6.9a) $$

$$ \eta_{\text{sm,AIP}}^{j,k} := ||\mu^j||_r, \quad (6.9b) $$
<figure><img src="image_10.png" /><figcaption>Figure 9: [Nonparametric interior-point method, stopping criterion (6.8)] Relative norm of the linearization residual vector (6.8) as a function of Newton iterations.</figcaption></figure>

where $H_{\mu^j}(\cdot)$ is defined in (5.1b), and the total estimator by $\eta_{\text{AIP}}^{j,k} := \eta_{\text{sm,AIP}}^{j,k} + \eta_{\text{lin,AIP}}^{j,k}$. Recall from (5.8) the definition of the total residual vector for $V \in \mathbb{R}^n$ as

$$ \mathbf{R}^{\text{AIP}}(\mathbf{V}) := \begin{bmatrix} \mathbf{F} - \mathbf{E}\mathbf{V} \\ -\mathbf{H}(\mathbf{V}) \end{bmatrix}, \qquad (6.10) $$

where $\mathbf{H}(\cdot)$ is defined in (5.7). The adaptive stopping criterion

$$ \eta_{\text{lin,AIP}}^{j,k} < \alpha_{\text{lin}} \eta_{\text{sm,AIP}}^{j,k} \qquad (6.11) $$

is used to stop the nonlinear solver and a criterion on the relative norm of the total residual vector is applied to stop the smoothing iterations

$$ \| \mathbf{R}^{\text{AIP}}(\mathbf{X}^{j, \bar{k}}) \|_{\mathrm{r}} < 10^{-8}. \qquad (6.12) $$

The initial smoothing vector is $\mu^1 = [1, \dots, 1]^T \in \mathbb{R}^N$ and $\alpha_{\text{lin}} = 1$. Concerning the update of the smoothing parameter $\mu$, we set $\alpha = 10^{-1}$. Table 3 summarizes the results. To achieve the stopping criterion (6.12), 11 smoothing iterations and 20 cumulated Newton iterations are needed (CPU time: 5.0s). In Figure 10, we plot the estimators (6.9) as a function of the cumulated Newton iterations (left), the smoothing iterations (middle), and the relative norm of the residual vector as a function of the smoothing iterations (right). The behavior of $\|\mathbf{R}^{\text{AIP}}(\mathbf{X}^{j, \bar{k}})\|_{\mathrm{r}}$ in Figure 10 appears a bit different from its behavior in Figure 4. This is related to the fact that the relative norm of the total residual given by (2.6) includes $\mathcal{C}(\mathbf{X})$ in the adaptive smoothing Newton method, whereas in this adaptive interior-point method, the relative norm of the total residual given by (6.10) includes $\mathbf{K}(\mathbf{X})\mathbf{G}(\mathbf{X})$.

<table><caption>Table 3: [Adaptive interior-point method] Number of Newton iterations, estimators (6.9), and relative norm of the total residual vector (6.10) at each smoothing step j, at convergence of the linearization solver.</caption><thead><tr><th>μ<sup>j</sup></th><th>Niter</th><th>η<sub>lin,AIP</sub><sup>j,bar</sup></th><th>η<sub>sm,AIP</sub><sup>j,bar</sup></th><th>||<b>R</b><sup>AIP</sup>(<b>X</b><sup>j,bar</sup>)||<sub>r</sub></th></tr></thead><tbody><tr><td>1e+00</td><td>2</td><td>1.11e+01</td><td>2.00e+01</td><td>3.00e+01</td></tr><tr><td>1e-01</td><td>2</td><td>1.24e+00</td><td>2.00e+00</td><td>3.20e+00</td></tr><tr><td>1e-02</td><td>2</td><td>1.15e-01</td><td>2.00e-01</td><td>3.11e-01</td></tr><tr><td>1e-03</td><td>2</td><td>6.51e-03</td><td>2.00e-02</td><td>2.43e-02</td></tr><tr><td>1e-04</td><td>2</td><td>3.38e-04</td><td>2.00e-03</td><td>2.14e-03</td></tr><tr><td>1e-05</td><td>1</td><td>1.58e-04</td><td>2.00e-04</td><td>2.82e-04</td></tr><tr><td>1e-06</td><td>2</td><td>3.67e-06</td><td>2.00e-05</td><td>2.10e-05</td></tr><tr><td>1e-07</td><td>2</td><td>1.00e-07</td><td>2.00e-06</td><td>2.02e-06</td></tr><tr><td>1e-08</td><td>1</td><td>1.86e-07</td><td>2.00e-07</td><td>3.84e-07</td></tr><tr><td>1e-09</td><td>2</td><td>9.33e-10</td><td>2.00e-08</td><td>2.01e-08</td></tr><tr><td>1e-10</td><td>2</td><td>2.55e-11</td><td>2.00e-09</td><td>2.00e-09</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>Figure 10: [Adaptive interior-point method] Estimators (6.9) as a function of cumulated Newton iterations (left). Estimators (6.9) (middle) and relative norm of the total residual vector (6.10) (right) as a function of smoothing iterations *j* at convergence of the linearization solver.</figcaption></figure>

## 6.8 Adaptive inexact interior-point method

Let us now present the numerical results of the adaptive inexact interior-point method, detailed in Section 5. We employ Algorithm 2 with the GMRES algebraic solver and an ILU preconditioner. The parameters in Algorithm 2 are set as $\mu^1 = [1, \dots, 1]^T \in \mathbb{R}^N$, $\varepsilon = 10^{-5}$, $\alpha_{alg} = 1$, $\alpha_{lin} = 1$, and $\alpha = 0.1$. The restart parameter of restarted GMRES is chosen equal to 300. From Table 4, we can see that the method converged after 8 smoothing iterations, 20 cumulated linearization iterations, and 760 cumulated GMRES iterations. Figure 11, left, displays the curves of the estimators (5.9) as a function of the smoothing iteration. One can see that the estimators satisfy the adaptive stopping criteria incorporated in Algorithm 2. In Figure 11, right, the estimators are shown as a function of cumulated Newton iterations, at convergence of the linear solver.

<table><caption>Table 4: [Adaptive inexact interior-point method, Algorithm 2] Number of cumulated Newton and GMRES iterations, estimators (5.9), and relative norm of the total residual vector (5.8) at each smoothing iteration *j*, at convergence of the algebraic and linearization solvers.</caption><thead><tr><th>μ<sup>j</sup></th><th>Niter</th><th>Giter</th><th>η<sup>j, ⃗k, ⃗i</sup><sub>lin, AIIP</sub></th><th>η<sup>j, ⃗k, ⃗i</sup><sub>sm, AIIP</sub></th><th>η<sup>j, ⃗k, ⃗i</sup><sub>alg, AIIP</sub></th><th>||R<sup>AIIP</sup>(X<sup>j, ⃗k, ⃗i</sup>)||<sub>r</sub></th></tr></thead><tbody><tr><td>1e+00</td><td>3</td><td>7</td><td>1.15e+01</td><td>2.00e+01</td><td>3.36e+00</td><td>5.59e+00</td></tr><tr><td>1e-01</td><td>2</td><td>12</td><td>5.44e-01</td><td>2.00e+00</td><td>1.78e-02</td><td>2.00e+00</td></tr><tr><td>1e-02</td><td>3</td><td>20</td><td>9.75e-02</td><td>2.00e-01</td><td>2.80e-02</td><td>2.04e-01</td></tr><tr><td>1e-03</td><td>3</td><td>29</td><td>4.82e-03</td><td>2.00e-02</td><td>1.74e-03</td><td>2.01e-02</td></tr><tr><td>1e-04</td><td>3</td><td>56</td><td>2.52e-04</td><td>2.00e-03</td><td>2.19e-04</td><td>2.01e-03</td></tr><tr><td>1e-05</td><td>2</td><td>62</td><td>1.77e-04</td><td>2.00e-04</td><td>1.08e-04</td><td>2.29e-04</td></tr><tr><td>1e-06</td><td>2</td><td>110</td><td>1.49e-05</td><td>2.00e-05</td><td>1.42e-05</td><td>2.46e-05</td></tr><tr><td>1e-07</td><td>2</td><td>464</td><td>1.34e-06</td><td>2.00e-06</td><td>1.16e-06</td><td>2.31e-06</td></tr></tbody></table>

## 6.9 Comparison of the methods

This section is devoted to compare the semismooth Newton method (SSN), semismooth Newton method with path- following (SSN-pf), nonparametric interior-point method (IP), adaptive smoothing Newton method (ASN), and adaptive interior-point method (AIP). For this purpose, we introduce a unified residual vector, for $V \in \mathbb{R}^n$

$$R_{\text{unif}}(V) := \begin{bmatrix} F - EV \\ \min(0, K(V)) \\ \min(0, G(V)) \\ K(V) \cdot G(V) \end{bmatrix}, \quad (6.13)$$
<figure><img src="image_12.png" /><figcaption>Figure 11: [Adaptive inexact interior-point method, Algorithm 2] Estimators (5.9) as a function of smoothing iterations $j$ at convergence of the algebraic and linearization solvers (left). Estimators as a function of cumulated Newton iterations $k$ at convergence of the algebraic solver (right).</figcaption></figure>

466

467

468

469

470

independent of the way the nonlinear complementarity constraints are reformulated. The stopping criterion of the nonlinear solver for the classical methods (SSN, SSN-pf, IP) is on the relative unified residual $\|\mathbf{R}_{\text{unif}}(\mathbf{X}^k)\|_r$ lying below $10^{-8}$. Regarding the adaptive methods (ASN, AIP), to stop the nonlinear solver, we use the adaptive stopping criteria given respectively in (6.5) and (6.11). To stop the smoothing iterations, $\|\mathbf{R}_{\text{unif}}(\mathbf{X}^{j,k})\|_r$ is requested to become smaller than $10^{-8}$.

<figure><img src="image_13.png" /><figcaption>Figure 12: [Semismooth Newton method (F-B function (2.2)), semismooth Newton method with a path-following strategy, nonparametric interior-point method, adaptive interior-point method, and adaptive smoothing Newton method (smoothed F-B function (3.2))] Number of cumulated Newton iterations (left) and CPU time (right) as a function of the number of mesh elements, employing a stopping criterion on the relative norm of the unified residual vector (6.13).</figcaption></figure>

471

472

473

474

475

476

477

478

In Figure 12, we plot the cumulated number of the Newton iterations (left) and the CPU time (right) required by each method, as a function of the number of mesh elements. It is clearly seen that the semismooth Newton method (green curve) is typically more costly, both in terms of the required number of iterations and the CPU time, in comparison with the other methods. Precisely, we can observe an important gain between the semismooth Newton method (green curve) and the adaptive smoothing Newton method (purple curve). Moreover, as we can remark from the red curve, the combination of a path-following strategy to the semismooth Newton method seems to be efficient. Finally, one does not see a remarkable difference between the results of the nonparametric interior-point method (cyan curve) and the adaptive interior-point method (black curve) in this test case.
# 7 Numerical experiments: Two-phase flow with phase transition

The second model problem that we consider in our numerical tests is a two-phase flow model (liquid-gas) with phase transition in porous media following [6, 13, 9]. Each of the liquid phase, denoted by l, and the gas phase, denoted by g, is composed of two components, water and hydrogen, denoted respectively by w and h.

## 7.1 Problem statement

The problem at hand can be formulated as a system of nonlinear partial differential equations with nonlinear complementarity constraints at each time step $\tau_v$. Let $\mathcal{T}_h$ be the spatial mesh, we denote respectively by $S_K^\nu, P_K^\nu$, and $\chi_K^\nu$ the discrete elementwise unknowns approximating the values of the saturation $S^1$, the pressure $P^1$, and the molar fraction of hydrogen in the liquid phase $\chi_h^1$ in the element $K \in \mathcal{T}_h$ and on time step $1 \le \nu \le N_t$. Let $N$ be the number of elements in the mesh $\mathcal{T}_h$. If one introduces the appropriate nonlinear function $H_{c,K}^\nu : \mathbb{R}^{3N} \to \mathbb{R}$, $c \in \{\text{w}, \text{h}\}$, and suitable functions $F_K : \mathbb{R}^3 \to \mathbb{R}$ and $G_K : \mathbb{R}^3 \to \mathbb{R}$, the discrete problem written elementwise consists in finding $\mathbf{X}^\nu := (\mathbf{X}_K^\nu)_{K \in \mathcal{T}_h} \in \mathbb{R}^n$, where $n = 3N$, and $\mathbf{X}_K^\nu := [S_K^\nu, P_K^\nu, \chi_K^\nu] \in \mathbb{R}^3$, such that for all $K \in \mathcal{T}_h$

$$H_{c,K}^{\nu}(\mathbf{X}^{\nu}) = 0, \quad c \in \{\text{w}, \text{h}\}, \qquad (7.1a)$$

$$F_K(\mathbf{X}_K^\nu) \ge 0, \quad G_K(\mathbf{X}_K^\nu) \ge 0, \quad F_K(\mathbf{X}_K^\nu)G_K(\mathbf{X}_K^\nu) = 0. \qquad (7.1b)$$

The formulation (7.1) allows to model the transition from a single-phase flow to a two-phase flow during the appearance and disappearance of the gas phase and vice versa. As an example, a detailed finite volume discretization can be found in [7, Section 3.2]. The first $2N$ lines of system (7.1) can be written globally as

$$\mathcal{H}^{\nu}(\mathbf{X}^{\nu}) = 0,$$

where $\mathcal{H}^n : \mathbb{R}^{3N} \to \mathbb{R}^{2N}$ is defined elementwise by (7.1a).

Considering a C-function $C^\nu$, for $1 \le \nu \le N_t$, we define a function $C^\nu : \mathbb{R}^{3N} \to \mathbb{R}^N$ as $C^\nu(\mathbf{X}^\nu) = C^\nu((F_K(\mathbf{X}_K^\nu))_{K \in \mathcal{T}_h}, (G_K(\mathbf{X}_K^\nu))_{K \in \mathcal{T}_h})$. This leads us to apply a semismooth Newton method to find a solution for problem (7.1) written as

$$\begin{aligned} \mathcal{H}^{\nu}(\mathbf{X}^{\nu}) &= 0, \\ C^{\nu}(\mathbf{X}^{\nu}) &= 0. \end{aligned} \qquad (7.2)$$

The total residual vector $\mathbf{R}(\mathbf{V})$ of problem (7.2) is thus given by

$$\mathbf{R}(\mathbf{V}) := \begin{bmatrix} -\mathcal{H}^{\nu}(\mathbf{V}) \\ -\mathcal{C}^{\nu}(\mathbf{V}) \end{bmatrix}, \quad \forall \mathbf{V} \in \mathbb{R}^n. \qquad (7.3)$$

## 7.2 Adaptive smoothing Newton method

We introduce a function $C_\mu^\nu : \mathbb{R}^{3N} \to \mathbb{R}^N$ defined as $C_\mu^\nu(\mathbf{X}^\nu) = C_\mu^\nu((F_K(\mathbf{X}_K^\nu))_{K \in \mathcal{T}_h}, (G_K(\mathbf{X}_K^\nu))_{K \in \mathcal{T}_h})$, for $1 \le \nu \le N_t$, where $C_\mu^\nu$ is a smoothed C-function. Line (7.1b) can be approximated as a smoothed nonlinear equation $C_\mu^\nu(\mathbf{X}^\nu) = 0$, making it possible to apply the standard Newton method to solve the resulting nonlinear system in the form: Find $\mathbf{X}^{\nu,j} \in \mathbb{R}^{3N}$ at each time step $\nu$, $1 \le \nu \le N_t$, satisfying

$$\begin{aligned} \mathcal{H}^{\nu}(\mathbf{X}^{\nu,j}) &= 0, \\ C_{\mu^{j\nu}}^{\nu}(\mathbf{X}^{\nu,j}) &= 0. \end{aligned} \qquad (7.4)$$

At each time step $1 \le \nu \le N_t$, each smoothing step $j \ge 1$, and each linearization step $k \ge 1$, fixing $\mathbf{X}^{\nu,j,0} \in \mathbb{R}^n$, we try to approach the solution of problem (7.4) by finding a solution $\mathbf{X}^{\nu,j,k} \in \mathbb{R}^n$ such that

$$A_{\mu^{j\nu}}^{\nu,j,k-1} \mathbf{X}^{\nu,j,k} = B_{\mu^{j\nu}}^{\nu,j,k-1}, \qquad (7.5)$$

where the Jacobian matrix $A_{\mu^{j\nu}}^{\nu,j,k-1} \in \mathbb{R}^{n,n}$ and the right-hand side vector $B_{\mu^{j\nu}}^{\nu,j,k-1} \in \mathbb{R}^n$ are defined by

$$A_{\mu^{j\nu}}^{\nu,j,k-1} := \begin{bmatrix} \mathbf{J}_{\mathcal{H}^{\nu}}(\mathbf{X}^{\nu,j,k-1}) \\ \mathbf{J}_{C_{\mu^{j\nu}}^{\nu}}(\mathbf{X}^{\nu,j,k-1}) \end{bmatrix}, \qquad (7.6a)$$

$$\mathbf{B}_{\mu^{j\nu}}^{\nu,j,k-1} := \begin{bmatrix} \mathbf{J}_{\mathcal{H}^{\nu}}(\mathbf{X}^{\nu,j,k-1})\mathbf{X}^{\nu,j,k-1} - \mathcal{H}^{\nu}(\mathbf{X}^{\nu,j,k-1}) \\ \mathbf{J}_{C_{\mu^{j\nu}}^{\nu}}(\mathbf{X}^{\nu,j,k-1})\mathbf{X}^{\nu,j,k-1} - C_{\mu^{j\nu}}^{\nu}(\mathbf{X}^{\nu,j,k-1}) \end{bmatrix}, \qquad (7.6b)$$

with $\mathbf{J}_{\mathcal{H}^{\nu}}(\mathbf{X}^{\nu,j,k-1})$ and $\mathbf{J}_{C_{\mu^{j\nu}}^{\nu}}(\mathbf{X}^{\nu,j,k-1})$ the Jacobian matrices of the function $\mathcal{H}^{\nu}$ and the smoothed function $C_{\mu^{j\nu}}^{\nu}$, respectively, at the point $\mathbf{X}^{\nu,j,k-1}$ obtained by a Newton linearization.
## 7.3 Adaptive smoothing Newton algorithm

522

523

524

**Algorithm 3** Unsteady adaptive smoothing Newton algorithm

**Initialization:** Fix $\epsilon > 0$, $\alpha \in ]0, 1[$, and $\alpha_{\text{lin}} \in ]0, 1]$. Set $\nu := 1$ and $t_{\nu} := 0$. Choose $\mathbf{X}^{\nu,0} \in \mathbb{R}^n$.

**Time loop**

1. Fix $\mu^{j\nu} > 0$ and set $j := 1$.

2. Smoothing loop

2.1 Set $\mathbf{X}^{\nu,j,0} := \mathbf{X}^{\nu,0}$ and $k := 1$.

2.2 Newton linearization loop

2.2.1 From $\mathbf{X}^{\nu,j,k-1}$ define $A_{\mu^{j\nu}}^{\nu,j,k-1} \in \mathbb{R}^{n,n}$ and $B_{\mu^{j\nu}}^{\nu,j,k-1} \in \mathbb{R}^n$ given by (7.6).

2.2.2 Find $\mathbf{X}^{\nu,j,k}$ solution to the linear system

$$A_{\mu^{j\nu}}^{\nu,j,k-1} \mathbf{X}^{\nu,j,k} = B_{\mu^{j\nu}}^{\nu,j,k-1}.$$

2.2.3 Compute the estimators and check the stopping criterion for the nonlinear solver

$$\left( \eta_{\text{lin,ASN}}^{\nu,j,k} < \alpha_{\text{lin}} \eta_{\text{sm,ASN}}^{\nu,j,k} \right) \quad \text{or} \quad \left( \eta_{\text{lin,ASN}}^{\nu,j,k} < \varepsilon \right). \qquad (7.7)$$

If satisfied, set $\bar{k} := k$ and stop. If not, set $k := k + 1$ and go to 2.2.1.

2.3 Check the stopping criterion for the smoothing iterations in the form:

$$\max \left\{ \eta_{\text{sm,ASN},}^{\nu,j,\bar{k}}, \left\| \mathbf{R}(\mathbf{X}^{\nu,j,\bar{k}}) \right\|_{\text{r}} \right\} < \varepsilon. \qquad (7.8)$$

If satisfied, set $\bar{j} := j$ and stop. If not, set $j := j + 1$, $\mathbf{X}^{\nu,j,0} := \mathbf{X}^{\nu,j-1,\bar{k}}$, and $\mu^{j\nu} := \alpha\mu^{(j-1)\nu}$. Then set $k := 1$ and go to 2.2.1.

If $\nu = N_t$, stop. If not, set $\nu := \nu + 1$, $j = 1$, $\mathbf{X}^{\nu,j,0} := \mathbf{X}^{\nu-1,\bar{j}}$, and $t_{\nu} := \tau_{\nu} + t_{\nu-1}$. Then set $\mu^{j\nu} = \mu^{\bar{j}\nu-1}$, $k = 1$, and go to 2.2.1.

525

526

527

528

529

530

**Description of Algorithm 3** For the first time step $\nu = 1$, starting with an initial approximation $\mathbf{X}^{\nu,0} \in \mathbb{R}^n$ and an initial smoothing parameter $\mu^{\nu,1} > 0$, we solve the smoothed nonlinear system (7.5) by the Newton linearization solver, and decrease the smoothing parameter $\mu^{\nu,j}$ at each smoothing step $j$, until the stopping criterion (7.8) on the smoothing estimator or the relative norm of the total residual vector is satisfied at step $\bar{j}$. Then, we continue the time loop, for $2 \le \nu \le N_t$, starting for $j = 1$ with $\mathbf{X}^{\nu,j,0} := \mathbf{X}^{\nu-1,\bar{j}}$ and $\mu^{j\nu} := \mu^{\bar{j}\nu-1}$, until satisfying the stopping criterion (7.8).

531

## 7.4 Numerical results

532

533

534

535

536

We consider a homogeneous porous medium in one dimension, supposed to be horizontal with length 2m, and a uniform spatial mesh with $N = 1000$ elements. The final time of simulation is $t_F = 100s$, and the time step is constant $\tau_{\nu} = 10s$. We assume that the medium is initially saturated with liquid, $S^l = 1$, and containing no hydrogen, $\chi_h^l = 0$, on which we impose an injection of gas (hydrogen), constant in time, in the first cell of the mesh. The initial conditions are $S^{l,\nu=0} = 1$, $P^{l,\nu=0} = 10^6 \text{Pa}$, and $\chi_h^{l,\nu=0} = 0$.

537

538

539

**Semismooth Newton method.** We begin by employing the semismooth Newton method presented in Section 2, with the min function (2.1) to solve the nonlinear system (7.2). On each time step $\nu \ge 1$, we request the relative norm of the total residual vector $\mathbf{R}(\mathbf{X}^{\nu,k})$ given by (7.3) to drop below $10^{-4}$.
In Figure 13, the evolution of $||R(X^{\nu,k})||_r$ is shown at each time step. 31 cumulated Newton iterations are needed.

<figure><img src="image_14.png" /><figcaption>Figure 13: [Semismooth Newton method, min function (2.1)] Relative norm of the total residual vector (7.3) as a function of cumulated Newton iterations along the time steps $\nu$.</figcaption></figure>

**Adaptive smoothing Newton method.** Next, we present the results obtained using the adaptive smoothing Newton method, summarized in Algorithm 3, with the smoothed min function (3.1) to solve the smoothed nonlinear problem (7.4) at each time step $\tau_\nu$, $1 \le \nu \le N_t$. The parameters are set as $\mu^{j1} = 10^{-1}$, $\varepsilon = 10^{-4}$, $\alpha_{\text{lin}} = 1$, and $\alpha = 0.1$.

<figure><img src="image_15.png" /><figcaption>Figure 14: [Adaptive smoothing Newton method, smoothed min function (3.1), Algorithm 3] Estimators (6.3) and relative norm of the total residual vector (7.3) at the first time step $\nu = 1$ as a function of smoothing iterations $j$, at convergence of the linearization solver ($\nu = 1$ fixed, $j$ varies, $k = \bar{k}$), left, and of cumulated Newton iterations, right, ($\nu = 1$ fixed, $j$ and $k$ vary).</figcaption></figure>
<table><caption>Table 5: [Adaptive smoothing Newton method, smoothed min function (3.1), Algorithm 3] Relative norm of the total residual vector (7.3) and estimators (6.3) at each time step $ν$, at convergence of the linearization solver.</caption><thead><tr><th>ν</th><th>μ<sup>jν</sup></th><th>Niter</th><th>η<sup>ν,j,bar</sup><sub>lin,ASN</sub></th><th>η<sup>ν,j,bar</sup><sub>sm,ASN</sub></th><th>R(X<sup>ν,j,bar</sup>)<sub>r</sub></th></tr></thead><tbody><tr><td>2</td><td>1e-05</td><td>3</td><td>2.15e-07</td><td>3.13e-07</td><td>2.86e-07</td></tr><tr><td>3</td><td>1e-05</td><td>3</td><td>3.13e-07</td><td>3.68e-07</td><td>4.24e-07</td></tr><tr><td>4</td><td>1e-05</td><td>3</td><td>3.93e-07</td><td>1.19e-07</td><td>3.47e-07</td></tr><tr><td>5</td><td>1e-05</td><td>3</td><td>4.62e-07</td><td>1.59e-07</td><td>4.01e-07</td></tr><tr><td>6</td><td>1e-05</td><td>3</td><td>5.04e-07</td><td>1.88e-06</td><td>1.87e-06</td></tr><tr><td>7</td><td>1e-05</td><td>3</td><td>5.58e-07</td><td>1.74e-07</td><td>4.94e-07</td></tr><tr><td>8</td><td>1e-05</td><td>3</td><td>5.94e-07</td><td>3.76e-07</td><td>7.08e-07</td></tr><tr><td>9</td><td>1e-05</td><td>3</td><td>6.64e-07</td><td>2.77e-07</td><td>7.50e-07</td></tr><tr><td>10</td><td>1e-05</td><td>3</td><td>7.01e-07</td><td>3.00e-07</td><td>7.89e-07</td></tr></tbody></table>

From Figure 14, one can see that at the first time step $ν = 1$ and at each smoothing step $j \le 4$, the linearization estimator decreases until lying below the smoothing estimator. The smoothing iterations are thus stopped in the first possibility according to the stopping criterion (7.7). On the other hand, at the 5th smoothing step, $η_{lin,ASN}^{ν,j,bar}$ is smaller than the fixed tolerance but not smaller than $η_{sm,ASN}^{ν,j,bar}$. Even after additional Newton iterations at this smoothing step, we will have the same observation. This justifies the modification applied in the adaptive stopping criterion (7.7). In Figure 14, right, we report the estimators and $||R(X^{1,j,k})||_r$ as a function of cumulated Newton iteration for $ν = 1$. The stopping criterion (7.8) is satisfied after 5 cumulated smoothing iterations, and 10 cumulated Newton iterations. Then, as presented in Table 5, starting at the second time step ($ν = 2$) with $μ^{jν} = 10^{-5}$, the smoothing parameter does not decrease since the stopping criterion (7.8) is satisfied at each time step after one smoothing step only. To reach the end of the simulation, 9 cumulated smoothing steps and 31 cumulated linearization steps are needed.

As a conclusion, the results confirm the expected behavior of Algorithm 3 featuring an adaptive stopping criterion for the nonlinear solver. In this case, though, the stopping criteria in the adaptive smoothing Newton method do not bring the number of iterations down since the semismooth Newton method already behaves very well here.

## 8 Conclusion and outlook

In this work, we have considered nonlinear algebraic systems with inequalities in a form of complementarity constraints. We have considered some existing methods, like the semismooth Newton method, possibly combined with a path-following strategy, or a nonparametric interior-point method. Our goal was to propose a systematic way to drive such methods with adaptive stopping criteria and possibly inexact algebraic solvers. We have achieved this by a reformulation of the complementarity constraints using a smoothed function and a posteriori error estimate that enables to distinguish the different error components. Numerical experiments confirmed that the proposed adaptive approaches yield significant computational savings compared to some standard approaches from literature. Moreover, their numerical performance seems to be notably good across a range of test problems. In [8], we also take into account the discretization error of the considered problem, enabling to adaptively stop the outer smoothing loop in Algorithm 1, and employ the method to solve more involved problems.

## References

[1] M. AGANAGIĆ, *Newton's method for linear complementarity problems*, Math. Programming, 28 (1984), pp. 349-362, https://doi.org/10.1007/BF02612339.

[2] M. AINSWORTH AND J. T. ODEN, *A posteriori error estimation in finite element analysis*, Pure and Applied Mathematics (New York), Wiley-Interscience [John Wiley & Sons], New York, 2000, https://doi.org/10.1002/9781118032824.

[3] S. BELLAVIA, *Inexact interior-point method*, J. Optim. Th. Appl., 96 (1998), pp. 109-121, https://doi.org/10.1023/A:1022663100715.
579

580

581

[4] S. BELLAVIA, M. MACCONI, AND B. MORINI, *An affine scaling trust-region approach to bound-constrained nonlinear systems*, Appl. Numer. Math., 44 (2003), pp. 257–280, https://doi.org/10.1016/S0168-9274(02)00170-8.

582

583

584

[5] F. BEN BELGACEM, C. BERNARDI, A. BLOUZA, AND M. VOHRALÍK, *On the unilateral contact between membranes. Part 2: a posteriori analysis and numerical experiments*, IMA J. Numer. Anal., 32 (2012), pp. 1147–1172, https://doi.org/10.1093/imanum/drr003.

585

586

587

[6] I. BEN GHARBIA, *Résolution de problèmes de complémentarité. Application à un écoulement diphasique dans un milieu poreux*, PhD thesis, Université Paris Dauphine (Paris IX), 2012, http://tel.archives-ouvertes.fr/tel-00776617.

588

589

590

[7] I. BEN GHARBIA, J. DABAGHI, V. MARTIN, AND M. VOHRALÍK, *A posteriori error estimates for a compositional two-phase flow with nonlinear complementarity constraints*, Comput. Geosci., 24 (2020), pp. 1031–1055, https://doi.org/10.1007/s10596-019-09909-5.

591

592

593

[8] I. BEN GHARBIA, J. FERZLY, M. VOHRALÍK, AND S. YOUSEF, *Adaptive inexact smoothing Newton method for a nonconforming discretization of a variational inequality*. HAL Preprint 03696024, submitted for publication, 2022, https://hal.inria.fr/hal-03696024.

594

595

[9] I. BEN GHARBIA AND E. FLAURAUD, *Study of compositional multiphase flow formulation using complementarity conditions*, Oil Gas Sci. Technol., 74 (2019), p. 43, https://doi.org/10.2516/ogst/2019012.

596

597

598

[10] I. BEN GHARBIA AND J. C. GILBERT, *Nonconvergence of the plain Newton-min algorithm for linear complementarity problems with a P-matrix*, Math. Prog., 134 (2012), pp. 349–364, https://doi.org/10.1007/s10107-010-0439-6.

599

600

[11] I. BEN GHARBIA AND J. C. GILBERT, *An algorithmic characterization of P-matricity*, SIAM J. Matrix Anal. Appl., 34 (2013), pp. 904–916, https://doi.org/10.1137/120883025.

601

602

603

[12] I. BEN GHARBIA AND J. C. GILBERT, *An algorithmic characterization of P-matricity II: adjustments, refinements, and validation*, SIAM J. Matrix Anal. Appl., 40 (2019), pp. 800–813, https://doi.org/10.1137/18M1168522.

604

605

[13] I. BEN GHARBIA AND J. JAFFRÉ, *Gas phase appearance and disappearance as a problem with complementarity constraints*, Math. Comput. Simul., 99 (2014), pp. 28–36, https://doi.org/10.1016/j.matcom.2013.04.021.

606

607

[14] J. F. BONNANS, J. C. GILBERT, C. LEMARÉCHAL, AND C. A. SAGASTIZÁBAL, *Numerical optimization*, Universitext, Springer-Verlag, Berlin, second ed., 2006, https://doi.org/10.1007/978-3-540-35447-5.

608

609

610

[15] M. BÜRG AND A. SCHRÖDER, *A posteriori error control of hp-finite elements for variational inequalities of the first and second kind*, Comput. Math. Appl., 70 (2015), pp. 2783–2802, https://doi.org/10.1016/j.camwa.2015.08.031.

611

612

[16] F. H. CLARKE, *Optimization and Nonsmooth Analysis*, Society for Industrial and Applied Mathematics, 1990, https://doi.org/10.1137/1.9781611971309.

613

614

615

[17] J. DABAGHI, V. MARTIN, AND M. VOHRALÍK, *Adaptive inexact semismooth Newton methods for the contact problem between two membranes*, J. Sci. Comput., 84 (2020), p. Paper No. 28, https://doi.org/10.1007/s10915-020-01264-3.

616

617

618

[18] D. A. DI PIETRO, M. VOHRALÍK, AND S. YOUSEF, *An a posteriori-based, fully adaptive algorithm with adaptive stopping criteria and mesh refinement for thermal multiphase compositional flows in porous media*, Comput. Math. Appl., 68 (2014), pp. 2331–2347, https://doi.org/10.1016/j.camwa.2014.08.008.

619

620

621

[19] J.-P. DUSSAULT, M. FRAPPIER, AND J. C. GILBERT, *A lower bound on the iterative complexity of the Harker and Pang globalization technique of the Newton-min algorithm for solving the linear complementarity problem*, EURO J. Comput. Optim., 7 (2019), pp. 359–380, https://doi.org/10.1007/s13675-019-00116-6.

622

623

624

[20] J.-P. DUSSAULT, M. FRAPPIER, AND J. C. GILBERT, *Polyhedral Newton-min algorithms for complementarity problems*, research report, Inria Paris, France, Université de Sherbrooke, Canada, HAL Preprint 02306526, submitted for publication, 2019, https://hal.archives-ouvertes.fr/hal-02306526.
625

626

[21] A. ERN AND M. VOHRALÍK, *Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs*, SIAM J. Sci. Comput., 35 (2013), pp. A1761–A1791, https://doi.org/10.1137/120896918.

627

628

629

[22] F. FACCHINEI AND C. KANZOW, *A nonsmooth inexact Newton method for the solution of large-scale non-linear complementarity problems*, Math. Program., 76 (1997), pp. 493–512, https://doi.org/10.1016/S0025-5610(96)00058-5.

630

631

632

[23] F. FACCHINEI AND J.-S. PANG, *Finite-dimensional variational inequalities and complementarity problems. Vol. I*, Springer Series in Operations Research, Springer-Verlag, New York, 2003, https://doi.org/10.1007/b97543.

633

634

635

[24] F. FACCHINEI AND J.-S. PANG, *Finite-dimensional variational inequalities and complementarity problems. Vol. II*, Springer Series in Operations Research, Springer-Verlag, New York, 2003, https://doi.org/10.1007/b97544.

636

637

[25] A. GALÁNTAI, *Properties and construction of NCP functions*, Comput. Optim. Appl., 52 (2012), pp. 805–824, https://doi.org/10.1007/s10589-011-9428-9.

638

639

640

[26] Z. GE, Q. NI, AND X. ZHANG, *A smoothing inexact Newton method for variational inequalities with nonlinear constraints*, J. Inequal. Appl., (2017), pp. Paper No. 160, 12., https://doi.org/10.1186/s13660-017-1433-9.

641

642

[27] M. HINTERMÜLLER, K. ITO, AND K. KUNISCH, *The primal-dual active set strategy as a semismooth Newton method*, SIAM J. Optim., 13 (2002), pp. 865–888 (2003), https://doi.org/10.1137/S1052623401383558.

643

644

[28] M. HINTERMÜLLER AND K. KUNISCH, *Path-following methods for a class of constrained minimization problems in function space*, SIAM J. Optim., 17 (2006), pp. 159–187, https://doi.org/10.1137/040611598.

645

646

647

[29] K. ITO AND K. KUNISCH, *Augmented lagrangian methods for nonsmooth, convex optimization in hilbert spaces*, Nonlinear Analysis: Theory, Methods & Applications, 41 (2000), pp. 591–616, https://doi.org/10.1016/S0362-546X(98)00299-5.

648

649

650

[30] K. ITO AND K. KUNISCH, *Lagrange multiplier approach to variational problems and applications*, vol. 15 of Advances in Design and Control, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008, https://doi.org/10.1137/1.9780898718614.

651

652

653

[31] S. LACROIX, Y. VASSILEVSKI, J. WHEELER, AND M. WHEELER, *Iterative solution methods for modeling multiphase flow in porous media fully implicitly*, SIAM J. Sci. Comput., 25 (2003), pp. 905–926, https://doi.org/10.1137/S106482750240443X.

654

655

656

[32] J. M. MARTÍNEZ AND L. Q. QI, *Inexact Newton methods for solving nonsmooth equations*, J. Comput. Appl. Math., 60 (1995), pp. 127–145, https://doi.org/10.1016/0377-0427(94)00088-I. Linear/nonlinear iterative methods and verification of solution (Matsuyama, 1993).

657

658

659

[33] T. S. MUNSON, F. FACCHINEI, M. C. FERRIS, A. FISCHER, AND C. KANZOW, *The semismooth algorithm for large scale complementarity problems*, INFORMS J. Comput., 13 (2001), pp. 294–311, https://doi.org/10.1287/ijoc.13.4.294.9734.

660

661

662

[34] R. H. NOCHETTO, K. G. SIEBERT, AND A. VEESEL, *Theory of adaptive finite element methods: an introduction*, in Multiscale, nonlinear and adaptive approximation, Springer, Berlin, 2009, pp. 409–542, https://doi.org/10.1007/978-3-642-03413-8_12.

663

664

[35] H.-D. QI AND L.-Z. LIAO, *A smoothing Newton method for general nonlinear complementarity problems*, Comput. Optim. Appl., 17 (2000), pp. 231–253, https://doi.org/10.1023/A:1026554432668.

665

666

667

[36] L. QI AND D. SUN, *Smoothing functions and smoothing Newton method for complementarity and variational inequality problems*, J. Optim. Th. Appl., 113 (2002), pp. 121–147, https://doi.org/10.1023/A:1014861331301.

668

669

670

[37] S. REPIN, *A posteriori estimates for partial differential equations*, vol. 4 of Radon Series on Computational and Applied Mathematics, Walter de Gruyter GmbH & Co. KG, Berlin, 2008, https://doi.org/10.1515/9783110203042.
671

672

[38] S. I. REPIN, Functional a posteriori estimates for elliptic variational inequalities, J. Math. Sci., 152 (2008), pp. 702–712, https://doi.org/10.1007/s10958-008-9093-4.

673

674

[39] S.-P. RUI AND C.-X. XU, A smoothing inexact Newton method for nonlinear complementarity problems, J. Comput. Appl. Math., 233 (2010), pp. 2332–2338, https://doi.org/10.1016/j.cam.2009.10.018.

675

676

[40] G. STADLER, Semismooth Newton and augmented Lagrangian methods for a simplified friction problem, SIAM J. Optim., 15 (2004), pp. 39–62, https://doi.org/10.1137/S1052623403420833.

677

678

[41] G. STADLER, Path-following and augmented Lagrangian methods for contact problems in linear elasticity, J. Comput. Appl. Math., 203 (2007), pp. 533–547, https://doi.org/10.1016/j.cam.2006.04.017.

679

680

[42] D. SUN AND L. QI, On NCP-functions, Comput. Optim. Appl., 13 (1999), pp. 201–220, https://doi.org/10.1023/A:1008669226453.

681

682

683

684

[43] M. ULBRICH, Semismooth Newton methods for variational inequalities and constrained optimization problems in function spaces, vol. 11 of MOS-SIAM Series on Optimization, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Optimization Society, Philadelphia, PA, 2011, https://doi.org/10.1137/1.9781611970692.

685

686

687

[44] D. T. S. VU, Numerical resolution of algebraic systems with complementarity conditions. Application to the thermodynamics of compositional multiphase mixtures, PhD thesis, Université Paris-Saclay, 2020, https://tel.archives-ouvertes.fr/tel-02987892.

688

689

690

691

[45] D. T. S. VU, I. BEN GHARBIA, M. HADDOU, AND Q. H. TRAN, A new approach for solving nonlinear algebraic systems with complementarity conditions. Application to compositional multiphase equilibrium problems, Mathematics and Computers in Simulation, 190 (2021), pp. 1243–1274, https://doi.org/10.1016/j.matcom.2021.07.015.

692

693

[46] M. H. WRIGHT, The interior-point revolution in optimization: history, recent developments, and lasting consequences, Bull. Amer. Math. Soc., 42 (2005), pp. 39–56, https://doi.org/10.1090/S0273-0979-04-01040-7.

694

695

696

[47] X. XIAO, Y. LI, Z. WEN, AND L. ZHANG, A regularized semi-smooth Newton method with projection steps for composite convex programs, J. Sci. Comput., (2018), pp. 364–389, https://doi.org/10.1007/s10915-017-0624-3.

697

698

699

[48] S. ZHANG, Y. YAN, AND R. RAN, Path-following and semismooth Newton methods for the variational inequality arising from two membranes problem, J. Inequal. Appl., (2019), pp. Paper No. 1, 13., https://doi.org/10.1186/s13660-019-1955-4.

700

701

[49] X.-Y. ZHAO, D. SUN, AND K.-C. TOH, A Newton-CG augmented Lagrangian method for semidefinite programming, SIAM J. Optim., 20 (2010), pp. 1737–1765, https://doi.org/10.1137/080718206.