Samir Adly

Radek Cibulka

Quantitative Stability of a Generalized Equation

Keywords: Generalized equation, Nonsmooth and variational analysis, Aubin/Lipschitz-like property, Calmness, Outer subdifferential, Non-regular circuits

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

There are many particular applications in sciences and engineering where the model can be formulated in the form of the generalized equation. In this paper, we will focus on the models coming from electrical engineering. The complexity of modern integrated circuits makes the computer-aided circuit analysis inevitable in order to provide information about the behavior, performance and robustness of electrical circuits, by using accurate methods. Many software packages are used in the literature (in both industry and academic institutions) for the simulation of electronic circuits, among them we can cite, for example, SPICE, PLECS, PSIM, QUCS, etc. S. Adly (B) • R. Cibulka XLIM UMR-CNRS 6172, Université de Limoges, 870 60 Limoges, France e-mail: samir.adly@unilim.fr

In this paper, we study a class of non-regular circuits involving electrical devices like ideal diodes, practical diodes, Zener diodes, silicon controlled rectifiers or DIACs (Diodes for Alternating Current). It is well known, in the engineering community, that the study and the simulation of such nonsmooth electrical devices is not an easy task. Classical (smooth) mathematical and numerical analysis is applicable only to few problems. Therefore, natural extensions to the nonsmooth case are needed. Our methodology is based on tools from nonsmooth and convex analysis.

In [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF], the authors considered the special case of the generalized equation ( 1) with the single-valued part being linear and the multivalued part being the Clarke subdifferential of the so-called Moreau-Panagiotopoulos super-potential. Our aim is to investigate the stability properties such as Aubin continuity, calmness and isolated calmness of the corresponding solution mapping. The verifiable conditions ensuring these properties are given in terms of the input data. The stability properties under parameter perturbations is a very important topic in engineering, not only for the determination of the behavior of a system with respect to perturbations but also for the construction of algorithms for the numerical simulation of the problem.

The originality of this work, compared to that in [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF], lies in the fact that the singlevalued term may be nonlinear (which allows us sometimes to simplify the multivalued part of the generalized equation), and that the multivalued term is not necessarily a subdifferential. The main contribution of this paper is the study of the calmness property. More precisely, we give a sufficient condition ensuring the calmness property at a reference point; it is based on the outer subdifferential of a particular distance function (Proposition 5.1). This condition is not easy to check in practice in general. Nevertheless, we succeeded to apply this condition for a meaningful example from electrical engineering (Example 6.2). As far as we know, this is the first time that this kind of properties is applied to some concrete models in electronics.

The paper is organized as follows. In Sect. 2, we pose the problem and recall several standard notions from variational analysis. The three subsequent ones are devoted to the criteria ensuring the Aubin, the isolated calmness and the calmness properties of the solution mapping, respectively. Section 6 deals with some particular cases arising in electrical circuit theory. A conclusion is presented in Sect. 7.

Preliminaries

Given matrices B ∈ R n×m , C ∈ R m×n with m ≤ n, a vector p ∈ R n and mappings f : R n → R n , F : R m ⇒ R m , we consider the problem of finding a solution z ∈ R n to the generalized equation

p ∈ f (z) + BF (Cz).
(1)

We will refer to several additional assumptions in the main results. Namely,

(A1) B is injective; (A2) f is continuously differentiable on R n ; (A3) F has closed graph; (A4) C is surjective; and
(A5) there are F j : R ⇒ R, j ∈ {1, . . . , m} such that

F (x) = m j =1 F j (x j ) whenever x = (x 1 , . . . , x m ) T ∈ R m .
In [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF], the authors considered the special case of the above inclusion with the linear single-valued part f , with B = C T , and F being the Clarke subdifferential of the super-potential J defined by

J (x) := j 1 (x 1 ) + j 2 (x 2 ) + • • • + j m (x m ) whenever x = (x 1 , . . . , x m ) T ∈ R m ,
with j i : R → R being a locally Lipschitz continuous function for each index i = 1, . . . , m.

The notation is fairly standard. In R d , the norm, the scalar product, the closed ball with the center x ∈ R d with the radius r ≥ 0, and the unit ball are denoted by • ,

•, • , B(x, r), and B, respectively. First, recall some well-known subdifferential concepts. Given a lower semicontinuous function g : R d → R := R∪{±∞} and a point x ∈ R d such that g(x) ∈ R, the Fréchet/regular subdifferential of g at x is the set

∂ F g(x) containing those ξ ∈ R d that satisfy lim inf 0 =h→0 g(x + h) -g(x) -ξ, h h ≥ 0;
the limiting (Mordukhovich) subdifferential of g at x is the set ∂g(x) containing all ξ ∈ R d such that there are sequences (x k ) k∈N and (ξ k ) k∈N converging to x and ξ , respectively, with

g x k → g(x) as k → ∞, and ξ k ∈ ∂ F g x k for each k ∈ N;
and finally, the outer subdifferential of g at x is the set ∂ > g(x) containing those ξ ∈ R d for which there are sequences (x k ) k∈N and (ξ k ) k∈N converging to x and ξ , respectively, with

g x k ↓ g(x) as k → ∞ and ξ k ∈ ∂ F g x k for each k ∈ N.
Note that nothing will change if for each k ∈ N one will take ξ k from ∂g(x k ) instead of ∂ F g(x k ) in the definition of the outer subdifferential. If g(x) is infinite, then all the above subdifferentials of g at x are defined to be empty sets. Furthermore, note that the choice of any other (equivalent) norm on R d (instead of the usual one induced by the scalar product) does not affect the above subdifferential constructions. Given a non-empty subset

Ω of R d , then d(x, Ω) := inf { y -x : y ∈ Ω} is the distance from x ∈ R d to Ω and δ Ω denotes the indicator function of Ω (i.e., δ Ω (y)
is equal to 0 if y ∈ Ω and to +∞ otherwise). It is well known (e.g., see [2, p. 111

]) that, if Ω is closed and does not contain x ∈ R d , then ∂d(•, Ω)(x) = (x -w)/d(x, Ω) : w ∈ Ω with x -w = d(x, Ω) , (2)
and

∂ F d(•, Ω)(x) = (x -w)/d(x, Ω), if ∃!w ∈ Ω with d(x, Ω) = x -w , ∅, otherwise. (3)
Fix any x ∈ Ω for a while; the Fréchet/regular normal cone to Ω at x is the set

N( x; Ω) := ∂ F δ Ω ( x) = ξ ∈ R d : lim sup Ω x→ x ξ, x -x x -x ≤ 0 ;
similarly, N( x; Ω) := ∂δ Ω ( x) is the Mordukhovich limiting normal cone to Ω at x (which contains all ξ ∈ R d for which there are sequences (x k ) k∈N in Ω and (ξ k ) k∈N in R d converging to x and ξ , respectively, such that ξ k ∈ N(x k ; Ω) for each k ∈ N); and, finally, the Bouligand-Severi tangent cone T ( x; Ω) to Ω at x contains those v ∈ R d for which there are sequences (t k ) k∈N in ]0, ∞[ and (v k ) k∈N in R d converging to 0 and v, respectively, such that x + t k v k ∈ Ω whenever k ∈ N.

By S : R d ⇒ R l we mean a multivalued mapping from R d into R l with the domain dom S, the graph gph S and the range rge S. Fix a point ( x, ȳ) in gph S, then the contingent (graphical) derivative of S at ( x, ȳ) is the mapping DS( x, ȳ) : R d ⇒ R l defined by

DS( x, ȳ)(u) := v ∈ R l : (u, v) ∈ T ( x, ȳ); gph S , u∈ R d ;
and the Mordukhovich coderivative of S at ( x, ȳ) is the set-valued mapping D * S( x, ȳ) : R l ⇒ R d , assigning to each ξ ∈ R l the set

D * S( x, ȳ)(ξ ) := η ∈ R d : (η, -ξ) ∈ N ( x, ȳ); gph S .
The mapping S is said to have Aubin/Lipschitz-like property around ( x, ȳ) iff there is a constant κ ≥ 0 along with neighborhoods U of x and V of ȳ such that

S(x) ∩ V ⊂ S x + κ x -x B for all x, x ∈ U. ( 4
)
The infimum of the set of values κ, for which there are neighborhoods U of x and V of ȳ such that (4) holds, is called Lipschitz modulus of S at ( x, ȳ) and is denoted by lip(S; ( x, ȳ)). It is well known (see [START_REF] Dontchev | Implicit Functions and Solution Mappings. A View from Variational Analysis[END_REF]Theorem 3E.6]) that the Aubin/Lipschitz-like property of S around ( x, ȳ) is equivalent to the metric regularity of Φ := S -1 around ( ȳ, x) and that the regularity modulus of Φ at ( ȳ, x) equals the Lipschitz modulus of S at ( x, ȳ).

