Franz Chouly
email: franz.chouly@univ-fcomte.fr

Patrick Hild
email: patrick.hild@math.univ-toulouse.fr

A Nitsche-based method for unilateral contact problems: numerical analysis

Keywords: unilateral contact, finite elements, Nitsche's method. AMS Subject Classification: 65N12, 65N30, 74M15

We introduce a Nitsche-based formulation for the finite element discretization of the unilateral contact problem in linear elasticity. It features a weak treatment of the non-linear contact conditions through a consistent penalty term. Without any additional assumption on the contact set, we can prove theoretically its fully optimal convergence rate in the H 1 (Ω)norm for linear finite elements in two dimensions, which is O(h 1 2 +ν ) when the solution lies in H 3 2 +ν (Ω), 0 < ν ≤ 1/2. An interest of the formulation is that, conversely to Lagrange multiplier-based methods, no other unknown is introduced and no discrete inf-sup condition needs to be satisfied.

Introduction

The numerical implementation of contact and impact problems in solid mechanics generally uses the Finite Element Method (FEM) (see [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Han | Quasistatic contact problems in viscoelasticity and viscoplasticity[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]). These problems involve non-linear contact conditions on a part of the boundary that lead naturally to a variational inequality (see, e.g., [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF]). In this paper, we limit ourselves to the case of frictionless unilateral contact in plane linear elasticity. Indeed, this case contains the key difficulties related to the non-linear contact conditions. We also limit our study to the approximation with continuous piecewise linear finite elements, which is the most standard approach. We introduce a special FEM inspired from Nitsche method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF][START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF][START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF][START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF]. Conversely to standard penalization techniques [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF], the resulting method is consistent. Moreover, unlike mixed methods (see, e.g., [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF]), no additional unknown (Lagrange multiplier) is needed. However, the most interesting in this method is that we are able to prove its optimal convergence in the H 1 (Ω)-norm of order O(h 1 2 +ν ) provided the solution has a regularity H 3 2 +ν (Ω), 0 < ν ≤ 1/2. 1 Furthermore, we do not need any additional assumption on the contact zone, such as an increased regularity of the contact stress or a finite number of transitions between contact and non-contact. Besides, the standard FEM for contact consists in a direct conforming approximation of the variational inequality, with the elastic displacement as the only unknown. For this standard FEM and also for all the other approaches such as mixed/hybrid methods (e.g., [START_REF] Belgacem | Hybrid finite element methods for the Signorini problem[END_REF]), stabilized mixed methods (e.g., [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]), penalty methods (e.g., [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF]), no such proof of optimal convergence has been established to the best of our knowledge in the case the solution u is in H 3 2 +ν (Ω) (0 < ν ≤ 1/2). The first analyzes in [START_REF] Scarpini | Error estimates for the approximation of some unilateral problems[END_REF][START_REF] Haslinger | Finite element analysis for unilateral problems with obstacles on the boundary[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF] were sub-optimal with a convergence in O(h

1 2 + ν
2 ) and these results were recently improved in [START_REF]An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF] to obtain a convergence rate of O(h 1 2 + ν 2 +ν 2 ) when (0 < ν < 1/2), of O(h | ln(h)|) when ν = 1/2 and in [START_REF] Renard | A quasi-optimal a priori error estimate for the two-dimensional Signorini problem approximated by linear finite elements[END_REF] where a convergence of O(h) is obtained when ν = 1/2 + ε (ε > 0 arbitrarily small). In fact, additional assumptions on the finiteness of transition points between contact and non-contact are needed to recover optimality (see [START_REF] Belgacem | Numerical simulation of some variational inequalities arisen from unilateral contact problems by the finite element methods[END_REF] for the case 0 < ν < 1/2 and [START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF] for the case ν = 1/2 in the case of the standard FEM). Our paper is outlined as follows. In Section 2, we recall the continuous (strong and weak) formulations for unilateral contact problems and introduce our Nitsche-based FEM. In Section 3, we carry out the numerical analysis of this method: we prove its consistency, the existence and uniqueness of solutions, and at least its optimal convergence. Conclusion and perspectives are drawn in Section 4. Let us introduce some useful notations. In what follows, bold letters like u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K . . .) represent functional sets involving vector fields. As usual, we denote by (H s (.)) d , s ∈ R, d = 1, 2 the Sobolev spaces in one or two space dimensions (see [START_REF] Adams | Sobolev spaces[END_REF]). The usual norm of (H s (D)) d is denoted by • s,D and we keep the same notation when d = 1 or d = 2. The letter C stands for a generic constant, independent of the discretization parameters.

Setting

The unilateral contact problem

We consider an elastic body Ω in R 2 where plane small strain assumptions are made. The boundary ∂Ω of Ω is polygonal and we suppose that ∂Ω consists in three nonoverlapping parts Γ D , Γ N and the contact boundary Γ C , with meas(Γ D ) > 0 and meas(Γ C ) > 0. The contact boundary is supposed to be a straight line segment to simplify ; the normal unit outward vector on ∂Ω is denoted n. In its initial stage, the body is in contact on Γ C with a rigid foundation (the extension to two elastic bodies in contact can be easily made, at least for small strain models) and we suppose that the unknown final contact zone after deformation will be included into Γ C . The body is clamped on Γ D for the sake of simplicity. It is subjected to volume forces f ∈ (L 2 (Ω)) 2 and to surface loads g ∈ (L 2 (Γ N )) 2 . The unilateral contact problem in linear elasticity consists in finding the displacement field u : Ω → R 2 verifying the equations and conditions (1)-( 2):

div σ(u) + f = 0 in Ω, σ(u) = A ε(u) in Ω, u = 0 on Γ D , σ(u)n = g on Γ N , (1)
where σ = (σ ij ), 1 ≤ i, j ≤ 2, stands for the stress tensor field and div denotes the divergence operator of tensor valued functions. The notation ε(v) = (∇v+∇v T )/2 represents the linearized strain tensor field and A is the fourth order symmetric elasticity tensor having the usual uniform ellipticity and boundedness property. For any displacement field v and for any density of surface forces σ(v)n defined on ∂Ω we adopt the following notation

v = v n n + v t and σ(v)n = σ n (v)n + σ t (v),
where v t (resp. σ t (v)) are the tangential components of v (resp. σ(v)n). The conditions describing unilateral contact without friction on Γ C are:

u n ≤ 0 (i) σ n (u) ≤ 0 (ii) σ n (u) u n = 0 (iii) σ t (u) = 0 (iv) (2)
We introduce the Hilbert space V and the convex cone K of admissible displacements which satisfy the noninterpenetration on the contact zone Γ C :

V := v ∈ H 1 (Ω) 2 : v = 0 on Γ D , K := {v ∈ V : v n = v • n ≤ 0 on Γ C } . Define a(u, v) := Ω σ(u) : ε(v) dΩ, L(v) := Ω f • v dΩ + Γ N g • v dΓ,
for any u and v in V. From the previous assumptions, we know that a(•, •) is bilinear, symmetric, V-elliptic and continuous on V × V. We know also that L(•) is a continuous linear form on V.

The weak formulation of Problem (1)-( 2), as a variational inequality (see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]), is:

Find u ∈ K such that: a(u, v -u) ≥ L(v -u), ∀ v ∈ K. (3)
Stampacchia's Theorem ensures that Problem (3) admits a unique solution.

The Nitsche-based finite element method

Let V h ⊂ V be a family of finite dimensional vector spaces (see [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF]) indexed by h coming from a family T h of triangulations of the domain Ω (h = max T ∈T h h T where h T is the diameter of T ). The family of triangulations is supposed regular, i.e., there exists σ > 0 such that ∀T ∈ T h , h T /ρ T ≤ σ where ρ T denotes the radius of the inscribed circle in T . We choose standard continuous and piecewise affine functions, i.e.:

V h = v h ∈ (C(Ω)) 2 : v h |T ∈ (P 1 (T )) 2 , ∀T ∈ T h , v h = 0 on Γ D .
Let us introduce the notation [•] + for the positive part of a scalar quantity a ∈ R:

[a] + = a if a > 0, 0 otherwise.
In the rest of this paper, we will make an extensive use of the following properties:

a ≤ [a] + , a[a] + = [a] 2 + , ∀a ∈ R. (4)
Note that conditions (4) can be straightforwardly extended to real valued functions. The derivation of a Nitsche-based method comes from the observation that the contact conditions (2) can be reformulated as follows:

Proposition 2.1. Let γ > 0. The contact conditions (2) (i)-(iii) on Γ C are equivalent to:

σ n (u) = - 1 γ [u n -γσ n (u)] + . (5)
Proof: Let u be a regular vector field on Ω such that (2) (i)-(iii) holds. The condition ( 2)-(ii) yields either σ n (u) < 0 or σ n (u) = 0. Suppose first that σ n (u) < 0. Then ( 2)-(iii) implies that u n = 0. In this case, it holds:

- 1 γ [u n -γσ n (u)] + = - 1 γ [-γσ n (u)] + = σ n (u).
Then, suppose that σ n (u) = 0. The condition ( 2)-(i) can also be expressed as [u n ] + = 0, so:

- 1 γ [u n -γσ n (u)] + = - 1 γ [u n ] + = 0 = σ n (u).
Reciprocally, let u such that (5) holds. Note first that it implies σ n (u) ≤ 0, so that (2)-(ii) holds. Take first the case σ n (u) = 0, then (5) can be rewritten as:

0 = - 1 γ [u n ] + , which is equivalent to condition (2)-(i). Since σ n (u) = 0, then condition (2)-(iii) (σ n (u)u n = 0) also holds.
We now consider the case σ n (u) < 0. From ( 5), [u n -γσ n (u)] + > 0, so that in this case

-γσ n (u) = [u n -γσ n (u)] + = u n -γσ n (u),
from which comes u n = 0, so that both (2) (i) and (iii) hold.

Remark 2.2. Note that condition ( 5) is still equivalent to (2) (i)-(iii) on Γ C when γ is a positive function defined on Γ C instead of a positive constant.

Let now u be the solution of the unilateral contact problem in its strong form ( 1)-( 2), sufficiently regular so that all the following calculations make sense. From the Green formula and equations ( 1), ( 2)-(iv), it holds for every v ∈ V:

a(u, v) - Γ C σ n (u) v n dΓ = L(v). Note that v n = (v n -γσ n (v)) + γσ n (v), so that: a(u, v) - Γ C γ σ n (u) σ n (v) dΓ - Γ C σ n (u) (v n -γσ n (v)) dΓ = L(v).
Finally, using condition (5), we obtain:

a(u, v) - Γ C γ σ n (u) σ n (v) dΓ + Γ C 1 γ [u n -γσ n (u)] + (v n -γσ n (v)) dΓ = L(v). (6)
Formula ( 6) is the starting point of our Nitsche-based formulation. Remark that it may have no sense at the continuous level if u lacks of regularity (u ∈ V only is not sufficient for instance to justify the above calculations). Nevertheless, and as in the stabilized Lagrange multiplier method [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF], we consider in what follows that γ is a positive piecewise constant function on the contact interface Γ C : for any x ∈ Γ C , let T be an element such that x ∈ T and set:

γ(x) = γ 0 h T (7)
where γ 0 is a positive constant. This allows to define a discrete counterpart of (6). Let us introduce for this purpose the discrete linear operator

P γ : V h → L 2 (Γ C ) v h → v h n -γ σ n (v h )
, and also the bilinear form:

A γ (u h , v h ) = a(u h , v h ) - Γ C γ σ n (u h )σ n (v h ) dΓ.
Remark 2.3. With the previous notations, we see that problem (3) could be formally written as follows :

    
Find a sufficiently regular u ∈ V such that:

A γ (u, v) + Γ C 1 γ [P γ (u)] + P γ (v) dΓ = L(v)
, for all sufficiently regular v ∈ V.

Our Nitsche-based method then reads:

     Find u h ∈ V h such that: A γ (u h , v h ) + Γ C 1 γ [P γ (u h )] + P γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h . ( 8
)

A reinterpretation of the method

The aim of this subsection is to show some links between the proposed method and the stabilized mixed/hybrid finite element methods (which do not require inf-sup conditions) introduced in [START_REF] Hughes | A new finite element formulation for computational fluid dynamics. VII. The Stokes problem with various well-posed boundary conditions: symmetric formulations that converge for all velocity/pressure spaces[END_REF]. We define the auxiliary variable

λ h := - 1 γ [P γ (u h )] + ,
which can be interpreted as the L 2 -projection of the quantity -1 γ P γ (u h ) on the convex cone:

L 2 -(Γ C ) := {µ ∈ L 2 (Γ C ) | µ ≤ 0 a.e. on Γ C }. This means in particular that λ h ∈ L 2 -(Γ C ) verifies the following inequality, for all µ ∈ L 2 -(Γ C ): Γ C (µ -λ h )(- 1 γ P γ (u h ) -λ h ) dΓ ≤ 0, (9)
that is the standard characterization of the projection onto a closed convex set. We then use the definition of the operator P γ to obtain

Γ C (µ -λ h )(- 1 γ P γ (u h ) -λ h ) dΓ = Γ C (µ -λ h )(- 1 γ u h n + σ n (u h ) -λ h ) dΓ = - Γ C 1 γ (µ -λ h )u h n + γ(µ -λ h )(λ h -σ n (u h )) dΓ.
Since γ is a positive piecewise constant function on Γ C (see [START_REF] Belgacem | Hybrid finite element methods for the Signorini problem[END_REF]), the above equalities prove that the inequality ( 9) is equivalent to:

Γ C (µ -λ h )u h n dΓ + Γ C γ(µ -λ h )(λ h -σ n (u h )) dΓ ≥ 0. ( 10
)
The next step is to rewrite the left hand-side of the Nitsche-based method (8) as follows:

A γ (u h , v h ) + Γ C 1 γ [P γ (u h )] + P γ (v h ) dΓ = a(u h , v h ) - Γ C γ σ n (u h )σ n (v h ) dΓ - Γ C λ h (v h n -γ σ n (v h )) dΓ = a(u h , v h ) - Γ C λ h v h n dΓ + Γ C γ(λ h -σ n (u h ))σ n (v h ) dΓ.
We combine this last result to the inequality [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF] to obtain an equivalent formulation of the Nitsche-based method:

             Find (u h , λ h ) ∈ V h × L 2 -(Γ C ) such that: a(u h , v h ) - Γ C λ h v h n dΓ + Γ C γ(λ h -σ n (u h )) σ n (v h ) dΓ = L(v h ), ∀ v h ∈ V h , Γ C (µ -λ h )u h n dΓ + Γ C γ(µ -λ h )(λ h -σ n (u h )) dΓ ≥ 0, ∀ µ ∈ L 2 -(Γ C ).
We finally observe that the Nitsche-based method can be regarded as a mixed method with a stabilization term (see [START_REF] Barbosa | Circumventing the Babuška-Brezzi condition in mixed finite element approximations of elliptic variational inequalities[END_REF]). This makes it close to the stabilized method proposed and analyzed in [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF], the difference being that in our case, the space for the Lagrange multiplier λ h is L 2 -(Γ C ). A similar analogy has already been noticed in [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF] in the case of elliptic problems with Dirichlet boundary conditions.

Analysis of the Nitsche-based method

In this section, we carry out the mathematical analysis of the method [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF]. A difference between Nitsche and classical penalty methods [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF] is that this former is consistent, which we first show in §3.1. The proof of well-posedness of the (non-linear) discrete problem ( 8) is carried out in §3.2. The error analysis is finally detailed in §3.3. We show that the method converges optimally for a fixed value of γ 0 and when the mesh size h vanishes.

Consistency