Fixing one of the points in U , i.e. the request for the existence of a constant κ ≥ 0 along with neighborhoods U of x and V of ȳ such that

S(x) ∩ V ⊂ S( x) + κ x -x B whenever x ∈ U, ( 5
)
we get the definition of the calmness property of S at ( x, ȳ). Again, the infimum of the set of values κ, for which there are neighborhoods U of x and V of ȳ such that (5) holds, is called calmness modulus of S at ( x, ȳ) and is denoted by clm(S; ( x, ȳ)).

Remember that ref. [START_REF] Dontchev | Implicit Functions and Solution Mappings. A View from Variational Analysis[END_REF]Theorem 3H.3] establishes the equivalence between calmness property of S at ( x, ȳ) and the metric subregularity of Φ at ( ȳ, x). Moreover, the subregularity modulus of Φ at ( ȳ, x) and the calmness modulus of S at ( x, ȳ) are the same.

Further, the mapping S has the isolated calmness property at ( x, ȳ) provided that it has calmness property and ȳ is an isolated point of S( x), i.e. S( x) ∩ B( ȳ, r) = { ȳ} for some r > 0. Note that the isolated calmness property of S at ( x, ȳ) amounts to the existence of a constant κ ≥ 0 along with neighborhoods U of x and V of ȳ such that

S(x) ∩ V ⊂ ȳ + κ x -x B for each x ∈ U. ( 6
)
Theorem 3I.2 in [START_REF] Dontchev | Implicit Functions and Solution Mappings. A View from Variational Analysis[END_REF] says that the isolated calmness property of S at ( x, ȳ) is equivalent to the strong metric subregularity of Φ at ( ȳ, x).

From now on, we are going to use the following notations.

Standing Assumptions Denote by Φ the set-valued mapping from R n into itself defined by Φ(z)

:= f (z) + BF (Cz) whenever z ∈ R n . Set S := Φ -1 . Let us define the mappings Q : R n ⇒ R n and F C : R m ⇒ R m by Q(z) := BF (Cz), z ∈ R n , and F C (u) := F (u) if u = Cz for some z ∈ R n and F C (u) := ∅ other- wise.
We also suppose that we have in hand a point (z, p) ∈ gph Φ. Finally, put

v := (B T B) -1 B T ( p -f (z)).
3 Aubin/Lipschitz-Like Property of S at ( p, z)

The following result provides the formula for the coderivative of the mapping Φ at the reference point.

Proposition 3.1 Under the assumptions (A1)-(A3), for any

ξ ∈ R n one has D * Φ(z, p)(ξ ) = ∇f (z) T ξ + C T D * F C (C z, v) B T ξ .
Proof Fix any ξ ∈ R n . As (A2) ensures the strict differentiability of f at z, ref.

[2, Theorem 1.62(ii)] implies that

D * Φ(z, p)(ξ ) = ∇f (z) T ξ + D * Q z, p -f (z) (ξ ). Define h : R n × R m → R n and G : R n ⇒ R n for each z ∈ R n by h(z, v) := Bv, v ∈ R m , and G(z) := F (Cz). Thus h(z, G(z)) = {Bv : v ∈ F (Cz)} = Q(z) whenever z ∈ R n . Moreover, v ∈ G(z) and h(z, v) = B v = p -f (z). Because of (A1), ref. [2, Lemma 1.126] in the finite-dimensional setting reveals that D * Q z, p -f (z) (ξ ) = D * G(z, v) B T ξ .
Finally, ref. [START_REF] Mordukhovich | Restrictive metric regularity and generalized differential calculus in Banach spaces[END_REF]Theorem 3.10] yields that

D * G(z, v) B T ξ = C T D * F C (C z, v) B T ξ .
Combine the above expressions to conclude the proof.

Remark 3.1 This proposition is proved using the first order calculus rules instead of a second order chain rule used in [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF], as it is valid for the particular case of B = C T only.

The above statement allows us to characterize the Aubin/Lipschitz-like property of S at the reference point.

Theorem 3.1 Under the assumptions (A1)-(A3), S has the Aubin/Lipschitz-like property at ( p, z) if and only if

0 ∈ ∇f (z) T ξ + C T D * F C (C z, v) B T ξ ⇒ ξ = 0. ( 7
)
Moreover, its Lipschitz modulus is given by

lip S; ( p, z) = sup ξ : ∇f (z) T ξ + C T D * F C (C z, v) B T ξ ∩ B = ∅ . Proof Note that η ∈ D * Φ(z, p)(ξ ) if and only if -ξ ∈ D * S( p, z)(-η).
Having Proposition 3.1 in hand, one applies the well-known Mordukhovich criterion (see [START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation. I. Basic Theory[END_REF]Theorem 4.10]) in the finite-dimensional setting to conclude the proof.

Employing more assumptions, we may get the following corollaries.

((CC T ) -1 C∇f (z) T ξ, B T ξ) ∈ -N ((C z, v); gph F ) ∇f (z) T ξ ∈ rge C T =⇒ ξ = 0. (8) Proof Indeed, if C is surjective, then F C = F and CC T ∈ R m×m is non-singular. First, let ξ ∈ R n be such that 0 ∈ ∇f (z) T ξ + C T D * F (C z, v)(B T ξ). Find w ∈ D * F (C z, v)(B T ξ) with ∇f (z) T ξ + C T w = 0. Thus -(CC T ) -1 C∇f (z) T ξ is in D * F (C z, v)(B T ξ).
Clearly, we have ∇f (z) T ξ ∈ rge C T and the definition of the coderivative of F yields the rest.

On the other hand, pick any

ξ ∈ R n with ((CC T ) -1 C∇f (z) T ξ, B T ξ) in -N ((C z, v); gph F ) and ∇f (z) T ξ ∈ rge C T . The definition of the coderivative of F says that w := -(CC T ) -1 C∇f (z) T ξ ∈ D * F (C z, v)(B T ξ). Thus CC T w = -C∇f (z) T ξ . This implies that C T w + ∇f (z) T ξ ∈ ker C ∩ rge C T . Therefore 0 ∈ ∇f (z) T ξ + C T D * F (C z, v)(B T ξ). Remark 3.2 If (A5) also holds, then gph F = m j =1 gph F j , hence Proposition 6.41 in [5] implies that N ((C z, v); gph F ) = m j =1 N (((C z) j , vj ); gph F j ).
Thus the first condition in ( 8) can be checked coordinate-wise. [START_REF] Dontchev | Metric regularity under approximations[END_REF], the Aubin/Lipschitz-like property is stable with respect to a small perturbation of the function f . To be more precise, suppose that S has the Aubin/Lipschitz-like property at ( p, z). Then there exists > 0 such that for every function f :

Remark 3.3 In view of Theorem 1 in

R n → R n , every point ( p, z) ∈ B( p, ) × B(z, ) with p ∈ f (z) + BF (C z), f (z) -f (z) ≤ , and f z -f z -f z -f z ≤ z -z whenever z , z ∈ B(z, ),
the mapping S := ( f + BF (C•)) -1 has the Aubin/Lipschitz-like property at ( p, z).

Isolated Calmness of S at ( p, z)

In this case, we are going to compute the graphical derivative of Φ at the reference point. Let us start with the following geometric lemma which is only a slight generalization of Lemma 4.1 in [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF] (where the first matrix in question was supposed to be surjective and Γ was a closed set).

Lemma 4.1 Let E ∈ R k×d be any matrix, let G ∈ R l×d be injective, and let Γ be a subset of rge E. Put Ξ := E -1 (Γ ) and Λ := G(Ξ ). For x ∈ Λ denote by ȳ the

(unique) point in Ξ with G ȳ = x. Then T ( x; Λ) = u ∈ R l : ∃w ∈ R d such that u = Gw and Ew ∈ T (E ȳ; Γ ) . Proof We claim that T ( ȳ; Ξ) = {w ∈ R d : Ew ∈ T (E ȳ, Γ )}. First, take any w ∈ T ( ȳ; Ξ). Find (t n ) n∈N in ]0, ∞[ and (w n ) n∈N in R d converging
( ȳ + t n w n ) ∈ Γ , hence ȳ + t n w n ∈ E -1 (Γ ) = Ξ . So w ∈ T ( ȳ; Ξ). The claim is proved.
Using exactly the same steps as in the proof of Lemma 4.1 in [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF] one gets that T ( x; Λ) = {Gw : w ∈ T ( ȳ, Ξ )}. This and the claim yield the assertion.

Proposition 4.1 Under the assumptions (A1)-(A2), for any

b ∈ R n one has DΦ(z, p)(b) = ∇f (z)b + BDF C (C z, v)(Cb). Proof Fix any b ∈ R n . By Proposition 4A.2 in [3], we have DΦ(z, p)(b) = ∇f (z)b + DQ z, p -f (z) (b). Observe that gph Q = u v ∈ R 2n : ∃ b c ∈ R n × R m : u v = G b c and E b c ∈ gph F C with G := I n 0 0 B and E := C 0 0 I m . As B is injective, so is G. Lemma 4.1 (with k := 2m, l := 2n, d := n + m, Γ := gph F C , x := (z, p -f (z)) T , and ȳ := (z, v) T ) reveals that T z, p -f (z) ; gph Q = b Bc : Cb c ∈ T (C z, v); gph F C .
This means that DQ(z, p -

f (z))(b) = BDF C (C z, v)(Cb).
The assertion is proved.

0 ∈ ∇f (z)b + BDF C (C z, v)(Cb) ⇒ b = 0. ( 9
)
Moreover, its calmness modulus is given by

clm S; ( p, z) = sup b : ∇f (z)b + BDF C (C z, v)(Cb) ∩ B = ∅ . Proof Note that x ∈ DΦ(z, p)(y) if
(Cb, -(B T B) -1 B T ∇f (z)b) ∈ T ((C z, v); gph F ) ∇f (z)b ∈ rge B =⇒ b = 0. ( 10
) Proof Indeed, if C is surjective, then F C = F . Note that (A1) ensures that B T B ∈ R m×m is non-singular. First, let b ∈ R n be such that 0 ∈ ∇f (z)b + BDF (C z, v)(Cb). Find a point w ∈ DF (C z, v)(Cb) with ∇f (z)b + Bw = 0. Thus -(B T B) -1 B T ∇f (z)b is in DF (C z, v)(Cb)
. Clearly, we have ∇f (z)b ∈ rge B and the definition of the contingent derivative of F yields the rest.

On the other hand, pick any b ∈ R n with (Cb, -(

B T B) -1 B T ∇f (z)b) in T ((C z, v); gph F ) and ∇f (z)b ∈ rge B. The definition of the contingent derivative says that w := -(B T B) -1 B T ∇f (z)b ∈ DF (C z, v)(Cb). Thus we have B T Bw = -B T ∇f (z)b. So Bw + ∇f (z)b ∈ ker B T ∩ rge B = {0}. Therefore 0 ∈ ∇f (z)b + BDF (C z, v)(Cb).
Remark 4.1 If (A5) is also satisfied and gph F j is normally (Clarke) regular at ((C z) j , vj ) for each j ∈ {1, . . . , m} in sense of Definition 6.4 in [START_REF] Rockafellar | Variational Analysis[END_REF], then Proposition 6.41 in [START_REF] Rockafellar | Variational Analysis[END_REF] implies that T ((C z, v); gph F ) = m j =1 T (((C z) j , vj ); gph F j ). Hence the first condition in (10) can be checked coordinate-wise.

In view of Theorem 3 in [START_REF] Dontchev | Metric regularity under approximations[END_REF], also the isolated calmness property is stable with respect to the perturbation of f in the sense discussed in Remark 3.3.

Calmness of S at ( p, z)

Let us point out that the criteria guaranteeing the calmness are more complicated than those for the isolated calmness. Hence, one should always consider the latter property first, and employ more sophisticated tools of this section only in case of their failure. One can derive the following sufficient condition for calmness.

0 / ∈ ∂ > h(z). ( 11
)
Proof The inclusion (1) says that, for each p ∈ R n , one has

S(p) = z ∈ R n : Cz, B T B -1 B T p -f (z) , p -f (z) ∈ Λ .
Indeed, fix any p ∈ R n . First, pick any z ∈ R n with (z, p) ∈ gph Φ. Inclusion (1) reveals that pf (z) ∈ rge B. We have already mentioned that B T B ∈ R m×m is nonsingular. Hence, (1) implies that (Cz, (

B T B) -1 B T (p -f (z))) ∈ gph F . On the other hand, fix an arbitrary z ∈ R n with Cz, B T B -1 B T p -f (z) , p -f (z) ∈ Λ. Then B T (p -f (z)) = B T Bw for some w ∈ F (Cz). As p -f (z) ∈ rge B, we have p -f (z) -Bw ∈ ker B T ∩ rge B = {0}. Therefore (z, p) ∈ gph Φ.
Define the mapping M : R m × R m × R n ⇒ R n as follows:

M(y) := z ∈ R n : g(z) + y ∈ Λ , y ∈ R m × R m × R n .
As Λ is closed and the mapping

R n × R n (p, z) → Cz, B T B -1 B T p -f (z) , p -f (z) ∈ R m × R m × R n (12)
is continuously differentiable (hence strictly differentiable and therefore locally Lipschitz continuous), Lemma 1 in [START_REF] Henrion | Calmness of constraint systems with applications[END_REF] says that S is calm at ( p, z), provided that so is M at (0, z). As observed in [7, p. 438], M is calm at (0, z) if and only if there are L > 0 and ε > 0 such that

d z, g -1 (Λ) ≤ Ld g(z), Λ whenever z ∈ B(z, ε). Since g is locally Lipschitz continuous, so is h. Moreover, [h ≤ 0] = g -1 (Λ) where [h ≤ 0] := {u ∈ R n : h(u) ≤ 0}. Theorem 2.
1 [START_REF] Ioffe | On metric and calmness qualification conditions in subdifferential calculus[END_REF] reveals that M is calm at (0, z), provided that (11) holds. The proposition is proved.

Remark 5.1 Clearly, if the matrix B is surjective, then it suffices to consider Λ = gph F and g(z) = (Cz, B -1 ( pf (z))), z ∈ R n . Further, Theorem 2.1 in [START_REF] Ioffe | On metric and calmness qualification conditions in subdifferential calculus[END_REF] also gives an upper estimate of the corresponding calmness modulus. Namely, if there is γ > 0 such that ξ ≥ γ for each ξ ∈ ∂ > h(z), then the set-valued mapping M from the proof of the above proposition is calm at (0, z) with modulus not exceeding 1/γ . Denote by ϕ the mapping defined in (12). A closer look at the proof of Lemma 1 in [START_REF] Henrion | Calmness of constraint systems with applications[END_REF] reveals that S is calm at ( p, z) with modulus not exceeding K/γ , where K > 0 is the Lipschitz constant of ϕ at ( p, z). To be more precise, there is ε > 0 such that ϕ p , zϕ p , z ≤ K pp whenever p , p ∈ B( p, ε), z ∈ B(z, ε).

Remark 5.2 Similar steps, as in the proof of the above theorem, show that the condition 0 / ∈ ∂h(z) ensures the (stronger) Aubin/Lipschitz-like property of S at the reference point (as observed by one of the referees). However, we prefer to use the well-known Mordukhovich criterion.

The function h is a composite function, so one can try to use some chain rule. For example, since g(z) ∈ Λ, one may use Proposition 6.4 in [START_REF] Hoheisel | Exact penalty results for mathematical programs with vanishing constraints[END_REF] (with the closed set C := Λ, with the Lipschitz continuous non-negative function ϑ := d(•, Λ), and the continuously differentiable mapping F := g) to get

∂ > h(z) ⊂ ∇g(z) T ∂ > d(•, Λ) g(z) .
However, in our setting this does not seem to be the right way to proceed, since the resulting set can be too large.

f (z) := z, z ∈ R, F(z):= ⎧ ⎪ ⎨ ⎪ ⎩ -1 -z, z < 0, [-1, 1], z = 0, 1 -z, z > 0.
Let z := 0 and p := -1. It is easy to verify that both (8) and (10) are violated. So S has neither Aubin nor the isolated calmness property at the reference point. Note that

h(z) = d (z, -1 -z, -1 -z), gph F × R = d (z, -1 -z), gph F , z ∈ R.

Let us compute ∂ > h(0).

Recall that this set contains those ξ ∈ R such that there are sequences (z k ) k∈N and (ξ k ) k∈N converging to 0 and ξ , respectively, with

h z k ↓ h(0) = 0 as k → ∞ and ξ k ∈ ∂ F h z k for each k ∈ N.
If z < 0, then the point (z, -1z) T ∈ gph F , thus h(z) = 0, hence we have to consider (z k ) k∈N in ]0, ∞[ only. Moreover, there is δ > 0 such that, for each z ∈]0, δ[, the point (z, -1z) has the unique nearest point (0, -1) in gph F . Therefore

h(z) = z 2 + (-1 -z + 1) 2 = √ 2z, z ∈]0, δ[.

Hence, ∂ F h(z) = { √

2} for such a point z. Let (z k ) k∈N be any sequence converging to 0 such that h(z k ) ↓ 0 as k → ∞. We may assume without any loss of generality that

z k ∈]0, δ[ for each k ∈ N. Thus ∂ F h(z k ) = { √ 2}, hence ∂ > h(0) = { √ 2}.
Therefore S is calm at the reference point with modulus not exceeding 1/ √ 2. To show the failure when one uses the chain rule, let g : R → R 2 be for each z ∈ R defined by g(z) := (z, -1z) T . Then ∇g(0) = (1, -1) T . As ȳ := g(0) = (0, -1) ∈ gph F , we arrive at the condition

0 / ∈ α -β : (α, β) T ∈ ∂ > d(•, gph F )( ȳ) .
Let z < 0 be arbitrary. Then the point (z, -1 + z) has the unique nearest point (0, -1) in gph F . Hence (3) implies that (-1, -1)

T = (z, z) T /(|z| √ 2) is in ∂ F d(•, gph F )(z, -1 + z). For each k ∈ N, put y k := (-1/k, -1 -1/k) and ξ k := (-1, -1) T . So we have y k → ȳ and d(y k , gph F ) ↓ d( ȳ, gph F ) = 0 as k → ∞ with ξ k ∈ ∂ F d(•, gph F )(y k ) for each k ∈ N. Thus ξ := (-1, -1) T is in ∂ > d(•, gph F )( ȳ).
So the above sufficient condition for calmness is not satisfied.

Let z < 0. Then p = -1 is the only point of Φ(z). Again, for each z ∈ R we have h(z) = d((z, -1z), gph F ). Then there is δ > 0 such that (z, -1z) is in gph F whenever z ∈]zδ, z + δ[. Thus h(z) = 0 for such a point z. This means that ∂ > h(z) = ∅ because no sequence (z k ) k∈N converging to z satisfies h(z k ) > 0 for each k ∈ N. Consequently, S is calm at the reference point. Note that any such a point z is an interior point of S(-1).

Application to Electrical Circuits

A nonlinear term f in (1) allows often to simplify the set-valued mapping F by splitting it into a single-valued part and (possibly simpler) multivalued one. Then the first mentioned part can be added to f . We illustrate this in the following example. Example 6.1 Consider the circuit in Fig. 1 involving a load resistance R > 0, a source E > 0, an input-signal source u and corresponding instantaneous current i. Suppose that the voltage-current characteristic of the DIAC is given by a set-valued mapping Ψ : R ⇒ R. Assume that Ψ (0) := [-1, 1], that Ψ is single-valued and continuously differentiable on R \ {0}, that its graph is symmetric with respect to the origin, that -1 < Ψ (z) < -1z whenever z ∈] -∞, 0[, and finally that Ψ (0-) = -1 and Ψ (0-) = -a for some a > 0 (see Fig. 2). Kirchhoff's voltage law reveals that

u -E ∈ Ri V R + Ψ (i) V D .
Assume that R = 1. Putting F = ∂| • |, one infers that there is a continuously differentiable odd function ψ : R → R such that Ψ = ψ + F . Moreover, ψ(0) = 0, ψ (0) = -a, and 0 < ψ(z) < -z whenever z < 0 (see Fig. 3). Putting z := i and Let (z, p) = (0, -1) ∈ gph Φ. In [1, Example 5.1] it is proved that S has Aubin/Lipschitz-like property (isolated calmness property) at ( p, z) if and only if 0 < a < 1 (a = 1). There it is also shown that lip(S; ( p, z)) = 1/(1a). First, to illustrate our approach, let us use Corollary 4.1 to check isolated calmness and Theorem 4.1 to compute the calmness modulus (which was not done in [START_REF] Adly | Qualitative stability of a class of non-monotone variational inclusions. Application in electronics[END_REF]).

Note that f (z) = 1a and v = -1. A simple calculation yields that

T (0, -1); gph F = R + 0 1 ∪ R + -1 0 .
Thus, the necessary and sufficient condition (10) for the isolated calmness property of S at ( p, z)

with b ∈ R is b, -(1 -a)b T ∈ T (0, -1); gph F =⇒ b = 0.
In xy-plane (see Fig. 4), this means that the intersection of the above tangent cone with the graph of the function y = -(1a)x, x ∈ R, is a singleton containing the origin only. Clearly, this is the case if and only if a = 1. To compute clm(S; ( p, z)) one has to find a maximum of |x| subject to the constraint -1 ≤ (1a)x + y ≤ 1 with (x, y) ∈ T ((0, -1); gph F ). Which immediately gives (see Fig. 4) that clm S; ( p, z) =

1 |a-1| , if a = 1, ∞, otherwise.
So, if a = 1, then S has neither Aubin nor isolated calmness property at ( p, z). Note that the sufficient condition (11) for calmness property is not satisfied in this case. Indeed, (z, p) ∈ gph F ∩ gph Φ and h : R → R + is for each z ∈ R defined by Hence, ∂ F h(z) = {-1ψ (z)} whenever z < 0. For each k ∈ N, put z k := -1 k and ξ k := -1ψ (z k ). Then both (z k ) k∈N and (ξ k ) k∈N converge to 0. For any k ∈ N, we have

h(z) := d z, -1 -z -ψ(z) , gph F . Note that -1 -z -ψ(z) < -1 -z and ψ(z) < -z whenever z < 0. Given z < 0, the point (z, -1 -z -ψ(z)) T has the unique nearest point (z, -1) T in gph F . Therefore h(z) = |z + ψ(z)| = -z -ψ(z).
ξ k ∈ ∂ F h(z k ), and ψ(z k ) < -z k implies that -1 -z k -ψ(z k ) + 1 is positive, so h(z k ) > 0. Thus 0 ∈ ∂ > h(0).
However, S is not calm at (-1, 0). Indeed, suppose on the contrary that it is. As S(-1) = {z ∈ R : (z, -1zψ(z)) T ∈ gph F } = {0}, the point 0 is an isolated point of S(-1), thus S would have the isolated calmness property at (-1, 0), which is not the case as seen before.

Example 6.2 Let us consider the circuit depicted in Fig. 5 involving a load resistance R > 0, an input-signal source u and corresponding instantaneous current i, an inductor with inductance L > 0, a capacitor with capacitance C 0 > 0, and a diode with the voltage-current characteristics given for some α 1 < 0 < α 2 by

F (x) := ⎧ ⎪ ⎨ ⎪ ⎩ α 1 , x<0, [α 1 , α 2 ], x = 0, α 2 , x>0.
Using the Kirchhoff voltage law, we have

u -U R -U L -U C 0 = U D ,
where

U R = Ri, U L = L di dt , U C 0 = 1 C 0 i(t)
dt and U D ∈ -F (-i) denote the difference of potential across the resistor, the inductor, the capacitor, and the diode, respectively. Setting

z 1 := i dt and z 2 := ż1 (= i),
we have

Lż 2 = - 1 C 0 z 1 -Rz 2 + u -U D .
Hence, dividing by L, we arrive at the dynamic system

ż1 ż2 = 0 1 -1 LC 0 -R L z 1 z 2 + 0 1 L u + 0 -1 L U D with U D ∈ -F (0 -1) z 1 z 2 .
From now on, assume that R = L = C 0 = 1, that α 1 = -100, and that α 2 = 1. The steady states of the above dynamic system may be viewed as particular solutions to (1) with n = 2, m = 1,

A := 0 1 -1 -1 , B:= 0 1 , C:= (0 -1),
and

f (z) := Az, z ∈ R 2 and F (x) := ⎧ ⎪ ⎨ ⎪ ⎩ -100, x<0, [-100, 1], x = 0, 1, x>0.
So, for z = (z 1 , z 2 ) T ∈ R 2 and p = (p 1 , p 2 ) T ∈ R 2 , (1) reads as

p 1 = z 2 , p 2 ∈ -z 1 -z 2 + F (-z 2 ).
Moreover, CC T = B T B = 1 and rge B = rge C T = {0} × R. Let (z, p) ∈ gph Φ be arbitrary. The necessary and sufficient condition [START_REF] Ioffe | On metric and calmness qualification conditions in subdifferential calculus[END_REF] for Aubin/Lipschitz-like property of S at ( p, z)

with ξ = (ξ 1 , ξ 2 ) T ∈ R 2 is (ξ 2 -ξ 1 , ξ 2 ) T ∈ -N ((C z, v); gph F ) (-ξ 2 , ξ 1 -ξ 2 ) T ∈ {0} × R. =⇒ ξ 1 = 0 and ξ 2 = 0. ( 13
)
The necessary and sufficient condition (10) for the isolated calmness property of S at

( p, z) with b = (b 1 , b 2 ) T ∈ R 2 is (-b 2 , b 1 + b 2 ) T ∈ T ((C z, v); gph F ) (b 2 , -b 1 -b 2 ) T ∈ {0} × R =⇒ b 1 = 0 and b 2 = 0. ( 14
)
Note that the latter inclusion in (13) and ( 14) immediately implies that ξ 2 = 0 and b 2 = 0, respectively. First, consider five particular cases (see Fig. 6). Hence C z = 1, v = 1, and

T (1, 1); gph F = R 1 0 and N (1, 1); gph F = R 0 1 .
Clearly, (13) is satisfied. So S has the Aubin/Lipschitz-like property (hence it has the calmness property) at ( p, z). Moreover, (14) holds as well. Thus S has the isolated calmness property at ( p, z). 3. Let z := (0, 0) T and p := (0, 1) T .

Then C z = 0 and v = 1. A simple calculation reveals that the tangent cone T ((0, 1); gph F ) and the normal cone N((0, 1); gph F ) are equal to

R + 0 -1 ∪ R + 1 0 and R -1 0 ∪ R 0 1 ∪ cone -1 0 , 0 1 ,
respectively. Clearly, (13) is violated (take ξ 1 := 1 and ξ 2 := 0). So S does not have the Aubin/Lipschitz-like property. Moreover, ( 14) is violated too (take b 1 := -1 and b 2 := 0). Thus S does not have the isolated calmness property at ( p, z).

To check the calmness property, we are going to apply Proposition 5.1. Note that

Λ = gph F × {0} × R. The functions g : R 2 → R 2 × R 2 and h : R 2 → R + are for each z = (z 1 , z 2 ) T ∈ R 2 defined by g(z) = (-z 2 , 1 + z 1 + z 2 , -z 2 , 1 + z 1 + z 2 ) T and h(z) = d g(z), Λ = d (-z 2 , 1 + z 1 + z 2 , -z 2 , 1 + z 1 + z 2 ) T , gph F × {0} × R ,
respectively. Hence, we need to prove that (0, 0) T / ∈ ∂ > h(z). Without any loss of generality, assume that the norm on R 2 × R 2 is for each its point (u, v) given by (u, v) = u 2 + v 2 , where • 2 is Euclidean norm on R 2 . First, we compute ∂d(•, gph F )(u) for each u = (x, y) T ∈ R 2 \ gph F sufficiently close to ū := (0, 1) T . In view of (2), it suffices to find nearest points to u in gph F . This yields r ∈]0, 1[ such that for each u = (x, y) T ∈ B 2 ( ū, r) one has

∂d(•, gph F )(u) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ {(0, 1) T }, if x ≥ 0, y > 1, {(-1, 0) T }, if x < 0, y ≤ 1, 1 √ x 2 +(y-1) 2 (x, y -1) T , if x < 0, y > 1, {(0, -1) T }, if x > 0, 1 > y > 1 -x, {(1, 0) T }, if x > 0, y < 1 -x, and {(1, 0) T , (0, -1) T }, if x > 0, y = 1 -x.
Second, we claim that given z = (z 1 , z 2 ) T ∈ R 2 and η ∈ R 2 such that h(z) > 0, that p -Az ∈ B 2 ( ū, r), and that η ∈ ∂h(z) one has η ≥ 1. Indeed,

h(z) = d (-z 2 , 1 + z 1 + z 2 ) T , gph F + d (-z 2 , 1 + z 1 + z 2 ) T , {0} × R = d (-z 2 , 1 + z 1 + z 2 ) T , gph F + |z 2 | = d( p -Az, gph F ) + |z 2 |.
Since, h is the sum of two Lipschitz functions, Theorem 3.36 in [START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation. I. Basic Theory[END_REF] implies that

∂h(z) ⊂ ∂d p -A(•), gph F (z) + {0} × ∂| • |(z 2 ).
As A is surjective, Proposition 1.112(i) in [START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation. I. Basic Theory[END_REF] reveals that

∂d p -A(•), gph F (z) = -A T ∂d(•, gph F )(-z 2 , 1 + z 1 + z 2 ).
Collecting the above partial results we are in position to compute an upper estimate of ∂h(z). Recall that h(z) > 0, thus it suffices to consider the following cases.

(i) If z 2 = 0, then z 1 > 0 because 1 + z 1 / ∈ F (0) (the case z 1 < -101 is impossible as p -Az = (0, 1 + z 1 ) T ∈ B 2 ( ū, r)). Thus ∂h(z) ⊂ 0 1 -1 1 0 1 + 0 [-1, 1] = 1 [0, 2] . (ii) If z 2 < 0 and 1 + z 1 + z 2 > 1, then ∂h(z) ⊂ 0 1 -1 1 0 1 + 0 -1 = 1 0 . (iii) If z 2 < 0 and 1 + z 1 + z 2 < 1, then ∂h(z) ⊂ 0 1 -1 1 0 -1 , 1 0 + 0 -1 = -1 -2 , 0 -2 .
(iv) If z 2 > 0 and 1

+ z 1 + z 2 ≤ 1, then ∂h(z) ⊂ 0 1 -1 1 -1 0 + 0 1 = 0 2 . (v) If z 2 > 0 and 1 + z 1 + z 2 > 1, then ∂h(z) ⊂ 1 z 2 2 + (z 1 + z 2 ) 2 0 1 -1 1 -z 2 z 1 + z 2 + 0 1 = ⎛ ⎝ z 1 +z 2 √ z 2 2 +(z 1 +z 2 ) 2 z 1 +2z 2 √ z 2 2 +(z 1 +z 2 ) 2 + 1 ⎞ ⎠ .
In any case, η ≥ 1 for each η ∈ ∂h(z), the claim is proved.

To conclude, pick an arbitrary ξ ∈ ∂ > h(z) if there is any. Find a sequence

z k = (z k 1 , z k 2 ) T in R 2 converging to z = (0, 0) T along with a sequence (ξ k ) k∈N converging to ξ such that h z k ↓ h(z) = 0 as k → ∞ and ξ k ∈ ∂h z k for each k ∈ N. As p = ū, there is k 0 ∈ N such that p -Az k ∈ B 2 ( ū, r) whenever k > k 0 . So, the previous claim implies that ξ ≥ 1. In particular, (0, 0) T / ∈ ∂ > h(z).
Thus S is calm at ( p, z) with the calmness modulus not exceeding 2. Similarly to the previous case, one infers that S has neither the Aubin/Lipschitz-like nor the isolated calmness property at ( p, z).

To check the calmness, note that for each z

= (z 1 , z 2 ) T ∈ R 2 we have g(z) = (-z 2 , z 1 + z 2 , -z 2 , z 1 + z 2 ) T and h(z) = d (-z 2 , z 1 + z 2 , -z 2 , z 1 + z 2 ) T , gph F × {0} × R . First, find again r ∈]0, 1[ such that ∂d(•, gph F )(u) = {( x |x| , 0) T } whenever u = (x, y) T ∈ B 2 ( ū, r) with ū := (0, 0) T .
Second, we claim that given z = (z 1 , z 2 ) T ∈ R 2 and η ∈ R 2 such that h(z) > 0, that -Az ∈ B 2 ( ū, r), and that η ∈ ∂h(z) one has η ≥ 2. Indeed,

h(z) = d (-z 2 , z 1 + z 2 ) T , gph F + |z 2 | = d(-Az, gph F ) + |z 2 |. and ∂h(z) ⊂ -A T ∂d(•, gph F )(-Az) + {0} × ∂| • |(z 2 )
. Suppose that z 2 = 0, then the assumption h(z) > 0 implies that z 1 / ∈ F (0), so either z 1 > 1 or z 1 < -100. Both the cases are impossible because |z 1 | ≤ r < 1. Therefore z 2 = 0. Moreover, for u := -Az we get either

∂h(z) ⊂ 0 1 -1 1 -1 0 + 0 1 = 0 2 when z 2 > 0, or ∂h(z) ⊂ 0 1 -1 1 1 0 + 0 -1 = 0 -2 when z 2 < 0.
In any case, η ≥ 2 for each η ∈ ∂h(z) which establishes the claim.

As in the previous case, the claim yields that ξ ≥ 2 whenever ξ ∈ ∂ > h(z). Thus S is calm at ( p, z) with the calmness modulus not exceeding 1.

5. Let z = (0, 0) T and p = (0, -100) T .

Then C z = 0 and v = -100. A simple calculation reveals that the tangent cone T ((0, -100); gph F ) and the normal cone N((0, -100); gph F ) are equal to

R + -1 0 ∪ R + 0 1 and R 0 -1 ∪ R 1 0 ∪ cone 0 -1 , 1 0 ,
respectively. Similarly to the previous cases, one infers that S has neither the Aubin/Lipschitz-like nor the isolated calmness property at ( p, z).

To check the calmness, note that for each z

= (z 1 , z 2 ) T ∈ R 2 we have g(z) = (-z 2 , -100 + z 1 + z 2 , -z 2 , -100 + z 1 + z 2 ) T and h(z) = d (-z 2 , -100 + z 1 + z 2 , -z 2 , -100 + z 1 + z 2 ) T , gph F × {0} × R .
We shall use the same steps as in the case 3. First, find r ∈]0, 1[ such that for each u = (x, y) T ∈ B 2 ( ū, r) with ū := (0, -100) T one has ∂d(

•, gph F )(u) equal to ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ {(1, 0) T }, if x > 0, y ≥ -100, {(0, -1) T }, if x ≤ 0, y < -100, 1 √ x 2 +(y+100) 2 (x, y + 100) T , if x > 0, y < -100, {(0, 1) T }, if x < 0, -100 < y < -100 -x, {(-1, 0) T }, if x < 0, y > -100 -x, and {(-1, 0) T , (0, 1) T }, if x < 0, y = -100 -x.
Second, we claim that given z = (z 1 , z 2 ) T ∈ R 2 and η ∈ R 2 such that h(z) > 0, that p -Az ∈ B 2 ( ū, r), and that η ∈ ∂h(z) one has η ≥ 1. Indeed, (iii) If z 2 < 0 and -100 + z 1 + z 2 < -100, then ∂h(z) ⊂ 1

z 2 2 + (z 1 + z 2 ) 2 0 1 -1 1 -z 2 z 1 + z 2 + 0 -1 = ⎛ ⎝ z 1 +z 2 √ z 2 2 +(z 1 +z 2 ) 2 z 1 +2z 2 √ z 2 2 +(z 1 +z 2 ) 2 -1 ⎞ ⎠ .
(iv) If z 2 > 0 and -100 + z 1 + z 2 < -100, then

∂h(z) ⊂ 0 1 -1 1 0 -1 + 0 1 = -1 0 .
(v) If z 2 > 0 and -100 + z 1 + z 2 > -100, then ∂h(z) ⊂ 0 1 -1 1

-1 0 , 0 1 + 0 1 = 0 2 , 1 2
.
In any case, η ≥ 1 for each η ∈ ∂h(z). Indeed, this needs an explanation in case (iii) only. If z 2 < 0, then z 1 + 2z 2 < z 1 + z 2 < 0. So we may estimate

(z 1 + z 2 ) 2 z 2 2 + (z 1 + z 2 ) 2 + (z 1 + 2z 2 ) 2 z 2 2 + (z 1 + z 2 ) 2 -2 z 1 + 2z 2 z 2 2 + (z 1 + z 2 ) 2 + 1 ≥ 1.
The claim is proved.

As in the case 3, the claim yields that ξ ≥ 1 whenever ξ ∈ ∂ > h(z). Thus S is calm at ( p, z) with the calmness modulus not exceeding 2.

A closer look at the above discussion reveals that any other (z, p) ∈ gph xΦ fits in one of the above categories and that this depends only on the position of the point (C z, v) in the graph of F (see Fig. 6). This is clear if the point in question belongs to the first two categories. In the remaining ones one has to remember that C z = z2 = p1 = 0 and (Cz, B T ( p -Az)) T = p -Az whenever z = (z 1 , z 2 ) T ∈ R 2 . Remark 6.1 Suppose that the diode in the above circuit is replaced with a DIAC. The computation in the Example 6.2 shows that one can use the approach from Example 6.1 to establish the calmness property.

Conclusions

In this paper, we studied the quantitative stability of a generalized equation involving a general set-valued map. The key point in our analysis is the coderivative criterion for the Aubin/Lipschitz-like property, the graphical derivative for the isolated calmness and the outer subdifferential sufficient condition for the calmness. We applied our theoretical results to the theory of non-regular electrical circuits involving electronic devices like ideal diode, practical diode, and DIAC. Some other nonlinear elements like SCR and transistor could be studied in the same way. In order to apply the calculus rules, the injectivity of the matrix B is assumed. It would be interesting to relax this condition. As observed in Example 6.2, the sufficient condition for the calmness is not easy to check even for one or two dimensional problems. It would be interesting to prove a point-based condition which will be easier to check in practice. Another interesting part would be the simulation of the circuits studied in this paper. This is out of the scope of the current paper and will be the subject of another work.

Corollary 3 . 1

 31 Suppose that the assumptions (A1)-(A4) hold true. Then S has the Aubin/Lipschitz-like property at ( p, z) if and only if

Theorem 4 . 1

 41 Under the assumptions (A1)-(A3), S has the isolated calmness property at ( p, z) if and only if

Proposition 5 . 1

 51 Suppose that the assumptions (A1)-(A3) are satisfied. Put Λ = gph F × rge B and define the functions g : R n → R m × R m × R n and h : R n → R + for each z ∈ R n by g(z) := (Cz, (B T B) -1 B T ( pf (z)), pf (z)) and h(z) := d(g(z), Λ). Then S has the calmness property at ( p, z), provided that

Example 5 . 1

 51 Let m = n = 1, B = C = 1.Suppose that f : R → R and F : R ⇒ R are defined as follows:

Fig. 1

 1 Fig. 1 Two equivalent circuits of Example 6.1

Fig. 2 IFig. 3

 23 Fig. 2 I-V characteristic Ψ of a DIAC

Fig. 4

 4 Fig. 4 Calmness modulus from Example 6.1 (Color figure online)

Fig. 5

 5 Fig. 5 Circuit of Example 6.2

Fig. 6 F

 6 Fig. 6 F from Example 6.2 (Color figure online)

2 .

 2 Let z := (0, 1) T and p := (1, -101) T .Hence C z = -1 and v = 100. Moreover, T ((-1, 100); gph F ) = T ((1, 1); gph F ) and N((-1, 100); gph F ) = N((1, 1); gph F ). In view of the previous case, we may conclude that S has the Aubin/Lipschitz-like property and the isolated calmness property at ( p, z).

4 .

 4 Let z := (0, 0) T and p := (0, 0) T . Therefore C z = 0 and v = 0. Moreover, T (0, 0); gph F = R 0 1 and N (0, 0); gph F = R 1 0 .

  h(z) = d (-z 2 , -100 + z 1 + z 2 ) T , gph F + |z 2 | = d( p -Az, gph F ) + |z 2 | and ∂h(z) ⊂ -A T ∂d(•, gph F )(-z 2 , -100 + z 1 + z 2 ) + {0} × ∂| • |(z 2 ). Consider the following cases: (i) If z 2 = 0, then z 1 < 0 because -100 + z 1 /∈ F (0) (the case z 1 > 101 is impossible as p -Az = (0, -100 + z 1 ) T ∈ B 2 ( ū, r)). Thus ∂h(z) If z 2 < 0 and -100 + z 1 + z 2 ≥ -100, then ∂h(z)

Acknowledgements

The authors would like to thank to Huynh Van Ngai for the fruitful discussion concerning the computation of the outer subdifferential in Example 6.2, and also to anonymous referees for their valuable comments. The second author was partially supported by the Research Plan MSM 4977751301.


Mohamed Elawady
email: mohamed.elawady@univ-st-etienne.fr

Cécile Barat

Christophe Ducottet

Philippe Colantoni

Global Bilateral Symmetry Detection Using Multiscale Mirror Histograms

Keywords: Symmetry Detection, Reflection Symmetry, Edge Features, Pairwise Similarity, Symmetry Histogram

In recent years, there has been renewed interest in bilateral symmetry detection in images. It consists in detecting the main bilateral symmetry axis inside artificial or natural images. State-of-the-art methods combine feature point detection, pairwise comparison and voting in Hough-like space. In spite of their good performance, they fail to give reliable results over challenging real-world and artistic images. In this paper, we propose a novel symmetry detection method using multi-scale edge features combined with local orientation histograms. An experimental evaluation is conducted on public datasets plus a new aesthetic-oriented dataset. The results show that our approach outperforms all other concurrent methods.

Introduction

Symmetry is a very important measure in some computer vision applications such as aesthetic analysis [START_REF] Zhao | Exploring principlesof-art features for image emotion recognition[END_REF], object detection and segmentation [START_REF] Teo | Detection and segmentation of 2d curved reflection symmetric structures[END_REF], depth estimation [START_REF] Yang | Depth from water reflection[END_REF], or medical image processing [START_REF] Bairagi | Symmetry-based biomedical image compression[END_REF]. This paper focuses on detecting a single bilateral symmetry axis by exploring the similarity between regions which are mirror reflected around an image axis.

Bilateral symmetry detection has been improved over decades since 1930's [START_REF] Liu | Computational symmetry in computer vision and computer graphics[END_REF]. In particular, Loy and Eklundh [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] approach was a great milestone in the domain, thanks to the use of SIFT descriptor [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], an appropriate pairwise comparison, and Hough-voting scheme [START_REF] Duda | Use of the hough transformation to detect lines and curves in pictures[END_REF]. Since then, many methods [START_REF] Cho | Bilateral symmetry detection via symmetry-growing[END_REF][START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF][START_REF] Kondra | Multi-scale kernel operators for reflection and rotation symmetry: further achievements[END_REF][START_REF] Michaelsen | Recognition of symmetry structure by use of gestalt algebra[END_REF][START_REF] Patraucean | Detection of mirror-symmetric image patches[END_REF][START_REF] Cai | An adaptive symmetry detection algorithm based on local features[END_REF] proposed some improvements and refinements of Loy's method, but still, they fail in cases where SIFT is not adapted (uniform texture, scale variation, etc.). Our idea is to exploit another type of local features, edges, which are the most important visual cues used by the human visual system to detect symmetry. Only few attempts have been done in this direction [START_REF] Ming | Symmetry detection via contour grouping[END_REF][START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF][START_REF] Wang | Reflection symmetry detection using locally affine invariant edge correspondence[END_REF][START_REF] Atadjanov | Bilateral symmetry detection based on scale invariant structure feature[END_REF]. The general framework for symmetry detection proposed by [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] provides promising results, but only edge fragments are used and it fails to account for texture information.

We propose a novel method for symmetry detection preserving the shape and texture information inside an image. This work contributes with: (1) introducing a local edge orientation histogram, (2) utilizing multi-scale edge extraction to exploit the full resolution image. In addition, we introduce groundtruth of a global symmetry axis inside images extracted from the large-scale Aesthetic Visual Analysis (AVA) dataset [START_REF] Murray | Ava: A large-scale database for aesthetic visual analysis[END_REF]. This dataset is composed of artistic photographs which involved in on-line professional photography competitions.

In section 2, we analyze related work of symmetry detection. Then we describe the details of the proposed algorithm in section 3. We present the experimental results on PSU and AVA datasets in section 4 and the conclusion in the last section.

Related Work

The baseline algorithm to detect a global bilateral symmetry axis was introduced by Loy and Eklundh [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF]. The general scheme is as follows: (1) extraction of local-feature points using Lowe's SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], (2) matching of pairs of those keypoints based on the similarity of corresponding mirror SIFT descriptors, and (3) use of the pairs in Hough-voting space to find the best symmetrical candidate. Many attempts proposed refinements to Loy's general scheme, especially Mo and Draper [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] in selecting all symmetry candidate pairs instead of finding closest matches at each point, and using less complex hough voting scheme. The improvement produced more particles than [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] for the symmetry axis, to elongate the output axis according to the global texture information of an image. These intensity-based methods depend mainly on the properties of SIFT features, and their capabilities for symmetry matching with respect to intensity and orientation in the neighboring pixels. In the presence of smooth objects with noisy background, the global symmetry information is lost due to a small number of extracted feature points inside an image.

Cicconet et al. [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] approach exploits the idea of using edge detection in a symmetry measure based on edge orientations. The method can be summarized in the following steps. First, a regular set of feature points with local edge amplitude and orientation is extracted. Second, given any pair of points, a symmetry axis candidate is defined (parametrized by angle and displacement variables) and a symmetry coefficient is associated with this axis. Third, a symmetry voting image is defined as the sum of the weighted contribution of all pairs of points for a given axis. Finally, the best candidate is extracted as the maximum over this voting image. This approach lacks neighborhood's texture information inside the feature representation, and it is totally dependable on the scale parameter of the edge detector. Because of these drawbacks, some local symmetric characteristics are neglected.

Our method investigates Cicconet's edge features [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] within Loy's scheme [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] by adding neighboring-pixel information resulting a better symmetry represen-tation in the voting space. Additionally, it overcomes the scale limitations by adding multi-scale information.

Algorithm Details

The main steps of the proposed symmetry detection algorithm is: 1) dense edge extraction, 2) multi-scale edge description with a local orientation histogram h p as a texture feature, 3) edge triangulation to find symmetry pairs, and 4) voting based on weighted pairs to detect a global bilateral axis.