As well as the Nitsche's method for second order elliptic problems with Dirichlet boundary condition or domain decomposition [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF], our Nitsche-based formulation (8) for unilateral contact is consistent: Lemma 3.1. The Nitsche-based method for contact is consistent: suppose that the solution u of (1)-( 2) is in (H 3 2 +ν (Ω)) 2 , with ν > 0, then u is also solution of

A γ (u, v h ) + Γ C 1 γ [P γ (u)] + P γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h .
Proof: Let u be the solution of ( 1)-( 2) and take

v h ∈ V h . Since u ∈ (H 3 2 +ν (Ω)) 2 and ν > 0, we have σ n (u) ∈ H ν (Γ C ) ⊂ L 2 (Γ C
). As a result, P γ (u) ∈ L 2 (Γ C ) and A γ (u, v h ) makes sense. On the one hand, we use the definition of P γ , of A γ (•, •) and the reformulation (5) of contact conditions to obtain:

A γ (u, v h ) + Γ C 1 γ [P γ (u)] + P γ (v h ) dΓ = a(u, v h ) - Γ C γ σ n (u)σ n (v h ) dΓ + Γ C 1 γ (-γσ n (u))(v h n -γ σ n (v h )) dΓ = a(u, v h ) - Γ C σ n (u)v h n dΓ.
On the other hand, with equations ( 1)-( 2) and integration by parts, it holds:

a(u, v h ) - Γ C σ n (u)v h n dΓ = L(v h ),
which ends the proof.

Well-posedness

To prove well-posedness of our Nitsche-based formulation, we need first the following classical lemma that we recall here for convenience of the reader.

Lemma 3.2. Let γ 0 be sufficiently small. Then the bilinear form A γ (•, •) is elliptic on V h : there exists C, independent of the parameter γ 0 and of the mesh size h, such that:

A γ (v h , v h ) ≥ C v h 2 1,Ω , for all v h ∈ V h .
Proof: Note that we can suppose without loss of generality that Γ C is a straight line segment parallel to the x-axis. Let E be an edge of a triangle on Γ C and let T ∈ T h be the element containing E. Consequently we deduce, for any v h ∈ V h :

σ n (v h ) 0,E = σ yy (v h ) 0,E = |E| 1 2 |T | 1 2 σ yy (v h ) 0,T ≤ Ch -1 2 T σ yy (v h ) 0,T = C γ γ 0 -1 2 σ yy (v h ) 0,T .
By summation on all the edges E ⊂ Γ C we get

γ 1 2 σ n (v h ) 2 0,Γ C ≤ Cγ 0 σ yy (v h ) 2 0,Ω ≤ Cγ 0 v h 2 1,Ω . (11)
Hence, from Korn inequality and (11), when γ 0 is small enough, there exists C > 0 such that for any v h ∈ V h :

a(v h , v h ) - Γ C γ(σ n (v h )) 2 dΓ ≥ C v h 2 1,Ω .
We then show that Problem ( 8) is well-posed using an argument proposed by Brezis for M type and pseudo-monotone operators [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] (see also [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] and [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF]):

Theorem 3.3. Suppose that γ 0 is sufficiently small. Then the problem (8) admits one unique solution u h in V h .

Proof: Using the Riesz representation theorem, we define a (nonlinear) operator B h : V h → V h , with help of the formula:

(B h v h , w h ) 1,Ω = A γ (v h , w h ) + Γ C 1 γ [P γ (v h )] + P γ (w h ) dΓ,
for all v h , w h ∈ V h , and where (•, •) 1,Ω stands for the scalar product in (H 1 (Ω)) 2 . Note that Problem ( 8) is well-posed if and only if B is a one-to-one operator. Due to the properties (4), we observe that, for all a, b ∈ R:

([a] + -[b] + )(a -b) = a[a] + + b[b] + -b[a] + -a[b] + ≥ [a] 2 + + [b] 2 + -2[a] + [b] + = ([a] + -[b] + ) 2 ≥ 0.
Let v h , w h ∈ V h , using the above property and then Lemma 3.2, it holds:

(B h v h -B h w h , v h -w h ) 1,Ω = A γ (v h -w h , v h -w h ) + Γ C 1 γ ([P γ (v h )] + -[P γ (w h )] + )P γ (v h -w h ) dΓ ≥ A γ (v h -w h , v h -w h ) ≥ C v h -w h 2 1,Ω . (12)
Next, let us show that B h is also hemicontinuous. Since V h is a vector space, it is sufficient to show that

[0, 1] t → ϕ(t) = (B h (v h -tw h ), w h ) 1,Ω ∈ R is a continuous real function, for all v h , w h ∈ V h . Let s, t ∈ [0, 1]
, we bound:

|ϕ(t) -ϕ(s)| = |(B h (v h -tw h ) -B h (v h -sw h ), w h ) 1,Ω | = A γ ((s -t)w h , w h ) + Γ C 1 γ [P γ (v h -tw h )] + -[P γ (v h -sw h )] + P γ (w h ) dΓ ≤ |s -t| A γ (w h , w h ) + Γ C 1 γ [P γ (v h -tw h )] + -[P γ (v h -sw h )] + |P γ (w h )| dΓ.
With help of the bound

|[a] + -[b] + | ≤ |a -b|
, for all a, b ∈ R, and using the linearity of P γ , we then remark that:

Γ C 1 γ [P γ (v h -tw h )] + -[P γ (v h -sw h )] + |P γ (w h )| dΓ ≤ Γ C 1 γ P γ (v h -tw h ) -P γ (v h -sw h ) |P γ (w h )| dΓ = Γ C 1 γ |(s -t)P γ (w h )||P γ (w h )| dΓ = |s -t| Γ C 1 γ P γ (w h ) 2 dΓ.
It results that:

|ϕ(t) -ϕ(s)| ≤ |s -t| A γ (w h , w h ) + γ -1 2 P γ (w h ) 2 0,Γ C ,
which means that ϕ is Lipschitz, so that B h is hemicontinuous. Since property ( 12) also holds, we finally apply the Corollary 15 (p.126) of [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] to conclude that B h is a one-to-one operator. This ends the proof.

Error analysis

Our Nitsche-based method (8) converges optimally when the mesh parameter h vanishes. This is the object of the following theorems. First, we establish an abstract error estimate.

Theorem 3.4. Suppose that the solution u of Problem (3) belongs to (H 3 2 +ν (Ω)) 2 with ν > 0. Suppose also that the parameter γ 0 is sufficiently small. The solution u h of Problem (8) satisfies the following error estimate:

u -u h 1,Ω + γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C ≤ C inf v h ∈V h ( u -v h 1,Ω + γ -1 2 (u n -v h n ) 0,Γ C + γ 1 2 σ n (u -v h ) 0,Γ C , (13)
with C > 0 a constant, independent of h and u.

Proof: Let v h ∈ V h . We first use the V-ellipticity and the continuity of a(•, •), as well as Young's inequality, to obtain:

α u -u h 2 1,Ω ≤ a(u -u h , u -u h ) = a(u -u h , (u -v h ) + (v h -u h )) ≤ C u -u h 1,Ω u -v h 1,Ω + a(u -u h , v h -u h ) ≤ α 2 u -u h 2 1,Ω + C 2 2α u -v h 2 1,Ω + a(u, v h -u h ) -a(u h , v h -u h ),
with α > 0 the ellipticity constant. We can transform the term a(u,

v h -u h ) -a(u h , v h -u h )
since u solves (3) and u h solves (8), which yields:

α 2 u -u h 2 1,Ω ≤ C 2 2α u -v h 2 1,Ω + Γ C σ n (u) v h n -u h n dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ + Γ C 1 γ [P γ (u h )] + P γ (v h -u h ) dΓ = C 2 2α u -v h 2 1,Ω + Γ C σ n (u) v h n -u h n dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ + Γ C 1 γ [P γ (u h )] + v h n -γσ n (v h ) dΓ - Γ C 1 γ [P γ (u h )] 2 + dΓ, (14)
where we use [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF]. We now write

Γ C σ n (u) v h n -u h n dΓ + Γ C 1 γ [P γ (u h )] + v h n dΓ = Γ C σ n (u) + 1 γ [P γ (u h )] + v h n -u n dΓ + Γ C σ n (u) u n -u h n dΓ + Γ C 1 γ [P γ (u h )] + u n dΓ ≤ Γ C σ n (u) + 1 γ [P γ (u h )] + v h n -u n dΓ + Γ C σ n (u) u n -u h n dΓ,
because of contact condition (2)-(i). We use this last inequality into the global estimate ( 14) to obtain:

α 2 u -u h 2 1,Ω ≤ C 2 2α u -v h 2 1,Ω + Γ C σ n (u) + 1 γ [P γ (u h )] + v h n -u n dΓ + Γ C σ n (u) u n -u h n dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ - Γ C 1 γ [P γ (u h )] + γσ n (v h ) dΓ - Γ C 1 γ [P γ (u h )] 2 + dΓ = C 2 2α u -v h 2 1,Ω + Γ C σ n (u) + 1 γ [P γ (u h )] + v h n -u n dΓ + T . (15)
The integral term in ( 15) is bounded thanks to Cauchy-Schwarz and Young inequalities:

Γ C σ n (u) + 1 γ [P γ (u h )] + v h n -u n dΓ ≤ γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C γ -1 2 (v h n -u n ) 0,Γ C ≤ 1 2β 1 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + β 1 2 γ -1 2 (v h n -u n ) 2 0,Γ C , (16)
with β 1 > 0. We then estimate the (four) integrals in (15) coming from T . Using (2)-(iii), we get

T = - Γ C σ n (u)u h n dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ - Γ C 1 γ [P γ (u h )] + γσ n (v h ) dΓ - Γ C 1 γ [P γ (u h )] 2 + dΓ = - Γ C σ n (u)(u h n -γσ n (u h )) dΓ - Γ C γσ n (u)σ n (u h ) dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ - Γ C 1 γ [P γ (u h )] + γσ n (v h ) dΓ - Γ C 1 γ [P γ (u h )] 2 + dΓ.
Once again we use (4) and contact condition ( 2)-(ii) to note that:

- Γ C σ n (u)(u h n -γσ n (u h )) dΓ = - Γ C σ n (u)P γ (u h ) dΓ ≤ - Γ C σ n (u)[P γ (u h )] + dΓ.
It results into:

T ≤ - Γ C σ n (u)[P γ (u h )] + dΓ - Γ C γσ n (u)σ n (u h ) dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ - Γ C 1 γ [P γ (u h )] + γσ n (v h ) dΓ - Γ C 1 γ [P γ (u h )] 2 + dΓ = - Γ C σ n (u) + 1 γ [P γ (u h )] + [P γ (u h )] + dΓ - Γ C γσ n (u)σ n (u h ) dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ - Γ C 1 γ [P γ (u h )] + + σ n (u) γσ n (v h ) dΓ + Γ C γσ n (u)σ n (v h ) dΓ = - Γ C γ σ n (u) + 1 γ [P γ (u h )] + 1 γ [P γ (u h )] + + σ n (u) dΓ - Γ C γσ n (u)σ n (u h ) dΓ - Γ C γ σ n (u h )σ n (v h -u h ) dΓ + Γ C γ σ n (u) + 1 γ [P γ (u h )] + σ n (u -v h ) dΓ + Γ C γσ n (u)σ n (v h ) dΓ.
As a consequence

T + γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C ≤ Γ C γ σ n (u) + 1 γ [P γ (u h )] + σ n (u -v h ) dΓ + Γ C γσ n (u -u h )σ n (v h -u h ) dΓ. (17)
With help of Cauchy-Schwarz and Young inequalities, we bound the first term in [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF] :

Γ C γ σ n (u) + 1 γ [P γ (u h )] + σ n (u -v h ) dΓ ≤ 1 2β 2 γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + β 2 2 γ 1 2 σ n (u -v h ) 2 0,Γ C , (18)
with β 2 > 0.

The second term in ( 17) is bounded as follows by using [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF] and Young inequality:

Γ C γσ n (u -u h )σ n (v h -u h ) dΓ ≤ γ 1 2 σ n (u -u h ) 0,Γ C γ 1 2 σ n (v h -u h ) 0,Γ C ≤ Cγ 1 2 0 v h -u h 1,Ω γ 1 2 σ n (u -v h ) 0,Γ C + γ 1 2 σ n (v h -u h ) 0,Γ C ≤ C γ 0 v h -u h 2 1,Ω + γ 1 2 σ n (u -v h ) 2 0,Γ C ≤ C γ 0 u -u h 2 1,Ω + γ 0 u -v h 2 1,Ω + γ 1 2 σ n (u -v h ) 2 0,Γ C . (19)
Putting together ( 18) and ( 19) in ( 17) gives:

T + (1 - 1 2β 2 ) γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C ≤C γ 0 u -u h 2 1,Ω + γ 0 u -v h 2 1,Ω + (C + β 2 2 ) γ 1 2 σ n (u -v h ) 2 0,Γ C . (20)
We finally combine estimates ( 16) and ( 20) into (15) to get:

( α 2 -Cγ 0 ) u -u h 2 1,Ω + (1 - 1 2β 1 - 1 2β 2 ) γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C ≤( C 2 2α + Cγ 0 ) u -v h 2 1,Ω + β 1 2 γ -1 2 (u n -v h n ) 2 0,Γ C + (C + β 2 2 ) γ 1 2 σ n (u -v h ) 2 0,Γ C .
To conclude, take γ 0 sufficiently small and for instance

β 1 = β 2 = 2.
The optimal convergence of the method is stated below.

Theorem 3.5. Suppose that the solution u of Problem (3) belongs to (H 3 2 +ν (Ω)) 2 with 0 < ν ≤ 1/2. Suppose also that the parameter γ 0 is sufficiently small. The solution u h of Problem (8) satisfies the following error estimate:

u -u h 1,Ω + γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (21)
with C > 0 a constant, independent of h and u.

Proof: We need to bound the right terms in estimate (13) and we choose v h = I h u where I h stands for the Lagrange interpolation operator mapping onto V h . The estimation of the Lagrange interpolation error in H 1 -norm on a domain Ω is classical (see, e.g., [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF]):

u -I h u 1,Ω ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (22)
for -1/2 < ν ≤ 1/2. The estimation of the term γ -1 2 (u n -(I h u) n ) 0,Γ C can be done in a very similar manner to [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]. Indeed, let E be an edge of a triangle T ∈ T h on Γ C :

γ -1 2 (u n -(I h u) n ) 0,E ≤ Ch -1 2 T h 1+ν T u n 1+ν,E = Ch 1 2 +ν u n 1+ν,E ,
(see [START_REF] Ern | Theory and practice of finite elements[END_REF] for instance). By summation on all the edges and the trace theorem, it results:

γ -1 2 (u n -(I h u) n ) 0,Γ C ≤ Ch 1 2 +ν u n 1+ν,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω . (23)
From [START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF] (see also the details in [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF]), the following estimate also holds:

γ 1 2 σ n (u -I h u) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω . (24)
We conclude by inserting the three estimates ( 22)-( 24) into ( 13).

Actually we are not able to obtain estimates for the displacements in the L 2 -norm uu h 0,Ω (and also uu h 0,Γ C ) by using the Aubin-Nitsche argument as it is achieved in the linear case (see [START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF]). Note also that the L 2 (Ω)-norm estimates for contact problems are not easy to prove and there are to our knowledge only few estimates (see [START_REF] Coorevits | Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF]). Nevertheless, we can easily obtain the following error estimate on the weighted L 2 (Γ C )-norm on the normal constraint γ 1 2 σ n (u -u h ) 0,Γ C (note that σ n (u h ) = -1 γ [P γ (u h )] + on Γ C contrary to the continuous case).

Corollary 3.6. Suppose that the solution u of Problem (3) belongs to (H 3 2 +ν (Ω)) 2 with 0 < ν ≤ 1/2. Suppose also that the parameter γ 0 is sufficiently small. The solution u h of Problem (8) satisfies the following error estimate:

γ 1 2 σ n (u -u h ) 0,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (25)
with C > 0 a constant, independent of h and u.

Proof: We use ( 24), ( 11), ( 22) and ( 21) to establish the bound: Remark 3.7. Although contact problems are known to be limited in regularity (the regularity H 5/2 (Ω) can not generally be passed beyond for such inequality problems), the use of quadratic finite element methods can be of interest in particular for the most regular solutions lying in H s (Ω), 2 < s < 5/2 (as it is considered in e.g., [START_REF] Belhachmi | Quadratic finite element approximation of the Signorini problem[END_REF][START_REF] Hild | Quadratic finite element methods for unilateral contact problems[END_REF][START_REF] Popp | Dual quadratic mortar finite element methods for 3D finite deformation contact[END_REF][START_REF] Wohlmuth | An abstract framework for a priori estimates for contact problems in 3D with quadratic finite elements[END_REF]). It is easy to check that the trace inequality [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF] is satisfied in the case of quadratic elements so that Theorem 3.4 still holds. Moreover, since estimates ( 22)-( 24) remain true in the case of quadratic finite elements (see [START_REF] Fritz | A comparison of mortar and Nitsche techniques for linear elasticity[END_REF] and [START_REF]A comparison of mortar and Nitsche techniques for linear elasticity[END_REF] for estimate ( 24)), we come to the conclusion that the optimal convergence rates ( 21) and ( 25) hold for 0 < ν ≤ 3/2 when quadratic finite elements are used.

γ 1 2 σ n (u -u h ) 0,Γ C ≤ γ 1 2 σ n (u -I h u) 0,Γ C + γ 1 2 σ n (I h u -u h ) 0,Γ C ≤ Ch

Conclusion and perspectives

In this paper we propose a new finite element approximation for the unilateral contact problem in two-dimensional linear elasticity and we achieve the corresponding numerical analysis. This analysis allows us to obtain the first optimal a priori error estimate under H s regularity (3/2 < s ≤ 2) on the solution without any additional assumption. The numerical experiments as well as the extension to the three-dimensional case should be considered in a forthcoming study.

1 2 +ν u 3 2 +ν,Ω + Cγ 1 2 0

 12 ( I h uu 1,Ω + uu h 1,Ω )

Acknowledgements.

We thank Yves Renard for helpful discussions and comments.


David Buchman
email: davidbuc@cs.ubc.ca

Mark Schmidt
email: mark.schmidt@inria.fr

Shakir Mohamed
email: shakirm@cs.ubc.ca

David Poole
email: poole@cs.ubc.ca

Nando De Freitas

On Sparse, Spectral and Other Parameterizations of Binary Probabilistic Models

This paper studies issues relating to the parameterization of probability distributions over binary data sets. Several such parameterizations of models for binary data are known, including the Ising, generalized Ising, canonical and full parameterizations. We also discuss a parameterization that we call the "spectral parameterization", which has received significantly less coverage in existing literature. We provide this parameterization with a spectral interpretation by casting loglinear models in terms of orthogonal Walsh-Hadamard harmonic expansions. Using various standard and group sparse regularizers for structural learning, we provide a comprehensive theoretical and empirical comparison of these parameterizations. We show that the spectral parameterization, along with the canonical, has the best performance and sparsity levels, while the spectral does not depend on any particular reference state. The spectral interpretation also provides a new starting point for analyzing the statistics of binary data sets; we measure the magnitude of higher order interactions in the underlying distributions for several data sets.

Introduction

Log-linear models are used as efficient parameterizations for probability distributions in a wide variety of applications [START_REF] Bishop | Discrete multivariate analysis: Theory and practice[END_REF][START_REF] Whittaker | Graphical models in applied multivariate analysis[END_REF][START_REF] Steffen | Graphical models[END_REF][START_REF] Wasserman | All of statistics[END_REF][START_REF] Koller | Probabilistic graphical models: Principles and techniques[END_REF]. Due to their relatively small number of parameters, pairwise log-linear models have sometimes been advocated in scenarios where limited data is available (Whittaker, 1990, §9.3). However, pairwise models only focus on unary and pairwise statistical properties of the data; a pairwise assumption can be restrictive if higher-order moments of the data are important and we have sufficient training examples available to reliably estimate these higher-order statistics. Despite this fact, almost all previous work on structure learning with 1 -regularization has made the pairwise assumption, with a few exceptions. [START_REF] Dahinden | Penalized likelihood for sparse contingency tables with an application to full-length cDNA libraries[END_REF] consider log-linear models of discrete data where all potentials up to a fixed order are considered with (group) 1 -regularization to learn the structure, [START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF] address hierarchical log-linear models with (overlapping-group) 1regularization, and [START_REF] Ding | Learning higher-order graph structure with features by structure penalty[END_REF] further consider hierarchical log-linear models with covariates. This paper makes two contributions. First, it develops a spectral interpretation of an existing parameterization of full log-linear models in terms of orthogonal Walsh-Hadamard bases [START_REF] Beer | Walsh transforms[END_REF]. Although this parameterization has been used occasionally, to the best of our knowledge, its interpretation in terms of orthogonal expansions is new. We should also point out that a spectral expansion of probability measures in terms of Walsh basis functions has been studied in the context of univariate and bivariate probability density functions [START_REF] Maqusi | Walsh series expansions of probability distributions[END_REF], and that orthogonal Walsh-Hadamard bases were used for analyzing binary factorial designs (Rockmore, 1997, §2.2). Our focus on log-linear models enables us to use this parameterization to study the "natural statistics", or spectrum, of several popular binary data sets, in the same fashion that researchers have investigated the natural statistics of images and other signals. Our results show that in this spectral domain, lower-order potentials tend to have much more weight than higher-order potentials. This result is very intuitive, but it is not obvious. For example, one can construct distributions which only have high-order potentials, such as the "parity distri-bution":

p parity (x) ∝ e i xi ,
x i ∈ {-1, +1}.

The parameterization we describe here, which we coin the spectral parameterization, was suggested previously by [START_REF] Bishop | Discrete multivariate analysis: Theory and practice[END_REF]. That work examined the 'full' parameterization of general discrete probabilistic models, and suggested adding constraints to the model parameters in order to obtain a minimal parameterization. For binary variables, the constraints leave a single degree of freedom for each potential, and can thus be modeled using a single parameter -essentially leading to the spectral parameterization. This initial work failed to notice the harmonic properties of the parameterization and does not make a connection to Walsh-Hadamard expansions.

The spectral parameterization has received minimal attention in the machine learning and statistics literature. Most recently, both the spectral representation and the Hadamard transform appeared in a paper in the field of haplotype inference [START_REF] Kato | A haplotype inference method based on sparsely connected multi-body Ising model[END_REF]. [START_REF] Kato | A haplotype inference method based on sparsely connected multi-body Ising model[END_REF], however, do not use the Hadamard transform to describe the harmonic structure of the parameterization, but rather use it as a computational tool to geometrically average approximated marginal distributions in a computationally efficient manner, and to find maximum-probability states given a model with a specific structure. This work, which is the one we believe is closest to ours, does not address the issue of learning. However, when used for learning, the spectral parameterization implicitly defines a new category of priors -priors over the spectral parameters of the distribution.

This spectral interpretation is important for several reasons. First, it allows us to conduct empirical analyses similar to the ones that are carried out for other types of data using the Fourier and wavelet transforms.