Local Orientation Histogram

Given an image I to analyze, feature points are sampled along a regular grid (stride and cell size are computed with respect to image size). A feature point p and its local edge characteristics (J p , τ p ) are extracted within each cell using a Morlet wavelet ψ k,σ of constant scale σ and varying orientation {τ k , k = 1 . . . n}. The wavelet parameters are set to detect edges: the wavelength is 4σ, and the elongation along the main orientation is half value of its orthogonal. Formally, The mother Morlet wavelet at point p = (p x , p y ) is defined as:

ψ(p) = e jk0p T e -1 2 |Ap| 2 (1)
where

k 0 = ( π 2 , 0) and A = diag(1, 1
2 ). The wavelet function of scale σ and orientation τ k is:

ψ k,σ (p) = 1 σ ψ( Ω k p σ ) (2)
Ω k = cos(φ k ) sin(φ k ) -sin(φ k ) cos(φ k ) (3)
where φ k is the angle associated with τ k . If J k (p) denote the modulus of wavelet coefficients at point p, local edge characteristics J p and τ p are obtained by seeking the maximum wavelet response and orientation over all orientations. Formally:

J k (p) = |I * ψ k,σ |(p) (4)
k p = arg max k J k (p) (5)
J p = J kp (p) (6)
τ p = τ kp (7)
where * denotes the 2D convolution operation.

Although local edge orientation is important to detect symmetry, the surrounding texture around an edge segment is also an important feature especially for natural images. Thus, we introduce a local orientation histogram h p as the weighted directional contribution of all neighboring points r in the corresponding grid cell. As described in figures (1b-1d), the histogram count at a given orientation τ k is:

h p (k) = r∈N (p) J r δ φ k -φr (8)
where φ k and φ r are angles associated with τ k and τ r , J r and τ r are defined the same way as for p (equations ( 5)-( 7)) and δ x is the Kronecker delta. h p is subsequently 1 normalized and circular shifted so as the first bin corresponds to τ p .

Multi-scale Information

In most images, relevant information about the visual content may appear at different scales. To consider this information, we propose to detect the feature points at multiple scales. For that purpose, we consider a set of regular grids at different scales and a corresponding set of wavelet scales {σ l , l = 1..m}. The whole set of feature points is augmented by extracting wavelet local maxima at each scale and by computing their corresponding local characteristics.

Pairwise Symmetry Coefficient and Mirror Symmetry Histogram

As shown in figure 1a, given two feature points p and q and their corresponding local edge orientations τ p and τ q , T pq is the line passing through p and q, c = (c x , c y ) denotes the center of pq segment and T ⊥ pq represents the candidate symmetry axis defined as the line perpendicular to pq passing through c. -→ o is the reference unit vector of x-axis of the Cartesian plane.