Second, it provides an important bridge between harmonic analysis and related fields such as compressed sensing, and the problem of learning the parameters and structure of discrete probabilistic graphical models. Although we do not explore this theoretical connection in this paper, we believe it could potentially lead to new theoretical results about the sample complexity of sparse undirected probabilistic models [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Abbeel | Learning factor graphs in polynomial time and sample complexity[END_REF][START_REF] Ravikumar | High-dimensional Ising model selection using 1 -regularized logistic regression[END_REF].

The second contribution of this paper is to present a comprehensive comparison of different parameterizations of undirected discrete probabilistic models, including the full, Ising, canonical and spectral parameterizations. This comparison is done in the context of learning with several types of standard and group tions. A comparative analysis of such a wide range has not been previously undertaken in the literature. We believe researchers will find it useful as it sheds light on the choice of parameterization and variant of 1 regularization that are used when learning discrete probabilistic models.

Log-linear Models

Given n binary random variables x ∈ {-1, +1} n , we can express a positive joint probability as a globally normalized product of potential functions exp(φ A (x A )) defined for each possible subset A of S {1, 2, . . . , n}:

p(x) 1 Z A⊆S exp(φ A (x A )).
The normalizing constant Z ensures that the distribution sums to one. When the logarithm of each potential is linear in the parameters of the potential, and the normalizing constant is encoded as -log(Z) = φ ∅ (x), we can express the model in the standard log-linear form [START_REF] Bishop | Discrete multivariate analysis: Theory and practice[END_REF]:

log p(x) = A⊆S φ A (x A ) = A⊆S w T A f A (x A ), (1)
where f A (x A ) is a feature vector derived from x A . The vector w is the set of log-linear parameters. We use the short-hand w A to refer to all the parameters associated with the function φ A (x A ), and use w to refer to the concatenation of all w A .

Undirected probabilistic graphical models can be derived from log-linear models by connecting node i to j if variables x i and x j co-occur in some φ term. More precisely, a log-linear model is graphical if there is a non-zero φ term for every clique in the graph [START_REF] Wasserman | All of statistics[END_REF]. Conversely, φ A (x) = 0 if {s, t} ⊆ A for t = s, and (s, t) is not an edge. Second, a log-linear model is hierarchical if φ A = 0 and A ⊂ B implies that φ B = 0. Typically, attention is restricted to the class of hierarchical log-linear models due to the interpretability of their sparsity pattern in terms of conditional independence [START_REF] Whittaker | Graphical models in applied multivariate analysis[END_REF].

In practice, it is typically not feasible to include a potential φ A (x A ) for all 2 n subsets. Removing the potential φ A (x A ) from the model is equivalent to setting it to zero for all values of x A , or equivalently setting all elements of w A to zero. For example, we obtain the class of pairwise models if we enforce w A = 0 for all A with a cardinality greater than two. This effectively nullifies the effects of higher-order statistics present in the data on the model.

Parameterizations

Full Parameterization

With the full parameterization of log-linear models, pairwise potentials have the form:

φ ij (x i , x j ) = s1 s2 I <s1,s2> (x i , x j )w ijs1s2 ,
where s 1 , s 2 ∈ {-1, +1} and the indicator I <s1,s2> (x i , x j ) is one if x i = s 1 and x j = s 2 , and zero otherwise. For three-way potentials we have:

φ ijk (x i , x j , x k ) = s1,s2,s3 I <s1,s2,s3> (x i , x j , x k )w ijks1s2s3
and similarly for higher-order potentials. In general, if A contains k elements that can each take 2 values, φ A (x A ) will have 2 k parameters w A .

Ising Parameterizations

The Ising parameterization allows us to reduce the model complexity and consider potentials with a single parameter:

φ ij (x i , x j ) = s I <s,s> (x i , x j )w ij .
Generalized Ising models allow potentials to take c parameters rather than a single one, where c is the number of values that variables can take. For binary variables, c = 2, which we assume throughout. For pairwise potentials we have:

φ ij (x i , x j ) = s I <s,s> (x i , x j )w ijs .
And similarly, for three-way potentials:

φ ijk (x i , x j , x k ) = s I <s,s,s> (x i , x j , x k )w ijks .

Canonical Parameterizations

Another strategy for decreasing the number of parameters is the canonical parameterization [START_REF] Koller | Probabilistic graphical models: Principles and techniques[END_REF][START_REF] Steffen | Graphical models[END_REF]. To understand this parameterization, consider the set of features:

f A (x A ) = 1 iff x A = 1 0 otherwise,
where 1 denotes a vector of ones. It is a simple exercise to verify that for this choice one has:

log p(x A = 1, x A = -1) = C⊆A w C ,
where A is the set complement of A, and all f A and w A are scalars. Solving the linear system of 2 n equations for w yields:

w A = C⊆A (-1) |A-C| log p(x C = 1, x C = -1). (2)
This is a special case of the Möbius inversion lemma.

We have chosen here the specific reference state 1, which may be replaced with an arbitrary reference state (see [START_REF] Koller | Probabilistic graphical models: Principles and techniques[END_REF] for discussion).

In this canonical parameterization we need only one parameter per potential, hence 2 n parameters are needed to represent the entire distribution, one of which is the normalizing constant. The Ising parameterization also requires 2 n parameters, but it is not complete (it can't represent all positive distributions). The other parameterizations require more parameters: the full parameterization requires 3 n , since each parameter corresponds to an assignment of one of {not in A, -1, +1} to each variable, and generalized Ising requires 2 n c = 2 n+1 (c = 2 for binary variables).

Spectral Parameterization

We consider an alternative parameterization that has received less attention in the literature, but which has important properties of theoretical and practical significance. We refer to this fourth parameterization as the spectral parameterization. The motivation for considering such a parameterization was made aptly by Koller and Friedman (2009, §4.4.2.1, pg. 130): "... canonical parameters are not very intuitive, highlighting yet again the difficulties of constructing a reasonable parameterization of a Markov network by hand."

With n = 3, the spectral parameterization of log p(x) is as follows:

A⊆S φ A (x A ) = w 0 + w 1 x 1 + w 2 x 2 + w 3 x 3 + w 12 x 1 x 2 +w 13 x 1 x 3 + w 23 x 2 x 3 + w 123 x 1 x 2 x 3 .
That is, the spectral parameterization has only one parameter per potential. In general, we have:

log p(x) = A⊆S w A i∈A x i . (3)
Again, we only need 2 n parameters to represent the entire distribution. We are required to assess 2 n -1 parameters for all w A , where A = ∅; w ∅ can be used to ensure values sum to 1.

Although this spectral parameterization is not widely discussed in existing literature, it has been known (for general discrete log-linear models) for several decades [START_REF] Bishop | Discrete multivariate analysis: Theory and practice[END_REF]. Here, we show that the spectral representation provides a map between log-linear models and orthogonal expansions in terms of Walsh-Hadamard bases. As an example, consider a log-linear model with two random variables: q(x 1 , x 2 ) = w 0 + w 1 x 1 +w 2 x 2 +w 12 x 1 x 2 , where q(x 1 , x 2 ) log p(x 1 , x 2 ). We can ground this expression for all realizations of the random variables to obtain:

q(+1, +1) = w 0 + w 1 + w 2 + w 12 q(-1, +1) = w 0 -w 1 + w 2 -w 12 q(+1, -1) = w 0 + w 1 -w 2 -w 12 q(-1, -1) = w 0 -w 1 -w 2 + w 12 ,
which written in matrix notation is: q = 2H 2 w. H 2 is the Hadamard matrix for two variables:

H 2 = 2 -1    1 1 1 1 1 -1 1 -1 1 1 -1 -1 1 -1 -1 1    (4)
For n variables, the Hadamard matrix has entries

(H n ) ij = 2 -n/2 (-1) i•j
, where i • j is the bitwise dot product of the binary representations of the numbers i and j indexing the 2 n possible realizations of x (e.g., [START_REF] Pratt | Hadamard transform image coding[END_REF]). The rows (and columns) of the Hadamard matrix are orthogonal Walsh functions h(•) [START_REF] Beer | Walsh transforms[END_REF]. Let x (t) denote the t-th realization of x for t = 0, . . . , 2 n -1. Then, we can rewrite our linear system in terms of these basis functions as follows:

q(x (t) ) = A⊆S h A (x (t) )w A .
(5) This is the forward Walsh-Hadamard transform.

Since the Walsh functions are orthogonal, that is

2 n -1 t=0 h A (x (t) )h B (x (t) ) = 2 n I A=B ,
we can multiply both sides of the forward transform by h A and sum over t to obtain the reverse transform:

w A = 2 -n 2 n -1 t=0 q(x (t) )h A (x (t) ). (6)
Equations 5 and 6 define the Walsh-Hadamard transform. H n is symmetrical, and as with the Fourier transform, the Walsh-Hadamard transform pair can be expressed as:

q = 2 n/2 H n w and w = 2 -n/2 H n q. (7)
Likewise, it can be computed using the Fast Walsh-

Hadamard transform (FWHT) in O(n2 n ) (that is, m log m instead of m 2
, where m = 2 n ).

Comparing different parameterizations

Table 1 summarizes the properties of the parameterizations under consideration. We consider the generalization of the parameterizations to variables with c values, however the spectral parameterization is considered only for c = 2. "Canonical" refers to the canonical parameterization with a general reference state, while C1 and C2 use reference states 1 and -1 respectively.

A parameterization is complete if it can represent any positive distribution. A parameterization is minimal if it has no redundant parameters. The next column corresponds to symmetry with respect to the values of individual variables. Assume we flip the values of a specific variable. How would the parameters change, in order to reflect the modification to the distribution? Symmetric parameterizations would only need to make trivial modifications to their parameters, such as exchanging some of them with each other, or flipping their signs. Non-symmetric parameterizations require much more complex calculations. In particular, a sparse w may then become dense. The next column offers a similar analysis, for exchanging the values of two different variables. Finally, "uniquely defined" means the distribution is uniquely defined as a function of the parameter vector w, without additional external information such as an identity of a reference state.

Generalized Ising is complete only for binary variables, for which it has two times the minimal number of parameters. Ising has the minimal number, but with redundancy, as it is not complete. The canonical parameterization requires the specification of a reference state. Each possible reference state can be seen as giving rise to a different parameterization, therefore the canonical parameterization is not uniquely defined. The spectral parameterization is the only one that is minimal and symmetric w.r.t. values. In addition, it is also complete, symmetric w.r.t. variables, and it is unique.

Sparse Regularization for Structure Learning

All complete parameterizations share the same maximum-likelihood (ML) estimate. Once regularization is introduced, each combination of parameterization and regularizer defines a different prior over the space of distributions. Aside from computational issues, this choice also impacts prediction performance.

In our experiments, we consider four different forms of sparse regularization. The first is the standard 1regularization, used in the context of log-linear models by [START_REF] Lee | Efficient structure learning of Markov networks using L 1 -regularization[END_REF]:

||w A || 1 j |w (j) A |, ||w A || 2 j (w (j) A ) 2 1/2 , max w n i=1 log p(x i |w) - A⊆S λ A ||w A || 1 . (8)
The 1 -regularization encourages the parameter vector w to be sparse. To illustrate why one has to devote careful thinking to the pairing of parameterization and

p(x 1 , x 2 ) = δ 1 δ 2 δ 3 1 -δ 1 -δ 2 -δ 3
. If we place an 1 regularizer on the parameters δ i , then while the first three cells are forced to shrink to zero, the fourth one is forced to grow to 1. Thus, a naive choice of parameterization with 1 -regularization can induce arbitrary and potentially unwanted behavior.

For parameterizations where each w A is a scalar (Ising, canonical, and spectral), Equation 8encourages the removal of factors from the model. However, for parameterizations where each w A has more than one element (generalized Ising, full), 1 -regularization does not encourage entire factors w A to be equal to the zero vector simultaneously. If this is desired, one can use group 1 -regularization:

max w n i=1 log p(x i |w) - A⊆S λ A ||w A || 2 , (9)
which was first used in the context of log-linear models by [START_REF] Dahinden | Penalized likelihood for sparse contingency tables with an application to full-length cDNA libraries[END_REF]. In cases where w A is a scalar, Equations 8 and 9 are equivalent. Although group 1 -regularization encourages the removal of entire factors even for parameterizations that have more than one parameter per factor, it still does not directly encourage conditional independencies in the learned distribution. This is because it may estimate a sparse but non-hierarchical model, where w A = 0 in the solution of Equation 8 or 9 but this is not true under a re-parameterization of the same distribution (for non-minimal parameterizations). To encourage the removal of factors and the underlying model to be hierarchical, we also use the hierarchical group We use this regularizer only for the spectral parameterization, and for p = 1, 2. For p = 1, this is identical to the standard 1 regularizer with a constant λ A = λ. This regularizer penalizes all parameters independently and identically, a fact we rely on in Section 6.

Experiment 1: Comparing Parameterizations

Here, we compare different parameterizations and regularizers on several data sets in terms of test-set negative log-likelihood. We followed the same training protocol as [START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF]. For crossvalidation of the hyperparameter, we divided the data into equal training, validation and testing sets. We repeated this with 10 different splits to control for variability in the results. We refer the reader to [START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF] for further details. We experimented with a wide range of benchmark data sets, including the Yeast [START_REF] Elisseeff | Kernel methods for multi-labelled classification and categorical regression problems[END_REF], USPS (Roweis), Jokes [START_REF] Goldberg | Eigentaste: A constant time collaborative filtering algorithm[END_REF], Flow [START_REF] Sachs | Causal protein-signaling networks derived from multiparameter single-cell data[END_REF], Rochdale [START_REF] Whittaker | Graphical models in applied multivariate analysis[END_REF], Czech [START_REF] Edwards | A fast procedure for model search in multidimensional contingency tables[END_REF] and NLTCS [START_REF] Erosheva | Describing disability through individuallevel mixture models for multivariate binary data[END_REF] data sets. We found the same trends in the results across these data sets, so for clarity of presentation we focus only on the NLTCS, Yeast (labels only), and USPS (central 16 pixels). We purposely chose data sets with a limited number of variables, for which it is possible to compute the partition function by enumeration, so as to avoid introducing artifacts in the comparison due to approximations.

We studied the Ising (I), canonical (C1, C2 and CR), spectral (S), generalized Ising (GI), and full (F) parameterizations. We use C1 and C2 to denote the parameterizations using reference states 1 and -1, respectively. CR uses a random reference state. Note that Schmidt and Murphy (2010) consider the GI and F parameterizations, while [START_REF] Ding | Learning higher-order graph structure with features by structure penalty[END_REF] consider the C1 parameterization. Each parameterization was assessed using the standard 1 regularization (Equation 8). The over-complete parameterizations (generalized Ising and full) were also tested with the group 1 regularization (Equation 9), where each group refers to the parameters w A associated with a particular set of variables A. For the other parameterizations group 1 regularization is identical to the standard 1 regularization. In addition, we tested all parameterizations using the hierarchical group 1 regularization (Equation 10), where each group refers to the set of parameters associated with a particular A and all its supersets. We set λ A = I |A|≥2 2 |A|-2 λ. The hyperparameter λ was chosen using the validation sets. We used the optimization software and hierarchical-search strategy of [START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF]. Finally, we examined the spectral parameterization with the flat regularizer of Equation 11 with p = 1 ("flat 1 ") and p = 2 ("flat 2 2 "), along with a pseudo-counts (PC) estimator. The PC estimator uses Dirichlet smoothing to estimate the log-probabilities from data.

Figure 1 shows a comparison of the predictive performance for several combinations of regularizers and parameterizations, in which lower values mean better performance. Figure 2 shows the corresponding sparsity levels of the learned models. In the figures, we use '-group' to indicate that the parameterization is subject to the group 1 regularizer and '-h' for the hierarchical regularizer. Due to space considerations, not all combinations are shown, but the ones shown do capture most of the important information.

The flat and PC models under-perform (Figures 1(b), 1(c)) and produce dense models (Figures 2(b), 2(c)). We do not consider them attractive for practical modeling purposes, but we added them to the comparison to facilitate the analysis of the next section. For all parameterizations, we have found no significant difference between the standard, group and hierarchical 1 regularizers in terms of predictive performance, as shown in Figure 1(a) for the NLTCS data set. An exception is the full parameterization, for which standard 1 under-performed relative to the other regularizers (Figures 1(a-c)). A possible explanation could be the much larger number of parameters in the full parameterization: 3 n , versus 2 n for the Ising, canonical and spectral parameterizations, and 2 n+1 for the generalized Ising parameterization. More importantly, a much smaller percentage of its parameters correspond to the lower-cardinality potential functions, compared to the other parameterizations. For all parameterizations, standard 1 produced sparser models than group 1 and hierarchical group 1 (see Figure 2). Figure 1 shows that the Ising parameterization performs poorly. This is not surprising, considering the fact that it is not complete. The canonical, spectral and generalized Ising parameterizations seem to have performance similar to each other. However, both the canonical and spectral parameterizations seem to do best in terms of both sparsity and predictive performance. These are also the only parameterizations that are both complete and minimal. The performance of the canonical parameterization may depend on the particular reference state. This is particularly evident in Figure 2(c).

For problems that can naturally be described as having a base state with rare deviations from this state, the canonical parameterization with this state as the reference state may be an appealing choice. In other cases, our experiments suggest either trying different reference states for the canonical parameterization, or using the spectral parameterization instead. Our experiments also suggest that standard 1 regularization is preferred due to its simplicity and sparser resulting models.

6 Experiment 2: The Statistics (Spectrum) of Binary Data Sets

The spectral decomposition of binary distributions enables us to answer the question: What are the statistics of commonly used binary data sets? Researchers have used tools such as the Fourier and wavelet transforms to understand the properties of natural data, such as images. By casting binary distributions in terms of the Walsh-Hadamard transform, we can now develop a method for estimating the spectral coefficients of the underlying distribution. This experiment reveals that low-order statistics are much more significant than higher-order ones for binary data sets routinely used in practice.

We trained the spectral parameterization with the flat 1 and 2 2 regularizers on all the data sets mentioned in the previous experiment. We also computed the pseudo-counts estimate, and used FWHT to transform the log-probabilities from the pseudo-counts model to the spectral parameterization. Readers may wonder, in light of the previous experiment, why we use the flat regularizers. Flat regularizers do not impose a prior bias for smaller-cardinality factors, so it is interesting to see that without this bias there is still a preference for smaller factors: indeed this is what our experiments indicate. This then provides a rationale for a prior bias towards considering only low-order factors. The parameter λ was estimated using 10-fold cross-validation. We used the size of the factors |A| to group parameters. For each group, we calculated the mean magnitude E(|w A |) of the parameters in the group (for example, the average magnitude of the parameters for all factors of size 3). For the flat 1 estimates, we also calculated the density (percentage of non-zero parameters). The results in Figure 3 show that the average magnitude of the spectral parameters diminishes with increasing factor size |A|. This behavior is especially pronounced for the 1 regularizer. According to our results in the previous experiment, the flat 1 regularizer allows for a much better fit than the flat We can see that parameters associated with a large |A| are sparse, especially for the data sets that have more parameters (in Figure 3(b), Jokes has 100% density for |A| = 10, however, there is only one such parameter, so this is not significant). In addition, the "spectra" (Figures 3(a), 3(c)) illustrate that pair-wise and three-way models can capture a large portion of the variability in these data sets. This probably accounts for their popularity. The idea of ignoring higher-order potentials was also recently considered by [START_REF] Jalali | On learning discrete graphical models using group-sparse regularization[END_REF], who analyzed the effect of fitting a pairwise model to data in terms of consistency of estimating the presence of higher-order interactions in hierarchical models. Nonetheless, our results also show that while adding higher order potentials may lead to an improvement in performance, this potential improvement comes with an additional computational cost that might not be justifiable in practical domains. Our spectral analysis of binary distributions provides justification for the standard machine learning and statistical practices of ignoring higher-order potentials. Importantly, it enables researchers to assess the cost of ignoring these terms.

I C 1 C 2 C R S G I F G I -g r o u p F -g r o u p I -h C 1 -h C 2 -h C R -h S -h G i- h F -h
3 I C 1 C 2 C R S G I F G I -g r o u p F -g r o u p I -h C 1 -h C 2 -h C R -h S -h G i- h F -h Number of non-zero parameters (a) NLTCS 10 2 10 3 10 4 I C 1 C 2 C R S G I F F -g r o u p F -h f la t L 1 f la t L 2 P C (b) Yeast 10 3 10 4 10 5 I C 1 C 2 C R S G I F F -g r o u p F -h f la t L 1 f la t L 2 P C (c) USPS

Conclusions and Future Work

We have presented a comparison of different parameterizations for discrete probabilistic log-linear models, when learning their parameters and structure with sparsity-promoting regularizers. We found that the spectral parameterization is one of the best performing, and were able to interpret it as a harmonic series expansion of orthogonal Walsh-Hadamard bases. Since this interpretation brings closer the fields of harmonic analysis and discrete probabilistic modeling, it opens many doors for future research. We have already used the spectral parameterization to study the statistics, or the spectrum, of popular binary data sets. However, we believe this theoretical connection could be exploited in connection with compressed sensing to provide sample complexity theorems for learning discrete probabilistic models. For example, if we could measure log-probabilities directly, then the theorems of compressed sensing could be easily adapted to this domain to estimate the number of measurements needed to reconstruct the sparse probabilistic model. However, since log-probabilities cannot typically be measured directly, some theoretical challenges lie ahead. In future work, we also plan to build on the work of [START_REF] Bishop | Discrete multivariate analysis: Theory and practice[END_REF] to extend the results to distributions over random variables with more than two values.

  x i |w) -λ j |w (j) | p . (11)

  10

Figure 2 :

 2 Figure 2: Model sparsity for different parameterizations and regularizers.

  a),(c) Mean parameter magnitude (log-scale), and (b),(d) percentage of non-zero parameters, for different factor sizes |A|. Data sets are split into two groups to allow clearer plots. Solid lines correspond to the flat 1 regularizer, dashed lines to the flat 22 regularizer and dotted lines for the pseudo-counts estimator. The difference between the pseudo-counts and 2 2 regularizer is hardly discernible.