The pairwise symmetry coefficient f (p, q) is defined as [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF]:

f (p, q) = |τ q S(T ⊥ pq )τ p | (9)
where S(T ⊥ pq ) is the reflection matrix with respect to the line perpendicular to (pq):

S(T ⊥ pq ) = -cos(2θ pq ) -sin(2θ pq ) -sin(2θ pq ) cos(2θ pq ) (10)
In addition to the geometric edge information, the symmetry degree of the two regions around p and q can be measured by comparing their corresponding local orientation histogram h p and h q . Due to symmetry considerations, one of the two histograms must be shifted and reversed before comparison. If d I is the histogram intersection distance, and h * q the reversed histogram, the texturebased symmetry measure is given by:

d I (h p , h * q ) = n k=1 min(h p (k), h * q (k)) (11)
Given a pair of feature points (p, q), the candidate axis T ⊥ pq perpendicular to (pq) is parametrized by the orientation of its normal θ pq and its distance to the origin ρ pq = c x cos θ pq + c y sin θ pq where (c x , c y ) are the coordinates of point c located at the middle of (p, q) segment (see figure 1a). A mirror symmetry histogram H S (ρ, θ) is defined as the sum of the contribution of all pairs of feature points such as:

H(c x , c y , θ) = p,q p =q J p J q f (p, q)d I (h p , h q )δ (cx,cy)-p+q 2 δ θ-θpq (12)
H S (ρ, θ) = cx,cy H(c x , c y , θ)δ ρ-ρpq (13)
The best bilateral symmetry axis inside the image is detected by searching for the global maximum of the mirror symmetry histogram H S . For that purpose, we solve the orientation discontinuity problem by extending the voting space in a circular way. The spatial extension of this axis can be determined as the convex hull of the set of feature points contributing to this axis. Figures (2a-2d) show the symmetry detection for some example in terms of two voting histograms.

The first H presents one major axis (A3 -red color) in the middle surrounded by four minor axes (A1, A2, A4, A5 -red color) vertically in figure 2b. In H S histogram (figure 2c) and its zoom-in version (figure 2d), these axes are represented by 5 local maxima aligned along the line of orientation θ = 90 • . Additionally, a horizontal axis (B -green color) appears near the boundaries of θ axis, it describes another symmetry detection (the top half of steel and its mirror in the bottom side). The main axis is the vertical one (A3 -red color) located in the center.

Experiments and Discussions

The aim of this part is to evaluate the results of our method and to compare it to concurrent methods from the state-of-the-art. The evaluation of symmetry detection had received re-surging interest in the main computer vision conferences. Indeed, symmetry detection challenges [START_REF] Rauschert | First symmetry detection competition: Summary and results[END_REF][START_REF] Liu | Symmetry detection from realworld images competition 2013: Summary and results[END_REF] were successively proposed with ECCV2010 1 , CVPR2011 2 and CVPR2013 3 . They all rely on a dataset proposed by Liu's group publicly available with an associated groundtruth and evaluation method. We then use these three general purpose datasets for our evaluation and in addition, we provide an additional aesthetic-oriented dataset. To build this fourth dataset, we extracted some symmetrical images from the Aesthetic Visual Analysis (AVA) collection and generated axis groundtruth among those images.