Table 1 :

 1 Properties of parameterizations.

	Parameterization #Parameters Complete Minimal Symm w.r.t. Symm w.r.t. Uniquely
					values	variables	defined
	Full	(c + 1) n	Yes	No	Yes	Yes	Yes
	Ising	2 n	No	No	No	Yes	Yes
	Generalized Ising	2 n c	For c = 2	No	No	Yes	Yes
	Canonical	c n	Yes	Yes	No	No	No
	Canonical (C1/C2)	c n	Yes	Yes	No	Yes	Yes
	Spectral (for c = 2)	2 n	Yes	Yes	Yes	Yes	Yes
	sparse regularizer let us look at a very simple exam-
	ple where undesirable behavior can easily arise. As-
	sume we have a model with two binary random vari-
	ables described by the following joint probability table:

  Test set negative log-likelihood for different parameterizations and regularizers.

	Test set negative log-likelihood	(a) NLTCS	3300 3350 3400 3450 3500 3550 3600 3650	I C 1 C 2 C R	S	G I (b) Yeast F F -g r o u p F -h f la t L 1 f la t L 2 P C	2.25 2.3 2.35 2.4 2.45 2.5 2.55	I C 1 C 2 C R x 10 4	S	G I (c) USPS F F -g r o u p F -h f la t L 1 f la t L 2 P C
		Figure 1:

Appearing in Proceedings of the 15 th International Conference on Artificial Intelligence and Statistics (AISTATS) 2012, La Palma, Canary Islands. Volume XX of JMLR: W&CP XX. Copyright 2012 by the authors.

regularizers. These experiments can be seen as an extension of the ones conducted by[START_REF] Schmidt | Convex structure learning in log-linear models: Beyond pairwise potentials[END_REF] to a much broader range of parameteriza-

Acknowledgements

This research was supported by NSERC grants to N. de Freitas and D. Poole. M. Schmidt was supported by the European Research Council (SIERRA-ERC-239993). S. Mohamed is supported by the Canadian Institute for Advanced Research (CIFAR).