Implementation and Evaluation Details

We compare our approach against three different methods (Loy2006 [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF], Mo2011 [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF], Cicconet2014 [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF]). We found or requested their source code and ran them using default parameter values, assigned by the authors for stable performance. As described in [START_REF] Liu | Symmetry detection from realworld images competition 2013: Summary and results[END_REF], a global symmetry detection is considered to be correct if it satisfies the following conditions: [START_REF] Atadjanov | Bilateral symmetry detection based on scale invariant structure feature[END_REF] The angle between the detected symmetry axis and the groundtruth axis is less than 10 degrees. [START_REF] Bairagi | Symmetry-based biomedical image compression[END_REF] The distance between centers of these two axes is less than 0.2 multiplied by minimum length of the axes. Accuracy represents the correct detection rate (the number of correct detection images divided by the total number of dataset images).

Comparison on PSU Datasets (2010, 2011 and 2013)

From Flickr photo sharing website, Liu's vision group provides a single-symmetry groundtruth for synthetic/real images (size of 2010, 2011 and 2013 datasets are 134, 14 and 75 respectively). The quantitative comparisons are presented in Figure 3 with respect to each version of PSU dataset. Our method is ranked first for all PSU version with a difference of more than 10% among the nearest Fig. 3: Comparison of accuracy rates for our method "Our2016" against the baseline algorithm "Loy2006" [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF], its improvement "Mo2011" [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] and one of the recent algorithms "Cic2014" [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF]. The used datasets are: different versions of PSU (2010, 2011, 2013) with provided groundtruth [START_REF] Rauschert | First symmetry detection competition: Summary and results[END_REF][START_REF] Liu | Symmetry detection from realworld images competition 2013: Summary and results[END_REF] and AVA with our groundtruth.

competent Loy2006 [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF], which achieves the best detection accuracy among 2011 and 2013 symmetry challenges [START_REF] Rauschert | First symmetry detection competition: Summary and results[END_REF][START_REF] Liu | Symmetry detection from realworld images competition 2013: Summary and results[END_REF]. Cic2014 [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] approach has ranked third due to the absence of local neighborhood information in the mirror symmetry histogram for candidate selection. Figures (4a-4f) gives some detection examples 4 to present the superior performance of the proposed algorithm over others. In figures (4a,4d), our method correctly detects the symmetry axis of the bird with an extra line fitting on axis endpoints. However, [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] finds very few SIFT features generating a tiny non-meaningful detection, and [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF][START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] output misalignment axes representing different components within an image. Figure 4b shows a valid axis detection by our method, while [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF][START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] attempt to find an inaccurate global axis due to presence of the non-exact shapes inside the centered object. Inside figure 4c with different changes in environment illumination, our method perfectly aligned with the full symmetry axis and [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF][START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] select different secondary global axes due to surface's mirror effect of the glass. A synthetic butterfly of figure 4e represents an easy detection of all methods except [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF]. Our method has the most accurate axis output compared to the axis groundtruth. With a camera flash-on, facial figure 4f results a challenging axis detection, in which [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] shortly aligns over the groundtruth according to insufficient SIFT features. But our method gives the most appropriate axis output with a horizontal shift, such that it covers all human parts (including frontal face, neck, and hair). In all figures except 4e, [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF]

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n)
Fig. 4: Some images in (a-f) PSU datasets and (g-n) our labeled AVA dataset with groundtruth (blue), where our method (cyan) produces better or similar results over Loy2006 [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] (green), Mo2011 [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] (magenta) and Cic2014 [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] (yellow). Each axis is shown in a straight line with squared endpoints and a starred midpoint. Best seen on screen.

is very sensitive to noise in limited contrast environments and in high-texture foreground/background objects.

Comparison on AVA Dataset

From DPChallenge photo contest website, Murray et al. [START_REF] Murray | Ava: A large-scale database for aesthetic visual analysis[END_REF] introduces different annotations5 (aesthetic, semantic and photographic style) for more than 250,000 images for Aesthetic Visual Analysis "AVA". From the following photography challenges, we labeled global-axis symmetry groundtruth6 for 253 out of 878 images: (1) five challenges of "Reflections Without Mirrors": images containing bilateral representation without using mirror, (2) three challenges of "Symmetry": photographs composing symmetrical balance. These images are selected to neglect unclear instances of ambiguity symmetry, and to represent many comparison cases (non-centering viewpoint, perspective view, blurring reflection, etc.) for detection algorithms. Figure 3 shows that our approach has better performance over other methods. Figures (4g-4n) presents some experimental results7 representing different types of artistic scene. Natural landscapes in figures (4g,4k) illustrate that [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] behaves abnormally with perspective view conditions, and [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF][START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] act incomplete with high shape occlusions. Abstract figures (4h,4l) demonstrate the superior performance of [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] compared to [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF], in order to fit the global symmetry axis. [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] works better among less texture objects with clear contours. In figures (4i,4j), our method easily detects global symmetry axis, while [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF][START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] stuck with the local ones if exists. Our method only finds a clear symmetry detection in figure 4m with a high-reflected object with non-homogeneous background. Contourbased approaches, as ours and [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] ones, execute successfully within the nonexact symmetry figure 4n. In summary, our method and [START_REF] Cicconet | Mirror symmetry histograms for capturing geometric properties in images[END_REF] ones align better in figures containing objects with clear edges, [START_REF] Loy | Detecting symmetry and symmetric constellations of features[END_REF] locally detects shapes with smooth background, and [START_REF] Mo | Detecting bilateral symmetry with feature mirroring[END_REF] provides incorrect and unclear symmetry over the most visual cases.

Conclusion

This paper introduces multi-scale edge features combined with local orientation histograms, in order to develop a reliable global symmetry detection among variants of visual cues. Quantitative and qualitative comparisons show a substantial advantage for our proposed method on different types of public datasets. Our model can be improved to handle complex images with non-centered viewpoints and large degrees of perspective view. It can be extended to solve the over-fitting of axis endpoints. The future work is to introduce a stable balance measure; describing the existence and degree of global axes in terms of both geometrical and symmetrical properties inside an image, and to integrate this measure within existing retrieval systems of visual arts.

Fig. 1 :

 1 Fig. 1: Pairwise mirror symmetry and local orientation histogram. (a) Pairwise mirror symmetry. (b) Edge segment distribution within a cell. (c) Angular distribution of edge segments for 16 main orientations over 180 degrees. (d) Normalized local orientation histogram. Best seen in screen.

Fig. 2 :

 2 Fig. 2: Symmetry detection process: (a) Input image with local symmetry axis candidates indicated with red and green lines. (b) The output of the voting histogram H(c x , c y ) marginalized over the orientation variable θ and the corresponding symmetry axis IDs. (c) The output of the voting histogram H S (ρ, θ). And (d) the zoomed-in region inside the red rectangle. Local maxima correspond to candidate symmetry axes (framed in red and green).

All PSU output with detected symmetry axes can be found in: http://perso.univst-etienne.fr/em68594h/SupplementalFilesPSU.zip

http://www.lucamarchesotti.com/

Source code to generate AVA images and their symmetry labels: http://perso.univst-etienne.fr/em68594h/SymAVA.zip

[START_REF] Kondra | Multi-scale kernel operators for reflection and rotation symmetry: further achievements[END_REF] All AVA output with the detected symmetry axes can be found in: http://perso.univst-etienne.fr/em68594h/SupplementalFilesAVA.zip